非线性分析作业3

1、非线性分析第三次作业学 院（系）： 电子信息与电气工程学部专 业： 信号与信息处理 学 生 姓 名： 代 菊 学 号： 11409013 任 课 教 师： 梅 建 琴 大连理工大学Dalian University of Technology1. Given the ODE: 1) Plot the bifurcation diagram and phase diagrams as F varies, and investigate the routes to chaos.2) Compute the Lyapunov exponents, and plot the value as a fu

2、nction of F.答：1）令，上述微分方程可以化为：Matlab 程序代码如下：% 定义ODE方程%function dx=ode(ignore,X)global F wd;r=1;x=X(1);v=X(2);psi=X(3);dx=zeros(3,1);dx(1)=v;dx(2)=-r*v+x-x3+F*cos(psi);dx(3)=wd;%分岔图绘制程序%function duffing_bifur_Fclear;clc;global F wd;wd=1.2;range=0.4:0.0001:0.47;%F的范围% range=0.4:0.001:0.47;%F的范围period=2

3、*pi/wd;k=0;YY1=;rangelength=length(range);YY1=ones(rangelength,3000)*NaN;step=2*pi/300/wd; %步长,由于wd=1，周期即为2*pi，此步长为1周期取100个点。for F=range y0=2 0 0; k=k+1; %除去前面60个周期的数据，并将最后的结果作为下一次积分的初值 tspan=0:step:60*period; ignore,Y=ode45(duffing,tspan,y0); y0=Y(end,:); j=1; kkk=300; for ii=20:59 for point=(ii-1)

4、*kkk+2:ii*kkk if Y(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; end end %取出每一个周期内的第一个解的最后一个值。 y0=Y(end,:); endendplot(range,bifdata,'k.','markersize',5);运行上述程序，并对结果进行分析：以F为自变量，运动幅度为因变量的分岔图如下：其混沌道路描述如下：(

5、a) 当时，微分方系统为单周期运动，此时的相图如下所示：(b)当时，单摆处于双周期运动状态，此时的相图如下所示：(c)当，单摆经历倍周期分岔，此时相图如下所示(d) 当时，单摆进入混沌运动区，此时的系统相图如下所示：由该相图可知，系统在数个周期内作运动。(e) 当时，系统恢复规则运动，此时相图如下： 由上图可知，系统从混沌中恢复，且做单周期运动。(2)wolf算法来计算李雅普诺夫指数的matlab程序如下：% 杜芬方程的参数%function f=duff_ext(t,X);global F;r=1;x=X(1);y=X(2);psi=X(3);dx=zeros(3,1);f(1)=y;f(2

6、)=-r*y+x-x3+F*cos(psi);f(3)=0.2;%Linearized system.Jac=0 , 1, 0; 1-3*x2, -r, -F*sin(psi); 0, 0, 0;f(4:12)=Jac*Y; %变量方程% 计算李雅普诺夫指数谱函数%function Texp,Lexp=lyapunov2();global F;n=3;rhs_ext_fcn=duff_ext;fcn_integrator=ode45;tstart=0;stept=0.5;tend=300;ystart=1 1 1;ioutp=10;n1=n; n2=n1*(n1+1);% Number of

7、steps.nit = round(tend-tstart)/stept);% Memory allocation.y=zeros(n2,1); cum=zeros(n1,1); y0=y;gsc=cum; znorm=cum;% Initial values.y(1:n)=ystart(:);for i=1:n1 y(n1+1)*i)=1.0; end;t=tstart;% Main loop.for ITERLYAP=1:nit% Solutuion of extended ODE system. T,Y = feval(fcn_integrator,rhs_ext_fcn,t t+ste

8、pt,y); t=t+stept; y=Y(size(Y,1),:); for i=1:n1 for j=1:n1 y0(n1*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; for j=2:n1 for k=1:(j-1) gs

9、c(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 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; end; znorm(j)=sqrt(znorm(j); for k=1:n1 y0(n1*k+j)=y0(n1*k+j)/znorm(j); end; end;% Update running ve

10、ctor magnitudes. for k=1:n1 cum(k)=cum(k)+log(znorm(k); end;% Normalize exponent. for k=1:n1 lp(k)=cum(k)/(t-tstart); end;% Output modification. if ITERLYAP=1 Lexp=lp; Texp=t; else Lexp=Lexp; lp; Texp=Texp; t; end; for i=1:n1 for j=1:n1 y(n1*j+i)=y0(n1*i+j); end; end;end;%主函数%clc;clear;global F;rang

11、e=0.4:0.001:0.6;k=1;for F=range; Texp,Lexp=lyapunov2(); record(k)=Lexp(end,1); k=k+1;enda=1;运行上述方程得到李雅普诺夫指数随的变化曲线如下： 由上图可见，李雅普诺夫指数在处大于0，系统进入混沌状态。2. For Henon map: 1) Investigate the bifurcation diagram for the henon map by plotting the values as a function of as and give the analysis of the routes t

12、o chaos. 2) Compute the Lyapunov exponent spectrum of the henon map when and .3) Use the OGY algorithm to stabilize the point of period one in the henon map when and .（1） 求Henon映射的不动点：假定是不动点，可以得到：将二式带入一式可得：分两种情况讨论：1） 当时，上述方程为线性方程，没有分岔现象。2） 当时，求解上述方程，得到不动点：所以当时，x有实数解。即当时，Henon映射的不动点为：（，）和（，）。Matlab程序

13、代码如下：%画出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);% hold on% box on;x=0;y=0;for a=0:0.001:1.4 for k=1:N xm=x; ym=y; x=1-a*xm.*xm+ym; y=b*xm; endxn(1)=x; for n=2:N xm=x; ym=y; x=1-a*xm.*xm+ym; y=b*xm; xn(n)=x; end plot(an*a,xn,'k.','m

14、arkersize',10); hold onendxlim(0,a);MATLAB运行分岔图结果如下：由分岔图可知，当之后，系统进入混沌状态。2）求解李雅普诺夫指数%计算henon映射的lyapunov指数谱%备注：b=0.5时，得到NaN的非数值解，这里取参数：a=1.15,b=0.5%clc;clear;close all;M=10000;N=10000;D2=1;D3=0.45;D4=0;L1=0;L2=0;q=1;for k=1:M; x=zeros(1,N); y=zeros(1,N); x(1)=rand; y(1)=rand; for L=1:N-1; x(L+1)=1

15、-1.15.*x(L)2+y(L); y(L+1)=0.45*x(L); end if abs(x(end)<2; D1=-2.3*x(end); JT=D1,D2;D3,D4;%Jaccob 矩阵 v,d=eig(JT); %特征向量和特征值 d=diag(d);% 取出特征值 L1=L1+log(abs(d(1); %第一李雅普诺夫指数 L2=L2+log(abs(d(2); %第二李雅普诺夫指数 Xp(q)=x(end); Yp(q)=y(end); q=q+1; end end% display the first and second Lyapunonv exponentL1=

16、L1/(q-1),L2=L2/(q-1),% Draw figure for Henon maping:figure; plot(Xp,Yp,'k.','markersize',2);运行上述程序，计算结果为：L1 =0.5837 L2 = -1.3822此时李雅普诺夫指数相图：3）OGA算法控制周期1的一个点Matlab代码：clearclcC=1.0;A=1.15;B=0.5;x=0.32; y=0.32;xF=(B-1+sqrt(1-B)2+4*A)/(2*A); %Fix-point g=(1-1).2+4*1.15).0.5*1;1; ju=(A*xF

17、+(xF.2)*(A2)+B).0.5)/B; hu=-B*(A*xF+(xF.2)*(A.2)+B).0.5)/(B+(A*xF+(xF.2)*(A.2)+B).0.5).2) B/(B+(A*xF+(xF.2)*(A.2)+B).0.5).2);z=zeros(1,140);p=zeros(1,140);for n=1:140 xpre=x; ypre=y; diag=x-xF;y-xF; if n<100 p(n)=0; else p(n)=(ju*hu*diag)/(ju-1)*(hu*g); end x=C+xpre*ypre-A*xpre.2+p(n); y=B*xpre;z

18、(n)=z(n)+x; endplot(z,'-k.')程序运行：初始条件为：（0.32,0.32），不动点为（0.8732,0.8732）3. For the Rossler equation:l Investigate the chaotic behavior by plotting the phase diagrams and the Poincare sections as vary.答：求Rossler映射的不动点：假定是不动点，可以得到： 解方程组可得：。所以当时，系统有，有实数解，对应的不动点分别为：和matlab程序代码如下% 定义rossler方程%funct

19、ion r=rossler(t,x)global a;global b;global c;r=-x(2)-x(3);x(1)+a*x(2);b+x(3)*(x(1)-c);% 绘制rossler方程相图和庞加莱截面图%clc;clear;global a; global b; global c;% a,b,c逐渐变化时，绘制rossler相图t0=0,200;f0=0,0,0;for c=2:0.02:4 for b=0:0.02:2 for a=0:0.01:0.1 t,x=ode45(rossler,t0,f0); t(1:length(t)-100)=; %取后面100个点 x(1:le

20、ngth(x)-100,:)=; % 绘制rossler相图 subplot(2,2,1); plot(t,x(:,1),'r',t,x(:,2),'g',t,x(:,3),'b'); title('x(红色),y(绿色),z(篮色)随t变化情况');xlabel('t'); subplot(2,2,2); plot3(x(:,1),x(:,2),x(:,3); title('rossler相图');xlabel('x');ylabel('y');zlabel(

21、9;z'); subplot(2,2,3); plot(x(:,1),x(:,2); title('x,y相图');xlabel('x');ylabel('y'); % 绘制rossler庞加莱截面图 z0=mean(x(:,3); % 选择z的均值所在的截面 j = 0; X1=; X2=; for k = 1:length(x(:,3)-1 dx = x(k,3)-z0; dy = x(k+1,3)-z0; if abs(dx)<1e-8 j = j+1; X1(j) = x(k,1); X2(j) = x(k,2); continue; end if sign(dx)*sign(dy)<0 j=j+1; Q=polyfit(x(k,3),x(k+1,3),x(k,1),x(k+1,1),1); X1(j)=polyval(Q,z0); Q=polyfit(x(k,3),x(k+1,3),x(k,2),x(k+1,2),1); X2(j)=polyval(Q,z0); end end sub