常微分方程數值方法 是用以尋找常微分方程 （ODE）解的數值 近似值的方法。其使用也稱作「數值積分 」，不過「數值積分」主要是指積分 的計算。
微分方程
y
′
=
y
,
y
(
0
)
=
1
{\displaystyle y'=y,\ y(0)=1}
藍：歐拉法
紅：精確解：
y
=
e
t
{\displaystyle y=e^{t}}
h
=
1.0
{\displaystyle h=1.0}
h
=
0.25
{\displaystyle h=0.25}
h
→
0
{\displaystyle h\to 0}
很多微分方程無法精確求解。但在工程學等領域的實際應用中，通常只需得到數值近似解。本文介紹的算法 可用於計算這種近似值，另一種方法是用微積分 技術得到解的級數展開 表達。
常微分方程出現於物理學 、化學 、生物學 、經濟學 等學科中。[ 1] 偏微分方程數值方法 中的一部分將偏微分方程 轉為常微分方程，然後可用本文所述方法求解。
一階微分方程是有如下形式的初值問題 （IVP）：[ 2] :533–655
y
′
(
t
)
=
f
(
t
,
y
(
t
)
)
,
y
(
t
0
)
=
y
0
,
{\displaystyle y'(t)=f(t,y(t)),\qquad y(t_{0})=y_{0},}
1
其中f 是函數
f
:
[
t
0
,
∞
)
×
R
d
→
R
d
{\displaystyle f:[t_{0},\infty )\times \mathbb {R} ^{d}\to \mathbb {R} ^{d}}
y
0
∈
R
d
{\displaystyle y_{0}\in \mathbb {R} ^{d}}
y 的一階導。
在不失高階系統一般性的前提下，限制(l)式為一階微分方程，因為高階ODE可通過引入額外變量轉換為更大的一階方程組。例如，二階方程
y
″
=
−
y
{\displaystyle y''=-y}
y
′
=
z
,
z
′
=
−
y
{\displaystyle y'=z,\ z'=-y}
本文介紹IVP的數值方法，並指出邊值問題 （BVP）需要一套不同的工具：需要在多個點上定義解y 的值或成分，因此要用不同方法求解BVP，如打靶法 （及其變體），或差分 、[ 3] 伽遼金法 、[ 4] 配置法
Picard–Lindelöf定理 指出，只要f 是利普希茨連續 的，就有唯一解。
解一階IVP的數值方法可分為兩大類：[ 5] 線性多步法 龍格-庫塔法 。還可進一步劃為顯式或隱式，例如隱式線性多步法包括亞當斯-莫爾頓法（Adams-Moulton methods）與向後微分公式 [ 6] :204–215 包括對角隱式龍格-庫塔法（DIRK）、[ 7] [ 8] [ 9] 高斯求積 [ 10] [ 11] 亞當斯-巴什福思法 ，布徹表（Butcher tableau）為下對角的龍格-庫塔法都是顯式方法。根據經驗，剛性 微分方程需要用隱式方案，而非剛性問題則可用顯示方案更有效地求解。
所謂一般線性方法 [ 12]
從曲線上任意一點出發，沿與曲線相切的直線移動一小段距離，就能得到曲線上鄰近點的近似值。
從微分方程(1 差分 近似代替導數y ′：
y
′
(
t
)
≈
y
(
t
+
h
)
−
y
(
t
)
h
,
{\displaystyle y'(t)\approx {\frac {y(t+h)-y(t)}{h}},}
2
重新排列後得到以下公式
y
(
t
+
h
)
≈
y
(
t
)
+
h
y
′
(
t
)
{\displaystyle y(t+h)\approx y(t)+hy'(t)}
利用(1
y
(
t
+
h
)
≈
y
(
t
)
+
h
f
(
t
,
y
(
t
)
)
.
{\displaystyle y(t+h)\approx y(t)+hf(t,y(t)).}
3
此式的應用通常如下：擇步長h ，構造序列
t
0
,
t
1
=
t
0
+
h
,
t
2
=
t
0
+
2
h
,
.
.
.
{\displaystyle t_{0},t_{1}=t_{0}+h,t_{2}=t_{0}+2h,...}
y
(
t
n
)
{\displaystyle y(t_{n})}
y
n
{\displaystyle y_{n}}
3 遞歸 方法計算估計值：
y
n
+
1
=
y
n
+
h
f
(
t
n
,
y
n
)
.
{\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n}).}
4
這就是（前向）歐拉法 ，是萊昂哈德·歐拉 （1768）描述的方法。
歐拉法是顯式 方法的一個例子，這是說新值
y
n
+
1
{\displaystyle y_{n+1}}
y
n
{\displaystyle y_{n}}
若不用(2
y
′
(
t
)
≈
y
(
t
)
−
y
(
t
−
h
)
h
,
{\displaystyle y'(t)\approx {\frac {y(t)-y(t-h)}{h}},}
5
則得到反向歐拉法 ：
y
n
+
1
=
y
n
+
h
f
(
t
n
+
1
,
y
n
+
1
)
.
{\displaystyle y_{n+1}=y_{n}+hf(t_{n+1},y_{n+1}).}
6
反向歐拉法是隱式 方法，這是說需要求解一個方程才能得到新值
y
n
+
1
{\displaystyle y_{n+1}}
定點迭代 或牛頓-拉弗森法 （的某種修改版）實現之。
隱式方法求解這方程比顯示方法直接代入要花更多時間，選擇方法時必須考慮這一成本。隱式方法（如(6 剛性方程 時通常更穩定，便可以使用更大的步長h 。
更多資訊：指數積分
指數積分 [ 13]
此處不用(1
y
′
(
t
)
=
−
A
y
+
N
(
y
)
,
{\displaystyle y'(t)=-A\,y+{\mathcal {N}}(y),}
7
或已被局部線性化，圍繞背景狀態產生線性項
−
A
y
{\displaystyle -Ay}
N
(
y
)
{\displaystyle {\mathcal {N}}(y)}
將(7
e
A
t
{\textstyle e^{At}}
[
t
n
,
t
n
+
1
=
t
n
+
h
]
{\displaystyle [t_{n},t_{n+1}=t_{n}+h]}
y
n
+
1
=
e
−
A
h
y
n
+
∫
0
h
e
−
(
h
−
τ
)
A
N
(
y
(
t
n
+
τ
)
)
d
τ
.
{\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 .}
此積分方程是精確的，但並沒有定義積分。
使
N
(
y
(
t
n
+
τ
)
)
{\displaystyle {\mathcal {N}}(y(t_{n}+\tau ))}
y
n
+
1
=
e
−
A
h
y
n
+
A
−
1
(
1
−
e
−
A
h
)
N
(
y
(
t
n
)
)
.
{\displaystyle y_{n+1}=e^{-Ah}y_{n}+A^{-1}(1-e^{-Ah}){\mathcal {N}}(y(t_{n}))\ .}
8
歐拉法往往不夠精確。更準確地說，只有一階（下面將介紹階的概念），這就促使數學家尋找高階方法。
一種方法是，用
y
n
{\displaystyle y_{n}}
y
n
+
1
{\displaystyle y_{n+1}}
多步法 便實現了這種想法。最簡單的可能是蛙跳積分法 ，其是二階方法，（粗略地說）依賴於2個時間值。
幾乎所有實用的多步法都屬於線性多步法
α
k
y
n
+
k
+
α
k
−
1
y
n
+
k
−
1
+
⋯
+
α
0
y
n
=
h
[
β
k
f
(
t
n
+
k
,
y
n
+
k
)
+
β
k
−
1
f
(
t
n
+
k
−
1
,
y
n
+
k
−
1
)
+
⋯
+
β
0
f
(
t
n
,
y
n
)
]
.
{\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}}}
另一種方法是在區間
[
t
n
,
t
n
+
1
]
{\displaystyle [t_{n},t_{n+1}]}
卡爾·龍格 與馬丁·威廉·庫塔 的龍格-庫塔法 ，其中一種4階方法尤為流行。
要用這些方法求解ODE，需要的不僅是時間步長公式。始終相同的步長效率不高，於是開發了可變步長方法 。通常，步長的選擇應使每步的（局部）誤差低於某容差水平，這意味着方法還要計算誤差指標（error indicator），即對局部誤差的估計。
這一思想的延伸是在不同階方法之間動態選擇（稱作可變階數方法 ）。基於理查德森外推法 [ 14] Bulirsch–Stoer算法 之類方法[ 15] [ 16]
其他理想特徵還有：
輸出稠密 ：不僅對點
t
0
,
t
1
,
t
2
,
…
{\displaystyle t_{0},\ t_{1},\ t_{2},\ \dots }
事件定位 ：找到事件發生的位置，例如某函數取0的時間。通常用求根算法 ；支持並行計算 ；
用於時間積分時，具有時間可逆性。
很多方法不在討論的框架內，如
多階導數法 ，不僅使用f ，還用其導數。這類方法包括Hermite–Obreschkoff法 、費爾貝格法 y 的泰勒級數 係數的Parker–Sochacki法 [ 17] 二階常微分方程組法 。前面說過，所有高階常微分方程都可化為形式如(1)的一階常微分方程組，這固然沒錯，但未必是最佳方法。尼斯特倫法 幾何積分 法[ 18] [ 19] 哈密頓力學 方程的辛積分 量化狀態系統方法
有些IVP要求積分具有很高的時間分辨率和/或很長的時間區間，傳統的序列時間步進法無法實時計算（如數值天氣預報、等離子體建模與分子動力學中的IVP）。針對這些問題，人們開發了實時並行（PinT）法以便用並行計算 縮短運行時間。
早期的PinT法（最早提出於20世紀60年代）[ 20] Parareal算法 百億億次級運算 超級計算機 的算法開發。截至2023年，最流行的方法有Parareal、PFASST、ParaDiag、MGRIT等，[ 21]
數值分析 不僅包括數值方法的設計，還包括其分析。分析中的3個核心概念是：
若數值解隨着步長h 趨近於0而逼近精確解，則稱此數值方法具有收斂性 （convergent）。更確切地說，要求對利普希茨 函數f 與每個
t
∗
>
0
{\displaystyle t^{*}>0}
lim
h
→
0
+
max
n
=
0
,
1
,
…
,
⌊
t
∗
/
h
⌋
‖
y
n
,
h
−
y
(
t
n
)
‖
=
0.
{\displaystyle \lim _{h\to 0^{+}}\max _{n=0,1,\dots ,\lfloor t^{*}/h\rfloor }\left\|y_{n,h}-y(t_{n})\right\|=0.}
上述所有方法都收斂。
更多資訊：截斷誤差
設數值方法
y
n
+
k
=
Ψ
(
t
n
+
k
;
y
n
,
y
n
+
1
,
…
,
y
n
+
k
−
1
;
h
)
.
{\displaystyle y_{n+k}=\Psi (t_{n+k};y_{n},y_{n+1},\dots ,y_{n+k-1};h).\,}
方法的局部（截斷）誤差 定義為迭代一步產生的誤差。即，假設之前的迭代無誤差，則是此方法結果與精確解之間的差
δ
n
+
k
h
=
Ψ
(
t
n
+
k
;
y
(
t
n
)
,
y
(
t
n
+
1
)
,
…
,
y
(
t
n
+
k
−
1
)
;
h
)
−
y
(
t
n
+
k
)
.
{\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}).}
若
lim
h
→
0
δ
n
+
k
h
h
=
0.
{\displaystyle \lim _{h\to 0}{\frac {\delta _{n+k}^{h}}{h}}=0.}
則稱此方法一致（consistent）。若
δ
n
+
k
h
=
O
(
h
p
+
1
)
as
h
→
0.
{\displaystyle \delta _{n+k}^{h}=O(h^{p+1})\quad {\mbox{as }}h\to 0.}
則稱此方法階數為p 。因此，階數不為0的方法是一致的。上文介紹的兩種歐拉法(4、6)都是1階的，因此是一致的。實踐中使用的大多數方法階數更高。一致性是收斂的必要條件[來源請求]
一個相關概念是全局（截斷）誤差 ，即達到固定時間t 所需所有步驟中持續存在的誤差。明確地說，t 時刻的全局誤差是
y
N
−
y
(
t
)
{\displaystyle y_{N}-y(t)}
N
=
(
t
−
t
0
)
/
h
{\displaystyle N=(t-t_{0})/h}
p 階一步法是收斂的，其全局誤差是
O
(
h
p
)
{\displaystyle O(h^{p})}
對某些微分方程，歐拉法、顯式龍格-庫塔法 、多步法（如亞當斯-巴什福思法）之類標準方法會表現出解的不穩定性，其他方法則可能產生穩定的解。方程中的這種「困難行為」（本身不一定複雜）稱作「剛性」（stiffness），通常是由於底問題中存在不同時間尺度造成的。[ 23]
剛性問題在化學動力學 、控制論 、固體力學 、天氣預報 、生物學 、等離子體 物理、電子學 等領域中無處不在。克服剛性的一種方法是將微分方程推廣到微分包含式 ，從而允許非光滑性並建立其模型。[ 24] [ 25]
下面是該領域一些重要進展的時間線。[ 26] [ 27]
邊值問題（BVPs）通常要離散化，得到近似相等的矩陣問題再數值求解。[ 28] 有限差分法 。[ 3] 有限差分係數 ，例如一階導數的二階中心差分近似為：
u
i
+
1
−
u
i
−
1
2
h
=
u
′
(
x
i
)
+
O
(
h
2
)
,
{\displaystyle {\frac {u_{i+1}-u_{i-1}}{2h}}=u'(x_{i})+{\mathcal {O}}(h^{2}),}
二階導數的二階中心差分近似為：
u
i
+
1
−
2
u
i
+
u
i
−
1
h
2
=
u
″
(
x
i
)
+
O
(
h
2
)
.
{\displaystyle {\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}=u''(x_{i})+{\mathcal {O}}(h^{2}).}
兩式中，
h
=
x
i
−
x
i
−
1
{\displaystyle h=x_{i}-x_{i-1}}
x 值間的距離。這樣就構建了線性系統，然後可用標準矩陣法 求解。例如，待解方程：
d
2
u
d
x
2
−
u
=
0
,
u
(
0
)
=
0
,
u
(
1
)
=
1.
{\displaystyle {\begin{aligned}&{}{\frac {d^{2}u}{dx^{2}}}-u=0,\\&{}u(0)=0,\\&{}u(1)=1.\end{aligned}}}
下一步是將問題離散化，用線性導數近似，如
u
i
″
=
u
i
+
1
−
2
u
i
+
u
i
−
1
h
2
{\displaystyle u''_{i}={\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}}
並解所得線性方程組。將有如下方程
u
i
+
1
−
2
u
i
+
u
i
−
1
h
2
−
u
i
=
0
,
∀
i
=
1
,
2
,
3
,
.
.
.
,
n
−
1
.
{\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時，有一項涉及邊值
u
(
0
)
=
u
0
{\displaystyle u(0)=u_{0}}
u
(
1
)
=
u
n
{\displaystyle u(1)=u_{n}}
線性方程組 。
CFL條件 能量漂移 一般線性方法 Modelica 語言及OpenModelica
Chicone, C. (2006). Ordinary differential equations with applications (Vol. 34). Springer Science & Business Media.
LeVeque, R. J. (2007). Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems (Vol. 98). SIAM.
Slimane Adjerid and Mahboub Baccouch (2010) Galerkin methods. Scholarpedia, 5(10):10056.
Griffiths, D. F., & Higham, D. J. (2010). Numerical methods for ordinary differential equations: initial value problems. Springer Science & Business Media.
Alexander, R. (1977). Diagonally implicit Runge–Kutta methods for stiff ODE’s. SIAM Journal on Numerical Analysis, 14(6), 1006-1021.
Cash, J. R. (1979). Diagonally implicit Runge-Kutta formulae with error estimates. IMA Journal of Applied Mathematics, 24(3), 293-301.
Ferracina, L., & Spijker, M. N. (2008). Strong stability of singly-diagonally-implicit Runge–Kutta methods. Applied Numerical Mathematics, 58(11), 1675-1686.
Everhart, E. (1985). An efficient integrator that uses Gauss-Radau spacings. In International Astronomical Union Colloquium (Vol. 83, pp. 185-202). Cambridge University Press.
Butcher, J. C. (1987). The numerical analysis of ordinary differential equations: Runge-Kutta and general linear methods. Wiley-Interscience.
Hochbruck (2010 ，第209–286頁) harvtxt模板錯誤: 無指向目標: CITEREFHochbruck2010 (幫助 ) This is a modern and extensive review paper for exponential integratorsBrezinski, C., & Zaglia, M. R. (2013). Extrapolation methods: theory and practice. Elsevier.
Monroe, J. L. (2002). Extrapolation and the Bulirsch-Stoer algorithm. Physical Review E, 65(6), 066116.
Kirpekar, S. (2003). Implementation of the Bulirsch Stoer extrapolation method. Department of Mechanical Engineering, UC Berkeley/California.
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.
Hairer, E., Lubich, C., & Wanner, G. (2006). Geometric numerical integration: structure-preserving algorithms for ordinary differential equations (Vol. 31). Springer Science & Business Media.
Hairer, E., Lubich, C., & Wanner, G. (2003). Geometric numerical integration illustrated by the Störmer–Verlet method. Acta Numerica, 12, 399-450.
Higham, N. J. (2002). Accuracy and stability of numerical algorithms (Vol. 80). SIAM.
Miranker, A. (2001). Numerical Methods for Stiff Equations and Singular Perturbation Problems: and singular perturbation problems (Vol. 5). Springer Science & Business Media.
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 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 Brezinski, C., & Wuytack, L. (2012). Numerical analysis: Historical developments in the 20th century. Elsevier.
Butcher, J. C. (1996). A history of Runge-Kutta methods. Applied numerical mathematics, 20(3), 247-260.
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
Dominique Tournès, L'intégration approchée des équations différentielles ordinaires (1671-1914)
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) 
![{\displaystyle [t_{n},t_{n+1}=t_{n}+h]}](http://wikimedia.org/api/rest_v1/media/math/render/svg/12509e2c621cd352867e237af8a03b79eddcc166)
![{\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}}}](http://wikimedia.org/api/rest_v1/media/math/render/svg/47bb8bad1854dddf30061661932c07d46850ceff)
![{\displaystyle [t_{n},t_{n+1}]}](http://wikimedia.org/api/rest_v1/media/math/render/svg/ea7558633c567b46e0e6a00c0284e61646613660)
