一、計算公式

對于形如以下的常微分方程：

采用四階龍格庫塔法的計算公式：

四階 龍格庫塔法精度為4，屬于單步遞推法。

單步遞推法的基本思想是從（x(i),y(i)）點出發，以某一斜率沿直線達到（x(i+1),y(i+1)）點。

二、實例計算

從上述定義可以看出，龍格庫塔實質上是求一階微分方程，但是如果將一階導看作變量，則二階導也不過是這個變量的一階導而已。

接下來就看實例吧：

對于下述二階方程：

1、基本思想：

令位移為 ??

q的一階導，即位移的一階導（速度）為 ?

q的二階導??

求解位移u(1)的龍格庫塔計算方法如下：

KK1=u(2); KK2=u(2)+h/2*KK1; KK3=u(2)+h/2*KK2; KK4=u(2)+h*KK3; u(1)=u(1)+h/6*(KK1+2*KK2+2*KK3+KK4);

求解速度u(2)的龍格庫塔計算方法如下：

K1=-2*eptheton*u(2)-u(1)+Fm+Fah*wh*wh*sin(wh*tao+fav_h); K2=-2*eptheton*(u(2)+h/2*K1)-(u(1)+h/2)+Fm+Fah*wh*wh*sin(wh*tao+fav_h); K3=-2*eptheton*(u(2)+h/2*K2)-(u(1)+h/2)+Fm+Fah*wh*wh*sin(wh*tao+fav_h); K4=-2*eptheton*(u(2)+h*K3)-(u(1)+h)+Fm+Fah*wh*wh*sin(wh*tao+fav_h); u(2)=u(2)+h/6*(K1+2*K2+2*K3+K4);

2、編程實現

參數設置：

h=0.001; % 步長 t0=0:h:400; % 總時長 w=5; ep=0.02; Fm=0.1; Fah=0.05; u(1)=0;u(2)=0; % 初值

總的程序實現

h=0.001; t0=0:h:400; w=5; ep=0.02; Fm=0.1; Fah=0.05; u(1)=0;u(2)=0; for i=1:length(t0) % 進行多次迭代tao=t0(i);u=RK(u,tao,h,ep,w,Fm,Fah); Result(i,:)=u; % 將每次迭代的位移和速度保存 end figure(1) subplot(2,1,1) plot(t0,Result(:,1)) % 繪制位移圖 xlabel('Time') ylabel('displacement') subplot(2,1,2) plot(t0,Result(:,2)) % 繪制速度圖 xlabel('Time') ylabel('velocity')function u=RK(u,tao,h,ep,w,Fm,Fah) KK1=u(2); KK2=u(2)+h/2*KK1; KK3=u(2)+h/2*KK2; KK4=u(2)+h*KK3; u(1)=u(1)+h/6*(KK1+2*KK2+2*KK3+KK4); K1=-2*ep*u(2)-u(1)+Fm+Fah*cos(w*tao); K2=-2*ep*(u(2)+h/2*K1)-u(1)-h/2+Fm+Fah*cos(w*tao); K3=-2*ep*(u(2)+h/2*K2)-u(1)-h/2+Fm+Fah*cos(w*tao); K4=-2*ep*(u(2)+h*K3)-u(1)-h+Fm+Fah*cos(w*tao); u(2)=u(2)+h/6*(K1+2*K2+2*K3+K4); end

結果圖如下所示

值得特別注意的是

u=RK(u,tao,h,ep,w,Fm,Fah);

輸入的u與輸出的u一定要符號一致，從而確保下一次迭代的初始值是上一次的值。確保迭代能一直進行下去。

三、輔助驗證

接下來用MATLAB自帶的ode45函數來進行驗證。

之前已經寫過ode45函數的用法，在此不再進行介紹。

源程序如下：

t0=0:0.001:400; w=5; ep=0.02; Fm=0.1; Fah=0.05; y0=[0 0]; [t,u]=ode45(@(t,u) odefun(t,u,w,ep,Fm,Fah),t0,y0); figure(1) subplot(2,1,1) plot(t,u(:,1)) xlabel('Time') ylabel('displacement') subplot(2,1,2) plot(t,u(:,2)) xlabel('Time') ylabel('velocity') function du=odefun(t,u,w,ep,Fm,Fah) du=[u(2);-2*ep*u(2)-u(1)+Fm+Fah*cos(w*t)]; end

運算結果圖如下所示

兩中方法計算的結果是一樣的。

創作不易，如有幫助，記得點個贊吶。

總結

以上是生活随笔為你收集整理的基于MATLAB实现四阶龙格库塔法求解一、二阶微分方程实例的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。