matlab动力学解析程序详解

1、1 .微分方程的定义对于duffing 方程x' +82x+x3 =0 ,先将方程写作x1 X2,23X2 = -3 xi -xifunction dy=duffing(t,x)omega=1;%E 义参数f1=x(2);f2=-omegaA2*x(1)-x(1)A3;dy=fi;f2;2 .微分方程的求解function solve (tstop) tstop=500;%定义时间长度 y0=0.01;0;%定义初始条件 t,y=ode45('duffing',tstop,y0,);function solve (tstop)step=0.01;% 定义步长y0=ran

2、d(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:e

3、nd ,1),y( end-1000:end ,2);% 绘制 y 的时间历 程 xlabel('y')% 横轴为 y ylabel('dy/dt')% 纵轴为 dy/dt grid;%显示网格线5 .Poincare映射的绘制对于不同白系统，Poincare截面的选取方法也不同 对于自治系统一般每过其对应线性系统的固有周期，截取一次 对于非自治系统，一般每过其激励的周期，截取一次例程：duffing 方程 x +82x+x3 =0 的 poincare 映射function poincare(tstop)global omega;omega=1;T=2*pi/

4、omega;%线性系统的周期或激励的周期step=T/100;% 定义步长为 T/100y0=0.01;0;% 初始条件tspan=0:step:100*T;% 定义时间范围 t,y=ode45('duffing',tspan,y0);fo门=5000:100:10000%稳态过程每个周期取一个点plot(y(i,1),y(i,2),'b.');hold on;% 保留上一次的图形endxlabel('y');ylabel('dy/dt');Poincare映射也可以通过取极值点得到function poincare(tstop)

5、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)+eps<y(i,1) && y(i,1)>y(i+1,1)+eps %简单的取出局部最大值plot(y(i,1),y(i,2),'.');hold onendendxlabel('y');ylabel('dy/dt')

6、;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 方程x +x +x3 =0的时间历程，相图，频谱和 poincare 映射。function dy=duffing(t,x)omega=1;%定义参数f1=x(2);f2=-omegaA2*x(1)-x(1)A3;dy=fi;f2;%function duffsim(ts

7、top)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(460 500 -Inf Inf)%显示范围设置%subplot(2,2,2)plot(y(end-1000:end,1),y(end-1000:end,2);% 绘制 y 的时间历程xlabel('y')%横轴为 yy

8、label('dy/dt')% 纵轴为 dy/dtgrid;%显示网格线%subplot(2,2,3)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')%subplot(2,2,4)count=find(t>100);%截取稳态过程y=y(count,:);n=length(y(:,1);% 计算点的总数for i=2:n-1if y(

9、i-1,1)+eps<y(i,1) && y(i,1)>y(i+1,1)+eps %简单的取出局部最大值plot(y(i,1),y(i,2),'.');hold on;endendxlabel('y');ylabel('dy/dt');D1QI8 .分岔图的绘制x +0.3x x +x3 = F cosl.2t 随 F 变化的分岔图。function dy=duffing(t,x)global c;omega=1;%E 义参数f1=x(2);f2=omegaA2*x(1)-x(1)A3-0.3*x(2)+c*cos(1

10、.2*t);dy=fi;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(t>200);Y=Y(count,:);j=1;n=length(Y(:,1);for i=2:n-1if Y(i-1,1)+eps<Y(i,1) && Y(i,1)>Y(i+1,1)+eps %简单的取出局部最大值。YY(k,j尸丫(i,1);j=j+1；endendif j>1plo