1、非线性分析第三次作业学 院（系）： 电子信息与电气工程学部专 业： 信号与信息处理 学 生 姓 名： 代 菊 学 号： 11409013 任 课 教 师： 梅 建 琴 大连理工大学Dalian University of Technology1. Given the ODE: 23dcos(0.2t)xxFt1) Plot the bifurcation diagram and phase diagrams as F varies, and investigate the routes to chaos.2) Compute the Lyapunov exponents, and plot t
2、he value as a function of F.答：1）令 ，上述微分方程可以化为：dxvt3cos(0.2t)dxvtFMatlab 程序代码如下：% 定义 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
3、的范围% 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 个点。for F=rangey0=2 0 0;k=k+1;%除去前面 60 个周期的数据，并将最后的结果作为下一次积分的初值tspan=0:step:60*period;ignore,Y=ode45(duffing,tspan,y0);y0=Y(end,:);j
4、=1;kkk=300;for ii=20:59for point=(ii-1)*kkk+2:ii*kkkif Y(point,1)Y(point-2,1)j=j+1;endend%取出每一个周期内的第一个解的最后一个值。y0=Y(end,:);endendplot(range,bifdata,k.,markersize,5);运行上述程序，并对结果进行分析：以 F 为自变量，运动幅度为因变量的分岔图如下：其混沌道路描述如下：(a) 当 时，微分方系统为单周期运动，此时的相图如下所示：0.435F(b)当 时，单摆处于双周期运动状态，此时的相图如下所示：0.435.F(c)当 ，单摆经历倍周期分
5、岔，此时相图如下所示0.45.608F(d) 当 时，单摆进入混沌运动区，此时的系统相图如下所示：0.4680.463F由该相图可知，系统在数个周期内作运动。(e) 当 时，系统恢复规则运动，此时相图如下：0.463F由上图可知，系统从混沌中恢复，且做单周期运动。(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)=-r*y+x-x3+F*cos(psi);f(3)=0.2;%Linear
6、ized 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 steps.nit = round(tend-tstart)/stept);% M
7、emory 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+stept,y); t=t+stept;y=Y(size(Y,1),:);for i=1
8、: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:n1for k=1:(j-1)gsc(k)=0.0;for l=1:n1 gsc(k)=gsc(k)+y0(n1*l+j)*y0(n1*
9、l+k); end;end;for k=1:n1for 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 vector magnitudes.for k=1:n1 cum(k)=cum(k)+log(znorm(k); end;% No
10、rmalize exponent.for k=1:n1 lp(k)=cum(k)/(t-tstart); end;% Output modification.if ITERLYAP=1Lexp=lp;Texp=t;elseLexp=Lexp; lp;Texp=Texp; t;end;for i=1:n1 for j=1:n1y(n1*j+i)=y0(n1*i+j);end;end;end;%主函数%clc;clear;global F;range=0.4:0.001:0.6;k=1;for F=range; Texp,Lexp=lyapunov2();record(k)=Lexp(end,1)
11、;k=k+1;enda=1;运行上述方程得到李雅普诺夫指数随 的变化曲线如下：F由上图可见，李雅普诺夫指数在 处大于 0，系统进入混沌状态。0.468.463F2. For Henon map: 211,nnnnxyx1) Investigate the bifurcation diagram for the henon map by plotting the values as _xna function of as and give the analysis of the routes to chaos. 0.52) Compute the Lyapunov exponent spectr
12、um of the henon map when and 1.5.053) Use the OGY algorithm to stabilize the point of period one in the henon map when and .1.（1 ） 求 Henon 映射的不动点：假定 是不动点，可以得到：xy21x将二式带入一式可得： 20x分两种情况讨论：1） 当 时，上述方程为线性方程，没有分岔现象。=02） 当 时，求解上述方程，得到不动点：21)4x（所以当 时，x 有实数解。21)40（即当 时，Henon 映射的不动点为：（ ， ）和21)4（ 21)4*（ ， ） 。
13、2-)（ 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);% hold on% box on;x=0;y=0;for a=0:0.001:1.4for k=1:Nxm=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;endplot(an*a,xn,k.,mark
14、ersize,10);hold onendxlim(0,a);MATLAB 运行分岔图结果如下：由分岔图可知，当 之后，系统进入混沌状态。a12）求解李雅普诺夫指数%计算 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-1.1
15、5.*x(L)2+y(L);y(L+1)=0.45*x(L);endif 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=L1/(q-1),L2=L2/(q-1)
16、,% 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+(xF.2)*(A2)+B).0.5)/B; hu=-B*(A*xF+
17、(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 n100p(n)=0; elsep(n)=(ju*hu*diag)/(ju-1)*(hu*g);endx=C+xpre*ypre-A*xpre.2+p(n); y=B*xpre;z(n)=z(n)+x; endplot(z,-k.)程序运行：初始条件为：（0.32,0.3
18、2） ，不动点为（0.8732,0.8732）3. For the Rossler equation: (xc)dxyztadzbt Investigate the chaotic behavior by plotting the phase diagrams and the Poincare sections as vary.,ac答：求 Rossler 映射的不动点：假定 是不动点，可以得到：xzy 0(xc)yzab解方程组可得： 。222444, ,cabcacxyza所以当 时，系统有， 有实数解，对应的不动点分别为：240cab,zx和2244(,)cabcab2224(,)cab
19、matlab 程序代码如下% 定义 rossler 方程%function 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:4for b=0:0.02:2for a=0:0.01:0.1t,x=ode45(rossler,t0,f0); t
20、(1:length(t)-100)=; %取后面 100 个点x(1:length(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(z); subplot(2,2,3); plot(x(:,1),x(:,2); titl
21、e(x,y 相图 );xlabel(x);ylabel(y);% 绘制 rossler 庞加莱截面图z0=mean(x(:,3); % 选择 z 的均值所在的截面 j = 0; X1=; X2=;for k = 1:length(x(:,3)-1dx = 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; endif sign(dx)*sign(dy)0j=j+1; Q=polyfit(x(k,3),x(k+1,3),x(k,1),x(k+1,1),1);
22、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); endendsubplot(2,2,4); plot(X1,X2,.);title(rossler 庞加莱截面图); xlabel(x,fontsize,14); ylabel(y,fontsize,14);pauseendendend运行上述程序，这里选取 b=2,c=4，a 在0,1之间变化时，画出 rossler 方程的相图和庞加莱截面图，并对结果进行分析。1)a=0.06, b=2,c=4 时，运行结果如下 2)
23、 a=0.12, b=2,c=4 时，运行结果如下3)a=0.18, b=2,c=4 时，运行结果如下 4)a=0.26, b=2,c=4 时，运行结果如下5)a=0.32, b=2,c=4 时，运行结果如下 6)a=0.4, b=2,c=4 时，运行结果如下7)a=0.46, b=2,c=4 时，运行结果如下 8)a=0.54, b=2,c=4 时，运行结果如下对上述 7 个图，可以做如下分析： 当 b，c 固定时，a 的值较小时，a=0.06 如图（1）所示， x，y，z 收敛于（0,0.5，-0.5） ，收敛速度很快； 随着 a 值变大,a=0.12 如图（2）所示，x，y ，z 是有收敛的趋势， ，收敛速度很慢，系统是一个极限环； a 继续增大时,a=0.18，如图（3）所示，x，y ，z 已经发散。但是 x，y，z 并不是发散于无穷大，而是周期性变化; 随着 a 的增大，x ，y，z 接近其极限环的速度加快， ，如图（4） 、 （5） 、 （6 ） 、 （7） ，进入倍周期运动 最后系统进入混沌状态（8）所示。
