chaoslyapunov.m

chaoslyapunov指数程序是基于非线性三维系统所得

clear;clc;
yinit = [0.82 0.29 0.48];
global a1 a2 b1 b2 b3 c1 c2 c3 M N K;
a1=0.09; a2=0.15;
b1=0.06; b3=0.07;
c1=0.2;  c2=0.5;   c3=0.4;
M=1.8;N=1;K=0.4;
k=0;
for b2=0.06:0.001:0.113   %参数变化取点
    m=0;
    k=k+1;    
y = zeros(12,1);

y(1:3) = yinit;
y(4:12) = eye(3);
initialtime = 0; 
tstep =0.01;  
finaltime =5000; 
steps =100; 
iterationtime=tstep*steps;
discarditeration=500;
discardtime=discarditeration*iterationtime; 
lp = [0 0 0];
for i=1:(finaltime-initialtime)/iterationtime-1
    T1=initialtime+(i-1)*iterationtime;
    T2=T1+iterationtime;
    tspan=[T1:tstep:T2];
    [T,Y] = ode45('chaossun_ly', tspan, y);
    
    y = Y(end,:);
    y0 = [y(4) y(7)  y(10);
          y(5) y(8)  y(11);
          y(6) y(9)  y(12)];
    [Q,R]=qr(y0); 
    if T2>discardtime
         y(4:12)=Q;
    else
         y(4:12)=eye(3);
    end
   permission=1;
   for j=1:3
      if R(j,j)==0
         permission=0;
         break;
      end
   end    
    if (T2>discardtime & permission)
       m=m+1;
       T=m*iterationtime;
       TT=i*iterationtime+initialtime;
       for j=1:3
           lp(j) = lp(j)+log(abs(R(j,j)));
       end
       
        Lyapunov1(k,i) = lp(1)/T;
        Lyapunov2(k,i) = lp(2)/T;
        Lyapunov3(k,i) = lp(3)/T;
       
    end
end
drawnow;
plot(b2,Lyapunov1(k,end),'r.');hold on;
plot(b2,Lyapunov2(k,end),'b.');hold on;
plot(b2,Lyapunov3(k,end),'k.');hold on;
end
plot([0 1],[0 0],'r:');
xlim([0.06 0.113]) %限定X轴的坐标范围
xlabel('parameter b_{2}');
ylabel('Lyapunov exponents');