matlab动力学分析报告程序详解.docx
《matlab动力学分析报告程序详解.docx》由会员分享，可在线阅读，更多相关《matlab动力学分析报告程序详解.docx（10页珍藏版）》请在冰豆网上搜索。
matlab动力学分析报告程序详解
1.微分方程的定义
对于duffing方程
，先将方程写作
functiondy=duffing（t,x）
omega=1;%定义参数
f1=x
（2）;
f2=-omega^2*x
（1）-x
（1）^3;
dy=[f1;f2];
2.微分方程的求解
functionsolve（tstop）
tstop=500;%定义时间长度
y0=[0.01;0];%定义初始条件
[t,y]=ode45（'duffing',tstop,y0,[]）;
functionsolve（tstop）
step=0.01;%定义步长
y0=rand（1,2）;%随机初始条件
tspan=[0:
step:
500];%定义时间范围
[t,y]=ode45（'duffing',tspan,y0）;
3.时间历程的绘制
时间历程横轴为t，纵轴为y，绘制时只取稳态部分。
plot（t,y（:
1））;%绘制y的时间历程
xlabel（'t'）%横轴为t
ylabel（'y'）%纵轴为y
grid;%显示网格线
axis（[460500-InfInf]）%图形显示范围设置
4.相图的绘制
相图的横轴为y，纵轴为dy/dt，绘制时也只取稳态部分。
红色部分表示只取最后1000个点。
plot（y（end-1000:
end,1）,y（end-1000:
end,2））;%绘制y的时间历程
xlabel（'y'）%横轴为y
ylabel（'dy/dt'）%纵轴为dy/dt
grid;%显示网格线
5.Poincare映射的绘制
对于不同的系统，Poincare截面的选取方法也不同
对于自治系统一般每过其对应线性系统的固有周期，截取一次
对于非自治系统，一般每过其激励的周期，截取一次
例程：
duffing方程
的poincare映射
functionpoincare（tstop）
globalomega;
omega=1;
T=2*pi/omega;%线性系统的周期或激励的周期
step=T/100;%定义步长为T/100
y0=[0.01;0];%初始条件
tspan=[0:
step:
100*T];%定义时间范围
[t,y]=ode45（'duffing',tspan,y0）;
fori=5000:
100:
10000%稳态过程每个周期取一个点
plot（y（i,1）,y（i,2）,'b.'）;
holdon;%保留上一次的图形
end
xlabel（'y'）;ylabel（'dy/dt'）;
Poincare映射也可以通过取极值点得到
functionpoincare（tstop）
y0=[0.01;0];
tspan=[0:
0.01:
500];
[t,y]=ode45（'duffing',tspan,y0）;
count=find（t>100）;%截取稳态过程
y=y（count,:
）;
n=length（y（:
1））;%计算点的总数
fori=2:
n-1
ify（i-1,1）+epsy（i+1,1）+eps%简单的取出局部最大值
plot（y（i,1）,y（i,2）,'.'）;
holdon
end
end
xlabel（'y'）;ylabel（'dy/dt'）;
6.频谱
yy=fft（y（end-1000:
end,1））;
N=length（yy）;
power=abs（yy）;
freq=（1:
N-1）*1/step/N;
plot（freq（1:
N/2）,power（1:
N/2））;
xlabel（'f（y）'）
ylabel（'y'）
7.算例
duffing方程
的时间历程，相图，频谱和poincare映射。
functiondy=duffing（t,x）
omega=1;%定义参数
f1=x
（2）;
f2=-omega^2*x
（1）-x
（1）^3;
dy=[f1;f2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functionduffsim（tstop）
step=0.01
y0=[0.1;0];
tspan=[0:
step:
500];
[t,y]=ode45（'duffing',tspan,y0）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot（2,2,1）
plot（t,y（:
1））;%绘制y的时间历程
xlabel（'t'）%横轴为t
ylabel（'y'）%纵轴为y
grid;%显示网格线
axis（[460500-InfInf]）%显示范围设置
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot（2,2,2）
plot（y（end-1000:
end,1）,y（end-1000:
end,2））;%绘制y的时间历程
xlabel（'y'）%横轴为y
ylabel（'dy/dt'）%纵轴为dy/dt
grid;%显示网格线
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot（2,2,3）
yy=fft（y（end-1000:
end,1））;
N=length（yy）;
power=abs（yy）;
freq=（1:
N-1）*1/step/N;
plot（freq（1:
N/2）,power（1:
N/2））;
xlabel（'f（y）'）
ylabel（'y'）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot（2,2,4）
count=find（t>100）;%截取稳态过程
y=y（count,:
）;
n=length（y（:
1））;%计算点的总数
fori=2:
n-1
ify（i-1,1）+epsy（i+1,1）+eps%简单的取出局部最大值
plot（y（i,1）,y（i,2）,'.'）;holdon;
end
end
xlabel（'y'）;ylabel（'dy/dt'）;
8.分岔图的绘制
随
变化的分岔图。
functiondy=duffing（t,x）
globalc;
omega=1;%定义参数
f1=x
（2）;
f2=omega^2*x
（1）-x
（1）^3-0.3*x
（2）+c*cos（1.2*t）;
dy=[f1;f2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
globalc;%定义全局变量
range=[0.1:
0.002:
0.9];%定义参数变化范围
k=0;
YY=[];%定义空数组
forc=range
y0=[0.1;0];%初始条件
k=k+1;
tspan=[0:
0.01:
400];
[t,Y]=ode45（'duffing',tspan,y0）;
count=find（t>200）;
Y=Y（count,:
）;
j=1;
n=length（Y（:
1））;
fori=2:
n-1
ifY（i-1,1）+epsY（i+1,1）+eps%简单的取出局部最大值。
YY（k,j）=Y（i,1）;
j=j+1;
end
end
ifj>1
plot（c,YY（k,[1:
j-1]）,'k.','markersize',3）;
end
holdon;
index（k）=j-1;
end
xlabel（'c'）;
ylabel（'y'）;
随F变化的分岔图
F=0.20
F=0.27
F=0.275
F=0.2875
F=0.32
F=0.36
F=0.4
F=0.652
F=0.8
升级会员 