Chaos Synchronization

lorenz_Synchronous.m

function lorenz_Synchronous % Lorenz 系统同步 % % Author:LDYU % Author's email: ustb03-07@yahoo.com.cn % [T,Y]=ode45(@lorenzf1,[0,3],[1;-10;1;10;2]); subplot(2,1,1) plot(T,Y(:,3),'k') hold on plot(T,Y(:,5),'b') [T,Y]=ode45(@lorenzf2,[0,3],[10;10;10;-20;2]); subplot(2,1,2) plot(T,Y(:,3),'k') hold on plot(T,Y(:,5),'b') %------------------------------------- function dx=lorenzf1(t,x) % dx= a*(y-x) % dy= -x*z+r*x-y % dz= x*y-b*z % dy1= -x*z1+r*x-y1 % dz1= x*y1-b*z1 a=16;b=4;r=45.92; dx=[ a*(x(2)-x(1)) -x(1)*x(3)+r*x(1)-x(2) x(1)*x(2)-b*x(3) -x(1)*x(5)+r*x(1)-x(4) x(1)*x(4)-b*x(5)]; %------------------------------------- function dx=lorenzf2(t,x) % dx= a*(y-x) % dy=-x*z+r*x-y % dz= x*y-b*z % dx1=a*(y-x1) % dz1=x1*y-b*z1 a=16;b=4;r=45.92; dx=[ a*(x(2)-x(1)) -x(1)*x(3)+r*x(1)-x(2) x(1)*x(2)-b*x(3) a*(x(2)-x(4)) x(4)*x(2)-b*x(5)];

用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

Chen Solution with delay

function dde_chen % 延迟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 % % Author's email: ustb03-07@yahoo.com.cn % history = [-3;-4;14]; tspan = [0,25]; opts = ddeset('RelTol',1e-5,'AbsTol',1e-8); a = 35;b = 3;c = 28; k = 2.4;tau = .2; % k = 2.8;tau = .15; % k = 4.6;tau = .3; % k = .8;tau = .5; % k = 2.65;tau = .16; % k = 4.2;tau = .1; % k = 3.29;tau = .13; % k = 2.3;tau = .18; sol = dde23(@dde_chenf,tau,history,tspan,opts,a,b,c,k); plot3(sol.y(1,4000:end),sol.y(2,4000:end),sol.y(3,4000:end)) title('Chen Solution with delay.') xlabel('x(t)') ylabel('y(t)') zlabel('z(t)') %-------------------------- function dydt = dde_chenf(t,y,Z,a,b,c,k) % Differential equations function for Chen. ylag = Z(:,1); dydt = [ a*(y(2)-y(1)) (c-a-y(3))*y(1)+c*y(2)+k*(ylag(2)-y(2)) y(1)*y(2)-b*y(3) ];

Rossler Solution with delay

function dde_rossler % 延迟Rossler系统(1) % % 方程如下: % dx/dt=-y-z+k*(x(t-tau)-x) % dy/dt=x+a*y % dz/dt=b+z*(x-c) % % % Author's email: ustb03-07@yahoo.com.cn % history = [-4;4;0]; tspan = [0,250]; opts = ddeset('RelTol',1e-5,'AbsTol',1e-8); a = 0.2; b = 0.2; c = 4.5; % 一周期 k = 0.2; % 二周期 % k = 0.08; % Solve the DDEs that arise when there is a delay of tau. tau = 5.691; sol = dde23(@dde_rosslerf,tau,history,tspan,opts,a,b,c,k); plot3(sol.y(1,4000:end),sol.y(2,4000:end),sol.y(3,4000:end)) title('Rossler Solution with delay.') xlabel('x(t)') ylabel('y(t)') zlabel('z(t)') %-------------------------- function dydt = dde_rosslerf(t,y,Z,a,b,c,k) % Differential equations function for Rossler. ylag = Z(:,1); dydt = [ -y(2)-y(3)+k*(ylag(1)-y(1)) y(1)+a*y(2) b+y(1)*y(3)-c*y(3) ];

function dde_rossler2 % 延迟Rossler系统(2) % % 方程如下: % dx/dt=-y-z % dy/dt=x+a*y+k*(x(t-tau)-x) % dz/dt=b+z*(x-c) % % % Author's email: ustb03-07@yahoo.com.cn % history = [-4;4;0]; tspan = [0,200]; opts = ddeset('RelTol',1e-5,'AbsTol',1e-8); a = 0.2; b = 0.2; c = 5.7; % Solve the DDEs that arise when there is a delay of tau. k = 0.025; tau = 17.4; sol = dde23(@dde_rosslerf,tau,history,tspan,opts,a,b,c,k); plot3(sol.y(1,3000:end),sol.y(2,3000:end),sol.y(3,3000:end)) title('Rossler Solution with delay.') xlabel('x(t)') ylabel('y(t)') zlabel('z(t)') %-------------------------- function dydt = dde_rosslerf(t,y,Z,a,b,c,k) % Differential equations function for Rossler. ylag = Z(:,1); dydt = [ -y(2)-y(3) y(1)+a*y(2)+k*(ylag(1)-y(1)) b+y(1)*y(3)-c*y(3) ];

Josephon系统(Poincare图象和Lyapunov指数)

Josephon.m

function dx=Josephon(t,x); % Josephon方程 % dx=Josephon(t,[x;y;b]) % t-时间，x,y-为自变量，b-为如下方程所示的参数 % eg: dx=Josephon(10,[0;0;.33803]) % % 方程如下: % θ''+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(Poincare图象): % [T,Y]=ode45('Josephon',[0,1000],[0;0;.33803]); % plot(mod(Y(:,1),6.283),Y(:,2),'.','markersize',2); % % Author's email: ustb03-07@yahoo.com.cn % G=.7;A=.4;w=.25; a=0.0125;I=.905; % b=.33803; b=x(3); dx(1,1)=x(2); dx(2,1)=-G*x(2)-sin(x(1))+I+A*sin(w*t)+a*sin(b*w*t); dx(3,1)=0;

jose_ly.m

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 % Author's email: ustb03-07@yahoo.com.cn % 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('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); end ly=lsum/k;

>> help jose_ly
 
  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
  Author's email: ustb03-07@yahoo.com.cn
  
 
>> ly=jose_ly(0,800)
 
ly =
 
    0.0522
 

josel.m

% the largest lyapunov exponents of josephson % % Author's email: ustb03-07@yahoo.com.cn % clear Z=[]; for b=linspace(0,1.5,50) Z=[Z b+i*jose_ly(b,300)]; end plot(Z,'-') title('Josephon最大Lyapunov指数图'),xlabel('b'),ylabel('lyapunov') grid on

>> josel

Duffing系统

Duffing.m

function dx=Duffing(t,x); % Duffing方程 % dx=Duffing(t,[x;y]) % t-时间，x,y-为自变量 % eg: dx=Duffing(10,[0;0;39]) % % 方程如下: % dx=y % dy=-r*y-x^3+b*cos(t) % r=0.3 b=39 % % Example(Duffing图象): % [T,Y]=ode45('Duffing',50,[0;0;39]); % plot(Y(:,1),Y(:,2)); % % Example(分岔图): % % Z=[]; % for p=linspace(6,13,280); % [T,Y]=ode45('Duffing',10,[3;0;p]); % [T,Y]=ode45('Duffing',100,Y(end,:)); % for k=2:length(Y) % f=k-1; % y=1; % if Y(k,1)<0 % if Y(f,1)>0 % y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1)); % end % else % if Y(f,1)<0 % y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1)); % end % end % if y<0 % Z=[Z p+y*i]; % end % end % end % plot(Z,'.','markersize',2) % title('Duffing映射分岔图(x=0,dx<0)'),xlabel('b'),ylabel('y') % % Author's email: ustb03-07@yahoo.com.cn % r=0.3; b=x(3); dx(1,1)=x(2); dx(2,1)=-r*x(2)-x(1)^3+b*cos(t); dx(3,1)=0;

基于奇异值分解的Lyapunov指数计算（3）

Rossler_euler.m

% Rossler图形(欧拉方法) % % Author：yujunjie % Author's email: ustb03-07@yahoo.com.cn % clear h=0.006;a=0.15;b=0.2;c=10; x=0;y=-12;z=0; Y=[]; for i=1:10000 x1=x+h*(-y-z); y1=y+h*(x+a*y); z1=z+h*(b+x*z-c*z); x=x1;y=y1;z=z1; Y(i,:)=[x y z]; end plot3(Y(:,1),Y(:,2),Y(:,3));

rossler_le_eu.m

% 奇异值分解求Lyapunov法 % 微分rossler系统 % % Author：yujunjie % Author's email: ustb03-07@yahoo.com.cn % h=0.005;a=0.15;b=0.2;c=10; x=0;y=-12;z=0; V=eye(3); S=V;b1=0; k=10000; for i=1:k x1=x+h*(-y-z); y1=y+h*(x+a*y); z1=z+h*(b+x*z-c*z); x=x1;y=y1;z=z1; J=[0 -1 -1 1 a 0 z 0 x-c]; J=eye(3)+h*J; B=J*V*S; [V,S,U]=svd(B); a_max=max(diag(S)); S=(1/a_max)*S; b1=b1+log(a_max); end Lyapunov=(log(diag(S))+b1)/(k*h)

>> rossler_le_eu
 
Lyapunov =
 
    0.0972
   -0.0152
   -9.8887

基于奇异值分解的Lyapunov指数计算（2）

Chen_euler.m

% Chen图形(欧拉方法) % % Author：yujunjie % Author's email: ustb03-07@yahoo.com.cn % clear h=0.004;a=35;b=3;c=28; x=10;y=10;z=20; Y=[]; for i=1:5000 x1=x+h*a*(y-x); y1=y+h*((c-a)*x-x*z+c*y); z1=z+h*(x*y-b*z); x=x1;y=y1;z=z1; Y(i,:)=[x y z]; end plot3(Y(:,1),Y(:,2),Y(:,3));

chen_le_eu.m

% 奇异值分解求Lyapunov法 % 微分chen系统 % % Author：yujunjie % Author's email: ustb03-07@yahoo.com.cn % h=0.002;a=35;b=3;c=28; x=10;y=10;z=20; V=eye(3); S=V;b1=0; k=6000; for i=1:k x1=x+h*a*(y-x); y1=y+h*((c-a)*x-x*z+c*y); z1=z+h*(x*y-b*z); x=x1;y=y1;z=z1; J=[-a a 0 c-a-z c -x y x -b]; J=eye(3)+h*J; B=J*V*S; [V,S,U]=svd(B); a_max=max(diag(S)); S=(1/a_max)*S; b1=b1+log(a_max); end Lyapunov=(log(diag(S))+b1)/(k*h)

>> chen_le_eu
 
Lyapunov =
 
    2.1531
    0.0051
  -11.8983