非线性分析作业3

文档来源为:从网络收集整理.word版本可编辑.欢迎下载支持.1文档来源为:从网络收集整理.word版本可编辑.非线性分析第三次作业学院（系）：电子信息与电气工程学部专业：信号与信息处理学生姓名：代菊学号：任课教师：梅建琴大连理工大学DalianUniversityofTechnology1.GiventheODE:1)PlotthebifurcationdiagramandphasediagramsasFvaries,andinvestigatetheroutestochaos.2)ComputetheLyapunovexponents,andplotthevalueasafunctionofF.答：1）令dxvdt，上述微分方程可以化为：Matlab程序代码如下：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%定义ODE方程%%%%%%%%%%%%%%%%%%%%%%%%%%%functiondx=ode(ignore,X)globalFwd;r=1;x=X(1);v=X(2);psi=X(3);dx=zeros(3,1);dx(1)=v;dx(2)=-r*v+x-x^3+F*cos(psi);dx(3)=wd;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%分岔图绘制程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionduffing_bifur_Fclear;clc;globalFwd;wd=1.2;文档来源为:从网络收集整理.word版本可编辑.欢迎下载支持.2文档来源为:从网络收集整理.word版本可编辑.range=0.4:0.0001:0.47;%F的范围%range=0.4:0.001:0.47;%F的范围period=2*pi/wd;k=0;YY1=[];rangelength=length(range);YY1=ones(rangelength,3000)*NaN;step=2*pi/300/wd;%步长,由于wd=1，周期即为2*pi，此步长为1周期取100个点。forF=rangey0=[200];k=k+1;%除去前面60个周期的数据，并将最后的结果作为下一次积分的初值tspan=0:step:60*period;[ignore,Y]=ode45(@duffing,tspan,y0);y0=Y(end,:);j=1;kkk=300;forii=20:59forpoint=(ii-1)*kkk+2:ii*kkkifY(point,1)Y(point-2,1)&&Y(point,1)Y(point+2,1)&&Y(point,1)Y(point-100,1)YY1(k,j)=Y(point,1);j=j+1;endend%取出每一个周期内的第一个解的最后一个值。y0=Y(end,:);endendplot(range,bifdata,'k.','markersize',5);运行上述程序，并对结果进行分析：以F为自变量，运动幅度为因变量的分岔图如下：其混沌道路描述如下：(a)当0.435F时，微分方系统为单周期运动，此时的相图如下所示：(b)当0.4350.455F时，单摆处于双周期运动状态，此时的相图如下所示：(c)当0.4550.4608F，单摆经历倍周期分岔，此时相图如下所示(d)当0.46080.463F时，单摆进入混沌运动区，此时的系统相图如下所示：由该相图可知，系统在数个周期内作运动。(e)当0.463F时，系统恢复规则运动，此时相图如下：由上图可知，系统从混沌中恢复，且做单周期运动。(2)wolf算法来计算李雅普诺夫指数的matlab程序如下：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%杜芬方程的参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%文档来源为:从网络收集整理.word版本可编辑.欢迎下载支持.3文档来源为:从网络收集整理.word版本可编辑.functionf=duff_ext(t,X);globalF;r=1;x=X(1);y=X(2);psi=X(3);dx=zeros(3,1);f(1)=y;f(2)=-r*y+x-x^3+F*cos(psi);f(3)=0.2;%Linearizedsystem.Jac=[0,1,0;1-3*x^2,-r,-F*sin(psi);0,0,0];f(4:12)=Jac*Y;%变量方程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算李雅普诺夫指数谱函数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[Texp,Lexp]=lyapunov2();globalF;n=3;rhs_ext_fcn=@duff_ext;fcn_integrator=@ode45;tstart=0;stept=0.5;tend=300;ystart=[111];ioutp=10;n1=n;n2=n1*(n1+1);%Numberofsteps.nit=round((tend-tstart)/stept);%Memoryallocation.y=zeros(n2,1);cum=zeros(n1,1);y0=y;gsc=cum;znorm=cum;%Initialvalues.y(1:n)=ystart(:);fori=1:n1y((n1+1)*i)=1.0;end;t=tstart;%Mainloop.forITERLYAP=1:nit%SolutuionofextendedODEsystem.[T,Y]=feval(fcn_integrator,rhs_ext_fcn,[tt+stept],y);t=t+stept;y=Y(size(Y,1),:);fori=1:n1forj=1:n1y0(n1*i+j)=y(n1*j+i);end;end;%ConstructneworthonormalbasisbyGram-Schmidt.znorm(1)=0.0;forj=1:n1znorm(1)=znorm(1)+y0(n1*j+1)^2;end;文档来源为:从网络收集整理.word版本可编辑.欢迎下载支持.4文档来源为:从网络收集整理.word版本可编辑.znorm(1)=sqrt(znorm(1));forj=1:n1y0(n1*j+1)=y0(n1*j+1)/znorm(1);end;forj=2:n1fork=1:(j-1)gsc(k)=0.0;forl=1:n1gsc(k)=gsc(k)+y0(n1*l+j)*y0(n1*l+k);end;end;fork=1:n1forl=1:(j-1)y0(n1*k+j)=y0(n1*k+j)-gsc(l)*y0(n1*k+l);end;end;znorm(j)=0.0;fork=1:n1znorm(j)=znorm(j)+y0(n1*k+j)^2;end;znorm(j)=sqrt(znorm(j));fork=1:n1y0(n1*k+j)=y0(n1*k+j)/znorm(j);end;end;%Updaterunningvectormagnitudes.fork=1:n1cum(k)=cum(k)+log(znorm(k));end;%Normalizeexponent.fork=1:n1lp(k)=cum(k)/(t-tstart);end;%Outputmodification.ifITERLYAP==1Lexp=lp;Texp=t;elseLexp=[Lexp;lp];Texp=[Texp;t];end;fori=1:n1forj=1:n1y(n1*j+i)=y0(n1*i+j);end;end;end;%%%%%%%%%%%%主函数%%%%%%%%%%%%clc;clear;globalF;range=0.4:0.001:0.6;k=1;forF=range;文档来源为:从网络收集整理.word版本可编辑.欢迎下载支持.5文档来源为:从网络收集整理.word版本可编辑.[Texp,Lexp]=lyapunov2();record(k)=Lexp(end,1);k=k+1;enda=1;运行上述方程得到李雅普诺夫指数随F的变化曲线如下：由上图可见，李雅普诺夫指数在0.46080.463F处大于0，系统进入混沌状态。2.ForHenonmap:1)Investigatethebifurcationdiagramforthehenonmapbyplottingthe_xnvaluesasafunctionofas0.5andgivetheanalysisoftheroutestochaos.2)ComputetheLyapunovexponentspectrumofthehenonmapwhen1.15and0.5.3)UsetheOGYalgorithmtostabilizethepointofperiodoneinthehenonmapwhen1.15and0.5.（1）求Henon映射的不动点：假定xy是不动点，可以得到：将二式带入一式可得：分两种情况讨论：1）当=0时，上述方程为线性方程，没有分岔现象。2）当0时，求解上述方程，得到不动点：所以当21)40（时，x有实数解。即当21)40（时，Henon映射的不动点为：（211)42（，211)4*2（）和（21-1)42（，211)4*2（）。Matlab程序代码如下：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%画出Henon映射在b=0.5时,a[0,1.4]，步长=0.001之间变化时的分岔图%设定x,y的初值为(0,0),%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%b=0.5;N=400;an=ones(1,N);xn=zeros(1,N);%holdon%boxon;文档来源为:从网络收集整理.word版本可编辑.欢迎下载支持.6文档来源为:从网络收集整理.word版本可编辑.x=0;y=0;fora=0:0.001:1.4fork=1:Nxm=x;ym=y;x=1-a*xm.*xm+ym;y=b*xm;endxn(1)=x;forn=2:Nxm=x;ym=y;x=1-a*xm.*xm+ym;y=b*xm;xn(n)=x;endplot(an*a,xn,'k.','markersize',10);holdonendxlim([0,a]);MATLAB运行分岔图结果如下：由分岔图可知，当a1之后，系统进入混沌状态。2）求解李雅普诺夫指数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算henon映射的lyapunov指数谱%备注：b=0.5时，得到NaN的非数值解，这里取参数：a=1.15,b=0.5%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clc;clear;closeall;M=10000;N=10000;D2=1;D3=0.45;D4=0;L1=0;L2=0;q=1;fork=1:M;x=zeros(1,N);y=zeros(1,N);x(1)=rand;y(1)=rand;forL=1:N-1;x(L+1)=1-1.15.*x(L)^2+y(L);文档来源为:从网络收集整理.word版本可编辑.欢迎下载支持.7文档来源为:从网络收集整理.word版本可编辑.y(L+1)=0.45*x(L);endifabs(x(end))2;D1=-2.3*x(end);JT=[D1,D2;D3,D4];%Jaccob矩阵[