Cat映射的周期性

function cattest2 % Cat 加密图象 clear all f=imread('cat.jpg'); % 读曲图象 A=double(f)/255; [m,n,l]=size(A); Cat=[321 40;8 1];N=124; m=m-rem(m,N); n=n-rem(n,N); A=A(1:m,1:n,:); subplot(2,3,1) imagesc(A,[0 1]); mk=m/N;nk=n/N; for k=1:5 % 将图象分割成124*124的数据块 for mki=0:mk-1 for nki=0:nk-1 % 处理N*N的数据块 for i=N*mki:N*mki+123 for j=N*nki:N*nki+123 c=[i j]'; c=Catmap(c,Cat,N); B(mki*N+c(1)+1,nki*N+c(2)+1,:)=A(i+1,j+1,:); end end end end subplot(2,3,k+1) imagesc(B,[0 1]); A=B; end

相关程序 Cat.m

Cat 映射

Cat.m

function newx=Catmap(x,A,N) % Cat映射 % x(n+1)=A*x(n) (mod N) % x是m*1向量，A是m*m矩阵,N是整数 % % Example: % global Cat; % global N; % A=[1 1;1 2];N=1; % newx=Catmap([.4 .5]') % % Author:LDYU % Author's email: ustb03-07@yahoo.com.cn % mewx=A*x; newx=mod(mewx,N);

混沌 Cat 映射

function newx=Catmap(x) % Cat映射 % x(n+1)=x(n)+y(n) (mod 1) % y(n+1)=x(n)+2*y(n) (mod 1) % % newx(1)=mod(x(1)+x(2),1); newx(2)=mod(x(1)+2*x(2),1);

function test clear all % 读曲图象（路径按自己情况改） f=imread('cat.jpg'); A=double(f)/255; [m,n,l]=size(A); subplot(3,2,1) imagesc(A,[0 1]); for k=1:5 S=zeros(m,n); for i=1:m for j=1:n c=[(i-1)/m (j-1)/n]; c=Catmap(c); ib=fix(m*c(1))+1; jb=fix(n*c(2))+1; if S(ib,jb)==0 B(ib,jb,:)=A(i,j,:); S(ib,jb)=1; % 如有变换后的点重叠，移到邻近单元，减少信息损失 elseif jb<n & S(ib,jb+1)==0 B(ib,jb+1,:)=A(i,j,:); S(ib,jb+1)=1; elseif ib<m & S(ib+1,jb)==0 B(ib+1,jb,:)=A(i,j,:); S(ib+1,jb)=1; end end end subplot(3,2,k+1) imagesc(B,[0 1]); A=B; end

变换的原图

庞加莱截面和倍周期分岔

以前的关于微分方程组的庞加莱截面和倍周期分岔程序有误，现在更正如下：

Z=[]; [T,Y]=ode45('Lorenz',[0,1000],[1;1;1;16;4;49.52]); Y(:,3)=Y(:,3)-30; for k=1:(length(Y)-1) if Y(k,3)<0 if Y(k+1,3)>0 % 用直线连接位于截面两侧的相邻离散点，近似得到穿过截面的点 x=Y(k,1)-Y(k,3)*(Y(k+1,1)-Y(k,1))/(Y(k+1,3)-Y(k,3)); y=Y(k,2)-Y(k,3)*(Y(k+1,2)-Y(k,2))/(Y(k+1,3)-Y(k,3)); Z=[Z x+i*y]; end % 如果只画从一侧穿过截面的图象，以下的if-end语句可删去 else if Y(k+1,3)<0 x=Y(k,1)-Y(k,3)*(Y(k+1,1)-Y(k,1))/(Y(k+1,3)-Y(k,3)); y=Y(k,2)-Y(k,3)*(Y(k+1,2)-Y(k,2))/(Y(k+1,3)-Y(k,3)); Z=[Z x+i*y]; end end end plot(Z,'.','markersize',2) title('Lorenz 系统的 Poincare 映像 z=30') xlabel('x'),ylabel('y')

Z=[]; [T,Y]=ode45('Chen',[0,1000],[1;1;1;35;3;28]); for k=1:(length(Y)-1) if Y(k,1)<0 if Y(k+1,1)>0 y=Y(k,2)-Y(k,1)*(Y(k+1,2)-Y(k,2))/(Y(k+1,1)-Y(k,1)); z=Y(k,3)-Y(k,1)*(Y(k+1,3)-Y(k,3))/(Y(k+1,1)-Y(k,1)); Z=[Z y+i*z]; end else if Y(k+1,1)<0 y=Y(k,2)-Y(k,1)*(Y(k+1,2)-Y(k,2))/(Y(k+1,1)-Y(k,1)); z=Y(k,3)-Y(k,1)*(Y(k+1,3)-Y(k,3))/(Y(k+1,1)-Y(k,1)); Z=[Z y+i*z]; end end end plot(Z,'.','markersize',2) title('Chen 系统的 Poincare 映像 x=0') xlabel('y'),ylabel('z')

Z=[]; [T,Y]=ode45('Lorenz3',[0,5000],[0.5;0.5;2;.9;.1;2;0]); Y(:,2)=Y(:,2)-.2; for k=1:(length(Y)-1) if Y(k,2)<0 if Y(k+1,2)>0 x=Y(k,1)-Y(k,2)*(Y(k+1,1)-Y(k,1))/(Y(k+1,2)-Y(k,2)); z=Y(k,3)-Y(k,2)*(Y(k+1,3)-Y(k,3))/(Y(k+1,2)-Y(k,2)); Z=[Z x+i*z]; end else if Y(k+1,2)<0 x=Y(k,1)-Y(k,2)*(Y(k+1,1)-Y(k,1))/(Y(k+1,2)-Y(k,2)); z=Y(k,3)-Y(k,2)*(Y(k+1,3)-Y(k,3))/(Y(k+1,2)-Y(k,2)); Z=[Z x+i*z]; end end end plot(Z,'.','markersize',2) title('模拟 Lorenz 系统的 Poincare 映像 y=0.2') xlabel('x'),ylabel('z')

Z=[]; for r=linspace(1,500,250); % 舍弃前面迭带的结果，用后面的结果画图 [T,Y]=ode45('Lorenz',1,[1;1;1;10;8/3;r]); [T,Y]=ode45('Lorenz',50,Y(length(Y),:)); Y(:,1)=Y(:,2)-Y(:,1); for k=2:length(Y) f=k-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)); Z=[Z r+abs(y)*i]; 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)); Z=[Z r+abs(y)*i]; end end end end plot(Z,'.','markersize',1) title('Lorenz映射分岔图') xlabel('r'),ylabel('|y| where x=y')

Chen 系统的 Poincare 映像

Z=[]; [T,Y]=ode45('Chen',[0,5000],[1;1;1;35;3;28]); for k=1:length(Y) if abs(Y(k,1))<1e-2 Z=[Z Y(k,2)+i*Y(k,3)]; end end plot(Z,'.','markersize',2) title('Chen 系统的 Poincare 映像 x=0') xlabel('y'),ylabel('z')

Z=[]; [T,Y]=ode45('Chen',[0,5000],[1;1;1;35;3;28]); for k=1:length(Y) if abs(Y(k,3)-28)<1e-2 Z=[Z Y(k,1)+i*Y(k,2)]; end end plot(Z,'.','markersize',2) title('Chen 系统的 Poincare 映像 z=c') xlabel('x'),ylabel('y')

部分放大



相关程序 Chen.m

模拟 Lorenz 系统的 Poincare 映像

Z=[]; [T,Y]=ode45('Lorenz3',[0,6000],[0.5;0.5;2;.9;.1;2;0]); for k=1:length(Y) if (abs(Y(k,2)-0.2)<1e-2) Z=[Z Y(k,1)+i*Y(k,3)]; end end plot(Z,'.','markersize',2) title('模拟 Lorenz 系统的 Poincare 映像 y=0.2') xlabel('x'),ylabel('z')

相关程序 Lorenz3.m

第一阶段个人工作小结

第一阶段个人工作小结 第一阶段的任务是查资料及上机实现Lyapunov指数的算法,在这些日子里我查阅 了大量的资料并且上机实现了一些求Lyapunov指数的算法。在查资料的过程中虽然这 方面的东西有很多,但却没有什么现成的程序,网上好不容易找到一些却不是我需要的, 有时找到可以用的程序,运行结果却和理论值有较大的偏差。书上和论文上几乎没给什 么程序，有伪代码的也是十分简略,增加了实现的难度。不过还是发现一些有价值的程 序，这些不全是用Matlab写的,有些还是用Basic写的，这些程序给了我很大的启发。 以前按算法写的程序总得不到正确的结果,现在回头看来很大程度是我理解的意思和书 上要表达的意思不一致，如果书上能提供这些算法的程序，无论是用什么语言实现,也 可以帮助读者理解，而不至于让人走弯路。 虽然遇到了些困难，但实现的程序结果和书上,论文上的结果比较吻合。总之是完 成了第一阶段的任务了，另外今天是中秋加十一放假期间，不熄灯，庆祝一下:)

Lyapunov exponents of Chen

Chen2.m

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*z global 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)]; end a=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

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 i>50 lsum=lsum+log(d1/d0); end end Z=[Z lsum/(i-50)]; end a=linspace(32,40,80); plot(a,Z,'-'); title('Chen 系统最大lyapunov指数') xlabel('parameter a'),ylabel('lyapunov exponents')

相连的程序有Chen.m

模拟 Lorenz 系统最大lyapunov指数谱

Z=[]; d0=1e-8; for a=linspace(0.3,2,70) lsum=0; x=1;y=1;z=1; x1=1;y1=1;z1=1+d0; for i=1:100 [T1,Y1]=ode45('Lorenz3',1,[x;y;z;a;.1;2;0]); [T2,Y2]=ode45('Lorenz3',1,[x1;y1;z1;a;.1;2;0]); 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 i>50 lsum=lsum+log(d1/d0); end end Z=[Z lsum/(i-50)]; end a=linspace(0.3,2,70); plot(a,Z,'-'); title('模拟 Lorenz 系统最大lyapunov指数谱') xlabel('parameter a') ylabel('lyapunov exponents')

相连的程序有Lorenz3.m

模拟洛仑滋系统

Lorenz3.m

function dx=Lorenz3(t,x); % 模拟洛仑滋系统[不显含时间t的自治系统] % dx=Lorenz(t,[x;y;z;a;b;c;m]) % t-可以取任何数字，x,y,z-为自变量，a,b,r-为如下方程所示的参数 % eg: dx=Lorenz(0,[1;1;1;.9;.1;2;0]) % % 方程如下: % dx/dt=a*(y-x) % dy/dt=sign(x)*(c-z)+m % dz/dt=|x|-b*z % % % Author：yujunjie % Author's email: ustb03-07@yahoo.com.cn % dx(1,1)=x(4)*(x(2)-x(1)); dx(2,1)=sign(x(1))*(x(6)-x(3))+x(7); dx(3,1)=abs(x(1))-x(5)*x(3); dx(4,1)=0; dx(5,1)=0; dx(6,1)=0; dx(7,1)=0;

>>    [T,Y]=ode45('Lorenz3',[0,10],[1;1;1;.9;.1;2;0]);
      x=Y(length(Y),:);
      [T,Y]=ode45('Lorenz3',[0,500],x);
      plot3(Y(:,1),Y(:,2),Y(:,3));
 

>>    [T,Y]=ode45('Lorenz3',[0,10],[1;1;1;.9;.1;2;1]);
      x=Y(length(Y),:);
      [T,Y]=ode45('Lorenz3',[0,500],x);
      plot3(Y(:,1),Y(:,2),Y(:,3));
 

>>    [T,Y]=ode45('Lorenz3',[0,10],[1;1;1;.9;.1;2;-1]);
      x=Y(length(Y),:);
      [T,Y]=ode45('Lorenz3',[0,500],x);
      plot3(Y(:,1),Y(:,2),Y(:,3));
 

系统对初始条件敏感依赖性

[T1,Y1]=ode45('Lorenz',20,[1;1;1;16;4;45.92]); plot(T1,Y1(:,1),'-'); hold on [T2,Y2]=ode45('Lorenz',20,[1+.001;1;1;16;4;45.92]); plot(T2,Y2(:,1),'r--');

相连的程序有Lorenz.m

Lorenz指数谱

Lorenlp.m

% Lyapunov exponents of Lorenz % % Author：yujunjie % Author's email: ustb03-07@yahoo.com.cn % function lorenzlp Z1=[]; Z2=[]; Z3=[]; global a; global b; global r; a=10;b=8/3; for r=linspace(0,500,80); y=[1;1;1;1;0;0;0;1;0;0;0;1]; lp=0; for k=1:100 [T,Y] = ode45('Lorenz2', 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/100; Z1=[Z1 lp(1)]; Z2=[Z2 lp(2)]; Z3=[Z3 lp(3)]; end r=linspace(0,500,80); plot(r,Z1,'-',r,Z2,'-',r,Z3,'-'); title('Lyapunov exponents of Lorenz') xlabel('parameter r'),ylabel('lyapunov exponents') grid on

>> lorenzlp
 

相连的程序有GS.m and Lorenz2.m

计算Lorenz系统的Lyapunov指数

Lorenz2.m

function dX = Lorenz2(t,X) % % Lorenz吸引子，用来计算Lyapunov指数 % dx/dt = a(y-x), % dy/dt = x(r-z)-y, % dz/dt = xy-bz, global a; % 变量不放入参数表中 global b; global r; 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) = x*(r-z)-y; dX(3) = x*y-b*z; % Lorenz吸引子的Jacobi矩阵 Jaco = [-a a 0; r-z -1 -x; y x -b]; dX(4:12) = Jaco*Y;

llp.m

function llp global a; global b; global r; a=16;b=4;r=45.92; y=[1;1;1;1;0;0;0;1;0;0;0;1]; lp=0; for k=1:300 [T,Y] = ode45('Lorenz2', 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/300

>> llp
 
lp =
 
    1.5096   -0.5995  -21.8105
 

相连的程序有GS.m

Lyapunov exponents of Henon

要用的相关程序 Henon.m 和 GS.m

% Lyapunov exponents of Henon % % Author：yujunjie % Author's email: ustb03-07@yahoo.com.cn % function henonlp Zmax=[]; Zmin=[]; for p=linspace(0,1.4,300); x=[1;0;p;.3]; e=[1 0 0 1]; lp=0; for k=1:500 x=Henon(x); J=[-2*p*x(1) 1 .3 0]; e=GS(J*e); mod(1)=norm(e(:,1)); mod(2)=norm(e(:,2)); lp = lp+log(abs(mod)); e(:,1)=e(:,1)/mod(1); e(:,2)=e(:,2)/mod(2); end lp=lp/500; Zmax=[Zmax lp(1)]; Zmin=[Zmin lp(2)]; end p=linspace(0,1.4,300); plot(p,Zmax,'-',p,Zmin,'-'); title('Lyapunov exponents of Henon') xlabel('parameter p'),ylabel('lyapunov exponents') grid on

>> henonlp
 

GS 正交化算法

GS.m

% G-S正交化 % A = GS(V) % V 为n*n向量 % A 为n*n正交向量 % % Example: % A=GS([1 0 1;2 2 0;3 1 1]) % % Author：yujunjie % Author's email: ustb03-07@yahoo.com.cn % function A = GS(V) n=length(V); for i=1:n k=0; for j=1:i-1 k=k+(A(:,j)'*V(:,i)/(A(:,j)'*A(:,j)))*A(:,j); end A(:,i)=V(:,i)-k; end

>> A=GS([1 0 1;2 2 0;3 1 1])
 
A =
 
    1.0000   -0.5000    0.3810
    2.0000    1.0000    0.0952
    3.0000   -0.5000   -0.1905
 

Logistic系统

一个网上的程序

% This file draws the bifurcation diagram and Lyapunov exponent % for the Logistic Iterated Map x_n+1 = r*x_n.*(1-x_n); 0 < x < 1. % The parameter r is in interval [0, 4] % adapted from a program by G Lindblad gli@theophys.kth.se % http://www.theophys.kth.se/5a1352/mat.html rect = [200 80 700 650]; set(0, 'defaultfigureposition',rect); n1=200; %% no of lattice points in coordinate and parameter r n2=500; %% no of iterations to reach attractor n3=250; %% no of iterations for bifurcation diagram n4=4000; %% no of iterations for Lyapunov exponent k1=[]; kk=[]; q1=[]; home, if isempty(kk) %> setting default kk=[3 4]; end q3=input('> Choose a r-interval [a b] or >>> '); q=isempty(q3); if q==0 disp('> The r-interval is set = '); disp(q3); kk = q3; end disp('> The calculations could take a few seconds .... '); k0=linspace(kk(1), kk(2) ,n1)'; seed=rand(1); x1=seed; ww=[]; for iter1=1:n2 %% this is just to reach the attractor x1=k0.*x1.*(1-x1); end for it3=1:n3 %% points on the attractor, we hope x1=k0.*x1.*(1-x1); ww=[ww, x1]; end % We pick a lattice of n1 values of the parameter k1 k0=linspace(kk(1), kk(2) ,n1)'; % We start the iteration at a random point seed=rand(1);x1=seed; % The first part of the iteration is discarded - % this is just to reach the attractor for it2=1:n1 % x1=k0.*x1.*(1-x1); end % Now we generate a sequence of approximations to the Lyapunov exponent % as a function of the vector k0! aa=log(abs(1-2*x1)); for it3=1:n4 %% x1=k0.*x1.*(1-x1); y1=log(abs(1-2*x1)); aa=(y1+it3*aa)/(1+it3); end aa = aa + log(k0+eps); % adding a term coming from the k0 factor subplot(211) plot(k0,ww, '. k','MarkerSize',4); axis([kk(1) kk(2) 0 1]); %axis tight; title('A bifurcation diagram of the logistic equation'); subplot(212) plot(k0,[aa, zeros(size(k0))]); axis([kk(1) kk(2) -1 1]), title('Lyapunov exponent of the logistic equation');

微分混沌方程组最大Lyapunov指数算法

求Lorenz系统
在参数a=16，b=4，r=46.92下的最大Lyapunov指数

d0=1e-8; le=0; lsum=0; x=1;y=1;z=1; x1=1;y1=1;z1=1+d0; for i=1:500 [T1,Y1]=ode45('Lorenz',[0,1],[x;y;z;16;4;46.92]); [T2,Y2]=ode45('Lorenz',[0,1],[x1;y1;z1;16;4;46.92]); 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 i>100 lsum=lsum+log(d1/d0); end end le=lsum/(i-100)

le =
 
    1.4954

Chen系统和Rossler系统也可以用以上的算法只不过将Lorenz换成Chen或Rossler即可。

Lorenz系统

Lorenz.m

function dx=Lorenz(t,x); % 洛仑滋方程[不显含时间t的自治系统] % dx=Lorenz(t,[x;y;z;a;b;r]) % t-可以取任何数字，x,y,z-为自变量，a,b,r-为如下方程所示的参数 % eg: dx=Lorenz(0,[1;1;1;16;4;49]) % % 方程如下: % dx/dt=a*(y-x) % dy/dt=x*(r-z)-y % dz/dt=x*y-b*z % % Example(函数图象): % [T,Y]=ode45('Lorenz',[0,50],[1;1;1;16;4;49]); % plot3(Y(:,1),Y(:,2),Y(:,3)); % % Example(分岔图): % % Z=[]; % for r=linspace(1,500,350); % [T,Y]=ode45('Lorenz',[0,50],[1;1;1;10;8/3;r]); % n=length(Y); % for k=round(n/2):n % if abs(Y(k,1))<1 % Z=[Z,r+abs(Y(k,2))*i]; % end % end % end % plot(Z,'.','markersize',1) % title('Lorenz映射分岔图') % xlabel('r'),ylabel('|y| where x=0') % % Author：yujunjie % Author's email: ustb03-07@yahoo.com.cn % dx(1,1)=x(4)*(x(2)-x(1)); dx(2,1)=x(1)*(x(6)-x(3))-x(2); dx(3,1)=x(1)*x(2)-x(5)*x(3); dx(4,1)=0; dx(5,1)=0; dx(6,1)=0;

Chen and Rossler

Chen.m

function dx=Chen(t,x); % Chen方程[不显含时间t的自治系统] % dx=Chen(t,[x;y;z;a;b;c]) % t-可以取任何数字，x,y,z-为自变量，a,b,c-为如下方程所示的参数 % eg: dx=Chen(0,[1;1;1;35;3;28]) % % 方程如下: % dx/dt=a*(y-x) % dy/dt=(c-a)*x+c*y-x*z % dz/dt=x*y-b*z % % Example(函数图象): % [T,Y]=ode45('Chen',[0,50],[1;1;1;35;3;28]); % plot3(Y(:,1),Y(:,2),Y(:,3)); % % Author's email: ustb03-07@yahoo.com.cn % dx(1,1)=x(4)*(x(2)-x(1)); dx(2,1)=(x(6)-x(4))*x(1)+x(6)*x(2)-x(1)*x(3); dx(3,1)=x(1)*x(2)-x(5)*x(3); dx(4,1)=0; dx(5,1)=0; dx(6,1)=0;

Rossler.m

function dx=Rossler(t,x); % Rossler方程[不显含时间t的自治系统] % dx=Rossler(t,[x;y;z;a;b;c]) % t-可以取任何数字，x,y,z-为自变量，a,b,c-为如下方程所示的参数 % eg: dx=Rossler(0,[1;1;1;.15;.2;10]) % % 方程如下: % dx/dt=-y-z % dy/dt=x+a*y % dz/dt=b+z*(x-c) % % Example(函数图象): % [T,Y]=ode45('Rossler',[0,400],[1;1;1;.15;.2;10]); % plot3(Y(:,1),Y(:,2),Y(:,3)); % % Author's email: ustb03-07@yahoo.com.cn % dx(1,1)=-x(2)-x(3); dx(2,1)=x(1)+x(4)*x(2); dx(3,1)=x(5)+x(3)*(x(1)-x(6)); dx(4,1)=0; dx(5,1)=0; dx(6,1)=0;

Henon系统 （2）

Henon.m

function newx=Henon(x); % Henon方程[差分方程] % newx=Henon([x;y;p;q]) % % 方程如下: % x(k+1)=-p*x(k)*x(k)+y(k)+1 % y(k+1)=q*x(k) % % Example(函数图象): % % X=[1;0;1.4;.3];Y=[]; % for i=1:1000 % X=feval(@Henon,X); % Y(i,:)=X(1:2,1); % end % plot(Y(:,1),Y(:,2),'.'); % title('Henon映射图'),xlabel('x'),ylabel('y') % % Example(分岔图): % % Z=[]; % for p=linspace(0,1.4,300); % x=[1;0;p;.3]; % for k=1:300; % x=Henon(x); % if k>60 % Z=[Z,p+x(1)*i]; % end % end % end % plot(Z,'.','markersize',1) % title('Henon映射分岔图'),xlabel('p'),ylabel('x') % % Example(最大Lyapunov指数谱图): % % d0=1e-8; % Z=[]; % for p=linspace(0,1.4,300) % le=0; % lsum=0; % x=[1;0;p;.3]; % x1=[1;d0;p;.3]; % for k=1:800 % x=Henon(x); % x1=Henon(x1); % d1=sqrt((x(1)-x1(1))^2+(x(2)-x1(2))^2); % x1=x+(d0/d1)*(x1-x); % if k>100 % lsum=lsum+log(d1/d0); % end % end % le=lsum/(k-100); % Z=[Z,p+le*i]; % end % plot(Z,'-') % title('Henon最大Lyapunov指数谱图'),xlabel('p'),ylabel('lyapunov') % grid on % % Author：yujunjie % Author's email: ustb03-07@yahoo.com.cn newx(1,1)=-x(3)*x(1)*x(1)+x(2)+1; newx(2,1)=x(4)*x(1); newx(3,1)=x(3); newx(4,1)=x(4);

Lozi系统

Lozi.m

function newx=Lozi(x); % Lozi方程[差分方程] % newx=Lozi([x;y;p;q]); % 方程如下: % x(k+1)=-p*x(k)+y(k)+1 % y(k+1)=q*x(k) % % Example(函数图象): % % X=[1;0;1.7;.5];Y=[]; % for i=1:1000 % X=feval(@Lozi,X); % Y(i,:)=X(1:2,1); % end % plot(Y(:,1),Y(:,2),'.'); % title('Lozi映射图'),xlabel('x'),ylabel('y') % % Example(分岔图): % % Z=[]; % for p=linspace(0,1.7,300); % x=[1;0;p;.5]; % for k=1:300; % x=Lozi(x); % if k>60 % Z=[Z,p+x(1)*i]; % end % end % end % plot(Z,'.','markersize',1) % title('Lozi映射分岔图'),xlabel('p'),ylabel('x') % % Example(最大Lyapunov指数谱图): % % d0=1e-8; % Z=[]; % for p=linspace(0,1.7,300) % le=0; % lsum=0; % x=[1;0;p;.5]; % x1=[1;d0;p;.5]; % for k=1:800 % x=Lozi(x); % x1=Lozi(x1); % d1=sqrt((x(1)-x1(1))^2+(x(2)-x1(2))^2); % x1=x+(d0/d1)*(x1-x); % if k>100 % lsum=lsum+log(d1/d0); % end % end % le=lsum/(k-100); % Z=[Z,p+le*i]; % end % plot(Z,'-') % title('Lozi最大Lyapunov指数图'),xlabel('p'),ylabel('lyapunov') % grid on % % Author：yujunjie % Author's email: ustb03-07@yahoo.com.cn newx(1,1)=-x(3)*abs(x(1))+x(2)+1; newx(2,1)=x(4)*x(1); newx(3,1)=x(3); newx(4,1)=x(4);