大三時候在跳蚤市場閑逛，從一位數學院的學長那里買了一些閑書，最近翻出來剛好有李榮華、劉播老師的《微分方程數值解法》和王仁宏老師的《數值逼近》，結合周善貴老師的《計算物理》課程，整理一下筆記。

本文整理常微分方程數值求解的歐拉法與龍格-庫塔法。

一般地，動力學系統的時間演化可以用常微分方程的初值問題來描述，例如設一維簡諧運動的回復力：

，有則運動方程： 。令 ，可以將二階微分方程轉化為一階微分方程組：

因此本文主要整理一階常微分方程初值問題的數值解法。

一階常微分方程初值問題

設

在區域 ： 上連續，對于一個給定的常微分方程 及初值 ，求解 。為了保證解 存在、唯一且連續依賴初值 ，要求 滿足Lipschitz條件：

存在常數L，使得

對所有 和 成立。

假設

總滿足上述條件。常用的近似解法有級數解法等近似解析方法，以及下文整理的數值方法：歐拉法與龍格-庫塔法。

歐拉法

將區間

作N等分，每一小區間長度 稱為步長， 稱為節點。根據初值 ，代入微分方程可直接解出 的導數值 。

推導

1、根據泰勒展開式：

略去二階小量，得：

以此類推，得到遞推公式：

2、數值積分推導

由

可得： ，使用左矩形積分得：

以此類推，可得到：

為了提高精度，可以使用梯形積分代替矩形積分，即：

以此類推，得到改進的歐拉法：

Python計算實例

以

為例，其精確解為 ，使用歐拉法求解的Python代碼如下：import math from matplotlib import pyplot as pltt_0 = 0 y_0 = 1 tau = 0.1 i = 1 solve = [] Euler = [] t = [] while i < 100:if i == 1:y_n = y_0t_n = t_0Euler.append(y_n)solve.append(math.exp(t_n))t.append(t_n)func = y_ny_n = y_n + tau * funct_n = t_n + taui += 1plt.plot(t, Euler, c='green', label=' Euler method') plt.plot(t, solve, c='red', label=' accuracy') plt.fill_between(t, solve, Euler, facecolor='blue', alpha=0.2) plt.title('Euler method', fontsize=19) plt.xlabel('t', fontsize=19) plt.ylabel('y', fontsize=19) plt.legend() plt.show()

作圖可以看到，當迭代步數較多后，歐拉法的結果逐漸落后于精確指數解的增長速度。下面分析歐拉法的誤差來源。

在

中，略去了高階小量 ，因此在每一步的遞推中，都有局部截斷誤差 ，其階為 。

在計算中，我們更關心精確解和數值解之間的誤差

，稱為整體誤差，其滿足

根據Lipschitz條件，可得：

由

且 ，可得：

局部截斷誤差

是 的二階量，設 ，得： ，整體誤差是 的一階量。同理可得，改進的歐拉法局部截斷誤差 是 的三階量 ，整體誤差是 的二階量 。

穩定性分析

如果計算的初值不能精確給定，例如存在測量、舍入誤差等，在計算過程中，每一步傳遞的誤差連續依賴于初始誤差，則稱算法穩定，否則該算法不穩定。

對于不同的初值

和 ，有

兩式相減，得：

根據Lipschitz條件，可得：

連續依賴于初始誤差，歐拉法穩定。同理，改進的歐拉法也穩定。

龍格-庫塔法

龍格庫塔法的主要思想：在

點的附近選取一些特定的點，然后把這些點的函數值進行線性組合，使用組合值代替泰勒展開中 點的導數值。

將

泰勒展開：

令

，根據多元函數求導法則有：

以此類推，可以得到：

同時，我們可以寫出泰勒展開的形式解:

其中：

通項為：

基本思路是，利用當前點的函數值

可以計算出 ，然后引入參數 和步長 可以計算出 ，之后使用 和步長 計算 ，以此類推，直到 。

現在把

展開：

將

代入得：

將

代入 可得：

與泰勒展開式

相比較，可知：

2個方程有3個未知數，因此有無窮多個解，可采用

， ， ，則：

令

， ，可以改寫為：

此即為二階龍格-庫塔法。

與上一節的歐拉法公式對比：

，因此二階龍格-庫塔法取參數 ， ， 時，即為改進的歐拉法。

Python計算實例

仍以

為例，其精確解為 ，使用二階龍格-庫塔法求解的Python代碼如下：import math from matplotlib import pyplot as pltt_0 = 0 y_0 = 1 z_0 = 1 tau = 0.1 i = 1 j = 1 solve = [] Euler = [] R_K = [] t = [] while i < 100:if i == 1:y_n = y_0t_n = t_0R_K.append(y_n)solve.append(math.exp(t_n))t.append(t_n)func_n = y_nfunc_m = y_n + tau * func_ny_n = y_n + 0.5 * tau * (func_n + func_m)t_n = t_n + taui += 1 t = [] while j < 100:if j == 1:z_n = z_0t_n = t_0Euler.append(z_n)t.append(t_n)func = z_nz_n = z_n + tau * funct_n = t_n + tauj += 1plt.scatter(t, R_K, marker='^', c='blue', s=70, label=' R-K method') plt.plot(t, Euler, c='green', label=' Euler method') plt.plot(t, solve, c='red', label=' accuracy') plt.fill_between(t, solve, Euler, facecolor='yellow', alpha=0.2) plt.title('Euler method & R-K method', fontsize=19) plt.xlabel('t', fontsize=19) plt.ylabel('y', fontsize=19) plt.legend() plt.show()

黃色部分表示數值解和精確解的偏離，可以看到，二階龍格-庫塔法（改進的歐拉法）精確度得到了很大的提升。

二階龍格-庫塔法中，泰勒展開到了

階，通過與泰勒展開系數進行對比，可以得到含3個未知數的2個方程。依次類推，如果泰勒展開到了 階，對比 可以得到 級 階龍格-庫塔法。常用經典四階龍格-庫塔法：

Reference：

1、周善貴，《計算物理》課程講義

2、李榮華，劉播，《微分方程數值解法》

3、王仁宏，《數值逼近》

創作挑戰賽新人創作獎勵來咯，堅持創作打卡瓜分現金大獎

總結

以上是生活随笔為你收集整理的四阶龙格库塔法的基本思想_数值常微分方程-欧拉法与龙格-库塔法的全部內容，希望文章能夠幫你解決所遇到的問題。

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