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’)％横轴为 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 的时间历程xlabel(’y’）％横轴为 yylabel('dy/dt'）%纵轴为 dy/dtgrid;%显示网格线5.Poincare 映射的绘制对于不同的系统,Poincare 截面的选取方法也不同对于自治系统一般每过其对应线性系统的固有周期，截取一次对于非自治系统，一般每过其激励的周期，截取一次例程：duffing 方程的 poincare 映射function poincare(tstop）global omega;omega=1;T=2＊pi/omega；％线性系统的周期或激励的周期step=T/100;％定义步长为 T/100y0=［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;% 保留上一次的图形endxlabel（'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—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(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 方程...
