matlab编写的Lyapunov指数计算程序Word格式.docx

1、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 ioutp0 then each ioutp-th point will be

2、 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%Th

3、e 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）,

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

5、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%options=odeset（RelTol,

6、DS（1）.rel_error,AbsTol,DS（1）.abs_error,MaxStep,DS（1）.max_step,.OutputFcn,odeoutp,Refine,0,InitialStep,0.001）;n_exp=DS（1）.n_lyapunov;n1=n;n2=n1*（n_exp+1）;neq=n2;%Number of steps nit=round（tend-tstart）/stept）+1;%Memory allocation y=zeros（n2,1）;cum=zeros（n2,1）;y0=y;gsc=cum;znorm=cum;%Initial values y（1

7、:n）=ystart（:）;for i=1:n_exp y（n1+1）*i）=1.0;end;t=tstart;Fig_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）;ylabel（Lyapunov expon

8、ents）;Fig_Lyap_Axes=findobj（Fig_Lyap,Type,axes）;n_exp PlotLyapi=plot（tstart,0）;uu=findobj（Fig_Lyap,Type,line）;size（uu,1）set（uu（i）,EraseMode,none）;set（uu（i）,XData,YData,）;set（uu（i）,Color,0 0 rand）;end ITERLYAP=0;%Main loop calculation_progress=1;while ttend,tt=tend;%Solutuion of extended ODE system%T

9、,Y=feval（fcn_integrator,rhs_ext_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;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（bufe

10、r_i）;calculation_progress=1;else calculation_progress=0;if（calculation_progress=99）break;else calculation_progress=1;t=tt;y=Y（size（Y,1）,:%construct new orthonormal basis by gram-schmidt%znorm（1）=0.0;for j=1:n1 znorm（1）=znorm（1）+y（n1+j）2;znorm（1）=sqrt（znorm（1）;n1 y（n1+j）=y（n1+j）/znorm（1）;for j=2:n_ex

11、p 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）;for k=1:n1 for l=1:（j-1）y（n1*j+k）=y（n1*j+k）-gsc（l）*y（n1*l+k）;znorm（j）=0.0;n1 znorm（j）=znorm（j）+y（n1*j+k）2;znorm（j）=sqrt（znorm（j）;n1 y（n1*j+k）=y（n1*j+k）/znorm（j）;%update running vector magnitudes%for k=1:n_exp cum（k）=cum（k）+log（zn

12、orm（k）;%normalize exponent%rescale=0;u1=1.e10;u2=-1.e10;n_exp lp（k）=cum（k）/（t-tstart）;%Plot Xd=get（PlotLyapk,Xdata）;Yd=get（PlotLyapk,Ydata）;if timeplott u1=min（u1,min（Yd）;u2=max（u2,max（Yd）;Xd=Xd t;Yd=Yd lp（k）;set（PlotLyapk,Xdata,Xd,Ydata,Yd）;if timeplot0 then each ioutp-th point will be print.%Outpu

13、t parameters:/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）;%Number of steps nit=round（tend-tstart）/stept）;cum=zeros（n1,

14、1）;n1 y（n1+1）*i）=1.0;%Main loop for ITERLYAP=1:nit%Solutuion of extended ODE system T,Y=unit（t,stept,y,d）;t=t+stept;n1 for j=1:n1 y0（n1*i+j）=y（n1*j+i）;n1 znorm（1）=znorm（1）+y0（n1*j+1）2;n1 y0（n1*j+1）=y0（n1*j+1）/znorm（1）;n1 for k=1:n1 gsc（k）=gsc（k）+y0（n1*l+j）*y0（n1*l+k）;（j-1）y0（n1*k+j）=y0（n1*k+j）-gsc（l

15、）*y0（n1*k+l）;n1 znorm（j）=znorm（j）+y0（n1*k+j）2;n1 y0（n1*k+j）=y0（n1*k+j）/znorm（j）;n1 cum（k）=cum（k）+log（znorm（k）;%normalize exponent%for k=1:n1 lp（k）=cum（k）/（t-tstart）;%Output modification if ITERLYAP=1 Lexp=lp;Texp=t;else Lexp=lp;n1 y（n1*j+i）=y0（n1*i+j）;五、小数据量法计算 Lyapunov 指数的 Matlab 程序-（mex 函数，超快）http

16、:/ 六、计算 lyapunov 指数的程序 program lylorenz parameter（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.a=10.b=8./3.r=28.t=0.h=0.01!initial conditions do i=n+1,m y（i）=0.end do do i=1,n y（n+1）*i）=1.s（i）=0 end do open（1,file=lorenz1.dat）open（2,file=l

17、y lorenz.dat）do 100 k=1,st!st iterations call rgkt（m,h,t,y,f,yc,y1）!normarize vector 123 z=0.do i=1,n 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）enddo 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+

18、2）k3=k3+y（3*i+1）*y（3*i+3）end do !conduct new vector2 and 3 do 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 new vector2 and 3,normarize them do 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）/z（i）end d

19、o end do !update lyapunov exponent do i=1,n if（z（i）0）s（i）=s（i）+log（z（i）enddo 100 continue do i=1,n s（i）=s（i）/（1.*st*h*1000.）write（2,*）s（i）enddo end!subroutine rgkt（m,h,t,y,f,yc,y1）real y（m）,f（m）,y1（m）,yc（m）,a,b,r integer:i,j do j=1,1000 call df（m,t,y,f）t=t+h/2.0 do i=1,m yc（i）=y（i）+h*f（i）/2.0 y1（i）=

20、y（i）+h*f（i）/6.0 end do call df（m,t,yc,f）do i=1,m yc（i）=y（i）+h*f（i）/2.0 y1（i）=y1（i）+h*f（i）/3.0 end do call df（m,t,yc,f）t=t+h/2.0 do i=1,m yc（i）=y（i）+h*f（i）y1（i）=y1（i）+h*f（i）/3.0 end do call df（m,t,yc,f）do i=1,m y（i）=y1（i）+h*f（i）/6.0 end do if（j500）write（1,*）t,y（1）,y（2）,y（3）end do return end!subroutine df（m,t,y,f）real y（m）,a,b,r,f（m）common a,b,r a=10.b=8./3.r=28.f（1）=a*（y（2）-y（1）f（2）=y（1）*（r-y（3）-y（2）f（3）=y（1）*y（2）-b*y（3）do i=0,2 f（4+i）=a*y（7+i）-y（4+i）f（7+i）=y（4+i）*（r-y（3）-y（7+i）-y（1）*y（10+i）f（10+i）=y（2）*y（4+i）-b*y（10+i）+y（1）*y（7+i）enddo return end 七、计算各种混沌系统李雅普洛夫指数的 MATLAB 源程序 http:/