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