matlab動力学分析程序详解.doc

1.微分方程的定义 对于duffing方程，先将方程写作 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方程的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(t100);%截取稳态过程 y=y(count,:); n=length(y(:,1));%计算点的总数 for i=2:n-1 if y(i-1,1)+epsy(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方程的时间历程，相图，频谱和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