matlab动力学解析程序详解

··
1.微分方程的定义
对于duffing 方程03
2=++x x x ω ，先将方程写作⎩⎨⎧--==31122
21x x x x x ω function dy=duffing(t,x)
omega=1;%定义参数
f1=x(2);
f2=-omega^2*x(1)-x(1)^3;
dy=[f1;f2];
2.微分方程的求解
function solve (tstop)
tstop=500;%定义时间长度
y0=[0.01;0];%定义初始条件
[t,y]=ode45('duffing',tstop,y0,[]);
function solve (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([460 500 -Inf Inf])%图形显示范围设置
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 方程03
2=++x x x ω
的poincare 映射 function poincare(tstop)
global omega;
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);
for i=5000:100:10000%稳态过程每个周期取一个点
plot(y(i,1),y(i,2),'b.');
hold on;% 保留上一次的图形
end
·· xlabel('y');ylabel('dy/dt');
Poincare 映射也可以通过取极值点得到
function poincare(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));%计算点的总数
for i=2:n-1
if y(i-1,1)+eps<y(i,1) && y(i,1)>y(i+1,1)+eps % 简单的取出局部最大值
plot(y(i,1),y(i,2),'.');
hold on
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 方程03
=++x x x
的时间历程，相图，频谱和poincare
映射。

function dy=duffing(t,x)
omega=1;%定义参数
f1=x(2);
f2=-omega^2*x(1)-x(1)^3;
dy=[f1;f2]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function duffsim(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([460 500 -Inf Inf])%显示范围设置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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));%计算点的总数
for i=2:n-1
if y(i-1,1)+eps<y(i,1) && y(i,1)>y(i+1,1)+eps % 简单的取出局部最大值
plot(y(i,1),y(i,2),'.');hold on;
end
end
xlabel('y');ylabel('dy/dt');
··
8.分岔图的绘制
随F变化的分岔图。

-
+
F
x
+
3.03=
x
x
x2.1
t
cos
function dy=duffing(t,x)
global c;
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;
global c; %定义全局变量
range=[0.1:0.002:0.9];%定义参数变化范围
··
k=0;
YY=[];%定义空数组
for c=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));
for i=2:n-1
if Y(i-1,1)+eps<Y(i,1) && Y(i,1)>Y(i+1,1)+eps % 简单的取出局部最大值。

YY(k,j)=Y(i,1);
j=j+1;
end
end
if j>1
plot(c,YY(k,[1:j-1]),'k.','markersize',3);
end
hold on;
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。