matlab编写的Lyapunov指数计算程序

1、matlab编写的Lyapunov指数计算程序 一、计算连续方程Lyapunov指数的程序其中给出了两个例子：! e  % u9 j* * b" m$ j$ I计算Rossler吸引子的Lyapunov指数  _9 : B" T* t! U4 J2 Q9 / n计算超混沌Rossler吸引子的Lyapunov指数二、recnstitution重构相空间，在非线性系统分析中有重要的作用function Texp,Lexp=lyapunov(n,tstart,stept,tend,ystart,ioutp);global DS;globa

2、l P;global calculation_progress first_call;global driver_window;global TRJ_bufer Time_bufer bufer_i;%    Lyapunov exponent calcullation for ODE-system.%    The alogrithm employed in this m-file for determining Lyapunov%    exponents was proposed in%     &

3、#160;   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

4、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%     

5、0;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

6、linearized systems, see Example).                   %      fcn_integrator - handle of ODE integrator function, for example: ode45                  %   

7、;   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.% 

8、    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%    

9、  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:%    

10、60;         dx/dt = sigma*(y - x)     = f1%               dy/dt = r*x - y - x*z = f2%               dz/dt = x*y - b*z     = f3%

11、60;   The Jacobian of system: %        | -sigma  sigma  0 |%    J = |   r-z    -1   -x |%        |    y      x   -b |%    The

12、n, 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

13、= 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);%      

14、;   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 expon

15、ent 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 w

16、as 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

17、odes%options = odeset('RelTol',DS(1).rel_error,'AbsTol',DS(1).abs_error,'MaxStep',DS(1).max_step,.                 'OutputFcn',odeoutp,'Refine',0,'InitialStep',0.001);      

18、;       n_exp = DS(1).n_lyapunov;n1=n; n2=n1*(n_exp+1);neq=n2;%  Number of stepsnit = round(tend-tstart)/stept)+1;% Memory allocation y=zeros(n2,1); cum=zeros(n2,1); y0=y;gsc=cum; znorm=cum;% Initial valuesy(1:n)=ystart(:);for i=1:n_exp y(n1+1)*i)=1.0; end;t=tstart;Fig_

19、Lyap = figure;set(Fig_Lyap,'Name','Lyapunov exponents','NumberTitle','off');set(Fig_Lyap,'CloseRequestFcn','');hold on;box on;timeplot = tstart+(tend-tstart)/10;axis(tstart timeplot -1 1);title('Dynamics of Lyapunov exponents');xlabel('t

20、9;);ylabel('Lyapunov exponents');Fig_Lyap_Axes = findobj(Fig_Lyap,'Type','axes');for i=1:n_exp       PlotLyapi=plot(tstart,0);        end;       uu=findobj(Fig_Lyap,'Type','line');   &#

21、160;    for i=1:size(uu,1)           set(uu(i),'EraseMode','none') ;           set(uu(i),'XData','YData',);           set(uu(i),'Color

22、',0 0 rand);       endITERLYAP = 0;% Main loopcalculation_progress = 1;while t<tend    tt = t + stept;    ITERLYAP = ITERLYAP+1;    if tt>tend, tt = tend; end;% Solutuion of extended ODE system %  T,Y = feval(fcn_integrator,rhs_ex

23、t_fcn,t t+stept,y);        while calculation_progress = 1       T,Y = integrator(DS(1).method_int,ode_lin,t tt,y,options,P,n,neq,n_exp);       first_call = 0;       if calculation_progress = 99, break; end;

24、0;      if ( T(size(T,1)<tt ) & (calculation_progress=0)         y=Y(size(Y,1),:);         y(1,1:n)=TRJ_bufer(bufer_i,1:n);         t = Time_bufer(bufer_i);       

25、  calculation_progress = 1;       else         calculation_progress = 0;       end;   end;   if (calculation_progress = 99)        break;    else       cal

26、culation_progress = 1;   end;     t=tt;  y=Y(size(Y,1),:);    first_call = 0;%       construct new orthonormal basis by gram-schmidt%  znorm(1)=0.0;  for j=1:n1 znorm(1)=znorm(1)+y(n1+j)2; end; &

27、#160;znorm(1)=sqrt(znorm(1);  for j=1:n1 y(n1+j)=y(n1+j)/znorm(1); end;  for j=2:n_exp      for k=1:(j-1)          gsc(k)=0.0;          for l=1:n1 gsc(k)=gsc(k)+y(n1*j+l)*y(n1*k+l); end;   &

28、#160;  end;      for k=1:n1          for l=1:(j-1)              y(n1*j+k)=y(n1*j+k)-gsc(l)*y(n1*l+k);          end;      end;    

29、 znorm(j)=0.0;      for k=1:n1 znorm(j)=znorm(j)+y(n1*j+k)2; end;      znorm(j)=sqrt(znorm(j);      for k=1:n1 y(n1*j+k)=y(n1*j+k)/znorm(j); end;  end;%       update running vector magnitudes%  f

30、or k=1:n_exp cum(k)=cum(k)+log(znorm(k); end;%       normalize exponent%  rescale = 0;  u1 =1.e10;  u2 =-1.e10;  for k=1:n_exp       lp(k)=cum(k)/(t-tstart); %  Plot       Xd=get(PlotLyapk,&

31、#39;Xdata');      Yd=get(PlotLyapk,'Ydata');      if timeplot<t            u1=min(u1,min(Yd);            u2=max(u2,max(Yd);      end;   

32、;   Xd=Xd t; Yd=Yd lp(k);      set(PlotLyapk,'Xdata',Xd,'Ydata',Yd);  end;  if timeplot<t     timeplot = timeplot+(tend-tstart)/20;     figure(Fig_Lyap);     axis(tstart t

33、imeplot u1 u2); end;  drawnow;% Output modification  if ITERLYAP=1     Lexp=lp;     Texp=t;  else     Lexp=Lexp; lp;     Texp=Texp; t;  end;    if (mod(ITERLYAP

34、,ioutp)=0)     for k=1:n_exp         txtstringk='lambda_' int2str(k) '=' num2str(lp(k);     end     legend(Fig_Lyap_Axes,txtstring,3);  end;end;   ss=warndlg('Attention!

35、Plot of lyapunov exponents will be closed!','Press OK to continue!');   uiwait(ss);   delete(Fig_Lyap);         fprintf('n n Results of Lyapunov exponents calculation: n t=%6.4f',t);     for k=1:n_exp fprintf(&

36、#39; L%d=%f; ',k,lp(k); end;     fprintf('n');            if isempty(driver_window)             if ishandle(driver_window)               &

37、#160; delete(driver_window);                 driver_window = ;             end;                  end;        &

38、#160;   calculation_progress = 0;   update_ds;三、wolf 方法计算李雅普诺夫指数 四、给出了分形计算的源代码的matlab程序，可以迅速帮助大家进行分形的分析与计算，参数容易设置function Texp,Lexp=new1lyapunov(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp,d);%    Lyapunov exponent calcullation for ODE-system.% 

39、60; 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,&quo

40、t; 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.rs

41、u.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 R

42、HS of ODE-system coupled with %              variational equation (n items of linearized systems, see Example).                   %      fcn_integrator - handle of ODE in

43、tegrator function, for example: ode45                  %      tstart - start values of independent value (time t)%      stept - step on t-variable for Gram-Schmidt renormalization procedure.%  &#

44、160;   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

45、 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 t

46、his function as rhs_ext_fcn - parameter.      %    Example. Lorenz system:%               dx/dt = sigma*(y - x)     = f1%               dy/dt = r*x -

47、y - x*z = f2%               dz/dt = x*y - b*z     = f3%    The Jacobian of system: %        | -sigma  sigma  0 |%    J = |   r-z    -1 

48、  -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:%&

49、#160;       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);%         &#

50、160;   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,SIG

51、MA,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_e

52、xt, 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 WIT

53、HOUT ANY WARRANTY. %       n=number of nonlinear odes%       n2=n*(n+1)=total number of odes%n1=n; n2=n1*(n1+1);%  Number of stepsnit = round(tend-tstart)/stept);% Memory allocation y=zeros(n2,1); cum=zeros(n1,1); y0=y;gsc=cum; znorm=cum;% Initial va

54、luesy(1:n)=ystart(:);for i=1:n1 y(n1+1)*i)=1.0; end;t=tstart;% Main loopfor ITERLYAP=1:nit% Solutuion of extended ODE system   T,Y = unit(t,stept,y,d);      t=t+stept;  y=Y(size(Y,1),:);  for i=1:n1       for j=1:n1 y0(n

55、1*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; 

56、60;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         

57、 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; e

58、nd;      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 expon

59、ent%  for k=1:n1       lp(k)=cum(k)/(t-tstart);   end;% Output modification  if ITERLYAP=1     Lexp=lp;     Texp=t;  else     Lexp=lp;     Texp=t;  en

60、d;  for i=1:n1       for j=1:n1          y(n1*j+i)=y0(n1*i+j);      end;  end;end;五、小数据量法计算 Lyapunov 指数的 Matlab 程序 - (mex 函数，超快) 六、计算lyapunov指数的程序program lylorenz        parameter(

61、n=3,m=12,st=100)        integer:i,j,k    real y(m),z(n),s(n),yc(m),h,y1(m),a,b,r,f(m),k1,k2,k3        y(1)=10.        y(2)=1.        y(3)=0.      

62、0; a=10.    b=8./3.    r=28.        t=0.        h=0.01!initial conditionsdo i=n+1,m        y(i)=0.end dodo  i=1,n        y(n+1)*i)=1.   

63、; s(i)=0end doopen(1,file='lorenz1.dat')open(2,file='ly lorenz.dat')do 100 k=1,st   !st iterationscall rgkt(m,h,t,y,f,yc,y1)    !normarize vector 123        z=0.       do i=1,n         

64、;do j=1,n        z(i)=z(i)+y(n*j+i)*2        enddo        if(z(i)>0.)z(i)=sqrt(z(i)        do j=1,n    y(n*j+i)=y(n*j+i)/z(i)        endd

65、o    end do   !generate gsr coefficient   k1=0.   k2=0.   k3=0.do i=1,n                k1=k1+y(3*i+1)*y(3*i+2)        k2=k2+y(3*i+3)*y(3*i+2)    

66、    k3=k3+y(3*i+1)*y(3*i+3)end do   !conduct new vector2 and 3do i=1,n        y(3*i+2)=y(3*i+2)-k1*y(3*i+1)            y(3*i+3)=y(3*i+3)-k2*y(3*i+2)-k3*y(3*i+1)end do       !generate ne

67、w vector2 and 3,normarize themdo i=2,n    z(i)=0.    do j=2,n        z(i)=z(i)+y(n*j+i)*2        enddo    if(z(i)>0.)z(i)=sqrt(z(i)        do j=2,n    y(n*j+i)=y(n*j+i

68、)/z(i)        end doend do    !update lyapunov exponentdo i=1,n   if(z(i)>0)s(i)=s(i)+log(z(i)enddo100 continuedo i=1,n  s(i)=s(i)/(1.*st*h*1000.)write(2,*)s(i)enddoend!subroutine rgkt(m,h,t,y,f,yc,y1)real y(m),f(m),y1(m),yc(m),a,b,rinteger:i,jdo j=1,1000