非线性分析作业3.docx

上传人:b****5 文档编号:7254605 上传时间:2023-01-22 格式:DOCX 页数:20 大小:310.20KB
下载 相关 举报
非线性分析作业3.docx_第1页
第1页 / 共20页
非线性分析作业3.docx_第2页
第2页 / 共20页
非线性分析作业3.docx_第3页
第3页 / 共20页
非线性分析作业3.docx_第4页
第4页 / 共20页
非线性分析作业3.docx_第5页
第5页 / 共20页
点击查看更多>>
下载资源
资源描述

非线性分析作业3.docx

《非线性分析作业3.docx》由会员分享,可在线阅读,更多相关《非线性分析作业3.docx(20页珍藏版)》请在冰豆网上搜索。

非线性分析作业3.docx

非线性分析作业3

 

非线性分析第三次作业

 

学院(系):

电子信息与电气工程学部

专业:

信号与信息处理

学生姓名:

代菊

学号:

********

***********

 

大连理工大学

DalianUniversityofTechnology

 

1.GiventheODE:

1)PlotthebifurcationdiagramandphasediagramsasFvaries,andinvestigatetheroutestochaos.

2)ComputetheLyapunovexponents,andplotthevalueasafunctionofF.

答:

1)令

,上述微分方程可以化为:

Matlab程序代码如下:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%定义ODE方程

%%%%%%%%%%%%%%%%%%%%%%%%%%%

functiondx=ode(ignore,X)

globalFwd;

r=1;

x=X

(1);

v=X

(2);

psi=X(3);

dx=zeros(3,1);

dx

(1)=v;

dx

(2)=-r*v+x-x^3+F*cos(psi);

dx(3)=wd;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%分岔图绘制程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

functionduffing_bifur_F

clear;

clc;

globalFwd;

wd=1.2;

range=0.4:

0.0001:

0.47;%F的范围

%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个点。

forF=range

y0=[200];

k=k+1;

%除去前面60个周期的数据,并将最后的结果作为下一次积分的初值

tspan=0:

step:

60*period;

[ignore,Y]=ode45(@duffing,tspan,y0);

y0=Y(end,:

);

j=1;

kkk=300;

forii=20:

59

forpoint=(ii-1)*kkk+2:

ii*kkk

ifY(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,:

);

end

end

plot(range,bifdata,'k.','markersize',5);

运行上述程序,并对结果进行分析:

以F为自变量,运动幅度为因变量的分岔图如下:

其混沌道路描述如下:

(a)当

时,微分方系统为单周期运动,此时的相图如下所示:

(b)当

时,单摆处于双周期运动状态,此时的相图如下所示:

(c)当

,单摆经历倍周期分岔,此时相图如下所示

(d)当

时,单摆进入混沌运动区,此时的系统相图如下所示:

由该相图可知,系统在数个周期内作运动。

(e)当

时,系统恢复规则运动,此时相图如下:

由上图可知,系统从混沌中恢复,且做单周期运动。

(2)wolf算法来计算李雅普诺夫指数的matlab程序如下:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%杜芬方程的参数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

functionf=duff_ext(t,X);

globalF;

r=1;

x=X

(1);

y=X

(2);

psi=X(3);

dx=zeros(3,1);

f

(1)=y;

f

(2)=-r*y+x-x^3+F*cos(psi);

f(3)=0.2;

%Linearizedsystem.

Jac=[0,1,0;

1-3*x^2,-r,-F*sin(psi);

0,0,0];

f(4:

12)=Jac*Y;%变量方程

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%计算李雅普诺夫指数谱函数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function[Texp,Lexp]=lyapunov2();

globalF;

n=3;rhs_ext_fcn=@duff_ext;fcn_integrator=@ode45;

tstart=0;stept=0.5;tend=300;

ystart=[111];ioutp=10;

n1=n;n2=n1*(n1+1);

%Numberofsteps.

nit=round((tend-tstart)/stept);

%Memoryallocation.

y=zeros(n2,1);cum=zeros(n1,1);y0=y;

gsc=cum;znorm=cum;

%Initialvalues.

y(1:

n)=ystart(:

);

fori=1:

n1y((n1+1)*i)=1.0;end;

t=tstart;

%Mainloop.

forITERLYAP=1:

nit

%SolutuionofextendedODEsystem.

[T,Y]=feval(fcn_integrator,rhs_ext_fcn,[tt+stept],y);

t=t+stept;

y=Y(size(Y,1),:

);

fori=1:

n1

forj=1:

n1y0(n1*i+j)=y(n1*j+i);end;

end;

%ConstructneworthonormalbasisbyGram-Schmidt.

znorm

(1)=0.0;

forj=1:

n1znorm

(1)=znorm

(1)+y0(n1*j+1)^2;end;

znorm

(1)=sqrt(znorm

(1));

forj=1:

n1y0(n1*j+1)=y0(n1*j+1)/znorm

(1);end;

forj=2:

n1

fork=1:

(j-1)

gsc(k)=0.0;

forl=1:

n1gsc(k)=gsc(k)+y0(n1*l+j)*y0(n1*l+k);end;

end;

fork=1:

n1

forl=1:

(j-1)

y0(n1*k+j)=y0(n1*k+j)-gsc(l)*y0(n1*k+l);

end;

end;

znorm(j)=0.0;

fork=1:

n1znorm(j)=znorm(j)+y0(n1*k+j)^2;end;

znorm(j)=sqrt(znorm(j));

fork=1:

n1y0(n1*k+j)=y0(n1*k+j)/znorm(j);end;

end;

%Updaterunningvectormagnitudes.

fork=1:

n1cum(k)=cum(k)+log(znorm(k));end;

%Normalizeexponent.

fork=1:

n1

lp(k)=cum(k)/(t-tstart);

end;

%Outputmodification.

ifITERLYAP==1

Lexp=lp;

Texp=t;

else

Lexp=[Lexp;lp];

Texp=[Texp;t];

end;

fori=1:

n1

forj=1:

n1

y(n1*j+i)=y0(n1*i+j);

end;

end;

end;

%%%%%%%%%%%%主函数%%%%%%%%%%%%

clc;

clear;

globalF;

range=0.4:

0.001:

0.6;

k=1;

forF=range;

[Texp,Lexp]=lyapunov2();

record(k)=Lexp(end,1);

k=k+1;

end

a=1;

运行上述方程得到李雅普诺夫指数随

的变化曲线如下:

由上图可见,李雅普诺夫指数在

处大于0,系统进入混沌状态。

2.ForHenonmap:

1)Investigatethebifurcationdiagramforthehenonmapbyplottingthe

valuesasafunctionof

as

andgivetheanalysisoftheroutestochaos.

2)ComputetheLyapunovexponentspectrumofthehenonmapwhen

and

.

3)UsetheOGYalgorithmtostabilizethepointofperiodoneinthehenonmapwhen

and

.

(1)求Henon映射的不动点:

假定

是不动点,可以得到:

将二式带入一式可得:

分两种情况讨论:

1)当

时,上述方程为线性方程,没有分岔现象。

2)当

时,求解上述方程,得到不动点:

所以当

时,x有实数解。

即当

时,Henon映射的不动点为:

)和

)。

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);

%holdon

%boxon;

x=0;

y=0;

fora=0:

0.001:

1.4

fork=1:

N

xm=x;

ym=y;

x=1-a*xm.*xm+ym;

y=b*xm;

end

xn

(1)=x;

forn=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);

holdon

end

xlim([0,a]);

MATLAB运行分岔图结果如下:

由分岔图可知,当

之后,系统进入混沌状态。

2)求解李雅普诺夫指数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算henon映射的lyapunov指数谱

%备注:

b=0.5时,得到NaN的非数值解,这里取参数:

a=1.15,b=0.5

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;clear;closeall;

M=10000;

N=10000;

D2=1;

D3=0.45;

D4=0;

L1=0;

L2=0;

q=1;

fork=1:

M;

x=zeros(1,N);

y=zeros(1,N);

x

(1)=rand;

y

(1)=rand;

forL=1:

N-1;

x(L+1)=1-1.15.*x(L)^2+y(L);

y(L+1)=0.45*x(L);

end

ifabs(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

%displaythefirstandsecondLyapunonvexponent

L1=L1/(q-1),

L2=L2/(q-1),

%DrawfigureforHenonmaping:

figure;

plot(Xp,Yp,'k.','markersize',2);

运行上述程序,计算结果为:

L1=0.5837L2=-1.3822

此时李雅普诺夫指数相图:

3)OGA算法控制周期1的一个点

Matlab代码:

clear

clc

C=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)*(A^2)+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);

forn=1:

140

xpre=x;

ypre=y;

diag=[x-xF;y-xF];

ifn<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(n)=z(n)+x;

end

plot(z,'-k.')

程序运行:

初始条件为:

(0.32,0.32),不动点为(0.8732,0.8732)

 

3.FortheRosslerequation:

●InvestigatethechaoticbehaviorbyplottingthephasediagramsandthePoincaresectionsas

vary.

答:

求Rossler映射的不动点:

假定

是不动点,可以得到:

解方程组可得:

所以当

时,系统有,

有实数解,对应的不动点分别为:

matlab程序代码如下

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%定义rossler方程

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

functionr=rossler(t,x)

globala;

globalb;

globalc;

r=[-x

(2)-x(3);x

(1)+a*x

(2);b+x(3)*(x

(1)-c)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%绘制rossler方程相图和庞加莱截面图

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;clear;

globala;globalb;globalc;

%%a,b,c逐渐变化时,绘制rossler相图

t0=[0,200];f0=[0,0,0];

forc=2:

0.02:

4

forb=0:

0.02:

2

fora=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(篮色)随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=[];

fork=1:

length(x(:

3))-1

dx=x(k,3)-z0;

dy=x(k+1,3)-z0;

ifabs(dx)<1e-8

j=j+1;

X1(j)=x(k,1);

X2(j)=x(k,2);

continue;

end

ifsign(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,2,4);plot(X1,X2,'.');

title('rossler庞加莱截面图');

xlabel('x','fontsize',14);ylabel('y','fontsize',14);

pause

end

end

end

运行上述程序,这里选取b=2,c=4,a在[0,1]之间变化时,画出rossler方程的相图和庞加莱截面图,并对结果进行分析。

1)a=0.06,b=2,c=4时,运行结果如下2)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)所示。

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 农林牧渔 > 农学

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1