Matlab画Lorenz系统的最大李雅普诺夫指数图.docx

PAGE PAGE 1 / 3 Lorenz 系统 文档分两个文件方程m 文件和计算L 指数 m 文件分开写，复制粘贴即可运行 matlab2012a，改写方程文件和参数即可算自己的系统，其中最大 L 指数用的是经典的柏内庭（G.Benettin）计算方法，准确快速无误！附计算结果图！！ 方 程 m 文 件 ： function dX = Loren(t,X) global a; %变量不放入参数表中global b; global c; x=X(1); y=X(2); z=X(3); % Y 的三个列向量为相互正交的单位向量 %输出向量的初始化dX = zeros(6,1); % Lorenz 吸引子dX(1)=a*(y-x); dX(2)=x*(b-z)-y; dX(3)=x*y-c*z; end 计算最大L 指数文件Z=[]; global a; global b; global c; a=10; c=8/3; d0=1e-7; for b=linspace(0,500,500) lsum=0; x=1;y=1;z=1; x1=1;y1=1;z1=1+d0; for i=1:100 [T1,Y1]=ode45(Loren,1,[x;y;z;16;b;4]); [T2,Y2]=ode45(Loren,1,[x1;y1;z1;16;b;4]); n1=length(Y1);n2=length(Y2); x=Y1(n1,1);y=Y1(n1,2);z=Y1(n1,3); x1=Y2(n2,1);y1=Y2(n2,2);z1=Y2(n2,3); d1=sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2); x1=x+(d0/d1)*(x1-x); y1=y+(d0/d1)*(y1-y); z1=z+(d0/d1)*(z1-z); if i50 lsum=lsum+log(d1/d0); end end Z=[Z lsum/(i-50)]; end b=linspace(0,500,500); plot(b,Z,-); title(JD_{1}系统最大 lyapunov 指数) xlabel(parameter b),ylabel(The largest Lyapunov exponents); grid on; 结果图