非线性分析报告作业3.docx
《非线性分析报告作业3.docx》由会员分享，可在线阅读，更多相关《非线性分析报告作业3.docx（23页珍藏版）》请在冰点文库上搜索。
非线性分析报告作业3
非线性分析第三次作业
学院（系）：
电子信息与电气工程学部
专业：
信号与信息处理
学生姓名：
代菊
学号：
11409013
任课教师：
梅建琴
大连理工大学
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）所示。
✧
升级会员 