matlab编写的Lyapunov指数计算程序Word格式.docx
《matlab编写的Lyapunov指数计算程序Word格式.docx》由会员分享，可在线阅读，更多相关《matlab编写的Lyapunov指数计算程序Word格式.docx（18页珍藏版）》请在冰豆网上搜索。
ode45%tstart-startvaluesofindependentvalue（timet）%stept-stepont-variableforGram-Schmidtrenormalizationprocedure.%tend-finishvalueoftime%ystart-startpointoftrajectoryofODEsystem.%ioutp-stepofprinttoMATLABmainwindow.ioutp=0-noprint,%ifioutp0theneachioutp-thpointwillbeprint.%Outputparameters:
%Texp-timevalues%Lexp-Lyapunovexponentstoeachtimevalue.%UsershavetowritetheirownODEfunctionsfortheirspecified%systemsandusehandleofthisfunctionasrhs_ext_fcn-parameter.%Example.Lorenzsystem:
%dx/dt=sigma*（y-x）=f1%dy/dt=r*x-y-x*z=f2%dz/dt=x*y-b*z=f3%TheJacobianofsystem:
%|-sigmasigma0|%J=|r-z-1-x|%|yx-b|%Then,thevariationalequationhasaform:
%F=J*Y%whereYisasquarematrixwiththesamedimensionasJ.%Correspondingm-file:
%functionf=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;
%RunLyapunovexponentcalculation:
%T,Res=lyapunov（3,lorenz_ext,ode45,0,0.5,200,010,10）;
%Seefiles:
lorenz_ext,run_lyap.%-%Copyright（C）2004,GovorukhinV.N.%ThisfileisintendedforusewithMATLABandwasproducedforMATDS-program%http:
/www.math.rsu.ru/mexmat/kvm/matds/%lyapunov.misfreesoftware.lyapunov.misdistributedinthehopethatit%willbeuseful,butWITHOUTANYWARRANTY.%n=numberofnonlinearodes%n2=n*（n+1）=totalnumberofodes%options=odeset（RelTol,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;
%Numberofstepsnit=round（tend-tstart）/stept）+1;
%Memoryallocationy=zeros（n2,1）;
cum=zeros（n2,1）;
y0=y;
gsc=cum;
znorm=cum;
%Initialvaluesy（1:
n）=ystart（:
）;
fori=1:
n_expy（n1+1）*i）=1.0;
end;
t=tstart;
Fig_Lyap=figure;
set（Fig_Lyap,Name,Lyapunovexponents,NumberTitle,off）;
set（Fig_Lyap,CloseRequestFcn,）;
holdon;
boxon;
timeplot=tstart+（tend-tstart）/10;
axis（tstarttimeplot-11）;
title（DynamicsofLyapunovexponents）;
xlabel（t）;
ylabel（Lyapunovexponents）;
Fig_Lyap_Axes=findobj（Fig_Lyap,Type,axes）;
n_expPlotLyapi=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,00rand）;
endITERLYAP=0;
%Mainloopcalculation_progress=1;
whilettend,tt=tend;
%SolutuionofextendedODEsystem%T,Y=feval（fcn_integrator,rhs_ext_fcn,tt+stept,y）;
whilecalculation_progress=1T,Y=integrator（DS
（1）.method_int,ode_lin,ttt,y,options,P,n,neq,n_exp）;
first_call=0;
ifcalculation_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（bufer_i）;
calculation_progress=1;
elsecalculation_progress=0;
if（calculation_progress=99）break;
elsecalculation_progress=1;
t=tt;
y=Y（size（Y,1）,:
%constructneworthonormalbasisbygram-schmidt%znorm
（1）=0.0;
forj=1:
n1znorm
（1）=znorm
（1）+y（n1+j）2;
znorm
（1）=sqrt（znorm
（1）;
n1y（n1+j）=y（n1+j）/znorm
（1）;
forj=2:
n_expfork=1:
（j-1）gsc（k）=0.0;
forl=1:
n1gsc（k）=gsc（k）+y（n1*j+l）*y（n1*k+l）;
fork=1:
n1forl=1:
（j-1）y（n1*j+k）=y（n1*j+k）-gsc（l）*y（n1*l+k）;
znorm（j）=0.0;
n1znorm（j）=znorm（j）+y（n1*j+k）2;
znorm（j）=sqrt（znorm（j）;
n1y（n1*j+k）=y（n1*j+k）/znorm（j）;
%updaterunningvectormagnitudes%fork=1:
n_expcum（k）=cum（k）+log（znorm（k）;
%normalizeexponent%rescale=0;
u1=1.e10;
u2=-1.e10;
n_explp（k）=cum（k）/（t-tstart）;
%PlotXd=get（PlotLyapk,Xdata）;
Yd=get（PlotLyapk,Ydata）;
iftimeplottu1=min（u1,min（Yd）;
u2=max（u2,max（Yd）;
Xd=Xdt;
Yd=Ydlp（k）;
set（PlotLyapk,Xdata,Xd,Ydata,Yd）;
iftimeplot0theneachioutp-thpointwillbeprint.%Outputparameters:
/www.math.rsu.ru/mexmat/kvm/matds/%lyapunov.misfreesoftware.lyapunov.misdistributedinthehopethatit%willbeuseful,butWITHOUTANYWARRANTY.%n=numberofnonlinearodes%n2=n*（n+1）=totalnumberofodes%n1=n;
n2=n1*（n1+1）;
%Numberofstepsnit=round（tend-tstart）/stept）;
cum=zeros（n1,1）;
n1y（n1+1）*i）=1.0;
%MainloopforITERLYAP=1:
nit%SolutuionofextendedODEsystemT,Y=unit（t,stept,y,d）;
t=t+stept;
n1forj=1:
n1y0（n1*i+j）=y（n1*j+i）;
n1znorm
（1）=znorm
（1）+y0（n1*j+1）2;
n1y0（n1*j+1）=y0（n1*j+1）/znorm
（1）;
n1fork=1:
n1gsc（k）=gsc（k）+y0（n1*l+j）*y0（n1*l+k）;
（j-1）y0（n1*k+j）=y0（n1*k+j）-gsc（l）*y0（n1*k+l）;
n1znorm（j）=znorm（j）+y0（n1*k+j）2;
n1y0（n1*k+j）=y0（n1*k+j）/znorm（j）;
n1cum（k）=cum（k）+log（znorm（k）;
%normalizeexponent%fork=1:
n1lp（k）=cum（k）/（t-tstart）;
%OutputmodificationifITERLYAP=1Lexp=lp;
Texp=t;
elseLexp=lp;
n1y（n1*j+i）=y0（n1*i+j）;
五、小数据量法计算Lyapunov指数的Matlab程序-（mex函数，超快）http:
/六、计算lyapunov指数的程序programlylorenzparameter（n=3,m=12,st=100）integer:
i,j,krealy（m）,z（n）,s（n）,yc（m）,h,y1（m）,a,b,r,f（m）,k1,k2,k3y
（1）=10.y
（2）=1.y（3）=0.a=10.b=8./3.r=28.t=0.h=0.01!
initialconditionsdoi=n+1,my（i）=0.enddodoi=1,ny（n+1）*i）=1.s（i）=0enddoopen（1,file=lorenz1.dat）open（2,file=lylorenz.dat）do100k=1,st!
stiterationscallrgkt（m,h,t,y,f,yc,y1）!
normarizevector123z=0.doi=1,ndoj=1,nz（i）=z（i）+y（n*j+i）*2enddoif（z（i）0.）z（i）=sqrt（z（i）doj=1,ny（n*j+i）=y（n*j+i）/z（i）enddoenddo!
generategsrcoefficientk1=0.k2=0.k3=0.doi=1,nk1=k1+y（3*i+1）*y（3*i+2）k2=k2+y（3*i+3）*y（3*i+2）k3=k3+y（3*i+1）*y（3*i+3）enddo!
conductnewvector2and3doi=1,ny（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）enddo!
generatenewvector2and3,normarizethemdoi=2,nz（i）=0.doj=2,nz（i）=z（i）+y（n*j+i）*2enddoif（z（i）0.）z（i）=sqrt（z（i）doj=2,ny（n*j+i）=y（n*j+i）/z（i）enddoenddo!
updatelyapunovexponentdoi=1,nif（z（i）0）s（i）=s（i）+log（z（i）enddo100continuedoi=1,ns（i）=s（i）/（1.*st*h*1000.）write（2,*）s（i）enddoend!
subroutinergkt（m,h,t,y,f,yc,y1）realy（m）,f（m）,y1（m）,yc（m）,a,b,rinteger:
i,jdoj=1,1000calldf（m,t,y,f）t=t+h/2.0doi=1,myc（i）=y（i）+h*f（i）/2.0y1（i）=y（i）+h*f（i）/6.0enddocalldf（m,t,yc,f）doi=1,myc（i）=y（i）+h*f（i）/2.0y1（i）=y1（i）+h*f（i）/3.0enddocalldf（m,t,yc,f）t=t+h/2.0doi=1,myc（i）=y（i）+h*f（i）y1（i）=y1（i）+h*f（i）/3.0enddocalldf（m,t,yc,f）doi=1,my（i）=y1（i）+h*f（i）/6.0enddoif（j500）write（1,*）t,y
（1）,y
（2）,y（3）enddoreturnend!
subroutinedf（m,t,y,f）realy（m）,a,b,r,f（m）commona,b,ra=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）doi=0,2f（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）enddoreturnend七、计算各种混沌系统李雅普洛夫指数的MATLAB源程序http:
/
升级会员 