matlab编写的Lyapunov指数计算程序Word格式.docx
- 文档编号:21528292
- 上传时间:2023-01-31
- 格式:DOCX
- 页数:18
- 大小:20.30KB
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:
/
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 编写 Lyapunov 指数 计算 程序
![提示](https://static.bdocx.com/images/bang_tan.gif)