matlab动力学解析程序详解.docx

1、matlab动力学解析程序详解1微分方程的定义对于duffing方程x 2x x3 -0，先将方程写作% =x2- 2 3X2 = 花fun cti on dy=duffi ng(t,x)omega=1;%定义参数f1=x(2);f2=-omegaA2*x(1)-x(1)A3;dy=f1;f2;2微分方程的求解fun cti on solve (tstop)tstop=500;%定义时间长度y0=0.01;0;%定义初始条件t,y=ode45(duffi ng,tstop,y0,);fun cti on solve (tstop)step=0.01;% 定义步长y0=ra nd(1,2);%随

2、机初始条件tspa n=0:step:500;% 定义时间范围t,y=ode45(duffi ng,tspa n,y0);3.时间历程的绘制时间历程横轴为t，纵轴为y,绘制时只取稳态部分。plot(t,y(:,1);% 绘制y的时间历程 xlabel(t)% 横轴为 tylabel(y)% 纵轴为 ygrid;%显示网格线axis(460 500 -Inf Inf)% 图形显示范围设置4.相图的绘制相图的横轴为y，纵轴为dy/dt，绘制时也只取稳态部分。红色部分 表示只取最后1000个点。plot(y( end-1000:end ,1),y( end-1000:end ,2);% 绘制 y 的

3、时间历 程xlabel(y)% 横轴为 yylabel(dy/dt)% 纵轴为 dy/dtgrid;%显示网格线5.Poi ncare映射的绘制对于不同的系统，Poincare截面的选取方法也不同对于自治系统一般每过其对应线性系统的固有周期，截取一次对于非自治系统，一般每过其激励的周期，截取一次2 3例程：duffing 方程 x x = 0的poincare映射fun cti on poin care(tstop)global omega; omega=1;T=2*pi/omega;%线性系统的周期或激励的周期step=T/100;% 定义步长为 T/100y0=0.01;0;% 初始条件t

4、spa n=0:step:100*T;% 定义时间范围t,y=ode45(duffi ng,tspa n,y0);for i=5000:100:10000% 稳态过程每个周期取一个点plot(y(i,1),y(i,2),b.);hold on;% 保留上一次的图形endxlabel(y);ylabel(dy/dt);Poincare映射也可以通过取极值点得到fun cti on poin care(tstop)y0=0.01;0;tspa n=0:0.01:500;t,y=ode45(duff ng,tspa n,yO);cou nt=fin d(t100);% 截取稳态过程y=y(co un

5、 t,:);n=le ngth(y(:,1);% 计算点的总数for i=2: n-1if y(i-1,1)+epsy(i+1,1)+eps %局部最大值plot(y(i,1),y(i,2),.);hold onendendxlabel(y);ylabel(dy/dt);6濒谱yy=fft(y(e nd-1000:e nd,1);N=le ngth(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.算例简单的取出poin careduffing 方程x x x

6、 0的时间历程，相图，频谱和 映射。fun cti on dy=duffi ng(t,x)omega=1;%定义参数f1=x(2);f2=-omegaA2*x(1)-x(1)A3;dy=f1;f2;% % fun cti on duffsim(tstop)step=0.01y0=0.1;0;tspa n=0:step:500;t,y=ode45(duffi ng,tspa n,yO);% subplot(2,2,1)plot(t,y(:,1);% 绘制y的时间历程xlabel(t)% 横轴为 tylabel(y)% 纵轴为 ygrid;%显示网格线axis(460 500 -Inf In f)

7、% 显示范围设置% subplot(2,2,2)plot(y(end-1000:end,1),y(end-1000:end,2);% 绘制 y 的时间历程xlabel(y)% 横轴为 yylabel(dy/dt)% 纵轴为 dy/dtgrid;%显示网格线% subplot(2,2,3)yy=fft(y(e nd-1000:e nd,1);N=le ngth(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)cou nt=fin d(t

8、100);% 截取稳态过程y=y(co un t,:);n=le ngth(y(:,1);% 计算点的总数for i=2: n-1if y(i-1,1)+epsy(i+1,1)+eps % 简单的取出局部最大值plot(y(i,1),y(i,2),.);hold on;endendxlabel(y);ylabel(dy/dt);D18.分岔图的绘制x 0.3x -x x3二F cos1.2t随F变化的分岔图。fun cti on dy=duffi ng(t,x)global c;omega=1;%定义参数f1=x(2);f2=omegaA2*x(1)-x(1)A3-0.3*x(2)+c*cos

9、(1.2*t);dy=f1；f2;% % clear;global c; % 定义全局变量ran ge=0.1:0.002:0.9;% 定义参数变化范围k=0;YY=；%定义空数组for c=ra ngey0=0.1;0;% 初始条件k=k+1;tspa n=0:0.01:400;t,Y=ode45(duffi ng,tspa n,yO);cou nt=fin d(t200);Y=Y(cou nt,:);j=1;n=le ngth(Y(:,1);for i=2: n-1简单的取出if Y(i-1,1)+epsY(i+1,1)+eps % 局部最大值。YY(k,j)=Y(i,1);j=j+1;endendif j1plot(c,YY(k,1:j-1),k.,markersize,3);endhold on;in dex(k)=j-1;endxlabel(c);ylabel(y);随F变化的分岔图F=0.20rfirna.3r.0-o *F=0.27F=0.275F=0.2875F=0.32-#县,F=0.36F=0.4F=0.652F=0.8