非線(xiàn)性分析作業(yè)(共15頁(yè))

1、精選優(yōu)質(zhì)文檔-傾情為你奉上非線(xiàn)性分析第三次作業(yè)學(xué) 院（系）： 電子信息與電氣工程學(xué)部專(zhuān) 業(yè)： 信號(hào)與信息處理 學(xué) 生 姓 名： 代 菊 學(xué) 號(hào)： 任 課 教 師： 梅 建 琴 大連理工大學(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

2、a function 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的范圍peri

3、od=2*pi/wd;k=0;YY1=;rangelength=length(range);YY1=ones(rangelength,3000)*NaN;step=2*pi/300/wd; %步長(zhǎng),由于wd=1，周期即為2*pi，此步長(zhǎng)為1周期取100個(gè)點(diǎn)。for F=range y0=2 0 0; k=k+1; %除去前面60個(gè)周期的數(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=(i

4、i-1)*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 %取出每一個(gè)周期內(nèi)的第一個(gè)解的最后一個(gè)值。 y0=Y(end,:); endendplot(range,bifdata,'k.','markersize',5);運(yùn)行上述程序，并對(duì)結(jié)果進(jìn)行分析：以F為自變量，運(yùn)動(dòng)幅度為因變量的分岔圖如下：其混沌道路描述

5、如下：(a) 當(dāng)時(shí)，微分方系統(tǒng)為單周期運(yùn)動(dòng)，此時(shí)的相圖如下所示：(b)當(dāng)時(shí)，單擺處于雙周期運(yùn)動(dòng)狀態(tài)，此時(shí)的相圖如下所示：(c)當(dāng)，單擺經(jīng)歷倍周期分岔，此時(shí)相圖如下所示(d) 當(dāng)時(shí)，單擺進(jìn)入混沌運(yùn)動(dòng)區(qū)，此時(shí)的系統(tǒng)相圖如下所示：由該相圖可知，系統(tǒng)在數(shù)個(gè)周期內(nèi)作運(yùn)動(dòng)。(e) 當(dāng)時(shí)，系統(tǒng)恢復(fù)規(guī)則運(yùn)動(dòng)，此時(shí)相圖如下： 由上圖可知，系統(tǒng)從混沌中恢復(fù)，且做單周期運(yùn)動(dòng)。(2)wolf算法來(lái)計(jì)算李雅普諾夫指數(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

6、;f(2)=-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; %變量方程% 計(jì)算李雅普諾夫指數(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

7、 of 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

8、+stept,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

9、) gsc(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 runnin

10、g vector 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;%主函數(shù)%clc;clear;global F;

11、range=0.4:0.001:0.6;k=1;for F=range; Texp,Lexp=lyapunov2(); record(k)=Lexp(end,1); k=k+1;enda=1;運(yùn)行上述方程得到李雅普諾夫指數(shù)隨的變化曲線(xiàn)如下： 由上圖可見(jiàn)，李雅普諾夫指數(shù)在處大于0，系統(tǒng)進(jìn)入混沌狀態(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 rout

12、es to 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映射的不動(dòng)點(diǎn)：假定是不動(dòng)點(diǎn)，可以得到：將二式帶入一式可得：分兩種情況討論：1） 當(dāng)時(shí)，上述方程為線(xiàn)性方程，沒(méi)有分岔現(xiàn)象。2） 當(dāng)時(shí)，求解上述方程，得到不動(dòng)點(diǎn)：所以當(dāng)時(shí)，x有實(shí)數(shù)解。即當(dāng)時(shí)，Henon映射的不動(dòng)點(diǎn)為：（，）和（，）。Matl

13、ab程序代碼如下：%畫(huà)出Henon映射在b=0.5時(shí), a 0,1.4，步長(zhǎng)=0.001之間變化時(shí)的分岔圖%設(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.',&#

14、39;markersize',10); hold onendxlim(0,a);MATLAB運(yùn)行分岔圖結(jié)果如下：由分岔圖可知，當(dāng)之后，系統(tǒng)進(jìn)入混沌狀態(tài)。2）求解李雅普諾夫指數(shù)%計(jì)算henon映射的lyapunov指數(shù)譜%備注：b=0.5時(shí)，得到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+

15、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(end); 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 exponen

16、tL1=L1/(q-1),L2=L2/(q-1),% Draw figure for Henon maping:figure; plot(Xp,Yp,'k.','markersize',2);運(yùn)行上述程序，計(jì)算結(jié)果為：L1 =0.5837 L2 = -1.3822此時(shí)李雅普諾夫指數(shù)相圖：3）OGA算法控制周期1的一個(gè)點(diǎn)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=(

17、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/(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*xp

18、re;z(n)=z(n)+x; endplot(z,'-k.')程序運(yùn)行：初始條件為：（0.32,0.32），不動(dòng)點(diǎn)為（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映射的不動(dòng)點(diǎn)：假定是不動(dòng)點(diǎn)，可以得到： 解方程組可得：。所以當(dāng)時(shí)，系統(tǒng)有，有實(shí)數(shù)解，對(duì)應(yīng)的不動(dòng)點(diǎn)分別為：和matlab程序代碼如下% 定義rossler方程%f

19、unction 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逐漸變化時(shí)，繪制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個(gè)點(diǎn) x(

20、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

21、('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 subplot(2