matlab动力学分析程序详解

··

1.微分方程的定义

对于duffing 方程03

2

=++x x x

ω ，先将方程写作⎩⎨⎧

--==3

1122

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)+epsy(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)+epsy(i+1,1)+eps % 简单的取出局部最大值

plot(y(i,1),y(i,2),'.');hold on;

end

end

··