洛伦兹系统ODE方程-MATLAB

1. 洛伦兹系统ODE方程

这是一个用于计算简单洛伦兹系统ODE方程的MATLAB函数。这个函数的输入参数包括一个时间向量t和一个状态向量y。y是一个3维向量，代表简单洛伦兹系统的三个状态变量：x、y和z。输出参数f也是一个3维向量，分别对应x、y和z的导数值。

函数首先使用全局变量c，来定义简单洛伦兹系统的一个参数值。
接着，函数将f向量的每个元素初始化为0。
然后，函数分别计算x、y和z的导数，将其分别保存到f向量的对应位置：

f(1)代表x的导数，其计算公式为：10*(y(2)-y(1))
f(2)代表y的导数，其计算公式为：(24-4*c)y(1)+cy(2)-y(1)*y(3)
f(3)代表z的导数，其计算公式为：y(1)y(2)-8/3y(3)
最后，函数返回计算得到的f向量。
请注意，这是一个函数，需要在其他代码中进行调用，并且需要提前声明全局变量c的值。

function f=chao_SimpleLorenz(t,y)
global c
f=zeros(3,1);
f(1)=10*(y(2)-y(1));
f(2)=(24-4*c)*y(1)+c*y(2)-y(1)*y(3);
f(3)=y(1)*y(2)-8/3*y(3);

2. 计算洛伦兹系统ODE方程

这是一个用于计算洛伦兹系统ODE方程的MATLAB函数。函数的输入参数包括一个时间向量t和一个状态向量X，其中X是一个12维向量，代表了洛伦兹系统的状态变量x、y、z及其对应的雅可比矩阵Y。
在函数体内，首先使用了全局变量c，来定义洛伦兹系统的参数值。
接下来，函数从输入的状态向量X中提取出x、y、z以及雅可比矩阵Y的值。
然后，函数初始化了一个12维的零向量f。
接着，函数计算了状态变量x、y、z的导数，即f(1)、f(2)和f(3)，计算公式参照了洛伦兹系统的ODE方程。
接下来，函数构建了雅可比矩阵Jac，并用它来计算状态变量X的雅可比矩阵Y的导数值，并将结果保存到f的4到12维。
最后，函数返回计算得到的f向量。
需要注意的是，这是一个函数，需要在其他代码中进行调用，并且需要提前声明全局变量c的值。

function f=SimLorenz_ly(t,X)
global c
x=X(1); y=X(2); z=X(3);
Y= [X(4), X(7), X(10);
    X(5), X(8), X(11);
    X(6), X(9), X(12)];
f=zeros(12,1);
f(1)=10*(y-x);
f(2)=-x*z+(24-4*c)*x+c*y;
f(3)=x*y-8/3*z;
 Jac=[-10,10,0;
     24-4*c-z,  c,  -x;
      y,    x,   -8/3];
  f(4:12)=Jac*Y;

3. 计算序列的C0复杂度
   这是一个用于计算序列的C0复杂度的MATLAB函数。这个函数的输入参数包括一个混沌序列x和一个容限度r。输出参数是C0，代表混沌序列的复杂度。函数首先对输入序列x进行FFT（快速傅里叶变换），然后计算FFT结果的平均幅度平方，保存为变量Gn。接下来，函数创建一个和输入序列x等长的零向量YY，并依次遍历输入序列的每个FFT结果元素。对于每个FFT结果元素，如果它的幅度平方超过容限度r乘以Gn，则将该元素保存到YY向量对应位置。然后，函数对YY进行IFFT（逆傅里叶变换），得到逆变换后的序列xx。接着，函数计算Ssum，即逐差平方和，通过对比原始序列x和逆变换后的序列xx。
   最后，函数将Ssum除以x的平方和，得到C0，即混沌序列的复杂度。
   该函数只提供了计算C0复杂度的功能，具体的使用方法和调用该函数的方式需要在其他代码中进行。

function CO=COFuZadu(x,r)
% 函数名称：COFuZadu
% 函数功能：计算序列的C0复杂度
% 输入参数x，r；r为容限度；x为混沌序列
% 输出参数：C0；C0为输出的混沌序列复杂度
Y=fft(x);
Gn=mean(abs(Y).^2);
YY=zeros(1,length(x));
for i=1:length(x)
    if abs(Y(i))^2>r*Gn
        YY(i)=Y(i);
    end
end
xx=ifft(YY);
Ssum=0;
for i=1:length(x)
    Ssum=Ssum+(xx(i)-x(i))^2;
end
CO=Ssum/sum(x.^2);

4. Lyapunov.m
   作为函数加入测试程序中，计算Lyapunov指数

function [Texp,Lexp]=lyapunov(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp)
%
%    Lyapunov exponent calcullation for ODE-system.
%
%    The alogrithm employed in this m-file for determining Lyapunov
%    exponents was proposed in
%
%         A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano,
%        "Determining Lyapunov Exponents from a Time Series," Physica D,
%        Vol. 16, pp. 285-317, 1985.
%
%    For integrating ODE system can be used any MATLAB ODE-suite methods. 
% This function is a part of MATDS program - toolbox for dynamical system investigation
%    See:    http://www.math.rsu.ru/mexmat/kvm/matds/
%
%    Input parameters:
%      n - number of equation
%      rhs_ext_fcn - handle of function with right hand side of extended ODE-system.
%              This function must include RHS of ODE-system coupled with 
%              variational equation (n items of linearized systems, see Example).                   
%      fcn_integrator - handle of ODE integrator function, for example: @ode45                  
%      tstart - start values of independent value (time t)
%      stept - step on t-variable for Gram-Schmidt renormalization procedure.
%      tend - finish value of time
%      ystart - start point of trajectory of ODE system.
%      ioutp - step of print to MATLAB main window. ioutp==0 - no print, 
%              if ioutp>0 then each ioutp-th point will be print.
%
%    Output parameters:
%      Texp - time values
%      Lexp - Lyapunov exponents to each time value.
%
%    Users have to write their own ODE functions for their specified
%    systems and use handle of this function as rhs_ext_fcn - parameter.      
%
%    Example. Lorenz system:
%               dx/dt = sigma*(y - x)     = f1
%               dy/dt = r*x - y - x*z = f2
%               dz/dt = x*y - b*z     = f3
%
%    The Jacobian of system: 
%        | -sigma  sigma  0 |
%    J = |   r-z    -1   -x |
%        |    y      x   -b |
%
%    Then, the variational equation has a form:
% 
%    F = J*Y
%    where Y is a square matrix with the same dimension as J.
%    Corresponding m-file:
%        function f=lorenz_ext(t,X)
%         SIGMA = 10; R = 28; BETA = 8/3;
%         x=X(1); y=X(2); z=X(3);
%
%         Y= [X(4), X(7), X(10);
%             X(5), X(8), X(11);
%             X(6), X(9), X(12)];
%         f=zeros(9,1);
%         f(1)=SIGMA*(y-x); f(2)=-x*z+R*x-y; f(3)=x*y-BETA*z;
%
%         Jac=[-SIGMA,SIGMA,0; R-z,-1,-x; y, x,-BETA];
%  
%         f(4:12)=Jac*Y;
%
%    Run Lyapunov exponent calculation:
%     
%    [T,Res]=lyapunov(3,@lorenz_ext,@ode45,0,0.5,200,[0 1 0],10);   
%   
%    See files: lorenz_ext, run_lyap.   
%  
% --------------------------------------------------------------------
% Copyright (C) 2004, Govorukhin V.N.
% This file is intended for use with MATLAB and was produced for MATDS-program
% http://www.math.rsu.ru/mexmat/kvm/matds/
% lyapunov.m is free software. lyapunov.m is distributed in the hope that it 
% will be useful, but WITHOUT ANY WARRANTY. 
%
%
%       n=number of nonlinear odes
%       n2=n*(n+1)=total number of odes
%


n1=n; n2=n1*(n1+1);%步数

nit = round((tend-tstart)/stept);

% Memory allocation 

y=zeros(n2,1); cum=zeros(n1,1); y0=y;
gsc=cum; znorm=cum;

% Initial values

y(1:n)=ystart(:);

for i=1:n1 y((n1+1)*i)=1.0; end;

t=tstart;

% Main loop

for ITERLYAP=1:nit

% Solutuion of extended ODE system 

  [T,Y] = feval(fcn_integrator,rhs_ext_fcn,[t t+stept],y);  
  
  t=t+stept;
  y=Y(size(Y,1),:);

  for i=1:n1 
      for j=1:n1 y0(n1*i+j)=y(n1*j+i); end;% 4 7 10 5 8 11  6  9  12
  end;                                     % 4 5  6 7 8 9  10 11  12

%
%       construct new orthonormal basis by gram-schmidt
%

  znorm(1)=0.0;
  for j=1:n1 znorm(1)=znorm(1)+y0(n1*j+1)^2; end;

  znorm(1)=sqrt(znorm(1));

  for j=1:n1 y0(n1*j+1)=y0(n1*j+1)/znorm(1); end;

  for j=2:n1
      for k=1:(j-1)
          gsc(k)=0.0;
          for l=1:n1 gsc(k)=gsc(k)+y0(n1*l+j)*y0(n1*l+k); end;
      end;
 
      for k=1:n1
          for l=1:(j-1)
              y0(n1*k+j)=y0(n1*k+j)-gsc(l)*y0(n1*k+l);
          end;
      end;

      znorm(j)=0.0;
      for k=1:n1 znorm(j)=znorm(j)+y0(n1*k+j)^2; end;
      znorm(j)=sqrt(znorm(j));

      for k=1:n1 y0(n1*k+j)=y0(n1*k+j)/znorm(j); end;
  end;

%
%       update running vector magnitudes
%

  for k=1:n1 cum(k)=cum(k)+log(znorm(k)); end;

%
%       normalize exponent
%

  for k=1:n1 
      lp(k)=cum(k)/(t-tstart); 
  end;

% Output modification

  if ITERLYAP==1
     Lexp=lp;
     Texp=t;
  else
     Lexp=[Lexp; lp];
     Texp=[Texp; t];
  end;

  if (mod(ITERLYAP,ioutp)==0)
     %fprintf('t=%6.4f',t);
     %for k=1:n1 fprintf(' %10.6f,',lp(k)); end;
     %fprintf('\n');
  end;

  for i=1:n1 
      for j=1:n1
          y(n1*j+i)=y0(n1*i+j);
      end;
  end;

end;

5. 测试程序TestSimpLorenz（全）

%绘制 Lorenz 吸引子在 x-z 平面上的投影图
clc;clear
global c
c=2;
[T,Y]=ode45(@chao_SimpleLorenz,0:0.01:500,[1 2 3 ]);
figure
plot(Y(40000:end,1),Y(40000:end,3),'b')
xlim([-15 15])
xlabel('\itx')
ylabel('\itz')
grid on;
set(gca,'linewidth',0.5,'fontsize',12,'fontname','Times');             %在时间范围内使用 “ode45” 函数求解 “chao_SimpleLorenz” 的微分方程，并绘制 Lorenz 吸引子的投影图。
%--------------------------------------------------------------------------
clc;clear
global c
c=2;                                                                   %使用 “lyapunov” 函数计算 Lorenz 系统的 Lyapunov 指数，并将结果存储在 “Res” 变量中。最后打印输出了最终的 Lyapunov 指数。

[T,Res]=lyapunov(3,@SimLorenz_ly,@ode45,0,0.5,500,[ 0.1 0 0.1],5);          %最后打印输出了最终的 Lyapunov 指数。
Ly=Res(end,:)
%----------------------Lyapunov 指数用于描述混沌系统中的指数敏感性-------------------------------
clc
clear
global c
% C=-2:0.05:8;L=length(C);

L=250;
C=linspace(-2,8,L);
ly=zeros(3,L);
Com=zeros(1,L);
figure
for i=1:L
    c=C(i);
    
    [T,Y]=ode45(@chao_SimpleLorenz,0:0.01:100,[1 2 3]);
    data=Y(5000:end,3);
    Com(i)=COFuZadu(data,15);
    for j=2:(length(data)-1)
        if data(j)>data(j-1)&&data(j)>data(j+1)
            plot(c,data(j),'.r','markersize',1);
            hold on;
            if j==20
                break;
            end
        end
    end
    
   [T,Res]=lyapunov(3,@SimLorenz_ly,@ode45,0,0.5,100,[ 0.1 0 0.1],5);
    ly(:,i)=Res(end,:);
    disp(i)				
end
xlabel('\itc')
ylabel('{\itz}_{max}')
grid on;
set(gca,'linewidth',0.5,'fontsize',12,'fontname','Times');
%------------------绘制参数 c 对系统行为的影响---------------------

figure
plot(C,ly(1,:),'r','linewidth',1)
hold on
plot(C,ly(2,:),'k','linewidth',1)
plot(C,ly(3,:),'b','linewidth',1)
ylim([-15,2])
xlabel('\itc')
ylabel('Lyapunov')
grid on;
set(gca,'linewidth',0.5,'fontsize',12,'fontname','Times');


figure
plot(C,Com,'r','linewidth',1)
xlabel('\itc')
ylabel('C_0')
grid on;
set(gca,'linewidth',0.5,'fontsize',12,'fontname','Times');