用WOLF方法求时间延迟CHEN系统最大LYAPUNOV指数

wolf_3Dlyapunov.m

function lyapunov=wolf_3Dlyapunov(t,Y) % 用wolf方法求三维实验数据最大李雅普洛夫指数 % % 输入：演化时间间隔t，3维数据Y(:,1:3) % 输出：lyapunov指数 % % Author：yujunjie % Author's email: ustb03-07@yahoo.com.cn % lyapunov=0; YL=length(Y); % 不同的系统以下两个参数可能要适当的改变 MaxFL=.5; MaxEL=2.5; fiducial1=1; fiducial2=3; for i=4:YL if (sum(abs(Y(fiducial1,:)-Y(i,:)))<MaxFL) fiducial2=i; break; end end evolveStep=1; L1=norm(Y(fiducial1,:)-Y(fiducial2,:)); L2=norm(Y(fiducial1+evolveStep,:)-Y(fiducial2+evolveStep,:)); while (1) if (fiducial1+evolveStep==YL)||(fiducial2+evolveStep==YL) lyapunov=lyapunov+log(L2/L1); fiducial1=fiducial1+evolveStep; break; elseif L2>MaxEL lyapunov=lyapunov+log(L2/L1); fiducial1=fiducial1+evolveStep; fiducial2=fiducial2+evolveStep; fiducial2=findreplace(Y,fiducial1,fiducial2,MaxFL); evolveStep=1; L1=norm(Y(fiducial1,:)-Y(fiducial2,:)); else evolveStep=evolveStep+1; end L2=norm(Y(fiducial1+evolveStep,:)-Y(fiducial2+evolveStep,:)); end fiducial1=fiducial1-1; lyapunov=lyapunov/(t*fiducial1); %--------------------------------------------------- function P=findreplace(Y,P1,P2,MaxL) % 寻找替代点 P=P2; YL=length(Y); HYL=fix(YL/2); k=1; % 优先选择距离较小的,与被替代的点处于以基点为原点的坐标下相同象限中 for i=1:YL if (sum(abs(Y(P1,:)-Y(i,:))) < MaxL ... && sign(Y(P2,1)-Y(P1,1))==sign(Y(i,1)-Y(P1,1))... && sign(Y(P2,2)-Y(P1,2))==sign(Y(i,2)-Y(P1,2))... && sign(Y(P2,3)-Y(P1,3))==sign(Y(i,3)-Y(P1,3))) temp(k)=i; k=k+1; end end % 如没有满足上述条件的点不进行替代,否则取尽可能保持方向不变的点 % 此处采用取外积模最小的点 if k>1 P=temp(1); PCross=norm(cross(Y(P1,:),Y(P,:))); for i=2:k-1 PCross1=norm(cross(Y(P1,:),Y(temp(i),:))); if PCross1 < PCross P=temp(i); PCross=PCross1; end end end

取延迟Chen系统，方程如下:

+----------------------------------------+
|   dx/dt=a*(y-x)                        |
|   dy/dt=(c-a)*x+c*y-x*z+k*(y(t-tau)-y) |
|   dz/dt=x*y-b*z                        |
+----------------------------------------+ 

其中 a = 35;b = 3;c = 28;
并且 k 和 tau 取以下值时都是非混沌

k = 2.4; tau = .2;
k = 2.8; tau = .15;
k = 4.6; tau = .3;
k = 0.8; tau = .5;
k = 2.65; tau = .16;
k = 4.2; tau = .1;
k = 3.29; tau = .13;
k = 2.3; tau = .18;

见dde_chen.m

如 k=0 即为不带延迟系统
 
>> history = [-3;-4;14];
tspan = [0,50];
a = 35;b = 3;c = 28;
k = 0;tau = .3;
opts = ddeset('MaxStep',1e-2);
sol = dde23(@dde_chenf,tau,history,tspan,opts,a,b,c,k);
Y=sol.y(1:3,2000:5000);
Y=Y';
plot3(Y(:,1),Y(:,2),Y(:,3));
l=wolf_3Dlyapunov(1e-2,Y)
 
l =
 
    2.0094

==========================================

取k = 2.4;tau = .2;
 
>> history = [-3;-4;14];
tspan = [0,50];
a = 35;b = 3;c = 28;
k = 2.4;tau = .2;
opts = ddeset('MaxStep',1e-2);
sol = dde23(@dde_chenf,tau,history,tspan,opts,a,b,c,k);
Y=sol.y(1:3,2000:5000);
Y=Y';
plot3(Y(:,1),Y(:,2),Y(:,3));
l=wolf_3Dlyapunov(1e-2,Y)
 
l =
 
   -0.0969

==========================================

取 k = 4.2;tau = .1; 
 
>> history = [-3;-4;14];
tspan = [0,50];
a = 35;b = 3;c = 28;
k = 4.2;tau = .1;
opts = ddeset('MaxStep',1e-2);
sol = dde23(@dde_chenf,tau,history,tspan,opts,a,b,c,k);
Y=sol.y(1:3,2000:5000);
Y=Y';
plot3(Y(:,1),Y(:,2),Y(:,3));
l=wolf_3Dlyapunov(1e-2,Y)
 
l =
 
   -0.2768

==========================================

取 k = .8;tau = .6; 是混沌系统

 
>> history = [-3;-4;14];
tspan = [0,50];
a = 35;b = 3;c = 28;
k = .8;tau = .6;
opts = ddeset('MaxStep',1e-2);
sol = dde23(@dde_chenf,tau,history,tspan,opts,a,b,c,k);
Y=sol.y(1:3,2000:5000);
Y=Y';
plot3(Y(:,1),Y(:,2),Y(:,3));
l=wolf_3Dlyapunov(1e-2,Y)
 
l =
 
    3.0002

参考程序
dde_chen.m
dde_rossler.m

参考文献：
[1]Alan WOLF,Jack B.SWIFT,Harry L.SWINNEY and John A.VASTANO DETERMINING LYAPUNOV EXPONENTS FROM A TIME SERIES Physica 16D(1985)285-317