常微分方程數值方法

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

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

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

## 問題

一階微分方程是有如下形式的初值問題（IVP）：

$$y'(t) = f(t,y(t)), \qquad y(t_0)=y_0,$$ 其中f是函數$f:[t_0, \infty) \times \R^d \to \R^d$，初始條件$y_0 \in \R^d$是已知向量。「一階」是說方程中只出現了y的一階導。

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

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

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

## 方法

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

所謂一般線性方法（GLM）是上述兩大類方法的概括。

### 歐拉法

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

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

$$y'(t) \approx \frac{y(t+h) - y(t)}{h},$$ 重新排列後得到以下公式

$$y(t+h) \approx y(t) + hy'(t)$$

利用(1)可得$$y(t+h) \approx y(t) + hf(t,y(t)).$$ 此式的應用通常如下：擇步長h，構造序列$t_0, t_1 = t_0 + h, t_2 = t_0 + 2h,...$記精確解$y(t_n)$的數值估計值為$y_n$。受(3)啟發，可用下面的遞歸方法計算估計值：

$$y_{n+1} = y_n + hf(t_n,y_n).$$ 這就是（前向）歐拉法，是萊昂哈德·歐拉（1768）描述的方法。

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

### 反向歐拉法

若不用(2)，而用近似值$$y'(t) \approx \frac{y(t) - y(t-h)}{h},$$ 則得到反向歐拉法：

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

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

### 一階指數積分法

指數積分 (微分方程)描述了一大類積分器，近來得到了長足發展。它們至少可以追溯到20世紀60年代。

此處不用(1)，假設微分方程形式為$$y'(t) = -A\, y+ \mathcal{N}(y),$$ 或已被局部線性化，圍繞背景狀態產生線性項$-Ay$與非線性項$\mathcal{N}(y)$。

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

$$y_{n+1} = e^{-A h } y_n + \int_{0}^{h} e^{ -(h-\tau) A } \mathcal{N}\left( y\left( t_n+\tau \right) \right)\, d\tau.$$

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

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

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

### 推廣

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

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

幾乎所有實用的多步法都屬於線性多步法，形式為$$\begin{align} &{} \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{align}$$

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

### 高級特徵

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

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

其他理想特徵還有：

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

### 其他方法

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

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

### 實時並行方法

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

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

## 分析

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

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

### 收斂性

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

$$\lim_{h\to0^+} \max_{n=0,1,\dots,\lfloor t^*/h\rfloor} \left\| y_{n,h} - y(t_n) \right\| = 0.$$

上述所有方法都收斂。

### 一致性與階數

設數值方法

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

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

$$\delta^h_{n+k} = \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\to 0} \frac{\delta^h_{n+k}}{h} = 0.$$

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

$$\delta^h_{n+k} = O(h^{p+1}) \quad\mbox{as } h\to0.$$

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

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

### 穩定性與剛性

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

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

## 歷史

下面是該領域一些重要進展的時間線。

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

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

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

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

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

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

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

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

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

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

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

$$\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$、$u(1)=u_n$，由於兩值已知，可以簡單地代入，就得到了具有非平凡解的非齊次線性方程組。

## 外部連結

- Joseph W. Rudmin, Application of the Parker–Sochacki Method to Celestial Mechanics , 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.)
- mskashi/kv (C++ library with rigorous ODE solvers)
- INTLAB (A library made by MATLAB/GNU Octave which includes rigorous ODE solvers)