请大家帮忙看看这个两个求最大Lyapunov指数的程序

wannete

有两种方法，都是根据指数工具箱来的，就是第二个降了一维，但是得到的结果就不一样！
function OUT=duffing(t,X)
k=0.5;
f=2.804;
wd=2*pi;
%Rearrange input data in desired format
%Note: the input data is a column vector
x=X(1); y=X(2);z=X(3);
Q=[X(4),X(7),X(10);
   X(5),X(8),X(11);
   X(6),X(9),X(12)];
%Duffing's equation
dx=wd*y;
dy=wd*(-k*y+x-x^3+f*cos(wd*z));
dz=wd*1;  %where z = t, this transformation is for changing
                %the non-autonomous system to a autonomous one
DX1=[dx;dy;dz]; %Output data
%Linearized system
J=[    0,           wd*1,             0;
   wd*(1-3*x^2),    -k*wd,      -f*wd*wd*sin(wd*z);
        0,             0,             0];
%Variational equation
F=J*Q;
%Put output data in a column vector
OUT=[DX1;F(:)];


function [Texp,Lexp]=lyapunov(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp);
n1=n; n2=n1*(n1+1);
%  Number of steps
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;
  end;
%
%       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;

  for k=1:n1 cum(k)=cum(k)+log(znorm(k)); end;
  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;



[T,Res]=lyapunov(3,@duffing,@ode45,0,0.5,200,[0 1 0],10);
plot(T,Res);
title('Dynamics of Lyapunov exponents');
xlabel('Time'); ylabel('Lyapunov exponents');
得到的结果是：

方法二：
%计算Duffing振子的lyapunov指数
clear all;
clc;
global X V
yinit=[1,1];
orthmatrix=[1 0;
            0 1];
y=zeros(6,1);
%初始化输入
IC(1:2)=yinit;
IC(3:6)=orthmatrix;
InitiaTime=0;  %时间初始值                                       0
TimeStep=1e-2;  %时间步长                                        0.02
wholetimes=10000;  %总的循环次数                                 10000
UpdateStepNum=10;%每次演化的步数                                 10
iteratetimes =wholetimes/UpdateStepNum; %演化的次数              1000  iterate 迭代的意思
DiscardItr=100;                  %                              100
Iteration=UpdateStepNum*TimeStep;   %                           0.2                           
DiscardTime=DiscardItr*Iteration+InitiaTime;      %             40
T1=InitiaTime;           %       T1=0
T2=T1+Iteration;         %       T2=0.2
Tspan=[T1:TimeStep:T2];  %       Tspan=[0:0.02:0.2]
options=odeset('RelTol',1e-5,'AbsTol',1e-6);
mod=zeros(1,2);
d=2;
Sum=zeros(1,d);
n=0;
k=0;
xData=[];
yData=[];

for i=1:iteratetimes                         %1:1000
    n=n+1
    [t,X]=ode45('twofun',Tspan,IC,options);
    %取积分得到的最后一个时刻的值
    [rX,cX]=size(X);
    L=cX-d*d;
   
    for i=1:d
        m1=L+1+(i-1)*d;
        m2=m1+d-1;
        A(:,i)=(X(rX,m1:m2))';
    end
   y0=TwoGS(A);
   
    %取两个向量的模
    mod(1)=sqrt(y0(:,1)'*y0(:,1));
    mod(2)=sqrt(y0(:,2)'*y0(:,2));
    y0(:,1)=y0(:,1)/mod(1);
    y0(:,2)=y0(:,2)/mod(2);
    if T2>DiscardTime
        Q0=y0;
    else
        Q0=eye(d);          %Y = eye(n) returns the n-by-n identity matrix.
    end
   
    if (T2>DiscardTime)
        k=k+1;
        T=k*Iteration;
        TT=n*Iteration+InitiaTime;
        Sum=Sum+log(abs(mod));
        lambda=Sum/T;
        
        Lambda=fliplr(sort(lambda));  %B = sort(A) sorts the elements along different dimensions of an array, and
                                      %arranges those elements in ascending order.
                               %B = fliplr(A) returns A with columns flipped in the left-right direction,
                               % that is, about a vertical axis.If A is a row vector, then fliplr(A) returns
                               %a vector of the same length with the order of its elements reversed. If A is
                               %a column vector, then fliplr(A) simply returns A.
         [rxD,cxD]=size(xData);
         [ryD,cyD]=size(yData);
         
         xData=[xData;TT];
         yData=[yData;lambda];
    end
    ic=X(rX,1:L);
    T1=T1+Iteration;
    T2=T2+Iteration;
    Tspan=[T1:TimeStep:T2];
   
    IC=[ic(:);y0(:)];
end
plot(xData,yData(:,1),xData,yData(:,2))


function dX=twofun(t,X)
k=0.5;
B=2.804;
wd=2*pi;
x=X(1);
y=X(2);
Y=[X(3) X(5);
   X(4) X(6)];
dX=zeros(6,1);
dX(1)=wd*y;
dX(2)=wd*(-k*y+x-x^3+B*cos(wd*t));
J=[ 0,  wd*1;
wd*(1-3*x^2),   wd*(-k)];
dX(3:6)=J*Y;


function A=TwoGS(V)  %V为3*3的向量
global a1 a2
v1=V(:,1);
v2=V(:,2);
a1=zeros(2,1);
a2=zeros(2,1);
a1=v1;
a2=v2-((a1'*v2)/(a1'*a1))*a1;
A=[a1,a2];
得到的结果是：

octopussheng

你犯了一个很典型的错误。

不知道你注意到没有，Lyapunov指数的定义，是基于自治系统，而不是非自治系统的

因此，我们在求LE的时候，都需要将原来的非自治系统转化为自治系统，在根据各种方法求解。

sizhiyuan2006

学习了，谢谢