---

序

沒有對比就沒有傷害，本文先給出很多時候直接采用的矩形法，然后與四階龍格庫塔法做比較，著重說明四階龍格庫塔法。

---

一、矩形法

---

1.1 原理

設微分方程

$$\dot{y}=f(y)\text{\hspace{1em}(1.1)}$$

求$y$。

使用數值方法，離散化得每一步的增量

$$\Delta y=f(y)\Delta t\text{\hspace{1em}(1.2)}$$

易得

$$y_{n+1}=y_n+f(y_n)\Delta t\text{\hspace{1em}(1.3)}$$

實際上，這就是矩形法計算積分。當 $\Delta t\to 0$時，可以得出很高精度的$y_n$，但實際工程中未必能夠取很小的$\Delta t$。

---

1.2 例子

以
$$\begin{array}{rl}\dot{y}& =e^{?y}\\ y_0& =0\end{array}\text{\hspace{1em}(1.4)}$$

為例，時間取1~10s，分別取 $\Delta t=0.1,\Delta t=1$，查看不同精度下的運算結果。式(1.4)可求出解析解為$y=\ln (t+1)$，用于比較求解精度。

%% 不同步長下的矩形法比較dt1 = 0.1; % 步長1 dt2 = 1.0; % 步長2t1 = 0:dt1:10; % 時間1 t2 = 0:dt2:10; % 時間2y1 = ode_rect(t1, 0); % 精度1下計算結果 y2 = ode_rect(t2, 0); % 精度2下計算結果plot(t1,log(t1+1),t1,y1,t2,y2); legend('理論值', 'dt=0.1', 'dt=1'); grid on;grid minor;xlabel 't';ylabel 'y'%% 導數方程 function dy=f(y)dy = exp(-y); end%% 矩形法 function y = ode_rect(t, y0)N = length(t);y = zeros(N,1);y(1) = y0;for n = 1:N - 1dt = t(n+1) - t(n); % 計算步長 dty(n+1) = y(n) + f(y(n)) * dt; % 累加計算 yend end

可見，當步長為0.1時，矩形法的精度較高，但步長為1時，矩形法誤差大。

---

二、龍格庫塔法

---

2.1 $\dot{y}=f(y)$ 形式

經過多年潛心研究，龍格庫塔站在前人的肩膀上，發現了一種高精度的方法。那就是把式(1.3)的計算改為

$$\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)\\ k_1& =f(y_n)\\ k_2& =f(y_n+k_1\frac{h}{2})\\ k_3& =f(y_n+k_2\frac{h}{2})\\ k_4& =f(y_n+k_{3}h)\end{array}\text{\hspace{1em}(2.1)}$$

此處，采用習慣上的符號$h$，上面的例子中，$h=\Delta t$。

依舊對于式(1.4)，步長取1，分別使用矩形法和四階龍格庫塔法求解，結果如下

%% 矩形法與龍格庫塔法比較 dt = 1.0; t = 0:dt:10; y1 = ode_rect(t, 0); % 矩形法計算 y2 = ode_rk(t, 0); % 龍格庫塔法計算plot(0:0.01:10,log([0:0.01:10]+1),t,y1,t,y2); legend('理論值', '矩形法', '龍格庫塔法'); grid on;grid minor;xlabel 't';ylabel 'y'%% 導數方程 function dy=f(y)dy = exp(-y); end%% 矩形法 function y = ode_rect(t, y0)N = length(t);y = zeros(N,1);y(1) = y0;for n = 1:N - 1dt = t(n+1) - t(n); % 計算步長 dty(n+1) = y(n) + f(y(n)) * dt; % 累加計算 yend end%% 四階龍格庫塔法 function y = ode_rk(t, y0)N = length(t);y = zeros(N,1);y(1) = y0;for n = 1:N-1h = t(n+1) - t(n); % 步長（即時間間隔）k1 = f(y(n)); % k1k2 = f(y(n) + h/2*k1); % k2k3 = f(y(n) + h/2*k2); % k3k4 = f(y(n) + h*k3); % k4y(n+1) = y(n) + h/6 * (k1 + 2*k2 + 2*k3 + k4); % 累加計算yend end

可見，四階龍格庫塔法很好地接近真實值。

---

2.2 $\dot{y}=f(t)$ 形式

設

$$\dot{y}=f(t)\text{\hspace{1em}(2.2)}$$

從理論計算微分方程的角度，式(1.1) 和式(2.2)有著截然不同的求解方式。但是使用數值方法，只不過是把 $y$變為$t$。四階龍格庫塔法求解式(2.2)的方法如下

$$\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)\\ k_1& =f(t_n)\\ k_2& =f(t_n+\frac{h}{2})\\ k_3& =f(t_n+\frac{h}{2})\\ k_4& =f(t_n+h)\end{array}\text{\hspace{1em}(2.3)}$$

一個簡單的例子是

$$\begin{array}{rl}& \dot{y}=\sin (t)\\ & y_0=0\end{array}\text{\hspace{1em}(2.4)}$$

其解析解為$y=1?\cos (t)$。設步長分別為1,0.1，使用四階龍格庫塔法發求解如下

%% 龍格庫塔法步長差異比較 dt1 = 1; % 步長為1 dt2 = 0.1; % 步長為0.1 T = 10; % 總時間10st1 = 0:dt1:T; % 時間t1 y1 = ode_rk(t1, 0); % 解y1 t2 = 0:dt2:T; % 時間t2 y2 = ode_rk(t2, 0); % 解y2figure(1) plot(0:0.01:T, 1-cos(0:0.01:T));hold on % 解析解計算值 plot(t1,y1, '-o', t2,y2, '--');hold off % 數值解計算值 legend('理論值','dt=1', 'dt=0.1');title('龍格庫塔法'); grid on;grid minor;xlabel 't';ylabel 'y'%% 導數方程 function dy=f(t)dy = sin(t); end%% 四階龍格庫塔法 function y = ode_rk(t, y0)N = length(t);y = zeros(N,1);y(1) = y0;for n = 1:N-1h = t(n+1) - t(n); % 步長（即時間間隔）k1 = f(t(n)); % k1k2 = f(t(n) + h/2); % k2k3 = f(t(n) + h/2); % k3k4 = f(t(n) + h); % k4y(n+1) = y(n) + h/6 * (k1 + 2*k2 + 2*k3 + k4); % 累加計算yend end

從例子中可見，步長為1時，龍格庫塔法依舊得到精確的結果。

---

2.3 $\dot{y}=f(y,t)$ 形式

設
$$\dot{y}=f(y,t)$$

求解$y$。不過是多了一個變量，使用四階龍格庫塔法計算方法為

$$\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)\\ k_1& =f(y_n,t_n)\\ k_2& =f(y_n+k_1\frac{h}{2},t_n+\frac{h}{2})\\ k_3& =f(y_n+k_2\frac{h}{2},t_n+\frac{h}{2})\\ k_4& =f(y_n+k_{3}h,t_n+h)\end{array}\text{\hspace{1em}(2.5)}$$

更多變量，以此類推，不再贅述。

總結

以上是生活随笔為你收集整理的四阶龙格库塔法的计算例子的全部內容，希望文章能夠幫你解決所遇到的問題。

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