lyapunov.m

LE计算方法。这里给出了经典的分析方法。

function [x,Lce,ErrorQ]=lyapunov(F1,F2,Ic,Ntot,LceFlag,Nskip,OrthFlag);
% lyapunov Computes the trajectory and the Lyapunov characteristic       
% exponents for discrete systems via the Householder QR Based Method.
% For more information, see "An efficient QR based method for the 
% computation of Lyapunov Exponents", by Hubertus F. von Bremen, 
% Firdaus  E. Udwadia and Wlodek Proskurowski, Physica D, 1997, to appear.
%  Program written by H. von Bremen 10/10/96 
%                                                                    
% [x,Lce,ErrorQ]=lyapunov('xfile','tanmfile',Ic,Ntot,LceFlag,Nskip,OrthFlag)
%   computes the trajectory of a discrete dynamical system which is      
%   specified in the M-file xfile.m .  It also computes the Lyapunov     
%   characteristic exponents (Lce's), provided the tangent maps are      
%   specified in the M-file tanmfile.m.  The trajectory is computed for 
%   up to Ntot iterations.  The LCE's are computed starting at the Nskip+1
%   iteration of the trajectory, and are computed for Ntot-Nskip iterations.                             
%                                                                       
%
%   INPUT:                                                               
%   xfile       String containing the user-supplied function file that defines
%               the trajectory.                                         
%                       CALL: xmap=fun(x,iter)  where fun = 'xfile'     
%                       x       - state vector                          
%                       iter    - iteration number                      
%                       xmap    - returned next trajectory point               
%                               {x(iter+1)=xmap=fun(x(iter),iter)}      
%   tanmfile    String containing the user-supplied file that defines the       
%               tangent map to the trajectory.                          
%                       CALL: tanmap=fun(x,iter)  where fun = 'tanmfile'
%                       x       - state vector                          
%                       iter    - iteration number                      
%                       tanmap  - returned next tangent map               
%                       {Tangent_Map(iter+1)=tanmap=fun(x(iter),iter)}  
%   Ic          Initial conditions for the trajectory, given as a column
%               vector                                                                  
%   Ntot        Number of iterations (time steps) to be used to calculate
%               the trajectory of the system.                           
%   LceFlag:    Flag used to determine if Lce's are computed.           
%               -If LceFlag = 1, then Lce's are computed.               
%               -If LceFlag = 0, then LCE's are not computed and only   
%                the trajectory is computed.                    
%   Nskip:      Number of iterations to be skipped before the LCE's are 
%               to be computed.  First Nskip points along the trajectory are
%               skipped in order to allow the disturbances due to the 
%               initial conditions to diminish, before LCE's are computed.                      
%   OrthFlag:   Flag used to determine if error in orthogonality is to  
%               be computed, and if computed, then of what type.
%               -If OrthFlag = 0 (recommended value), then the error in 
%               orthogonality is not computed. 
%               -If OrthFlag = 1, then the determinant error in orthogonality        
%                is computed.  Error in Orthogonality = abs(1-abs(det(Q)))
%               -If OrthFlag = 2, then the Two-Norm error in orthogonality    
%                is computed.  Error in Orthogonality = norm(Q'*Q-I)    
%               NOTE: when using OrthFlag = 1 or 2, the computational cost is
%               increased significantly (using 2 is more expensive than 1).                             
%                                                                       
%   OUTPUT:                                                              
%       x       Trajectory of dynamical system given in matrix form, where the 
%               k-th column are the trajectory coordinates for the k-th 
%               iteration (the first column is Ic). Size(x)=[n,Ntot+1].
%       Lce     Computed Lyapunov exponents given in matrix form, where the 
%               k-th column are the LCE's for the k-th iteration of LCE's.  
%               Size(Lce)=[n,Ntot-Nskip].
%       ErrorQ  Error in othogonality given in vector form, where the k-th 
%               entry is the error in orthogonality of the k-th Lce iteration.  
%               Size(ErrorQ)=[1,Ntot-Nskip].   
%               -If OrthFlag = 1, ErrorQ = abs(1-abs(det(Q)))
%               -If OrthFlag = 2, ErrorQ = Norm(Q'*Q-I)                 
%                       
% 
% **SAMPLE:
%   *************
%   Sample INPUT: 
%   Ic=[0.2  0.3]'; Ntot= 1000; LceFlag=1; Nskip=100; OrthFlag=0;
% 
%   Sample function files for the trajectory and the tangent map:
%   ---------------------------
%   Sample for xfile='map'; coupled logistic maps (stored in m-file map.m)
%   function xmap=map(x,iter)
%   d=.2; r=3.6;
%   xm(1)=d*r*x(1)*(1-x(1))     +    (1-d)*r*x(2)*(1-x(2));
%   xm(2)=(1-d)*r*x(1)*(1-x(1)) +    d*r*x(2)*(1-x(2));
%   xmap=[xm(1) xm(2)]';
%   ---------------------------
%   Sample for tanmfile='tanmap'; Jacobian for coupled logistic maps (stored in 
%   m-file tanmap.m)
%   function [T]=tanmap(x)
%   d=.2; r=3.6;
%   T=  [  d*r*(1-2*x(1))        (1-d)*r*(1-2*x(2))
%         (1-d)*r*(1-2*x(1))       d*r*(1-2*x(2)) ];
%   ---------------------------
%   Note: the above two sample function files are provided in the files map.m 
%   and tanmap.m.
%   --------------------------------
%   *************  
%   Sample CALL
%   [x,Lce,ErrorQ]=lyapunov('map','tanmap',Ic,Ntot,LceFlag,Nskip,OrthFlag);
%   *************
%   Sample OUTPUT:
%   To plot the LCE's for each iteration you can use:  
%   subplot(2,1,1), plot(Lce(1,:))
%   subplot(2,1,2), plot(Lce(2,:))
%   To plot the trajectory coordinates for each iteration you can use: 
%   subplot(2,1,1), plot(x(1,:))
%   subplot(2,1,2), plot(x(2,:))
%   
%   Note: By setting OrthFlag=1 (or = 2), you can plot the error in 
%   orthogonality by using: plot(ErrorQ)
%  
%   In general we have:
%   To plot the k-th LCE for each iteration:  plot(Lce(k,:))
%   To plot the k-th coordinate of the trajectory:  plot(x(k,:))
%   If OrthFlag is 1 or 2,then to plot the error in orthogonality: plot(ErrorQ)
%

% INPUT ERROR HANDLING
if (LceFlag~=0)&(LceFlag~=1)
  disp('You have entered an invalid value for LceFlag')
  disp('You should choose LceFlag to be either 0 or 1 (see help file)')
  disp('The value of LceFlag was set to 0 ')
  disp(' ')
  LceFlag=0;
end

if (OrthFlag~=0)&(OrthFlag~=1)&(OrthFlag~=2)
  disp('You have entered an invalid value for OrthFlag')
  disp('You should choose OrthFlag to be either 0 or 1 or 2 (see help file)')
  disp('The value of OrthFlag was set to 0')
  disp(' ')
  OrthFlag=0;
end

if Ntot<Nskip
   disp('You have entered an invalid value for either Ntot or Nskip')
   disp('You must have Ntot>=Nskip, (see help file)')
   disp('The values of Nskip and Ntot have been switched')
   disp(' ')
   temp1=Nskip;
   Nskip=Ntot;
   Ntot=temp1;
end


%  INITIALIZATION PART
  m=Ntot-Nskip;
  n=max(size(Ic));
  x=zeros(n,Ntot+1);
% END INITIALIZATION

x(:,1)=Ic;
for j=2:Nskip+1
  Ic=feval(F1,Ic);
  x(:,j)=Ic;
end  % (j)

if LceFlag==0
   for j=1:m
     Ic=feval(F1,Ic);
   x(:,j+Nskip+1)=Ic;
   end  % (j)
else 
  r=zeros(1,n);gama=zeros(1,n);
  q=eye(n);sl=zeros(n,1);Lce=zeros(n,m);
  if OrthFlag==1|OrthFlag==2
     ErrorQ=zeros(1,m);
  end  % (if)
  A=feval(F2,Ic);
  if OrthFlag==0
     a=A;
  end;
  N2=Nskip+1;
  for i=1:m
     if OrthFlag~=0
        a=A*q;
        Ic=feval(F1,Ic);
        x(:,i+N2)=Ic;
        A=feval(F2,Ic);
        B=eye(n);
     elseif OrthFlag == 0 
       Ic=feval(F1,Ic);
       x(:,i+N2)=Ic;
       B=feval(F2,Ic);
     end % (if OrthFlag)
% Computation of the Factorization 
        for k=1:n-1
% computation of the reflector
        if a(k,k)<0
          b=-1;
          else
            b=1;
          end
          sig=sqrt(a(k:n,k)'*a(k:n,k));
          gama=sig*(sig+abs(a(k,k)));
          r(k)=-b*sig;
          a(k,k)=a(k,k)-r(k);
% end computation of the reflector
         p=k+1;
          for j=p:n
             b=a(k:n,k)'*a(k:n,j)/gama;
             a(p:n,j)=a(p:n,j)-b*a(p:n,k);
          end
% Computation of B*Q
          for j=1:n  
             b=B(j,k:n)*a(k:n,k)/gama;
             B(j,k:n)=B(j,k:n)-a(k:n,k)'*b;
           end
        end
        r(n)=a(n,n);
% End Computation of the factorization 
     if OrthFlag~=0
        q=B;
        if OrthFlag==1
           ErrorQ(i)=abs(1-abs(det(q)));
        elseif OrthFlag==2
           ErrorQ(i)=norm(q'*q-eye(n));
        end  % (if)
     elseif OrthFlag == 0 % (if OrthFlag)
        a=B;
     end  % (if OrthFlag)
     sl=sl+log(abs(r'));
     Lce(:,i)=sl/i;
  end  % (i)
end % (else)