非線性分析第三次作業(yè)

1、.非線性分析第三次作業(yè)學(xué) 院（系）： 電子信息與電氣工程學(xué)部專 業(yè)： 信號與信息處理 學(xué) 生 姓 名： 代 菊 學(xué) 號： 11409013 任 課 教 師： 梅 建 琴 大連理工大學(xué)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 f

2、unction 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=

3、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=range y0=2 0 0; k=k+1; %除去前面60個周期的數(shù)據(jù)，并將最后的結(jié)果作為下一次積分的初值 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 %取出每一個周期內(nèi)的第一個解的最后一個值。 y0=Y(end,:); endendplot(range,bifdata,k.,markersize,5);運行上述程序，并對結(jié)果進行分析：以F為自變量，運動幅度為因變量的分岔圖如下：其混沌道路描述如下：(a) 當時，微分方系統(tǒng)為單周期運動，此時的相圖如下所示：(b)當時，單擺處于雙周期運動狀態(tài)，此時的

5、相圖如下所示：(c)當，單擺經(jīng)歷倍周期分岔，此時相圖如下所示(d) 當時，單擺進入混沌運動區(qū)，此時的系統(tǒng)相圖如下所示：由該相圖可知，系統(tǒng)在數(shù)個周期內(nèi)作運動。(e) 當時，系統(tǒng)恢復(fù)規(guī)則運動，此時相圖如下： 由上圖可知，系統(tǒng)從混沌中恢復(fù)，且做單周期運動。(2)wolf算法來計算李雅普諾夫指數(shù)的matlab程序如下：% 杜芬方程的參數(shù)%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;%Linearized syste

6、m.Jac=0 , 1, 0; 1-3*x2, -r, -F*sin(psi); 0, 0, 0;f(4:12)=Jac*Y; %變量方程% 計算李雅普諾夫指數(shù)譜函數(shù)%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);% Memory al

7、location.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:n1 fo

8、r 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) gsc(k)=0.0; for l=1:n1 gsc(k)=gsc(k)+y0(n1*l+j)*y0(

9、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 vector magnitudes. for k=1:n1 cum(k)=cum(k)+log(zno

10、rm(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;%主函數(shù)%clc;clear;global F;range=0.4:0.001:0.6;k=1;for F=range; Texp,Lexp=lyapun

11、ov2(); record(k)=Lexp(end,1); k=k+1;enda=1;運行上述方程得到李雅普諾夫指數(shù)隨的變化曲線如下： 由上圖可見，李雅普諾夫指數(shù)在處大于0，系統(tǒng)進入混沌狀態(tài)。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 to chaos. 2) Compute the Lyapunov exponent spectru

12、m 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） 當時，上述方程為線性方程，沒有分岔現(xiàn)象。2） 當時，求解上述方程，得到不動點：所以當時，x有實數(shù)解。即當時，Henon映射的不動點為：（，）和（，）。Matlab程序代碼如下：%畫出Henon映射在b=0.5時, a 0,1.4，步長=0.001之間變化時的分岔圖

13、%設(shè)定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.,markersize,10); hold onendxlim(0,a);MATLAB運行分岔圖結(jié)果如下：由分岔圖可知，當之后，系統(tǒng)

14、進入混沌狀態(tài)。2）求解李雅普諾夫指數(shù)%計算henon映射的lyapunov指數(shù)譜%備注：b=0.5時，得到NaN的非數(shù)值解，這里取參數(shù)：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.15.*x(L)2+y(L); y(L+1)=0.45*x(L); end if abs(x(end)2; D1=-2.3*x(en

15、d); JT=D1,D2;D3,D4;%Jaccob 矩陣 v,d=eig(JT); %特征向量和特征值 d=diag(d);% 取出特征值 L1=L1+log(abs(d(1); %第一李雅普諾夫指數(shù) L2=L2+log(abs(d(2); %第二李雅普諾夫指數(shù) 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),% Draw figure for Henon maping:figure; plot(Xp,Yp,k.

16、,markersize,2);運行上述程序，計算結(jié)果為：L1 =0.5837 L2 = -1.3822此時李雅普諾夫指數(shù)相圖：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+(xF.2)*(A.2)+B).0.5)/(B+(A*xF+(xF.2)*(A.2)+B).0.5).2) B/(

17、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 n100 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(n)=z(n)+x; endplot(z,-k.)程序運行：初始條件為：（0.32,0.32），不動點為（0.8732,0.8732）3. For the Rossler equation:l

18、 Investigate the chaotic behavior by plotting the phase diagrams and the Poincare sections as vary.答：求Rossler映射的不動點：假定是不動點，可以得到： 解方程組可得：。所以當時，系統(tǒng)有，有實數(shù)解，對應(yīng)的不動點分別為：和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方程相圖和龐加萊截面

19、圖%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: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(

20、籃色)隨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); 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