《matlab动力学解析程序详解.doc》由会员分享，可在线阅读，更多相关《matlab动力学解析程序详解.doc（10页珍藏版）》请在三一办公上搜索。
1、1.微分方程的定义对于duffing方程，先将方程写作function dy=duffing(t,x)omega=1;%定义参数f1=x(2);f2=-omega2*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,tspa
2、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的时间历程xlabel(y)%横轴为yylabel(dy/dt)%纵轴为dy/dtgrid;%显示网格线5.Poincare映射的绘制对于不同
3、的系统，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,
4、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(t100);%截取稳态过程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 one
5、ndendxlabel(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=-omega2*x(1)-x(1)3;dy=f1;f2;%function duffsim(tstop)step=0.
6、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(460 500 -Inf Inf)%显示范围设置%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(en
7、d-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)%subplot(2,2,4)count=find(t100);%截取稳态过程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 on;endendxlabel(y);ylabel(dy/dt);8
8、.分岔图的绘制随变化的分岔图。function dy=duffing(t,x)global c;omega=1;%定义参数f1=x(2);f2=omega2*x(1)-x(1)3-0.3*x(2)+c*cos(1.2*t);dy=f1;f2;%clear;global c; %定义全局变量range=0.1:0.002:0.9;%定义参数变化范围k=0;YY=;%定义空数组for c=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);for i=2:n-1if 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;index(k)=j-1;endxlabel(c);ylabel(y);随F变化的分岔图F=0.20F=0.27F=0.275F=0.2875F=0.32F=0.36F=0.4F=0.652F=0.8
