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

求解系统的Lyapunov指数谱程序Lyapunov 指数是描述时序数据所生成的相空间中两个极其相近的初值所产生的轨道，随时间推移按指数方式分散或收敛的平均变化率。任何一个系统，只要有一个Lyapunov 大于零，就认为该系统为混沌系统。李雅普诺夫指数是指在相空间中相互靠近的两条轨线随着时间的推移，按指数分离或聚合的平均变化速率。一 chen系统的Lyapunov指数谱function dX = Chen2(t,X)% Chen吸引子，用来计算Lyapunov指数% dx/dt=a*(y-x)% dy/dt=(c-a)*x+c*y-x*z% dz/dt=x*y-b*zglobal a; % 变量不放入参数表中global b;global c;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 = -a a 0; c-a-z c -x; y x -b; dX(4:12) = Jaco*Y;Z1=;Z2=;Z3=;global a;global b;global c;b=3;c=28;for a=linspace(32,40,100); y=1;1;1;1;0;0;0;1;0;0;0;1; lp=0; for k=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; end lp=lp/200; Z1=Z1 lp(1); Z2=Z2 lp(2); Z3=Z3 lp(3); enda=linspace(32,40,100);plot(a,Z1,-,a,Z2,-,a,Z3,-);title(Lyapunov exponents of Chen)xlabel(b=3,c=28,parameter a),ylabel(lyapunov exponents)grid on以上是三个变量的Lyapunov指数谱，下面是最大的Lyapunov指数谱：Z=;d0=1e-8;for a=linspace(32,40,80)lsum=0;x=1;y=1;z=1;x1=1;y1=1;z1=1+d0;for i=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); if i50 lsum=lsum+log(d1/d0); endendZ=Z lsum/(i-50);enda=linspace(32,40,80);plot(a,Z,-);title(Chen 系统最大lyapunov指数)xlabel(parameter a),ylabel(lyapunov exponents)二 模拟 Lorenz 系统最大lyapunov指数谱 function ly=jose_ly(b,k)% the largest lyapunov exponent of josephson% 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% Authors email: %d0=1e-8;ly=0;lsum=0;x=0;2;b;x1=d0;2;b;for t=1:k T1,Y1=ode45(Josephon,t-1,t,x); T2,Y2=ode45(Joseph