matlab求最大李雅普诺夫(Lyapunov)指数程序

求解系统的Lyapunov指数谱程序Lyapunov指数是描述时序数据所生成的相空间中两个极其相近的初值所产生的轨道，随时间推移按指数方式分散或收敛的平均变化率。任何一个系统，只要有一个Lyapunov大于零，就认为该系统为混沌系统。李雅普诺夫指数是指在相空间中相互靠近的两条轨线随着时间的推移，按指数分离或聚合的平均变化速率。一chen系统的Lyapunov指数谱functiondX=Chen2(t,X)%Chen吸引子，用来计算Lyapunov指数%dx/dt=a*(y-x)%dy/dt=(c-a)*x+c*y-x*z%dz/dt=x*y-b*zglobala;%变量不放入参数表中globalb;globalc;x=X(1);y=X(2);z=X(3);%Y的三个列向量为相互正交的单位向量Y=[X(4),X(7),X(10);X(5),X(8),X(11);X(6),X(9),X(12)];%输出向量的初始化dX=zeros(12,1);%Lorenz吸引子dX(1)=a*(y-x);dX(2)=(c-a)*x+c*y-x*z;dX(3)=x*y-b*z;%Lorenz吸引子的Jacobi矩阵Jaco=[-aa0;c-a-zc-x;yx-b];dX(4:12)=Jaco*Y;Z1=[];Z2=[];Z3=[];globala;globalb;globalc;b=3;c=28;fora=linspace(32,40,100);y=[1;1;1;1;0;0;0;1;0;0;0;1];lp=0;fork=1:200[T,Y]=ode45('Chen2',1,y);y=Y(size(Y,1),:);y0=[y(4)y(7)y(10);y(5)y(8)y(11);y(6)y(9)y(12)];y0=GS(y0);mod(1)=norm(y0(:,1));mod(2)=norm(y0(:,2));mod(3)=norm(y0(:,3));lp=lp+log(abs(mod));y0(:,1)=y0(:,1)/mod(1);y0(:,2)=y0(:,2)/mod(2);y0(:,3)=y0(:,3)/mod(3);y(4:12)=y0';endlp=lp/200;Z1=[Z1lp(1)];Z2=[Z2lp(2)];Z3=[Z3lp(3)];enda=linspace(32,40,100);plot(a,Z1,'-',a,Z2,'-',a,Z3,'-');title('LyapunovexponentsofChen')xlabel('b=3,c=28,parametera'),ylabel('lyapunovexponents')gridon以上是三个变量的Lyapunov指数谱，下面是最大的Lyapunov指数谱：Z=[];d0=1e-8;fora=linspace(32,40,80)lsum=0;x=1;y=1;z=1;x1=1;y1=1;z1=1+d0;fori=1:100[T1,Y1]=ode45('Chen',1,[x;y;z;a;3;28]);[T2,Y2]=ode45('Chen',1,[x1;y1;z1;a;3;28]);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);ifi50lsum=lsum+log(d1/d0);endendZ=[Zlsum/(i-50)];enda=linspace(32,40,80);plot(a,Z,'-');title('Chen系统最大lyapunov指数')xlabel('parametera'),ylabel('lyapunovexponents')二模拟Lorenz系统最大lyapunov指数谱functionly=jose_ly(b,k)%thelargestlyapunovexponentofjosephson%k迭代步数，b参数%方程如下:%θ''+G*θ'+sinθ=I+A*sin(ωt)+αsin(βωt)%变化：%dx=y%dy=-G*y-sin(x)+I+A*sin(w*t)+a*sin(b*w*t)%%Example:%ly=jose_ly(0,800)%%Author:LDYU%Author'semail:[e-mail]%d0=1e-8;ly=0;lsum=0;x=[0;2;b];x1=[d0;2;b];fort=1:k[T1,Y1]=ode45('Josephon',[t-1,t],x);[T2,Y2]=ode45('Josephon',[t-1,t],x1);x=Y1(end,:);x1=Y2(end,:);d1=norm(x-x1);x1=x+(d0/d1)*(x1-x);lsum=lsum+log(d1/d0);endly=lsum/k;