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