精品文档
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 i>50 精品文档
精品文档
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; 结果图
精品文档