常微分方程數值方法

常微分方程數值方法是用以尋找常微分方程（ODE）解的數值近似值的方法。其使用也稱作「數值積分」，不過「數值積分」主要是指積分的計算。

微分方程${\displaystyle y'=y,\ y(0)=1}$的數值積分求解圖。

  藍：歐拉法

  綠：中點法

  紅：精確解：${\displaystyle y=e^{t}}$

步長${\displaystyle h=1.0}$。

${\displaystyle h=0.25}$時的同一圖。${\displaystyle h\to 0}$時，中點法比歐拉法收斂得更快。

很多微分方程無法精確求解。但在工程學等領域的實際應用中，通常只需得到數值近似解。本文介紹的算法可用於計算這種近似值，另一種方法是用微積分技術得到解的級數展開表達。

常微分方程出現於物理學、化學、生物學、經濟學等學科中。^[1]此外，偏微分方程數值方法中的一部分將偏微分方程轉為常微分方程，然後可用本文所述方法求解。

問題

一階微分方程是有如下形式的初值問題（IVP）：^{[2]:533–655}

${\displaystyle y'(t)=f(t,y(t)),\qquad y(t_{0})=y_{0},}$

1

其中f是函數${\displaystyle f:[t_{0},\infty )\times \mathbb {R} ^{d}\to \mathbb {R} ^{d}}$，初始條件${\displaystyle y_{0}\in \mathbb {R} ^{d}}$是已知向量。「一階」是說方程中只出現了y的一階導。

在不失高階系統一般性的前提下，限制(l)式為一階微分方程，因為高階ODE可通過引入額外變量轉換為更大的一階方程組。例如，二階方程${\displaystyle y''=-y}$可重寫為2個一階方程：${\displaystyle y'=z,\ z'=-y}$。

本文介紹IVP的數值方法，並指出邊值問題（BVP）需要一套不同的工具：需要在多個點上定義解y的值或成分，因此要用不同方法求解BVP，如打靶法（及其變體），或差分、^[3]伽遼金法、^[4]配置法（英語：Collocation method）之類全局方法，都適用於此類問題。

Picard–Lindelöf定理指出，只要f是利普希茨連續的，就有唯一解。

方法

解一階IVP的數值方法可分為兩大類：^[5]線性多步法（英語：Linear multistep method）與龍格-庫塔法。還可進一步劃為顯式或隱式，例如隱式線性多步法包括亞當斯-莫爾頓法（Adams-Moulton methods）與向後微分公式（英語：Backward differentiation formula）（BDF），隱式龍格-庫塔法^{[6]:204–215}包括對角隱式龍格-庫塔法（DIRK）、^[7][8]單對角隱式龍格-庫塔法（SDIRK）^[9]與基於高斯求積^[10]的高斯–拉道法^[11]等等。線性多步法中的顯式方法有亞當斯-巴什福思法，布徹表（Butcher tableau）為下對角的龍格-庫塔法都是顯式方法。根據經驗，剛性微分方程需要用隱式方案，而非剛性問題則可用顯示方案更有效地求解。

所謂一般線性方法（英語：General linear methods）（GLM）是上述兩大類方法的概括。^[12]

歐拉法

更多資訊：歐拉法

從曲線上任意一點出發，沿與曲線相切的直線移動一小段距離，就能得到曲線上鄰近點的近似值。

從微分方程(1)開始，可用差分近似代替導數y′：

${\displaystyle y'(t)\approx {\frac {y(t+h)-y(t)}{h}},}$

2

重新排列後得到以下公式

${\displaystyle y(t+h)\approx y(t)+hy'(t)}$

利用(1)可得

${\displaystyle y(t+h)\approx y(t)+hf(t,y(t)).}$

3

此式的應用通常如下：擇步長h，構造序列${\displaystyle t_{0},t_{1}=t_{0}+h,t_{2}=t_{0}+2h,...}$記精確解${\displaystyle y(t_{n})}$的數值估計值為${\displaystyle y_{n}}$。受(3)啟發，可用下面的遞歸方法計算估計值：

${\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n}).}$

4

這就是（前向）歐拉法，是萊昂哈德·歐拉（1768）描述的方法。

歐拉法是顯式方法的一個例子，這是說新值${\displaystyle y_{n+1}}$是根據已知值（如${\displaystyle y_{n}}$）定義的。

反向歐拉法

更多資訊：反向歐拉法

若不用(2)，而用近似值

${\displaystyle y'(t)\approx {\frac {y(t)-y(t-h)}{h}},}$

5

則得到反向歐拉法：

${\displaystyle y_{n+1}=y_{n}+hf(t_{n+1},y_{n+1}).}$

6

反向歐拉法是隱式方法，這是說需要求解一個方程才能得到新值${\displaystyle y_{n+1}}$。通常用定點迭代或牛頓-拉弗森法（的某種修改版）實現之。

隱式方法求解這方程比顯示方法直接代入要花更多時間，選擇方法時必須考慮這一成本。隱式方法（如(6)）的優點是求解剛性方程時通常更穩定，便可以使用更大的步長h。

一階指數積分法

更多資訊：指數積分（英語：Exponential integrator）

指數積分（英語：Exponential integrator）描述了一大類積分器，近來得到了長足發展。^[13]它們至少可以追溯到20世紀60年代。

此處不用(1)，假設微分方程形式為

${\displaystyle y'(t)=-A\,y+{\mathcal {N}}(y),}$

7

或已被局部線性化，圍繞背景狀態產生線性項${\displaystyle -Ay}$與非線性項${\displaystyle {\mathcal {N}}(y)}$。

將(7)乘以${\textstyle e^{At}}$，並在時間區間${\displaystyle [t_{n},t_{n+1}=t_{n}+h]}$內精確積分，便構造得到了指數積分：

${\displaystyle y_{n+1}=e^{-Ah}y_{n}+\int _{0}^{h}e^{-(h-\tau )A}{\mathcal {N}}\left(y\left(t_{n}+\tau \right)\right)\,d\tau .}$

此積分方程是精確的，但並沒有定義積分。

使${\displaystyle {\mathcal {N}}(y(t_{n}+\tau ))}$在整個區間內不變，可實現一階指數積分：

${\displaystyle y_{n+1}=e^{-Ah}y_{n}+A^{-1}(1-e^{-Ah}){\mathcal {N}}(y(t_{n}))\ .}$

8

推廣

歐拉法往往不夠精確。更準確地說，只有一階（下面將介紹階的概念），這就促使數學家尋找高階方法。

一種方法是，用${\displaystyle y_{n}}$以及更多之前的值確定${\displaystyle y_{n+1}}$，所謂多步法便實現了這種想法。最簡單的可能是蛙跳積分法，其是二階方法，（粗略地說）依賴於2個時間值。

幾乎所有實用的多步法都屬於線性多步法（英語：Linear multistep method），形式為 $${\displaystyle {\begin{aligned}&{}\alpha _{k}y_{n+k}+\alpha _{k-1}y_{n+k-1}+\cdots +\alpha _{0}y_{n}\\&{}\quad =h\left[\beta _{k}f(t_{n+k},y_{n+k})+\beta _{k-1}f(t_{n+k-1},y_{n+k-1})+\cdots +\beta _{0}f(t_{n},y_{n})\right].\end{aligned}}}$$

另一種方法是在區間${\displaystyle [t_{n},t_{n+1}]}$內取更多點。這產生了得名於卡爾·龍格與馬丁·威廉·庫塔的龍格-庫塔法，其中一種4階方法尤為流行。

高級特徵

要用這些方法求解ODE，需要的不僅是時間步長公式。始終相同的步長效率不高，於是開發了可變步長方法。通常，步長的選擇應使每步的（局部）誤差低於某容差水平，這意味着方法還要計算誤差指標（error indicator），即對局部誤差的估計。

這一思想的延伸是在不同階方法之間動態選擇（稱作可變階數方法）。基於理查德森外推法^[14]的Bulirsch–Stoer算法之類方法^[15][16]通常用於構建各種不同階的方法。

其他理想特徵還有：

• 輸出稠密：不僅對點${\displaystyle t_{0},\ t_{1},\ t_{2},\ \dots }$，還對整個積分區間進行低成本數值逼近；
• 事件定位：找到事件發生的位置，例如某函數取0的時間。通常用求根算法；
• 支持並行計算；
• 用於時間積分時，具有時間可逆性。

其他方法

很多方法不在討論的框架內，如

• 多階導數法，不僅使用f，還用其導數。這類方法包括Hermite–Obreschkoff法、費爾貝格法（英語：Runge–Kutta–Fehlberg method）與遞歸計算解y的泰勒級數係數的Parker–Sochacki法（英語：Parker–Sochacki method）^[17]、Bychkov–Scherbakov法等。
• 二階常微分方程組法。前面說過，所有高階常微分方程都可化為形式如(1)的一階常微分方程組，這固然沒錯，但未必是最佳方法。尼斯特倫法（英語：Nyström method）可直接用於二階方程。
• 幾何積分法^[18][19]專門針對幾類特殊的常微分方程（如解哈密頓力學方程的辛積分（英語：Symplectic integrator）法），這些方法在數值求解時會尊重這些類的底結構或幾何結構。
• 量化狀態系統方法（英語：Quantized state systems method）是基於狀態量化思想的一系列ODE積分方法，在模擬頻繁不連續的稀疏系統時非常有效。

實時並行方法

有些IVP要求積分具有很高的時間解像度和/或很長的時間區間，傳統的序列時間步進法無法實時計算（如數值天氣預報、等離子體建模與分子動力學中的IVP）。針對這些問題，人們開發了實時並行（PinT）法以便用並行計算縮短運行時間。

早期的PinT法（最早提出於20世紀60年代）^[20]最初被研究人員忽視，因為其所需的並行計算架構尚未普及。2000年代初，隨着算力的提高，靈活易用的PinT算法——Parareal算法（英語：Parareal）重新吸引了興趣，它適用於求解各種IVP。百億億次級運算（英語：Exascale computing）的出現使PinT算法獲得更多關注，並啟動了能用於世界上最強大的超級計算機的算法開發。截至2023年，最流行的方法有Parareal、PFASST、ParaDiag、MGRIT等，^[21]

分析

數值分析不僅包括數值方法的設計，還包括其分析。分析中的3個核心概念是：

• 收斂性：方法是否逼近解；
• 階數：近似解的程度；
• 數值穩定性：誤差是否能得到抑制。^[22]

收斂性

主條目：序列、極限 (數學)和數列極限

若數值解隨着步長h趨近於0而逼近精確解，則稱此數值方法具有收斂性（convergent）。更確切地說，要求對利普希茨函數f與每個${\displaystyle t^{*}>0}$，ODE (1)

${\displaystyle \lim _{h\to 0^{+}}\max _{n=0,1,\dots ,\lfloor t^{*}/h\rfloor }\left\|y_{n,h}-y(t_{n})\right\|=0.}$

上述所有方法都收斂。

一致性與階數

更多資訊：截斷誤差（英語：Truncation error (numerical integration)）

設數值方法

${\displaystyle y_{n+k}=\Psi (t_{n+k};y_{n},y_{n+1},\dots ,y_{n+k-1};h).\,}$

方法的局部（截斷）誤差定義為迭代一步產生的誤差。即，假設之前的迭代無誤差，則是此方法結果與精確解之間的差

${\displaystyle \delta _{n+k}^{h}=\Psi \left(t_{n+k};y(t_{n}),y(t_{n+1}),\dots ,y(t_{n+k-1});h\right)-y(t_{n+k}).}$

若

${\displaystyle \lim _{h\to 0}{\frac {\delta _{n+k}^{h}}{h}}=0.}$

則稱此方法一致（consistent）。若

${\displaystyle \delta _{n+k}^{h}=O(h^{p+1})\quad {\mbox{as }}h\to 0.}$

則稱此方法階數為p。因此，階數不為0的方法是一致的。上文介紹的兩種歐拉法(4、6)都是1階的，因此是一致的。實踐中使用的大多數方法階數更高。一致性是收斂的必要條件^{[來源請求]}，但不是充分條件；方法要收斂，必須同時具有一致性與零穩性（zero-stable）。

一個相關概念是全局（截斷）誤差，即達到固定時間t所需所有步驟中持續存在的誤差。明確地說，t時刻的全局誤差是${\displaystyle y_{N}-y(t)}$，其中${\displaystyle N=(t-t_{0})/h}$。第p階一步法是收斂的，其全局誤差是${\displaystyle O(h^{p})}$。對多階方法，這一說法不一定成立。

穩定性與剛性

更多資訊：剛性方程

對某些微分方程，歐拉法、顯式龍格-庫塔法、多步法（如亞當斯-巴什福思法）之類標準方法會表現出解的不穩定性，其他方法則可能產生穩定的解。方程中的這種「困難行為」（本身不一定複雜）稱作「剛性」（stiffness），通常是由於底問題中存在不同時間尺度造成的。^[23]例如，機械系統中的碰撞（如碰撞振子中的）發生的時間尺度通常比物體運動的時間尺度小得多，這種差異使狀態參數曲線變得非常陡峭。

剛性問題在化學動力學、控制論、固體力學、天氣預報、生物學、等離子體物理、電子學等領域中無處不在。克服剛性的一種方法是將微分方程推廣到微分包含式，從而允許非光滑性並建立其模型。^[24][25]

歷史

下面是該領域一些重要進展的時間線。^[26][27]

• 1768 - 萊昂哈德·歐拉發現歐拉法。
• 1824 - 奧古斯丁-路易·柯西證明了歐拉法的收斂性。此證明中，柯西使用了隱性歐拉法。
• 1855 - Francis Bashforth在一封信中首次提到約翰·柯西·亞當斯的多步法。
• 1895 - 卡爾·龍格發現第一個龍格-庫塔法。
• 1901 - 馬丁·威廉·庫塔描述了現在流行的4階龍格-庫塔法。
• 1910 - 路易斯·弗萊·理查德森公佈了他的外推法——理查德森外推法。
• 1952 - Charles F. Curtiss與Joseph Oakland Hirschfelder提出剛性方程概念。
• 1963 - Germund Dahlquist引入積分法的A穩定性。

二階一維邊值問題的數值解法

邊值問題（BVPs）通常要離散化，得到近似相等的矩陣問題再數值求解。^[28]最常用的一維BVP數值求解方法稱作有限差分法。^[3]這種方法用點值的線性組合構造描述函數導數的有限差分係數，例如一階導數的二階中心差分近似為：

${\displaystyle {\frac {u_{i+1}-u_{i-1}}{2h}}=u'(x_{i})+{\mathcal {O}}(h^{2}),}$

二階導數的二階中心差分近似為：

${\displaystyle {\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}=u''(x_{i})+{\mathcal {O}}(h^{2}).}$

兩式中，${\displaystyle h=x_{i}-x_{i-1}}$是離散域上相鄰x值間的距離。這樣就構建了線性系統，然後可用標準矩陣法求解。例如，待解方程：

${\displaystyle {\begin{aligned}&{}{\frac {d^{2}u}{dx^{2}}}-u=0,\\&{}u(0)=0,\\&{}u(1)=1.\end{aligned}}}$

下一步是將問題離散化，用線性導數近似，如

${\displaystyle u''_{i}={\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}}$

並解所得線性方程組。將有如下方程

${\displaystyle {\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}-u_{i}=0,\quad \forall i={1,2,3,...,n-1}.}$

初看之下，這個方程組似乎有困難，因為方程中沒有不與變量相乘的項，但事實上這是錯誤的。i = 1或n − 1時，有一項涉及邊值${\displaystyle u(0)=u_{0}}$、${\displaystyle u(1)=u_{n}}$，由於兩值已知，可以簡單地代入，就得到了具有非平凡解的非齊次線性方程組。

相關條目

• CFL條件（英語：Courant–Friedrichs–Lewy condition）
• 能量漂移（英語：Energy drift）
• 一般線性方法（英語：General linear methods）
• Modelica語言及OpenModelica（英語：OpenModelica）軟件

註釋

1. Chicone, C. (2006). Ordinary differential equations with applications (Vol. 34). Springer Science & Business Media.
2. Bradie (2006)
3. LeVeque, R. J. (2007). Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems (Vol. 98). SIAM.
4. Slimane Adjerid and Mahboub Baccouch (2010) Galerkin methods. Scholarpedia, 5(10):10056.
5. Griffiths, D. F., & Higham, D. J. (2010). Numerical methods for ordinary differential equations: initial value problems. Springer Science & Business Media.
6. Hairer, Nørsett & Wanner (1993) harvtxt模板錯誤: 無指向目標: CITEREFHairerNørsettWanner1993 (幫助)
7. Alexander, R. (1977). Diagonally implicit Runge–Kutta methods for stiff ODE’s. SIAM Journal on Numerical Analysis, 14(6), 1006-1021.
8. Cash, J. R. (1979). Diagonally implicit Runge-Kutta formulae with error estimates. IMA Journal of Applied Mathematics, 24(3), 293-301.
9. Ferracina, L., & Spijker, M. N. (2008). Strong stability of singly-diagonally-implicit Runge–Kutta methods. Applied Numerical Mathematics, 58(11), 1675-1686.
10. Weisstein, Eric W. "Gaussian Quadrature." From MathWorld--A Wolfram Web Resource. https://mathworld.wolfram.com/GaussianQuadrature.html （頁面存檔備份，存於互聯網檔案館）
11. Everhart, E. (1985). An efficient integrator that uses Gauss-Radau spacings. In International Astronomical Union Colloquium (Vol. 83, pp. 185-202). Cambridge University Press.
12. Butcher, J. C. (1987). The numerical analysis of ordinary differential equations: Runge-Kutta and general linear methods. Wiley-Interscience.
13. Hochbruck (2010，第209–286頁) harvtxt模板錯誤: 無指向目標: CITEREFHochbruck2010 (幫助) This is a modern and extensive review paper for exponential integrators
14. Brezinski, C., & Zaglia, M. R. (2013). Extrapolation methods: theory and practice. Elsevier.
15. Monroe, J. L. (2002). Extrapolation and the Bulirsch-Stoer algorithm. Physical Review E, 65(6), 066116.
16. Kirpekar, S. (2003). Implementation of the Bulirsch Stoer extrapolation method. Department of Mechanical Engineering, UC Berkeley/California.
17. Nurminskii, E. A., & Buryi, A. A. (2011). Parker-Sochacki method for solving systems of ordinary differential equations using graphics processors. Numerical Analysis and Applications, 4(3), 223.
18. Hairer, E., Lubich, C., & Wanner, G. (2006). Geometric numerical integration: structure-preserving algorithms for ordinary differential equations (Vol. 31). Springer Science & Business Media.
19. Hairer, E., Lubich, C., & Wanner, G. (2003). Geometric numerical integration illustrated by the Störmer–Verlet method. Acta Numerica, 12, 399-450.
20. Nievergelt, Jürg. Parallel methods for integrating ordinary differential equations. Communications of the ACM. 1964, 7 (12): 731–733. S2CID 6361754. doi:10.1145/355588.365137 .
21. Parallel-in-Time.org. Parallel-in-Time.org. [2023-11-15]. （原始內容存檔於2023-11-15）.
22. Higham, N. J. (2002). Accuracy and stability of numerical algorithms (Vol. 80). SIAM.
23. Miranker, A. (2001). Numerical Methods for Stiff Equations and Singular Perturbation Problems: and singular perturbation problems (Vol. 5). Springer Science & Business Media.
24. Markus Kunze; Tassilo Kupper. Non-smooth Dynamical Systems: An Overview. Bernold Fiedler (編). Ergodic Theory, Analysis, and Efficient Simulation of Dynamical Systems. Springer Science & Business Media. 2001: 431. ISBN 978-3-540-41290-8.
25. Thao Dang. Model-Based Testing of Hybrid Systems. Justyna Zander, Ina Schieferdecker and Pieter J. Mosterman (編). Model-Based Testing for Embedded Systems. CRC Press. 2011: 411. ISBN 978-1-4398-1845-9.
26. Brezinski, C., & Wuytack, L. (2012). Numerical analysis: Historical developments in the 20th century. Elsevier.
27. Butcher, J. C. (1996). A history of Runge-Kutta methods. Applied numerical mathematics, 20(3), 247-260.
28. Ascher, U. M., Mattheij, R. M., & Russell, R. D. (1995). Numerical solution of boundary value problems for ordinary differential equations. Society for Industrial and Applied Mathematics.

參考文獻

• Bradie, Brian. A Friendly Introduction to Numerical Analysis. Upper Saddle River, New Jersey: Pearson Prentice Hall. 2006. ISBN 978-0-13-013054-9.
• J. C. Butcher, Numerical methods for ordinary differential equations, ISBN 0-471-96758-0
• Ernst Hairer, Syvert Paul Nørsett and Gerhard Wanner, Solving ordinary differential equations I: Nonstiff problems, second edition, Springer Verlag, Berlin, 1993. ISBN 3-540-56670-8.
• Ernst Hairer and Gerhard Wanner, Solving ordinary differential equations II: Stiff and differential-algebraic problems, second edition, Springer Verlag, Berlin, 1996. ISBN 3-540-60452-9.
  (This two-volume monograph systematically covers all aspects of the field.)
• Hochbruck, Marlis; Ostermann, Alexander. Exponential integrators. Acta Numerica. May 2010, 19: 209–286. Bibcode:2010AcNum..19..209H. CiteSeerX 10.1.1.187.6794 . S2CID 4841957. doi:10.1017/S0962492910000048.
• Arieh Iserles, A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, 1996. ISBN 0-521-55376-8 (hardback), ISBN 0-521-55655-4 (paperback).
  (Textbook, targeting advanced undergraduate and postgraduate students in mathematics, which also discusses numerical partial differential equations.)
• John Denholm Lambert, Numerical Methods for Ordinary Differential Systems, John Wiley & Sons, Chichester, 1991. ISBN 0-471-92990-5.
  (Textbook, slightly more demanding than the book by Iserles.)

外部連結

• Joseph W. Rudmin, Application of the Parker–Sochacki Method to Celestial Mechanics Portuguese Web Archive的存檔，存檔日期2016-05-16, 1998.
• Dominique Tournès, L'intégration approchée des équations différentielles ordinaires (1671-1914), thèse de doctorat de l'université Paris 7 - Denis Diderot, juin 1996. Réimp. Villeneuve d'Ascq : Presses universitaires du Septentrion, 1997, 468 p. (Extensive online material on ODE numerical analysis history, for English-language material on the history of ODE numerical analysis, see, for example, the paper books by Chabert and Goldstine quoted by him.)
• Pchelintsev, A.N. An accurate numerical method and algorithm for constructing solutions of chaotic systems. Journal of Applied Nonlinear Dynamics. 2020, 9 (2): 207–221. S2CID 225853788. arXiv:2011.10664 . doi:10.5890/JAND.2020.06.004.
• GitHub上的kv頁面 (C++ library with rigorous ODE solvers)
• INTLAB （頁面存檔備份，存於互聯網檔案館） (A library made by MATLAB/GNU Octave which includes rigorous ODE solvers)