学习混沌: October 2006

Sunday, October 22, 2006

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 映射

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

变换的原图

Saturday, October 21, 2006

庞加莱截面和倍周期分岔

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


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')

Saturday, October 14, 2006

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')

部分放大

模拟 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')

Friday, October 06, 2006

第一阶段个人工作小结


    第一阶段的任务是查资料及上机实现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

Thursday, October 05, 2006

模拟 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

Wednesday, October 04, 2006

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

Tuesday, October 03, 2006

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');

Monday, October 02, 2006

微分混沌方程组最大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);

博客归档