1、1.微分方程的定義對于duffing方程，先將方程寫作function dy=duffing(t,x)omega=1;%定義參數(shù)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，繪制時只取穩(wěn)態(tài)部分。plot(t,y(:,1);%繪制y的時間歷程xlabel(t)%橫軸為tylabel(y)%縱軸為ygrid;%顯示網(wǎng)格線axis(460 500 -Inf Inf)%圖形顯示范圍設(shè)置4.相圖的繪制相圖的橫軸為y，縱軸為dy/dt，繪制時也只取穩(wěn)態(tài)部分。紅色部分表示只取最后1000個點。plot(y(end-1000:end,1),y(end-1000:end,2);%繪制y的時間歷程xlabel(y)%橫軸為yylabel(dy/dt)%縱軸為dy/dtgrid;%顯示網(wǎng)格線5.Poincare映射的繪制對于不同

3、的系統(tǒng)，Poincare截面的選取方法也不同對于自治系統(tǒng)一般每過其對應(yīng)線性系統(tǒng)的固有周期，截取一次對于非自治系統(tǒng)，一般每過其激勵的周期，截取一次例程：duffing方程的poincare映射function poincare(tstop)global omega;omega=1;T=2*pi/omega;%線性系統(tǒng)的周期或激勵的周期step=T/100;%定義步長為T/100y0=0.01;0;%初始條件tspan=0:step:100*T;%定義時間范圍t,y=ode45(duffing,tspan,y0);for i=5000:100:10000%穩(wěn)態(tài)過程每個周期取一個點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);%截取穩(wěn)態(tài)過程y=y(count,:);n=length(y(:,1);%計算點的總數(shù)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;%定義參數(shù)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;%顯示網(wǎng)格線axis(460 500 -Inf Inf)%顯示范圍設(shè)置%subplot(2,2,2)plot(y(end-1000:end,1),y(end-1000:end,2);%繪制y的時間歷程xlabel(y)%橫軸為yylabel(dy/dt)%縱軸為dy/dtgrid;%顯示網(wǎng)格線%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);%截取穩(wěn)態(tài)過程y=y(count,:);n=length(y(:,1);%計算點的總數(shù)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;%定義參數(shù)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;%定義參數(shù)變化范圍k=0;YY=;%定義空數(shù)組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