Z=[];
d0=1e-8;
for a=linspace(0.3,2,70)
lsum=0;
x=1;y=1;z=1;
x1=1;y1=1;z1=1+d0;
for i=1:100
[T1,Y1]=ode45('Lorenz3',1,[x;y;z;a;.1;2;0]);
[T2,Y2]=ode45('Lorenz3',1,[x1;y1;z1;a;.1;2;0]);
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
a=linspace(0.3,2,70);
plot(a,Z,'-');
title('模拟 Lorenz 系统最大lyapunov指数谱')
xlabel('parameter a')
ylabel('lyapunov exponents')
|
No comments:
Post a Comment