matlab动力学分析报告程序详解.docx
- 文档编号:10865137
- 上传时间:2023-02-23
- 格式:DOCX
- 页数:10
- 大小:149.81KB
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)+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)+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)+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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 动力学 分析 报告 程序 详解