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