matlab动力学分析程序详解

11.微分方程的定义对于duffing方程032xxx，先将方程写作3112221xxxxxfunctiondy=duffing(t,x)omega=1;%定义参数f1=x(2);f2=-omega^2*x(1)-x(1)^3;dy=[f1;f2];2.微分方程的求解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:500];%定义时间范围[t,y]=ode45('duffing',tspan,y0);3.时间历程的绘制时间历程横轴为t，纵轴为y，绘制时只取稳态部分。plot(t,y(:,1));%绘制y的时间历程xlabel('t')%横轴为tylabel('y')%纵轴为ygrid;%显示网格线2axis([460500-InfInf])%图形显示范围设置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方程032xxx的poincare映射functionpoincare(tstop)globalomega;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);fori=5000:100:10000%稳态过程每个周期取一个点plot(y(i,1),y(i,2),'b.');holdon;%保留上一次的图形endxlabel('y');ylabel('dy/dt');3Poincare映射也可以通过取极值点得到functionpoincare(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));%计算点的总数fori=2:n-1ify(i-1,1)+epsy(i,1)&&y(i,1)y(i+1,1)+eps%简单的取出局部最大值plot(y(i,1),y(i,2),'.');holdonendendxlabel('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xxx的时间历程，相图，频谱和poincare映射。4functiondy=duffing(t,x)omega=1;%定义参数f1=x(2);f2=-omega^2*x(1)-x(1)^3;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,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));5N=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(t100);%截取稳态过程y=y(count,:);n=length(y(:,1));%计算点的总数fori=2:n-1ify(i-1,1)+epsy(i,1)&&y(i,1)y(i+1,1)+eps%简单的取出局部最大值plot(y(i,1),y(i,2),'.');holdon;endendxlabel('y');ylabel('dy/dt');68.分岔图的绘制tFxxxx2.1cos3.03随F变化的分岔图。functiondy=duffing(t,x)globalc;omega=1;%定义参数f1=x(2);f2=omega^2*x(1)-x(1)^3-0.3*x(2)+c*cos(1.2*t);dy=[f1;f2];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;globalc;%定义全局变量range=[0.1:0.002:0.9];%定义参数变化范围k=0;7YY=[];%定义空数组forc=rangey0=[0.1;0];%初始条件k=k+1;tspan=[0:0.01:400];[t,Y]=ode45('duffing',tspan,y0);count=find(t200);Y=Y(count,:);j=1;n=length(Y(:,1));fori=2:n-1ifY(i-1,1)+epsY(i,1)&&Y(i,1)Y(i+1,1)+eps%简单的取出局部最大值。YY(k,j)=Y(i,1);j=j+1;endendifj1plot(c,YY(k,[1:j-1]),'k.','markersize',3);endholdon;index(k)=j-1;endxlabel('c');ylabel('y');8随F变化的分岔图F=0.209F=0.27F=0.275F=0.287510F=0.32F=0.36F=0.411F=0.652F=0.8