matlab动力学解析程序详解

1、微分方程的定义对于duffing方程x'+与2x+x，=0,先将方程写作x1=x2-23*=-财*一*functiondy=duffing(t,x)omega=1;%定义参数f1=x(2);f2=-omegaA2*x(1)-x(1)A3;dy=f1;f2;微分方程的求解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:

2、500;%定义时间范围t,y=ode45('duffing',tspan,y0);时间历程的绘制时间历程横轴为t,纵轴为y,绘制时只取稳态部分。plot(t,y(:,1);%绘制y的时间历程xlabel('t')%横轴为tylabel('y')%纵轴为ygrid;%显示网格线axis(460500-InfInf)%图形显示范围设置相图的绘制相图的横轴为y,纵轴为dy/dt,绘制时也只取稳态部分。红色部分表示只取最后1000个点。plot(y(end-1000:end,1),y(end-1000:end,2);%绘制y的时间历程xlabel(

3、9;y')%横轴为yylabel('dy/dt')%纵轴为dy/dtgrid;%显示网格线Poincare映射的绘制对于不同的系统，Poincare截面的选取方法也不同对于自治系统一般每过其对应线性系统的固有周期，截取一次对于非自治系统，一般每过其激励的周期，截取一次例程：duffing方程x+切2x+x3=0的poincare映射functionpoincare(tstop)globalomega;omega=1;T=2*pi/omega;%线性系统的周期或激励的周期step=T/100;%定义步长为T/100y0=0.01;0;%初始条件tspan=0:step:1

4、00*T;%定义时间范围t,y=ode45('duffing',tspan,y0);fori=5000:100:10000%稳态过程每个周期取一个点plot(y(i,1),y(i,2),'b.');holdon;%保留上一次的图形endxlabel('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&

5、gt;100);%截取稳态过程y=y(count,:);n=length(y(:,1);%计算点的总数fori=2:n-1ify(i-1,1)+eps<y(i,1)&&y(i,1)>y(i+1,1)+eps%简单的取出局部最大值plot(y(i,1),y(i,2),'.');holdonendendxlabel('y');ylabel('dy/dt');频谱yy=fft(y(end-1000:end,1);N=length(yy);power=abs(yy);freq=(1:N-1)*1/step/N;plot(fre

6、q(1:N/2),power(1:N/2);xlabel('f(y)')ylabel('y')算例duffing方程x+x+x，=0的时间历程，相图，频谱和poincare映射。functiondy=duffing(t,x)omega=1;%定义参数f1=x(2);f2=-omegaA2*x(1)-x(1)A3;dy=f1;f2;%functionduffsim(tstop)step=0.01y0=0.1;0;tspan=0:step:500;t,y=ode45('duffing',tspan,y0);%subplot(2,2,1)plot(t,

7、y(:,1);%绘制y的时间历程xlabel('t')%横轴为tylabel('y')%纵轴为ygrid;%显示网格线axis(460500-InfInf)%显示范围设置%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(end-1000:end,1);N=length(yy);power=abs(yy);

8、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-1ify(i-1,1)+eps<y(i,1)&&y(i,1)>y(i+1,1)+eps%简单的取出局部最大值plot(y(i,1),y(i,2),'.');holdon;endendxlabel

9、('y');ylabel('dy/dt');D1分岔图的绘制x+0.3xx+x3=Fcosl.2t随F变化的分岔图。functiondy=duffing(t,x)globalc;omega=1;%定义参数f1=x(2);f2=omegaA2*x(1)-x(1)A3-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=rangey0=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-1ifY(i-1,1)+eps<Y(i,1)&&Y(i,1)>Y(i+1,1)+eps%简单的取出局部最大值。YY(k,j)=Y(i,1);j=j+1;endendifj>1plot(c,YY(k,1:j-1)