單步法

數值分析裡的單步法（one-step methods）和多步法（multistep methods）是求解初值問題的數值方法的兩種分類。常微分方程和初值問題在自然科學、工程物理學中都很重要，在經濟學和社會科學中的重要性也慢慢提高。

單步法近似初值問題，給定起始值${\displaystyle A_{0}}$，可以找到${\displaystyle A_{1}}$、${\displaystyle A_{2}}$等

單步法背後的基本概念是從給定的起點出發，沿著想要的解，一步一步的計算近似解。單步法只使用最近的資訊計算下一個點的資訊，多步法則不同，多步法還會考慮更早的資訊。單步法粗略可以分為兩類：顯式法是直接計算新的近似值，隱式法是將新的近似值和目前已知的資訊形成方程，需求解方程才能知道新的近似值。後者更適用於剛性方程（stiff equation，其數值分析的解只有在時間間隔很小時才會穩定，只要時間間隔略大，其解就會不穩定）。

最簡單也最古老的單步法，是顯式歐拉法，由萊昂哈德·歐拉在1768年發表。後來，在1883年，出現了許多的多步法。而卡爾·龍格、Karl Heun（英語：Karl Heun）和馬丁·威廉·庫塔在1900年針對歐拉法進行大幅改進，之後出現了許多龍格-庫塔法的算法，也是單步法中最重要的一組算法。二十世紀的發展包括外推的概念、以及考慮步長的控制，也就是在每一步經計算後，選擇下一步的步長。這些概念組成了求解複雜初值問題的基礎，常出現在現有的應用中，利用電腦可以有效率的求解問題，而且可以達到要求的精度。

基本概念

考慮以下初值問題的通式，要找到函數${\displaystyle y}$滿足以下的方程

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

其中${\displaystyle f}$是已知的函數

洛倫茲吸引子微分方程的解，是三維空間中的複雜曲線

若${\displaystyle y'(t)}$和${\displaystyle y(t)}$等比例，則${\displaystyle y(t)}$會指數增長。因此，${\displaystyle y'(t)=\lambda y(t)}$，其成長率是${\displaystyle \lambda }$，也有初始條件${\displaystyle y(0)=y_{0}}$。此時，可以用初等的微積分配合指數函數求得結果：${\displaystyle y(t)=y_{0}e^{\lambda t}}$。

微分方程中所要求的函數${\displaystyle y}$其值也可以是向量，也就是針對每一個${\displaystyle t}$，${\displaystyle y(t)=(y_{1}(t),\dotsc ,y_{d}(t))}$是一個${\displaystyle d}$個元素的向量。這也是${\displaystyle d}$維微分方程的系統。以移動物體為例，${\displaystyle y(t)}$是其d維空間中的位置，而${\displaystyle y'(t)}$是其時間${\displaystyle t}$的速度。因此微分方程標示了物體軌跡上每一點的速度方向和大小。就可以由此計算軌跡。

在上述簡化版的指數成長問題中，可以直接求解。但大多數的微分方程問題更複雜，也無法求得解析解。在特定條件下，可以證明針對函數${\displaystyle f}$，其初始值問題的解存在，但無法用分析的解法求解（例如分離變數法、指數法或是變分法），因此需要利用數值方法，找到其近似解。

初值問題的數值方法要在${\displaystyle t_{0}<t_{1}<t_{2}<\ldots }$各時間，逐步計算待求解函數的值${\displaystyle y(t_{0}),y(t_{1}),y(t_{2}),\dotsc }$。單步法在計算${\displaystyle y_{j+1}}$近似值時，只需要${\displaystyle y_{j}}$的近似值，與其相對的，多步法在計算${\displaystyle y_{j+1}}$近似值時，除了${\displaystyle y_{j}}$近似值以外，還至少需要${\displaystyle y_{j-1}}$，甚至${\displaystyle y_{j-2}}$......的值。

顯式歐拉法的前二步

最簡單也最基本的單步法就是顯式歐拉法，是瑞士數學家和物理學家萊昂哈德·歐拉於1768年提出的^[1]。此方法的概念是由分段線性函數來近似要求的解，而從點${\displaystyle t_{j+1}}$到點${\displaystyle f(t_{j},y_{j})}$的斜率是由${\displaystyle t_{j}}$時的資訊所決定。更深入的說：問題定義只提供了此函數的初值：${\displaystyle y(t_{0})=y_{0}}$，而且，此點的斜率也是已知的，${\displaystyle y'(t_{0})=f(t_{0},y_{0})}$。因此，可以確定函數在此點的切線斜率，並以此近似函數在往後時間的特性。在點${\displaystyle t_{1}>t_{0}}$，可以得到以下的結果，而步長為${\displaystyle h_{0}:=t_{1}-t_{0}}$

${\displaystyle y(t_{1})\approx y_{0}+h_{0}f(t_{0},y_{0})=:y_{1}}$.

此作法可以繼續進行。最後，以下計算的結果就是顯式歐拉法

${\displaystyle y_{j+1}=y_{j}+h_{j}f(t_{j},y_{j}),\quad j=0,1,2,\dotsc }$

其遞增量${\displaystyle h_{j}=t_{j+1}-t_{j}}$.^[2]

顯式歐拉法是許多單步法的基礎，其斜率${\displaystyle f(t_{j},y_{j})}$由其他近似函數在點${\displaystyle t_{j}}$和點${\displaystyle t_{j+1}}$之間行為的運算代替。單步法的另一個概念是來自隱式歐拉法，用${\displaystyle f(t_{j+1},y_{j+1})}$為斜率。初看會認為這作法不一定合適，因為此時${\displaystyle y_{j+1}}$未知，不過，若以此往下計算，就得到以下的方程

${\displaystyle y_{j+1}=y_{j}+h_{j}f(t_{j+1},y_{j+1})}$

可以以此求解${\displaystyle y_{j+1}}$（若有需要的話，也可以用數值方法求解）。例如用顯式歐拉法和隱式歐拉法的算術平均數作為斜率，這就是隱式梯形法。若等式右側未知的${\displaystyle y_{j+1}}$用顯式歐拉法近似，也可以得到另一個顯式法，這稱為Heun方法^[3]。所有方法和擴展都有一步法的基本概念：

${\displaystyle y_{j+1}=y_{j}+h_{j}\Phi }$

其斜率${\displaystyle \Phi }$是由${\displaystyle t_{j}}$ , ${\displaystyle y_{j}}$, ${\displaystyle h_{j}}$以及（針對隱式歐拉法）${\displaystyle y_{j+1}}$所決定。

定義

單步法可以定義如下：令 ${\displaystyle y}$ 是以下初值問題的解

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

假設解

${\displaystyle y\colon I\to \mathbb {R} ^{d}}$

存在在一定區間${\displaystyle I=[t_{0},T]}$內，而且其值唯一。以下時間

${\displaystyle t_{0}<t_{1}<t_{2}<\ldots <t_{n}=T}$

是區間${\displaystyle I}$內滿足上述條件的時間點，且${\displaystyle h_{j}=t_{j+1}-t_{j}}$是對應的增量，則以下

${\displaystyle y_{j+1}=y_{j}+h_{j}\Phi (t_{j},y_{j},y_{j+1},h_{j})}$, ${\displaystyle \quad j=0,\dotsc ,n-1}$

就是單步法的公式，其斜率函數為${\displaystyle \Phi }$。若${\displaystyle \Phi }$和${\displaystyle y_{j+1}}$無關，這個單步法就是顯式的，不然，在每一個${\displaystyle j}$都要求解方程，才能得到解在該時間的值，這就是隱式的單步法^[4]。

一致性和收斂性

收斂階數

針對實際使用的單步法，所計算的${\displaystyle y_{j}}$應該是其精確解${\displaystyle y}$在點${\displaystyle t_{j}}$的值${\displaystyle y(t_{j})}$之良好近似。若其變數是通用的${\displaystyle d}$維向量，此近似的好壞會用向量範${\displaystyle \|y_{j}-y(t_{j})\|}$來量測，也就是在該點${\displaystyle t_{j}}$的誤差。若可以將其步階收斂到零時，會希望針對所有${\displaystyle j}$的誤差都快速的收斂到零。為了考慮非定值步階的情形，會將${\displaystyle h}$準確定義為所使用步階的最大值，在所有點${\displaystyle j}$裡最大誤差的特性會和${\displaystyle h}$的冪次比較。求解給定初值問題的單步法，其收斂階數（Convergence order）為${\displaystyle p\geq 1}$，若以下估測值

${\displaystyle \max _{j=0,\dotsc ,n}\|y_{j}-y(t_{j})\|\leq Ch^{p}}$

適用在所有夠小的${\displaystyle h}$上，其常數${\displaystyle C>0}$，和${\displaystyle h}$無關^[5]。要比較不同的單步法時，收斂階數是最重要的參數^[6]。較高收斂階${\displaystyle p}$的方法，在固定步階下的總誤差較小，若要到給定的精度，需要的步數也比較少。針對${\displaystyle p=1}$的方法，可以預期若步階變一半，其總誤差也只會變成一半。若收斂階數${\displaystyle p=4}$的方法，步階變一半，總誤差會是${\displaystyle {\bigl (}{\tfrac {1}{2}}{\bigr )}^{4}={\tfrac {1}{16}}}$。

全域誤差和局部誤差

收斂階數的${\displaystyle \|y_{j}-y(t_{j})\|}$是由兩個獨立的因素所組成，乍看之下會很複雜。一方面，誤差是由每一步用數值方法近似函數未知梯度的誤差有關，但是，也需要考慮每一步的起始點${\displaystyle (t_{j},y_{j})}$和真正起始點${\displaystyle (t_{j},y(t_{j}))}$之間的誤差，目前這一步的誤差也會和以往各步的誤差有關。因為單步法一致性的定義，不同方式的差異只會在程序函數${\displaystyle \Phi }$的選擇，不過可以證明，在${\displaystyle \Phi }$符合特定條件下，可以從單一步的誤差階數推斷收斂階數，這稱為一致階數（consistency order）。

一致的概念是現代數值分析裡，通用的中心概念。數值分析方式的收斂性和數值近似如何近似其精確解有關，一致性要問的是逆問題的簡化版：其精確解可以多符合其數值方式的規格？在通用理論中，一個方法若是一致和穩定的，就是收斂的。為了簡化說明，以下的說明中假設存在顯式單步法

${\displaystyle y_{j+1}=y_{j}+h\Phi (t_{j},y_{j},h)}$

而且其固定步階為${\displaystyle h}$。其精確解為${\displaystyle t\mapsto y(t)}$，局部截尾誤差（也稱為局部處理誤差）${\displaystyle \eta }$ 定義為^[7]

${\displaystyle \eta (t,h)=y(t)+h\Phi (t,y(t),h)-y(t+h)}$.

因此，假設其精確解已知。從${\displaystyle (t,y(t))}$開始單步法的步階，形成了在${\displaystyle t+h}$和精確解的誤差。接著定義：單步法有一致階數（consistency order）${\displaystyle p\geq 1}$若以下估測值的算式成立

${\displaystyle \|\eta (t,h)\|\leq Ch^{p+1}}$

針對所有夠小的${\displaystyle h}$，其常數${\displaystyle C>0}$和${\displaystyle h}$無關。

一致階數和收斂階數中很大的差異是一致階數是${\displaystyle h^{p+1}}$而不是${\displaystyle h^{p}}$。這可以解釋成在局部誤差轉換到全域誤差時，少了一層步階${\displaystyle h}$的影響。以下定理是單步法的核心：^[8]

若過程函數${\displaystyle \Phi }$為利普希茨連續，且其相關的單步法有一致階數${\displaystyle p}$，則其也有收斂階數${\displaystyle p}$。

過程函數的利普希茨連續是針對穩定性的額外要求，若微分方程本身的${\displaystyle f}$是利普希茨連續，其過程函數多半也會是利普希茨連續。不過在大部份的應用都要假設此要求，以確保其初始值問題有唯一可解性，此時已可以確定單步法的一致階數。理論上，這可以用${\displaystyle \eta (t,h)}$的泰勒級數求得。但實務上，高階的公式很複雜，不易理解，因此需要額外的概念以及表示法來處理此議題^[9]。

剛性和A穩定性

主條目：剛性方程

數值方法的收斂階數是漸進式的表示方式，描述步階收斂到零時，其近似的行為。不過，它並沒有說明在給定步階時，該方法是否可以產生可用的近似值。Charles Francis Curtiss和Joseph O. Hirschfelder（英語：Joseph O. Hirschfelder）最早在1952年描述後者的問題可能是一些初始值問題的真正困難點。他們觀察到一些化學反應動力微分方程的解無法用顯式數值方法求解，稱這類初始值問題為「剛性」（stiff）^[10]。有許多的數學準則可以判斷給定問題的剛性程度。剛性初始值問題多半是微分方程系統，其中一些變數很快收斂變成常數，有些變數則很慢收斂。這類特性一般是發生在化學反應問題建模時。不過，實務應用問題中，最有用的剛性定義是：初始值問題是剛性問題，在用顯式方式解問題時，為了要得到可用的解，需要選擇一個「非常小」的步階。這問問題只能用隱式法求解^[11]

此問題中，若要計算指數遞減的解（藍色），若步階解太大，顯式法會嚴重振蕩，因此無法使用（紅色）。隱式歐拉法（綠色）可以在任意步階下，其解不會振盪

此效果可以檢驗各方法克服指數衰減的效果來說明。依照瑞典數學家Germund Dahlquist（英語：Germund Dahlquist），測試方程

${\displaystyle y'(t)=\lambda y(t),\quad y(0)=1}$

${\displaystyle \lambda <0}$時有指數遞減解${\displaystyle y(t)=e^{\lambda t}}$。旁邊的圖可以說明顯式法和隱式法在這種看似簡單問題的典型特性：若顯式法的步階太大，其數值會強烈振盪，比原來計算的值更大，而且越來越偏離正確解。另一方面，隱式法可以在任意步階求解，其解都是正確的，稱為近似值的指數遞減數列^[12]。

若將上式推廣，上述測試方程的${\displaystyle \lambda }$可以是複數。在此例中，當${\displaystyle \operatorname {Re} (\lambda )\leq 0}$（也就是${\displaystyle \lambda }$的實部小於等於0）時，解會振盪，而其振幅仍受到精確限制。因此可以制定可以用在剛性初值問題的單步法的理想屬性，所謂的A穩定性。一方法為A穩定，若一方法可以計算一個近似值數列，${\displaystyle h>0}$對應任何用在測試方程的步階${\displaystyle y_{0},y_{1},y_{2},\dotsc }$，針所有${\displaystyle \lambda }$其${\displaystyle \operatorname {Re} (\lambda )\leq 0}$，而且其近似值仍有然有界（類似真正的解）。隱式歐拉法和隱式梯型法是A穩定單步法最簡單的例子。另一方面，可以證明顯式法不可能是A穩定^[13]。

特殊方法和方法類

一階和二階的簡單方法

單步法的比較

如同法國數學家奧古斯丁-路易·柯西1820年所證明的，歐拉法的收斂階數為1。若將顯式法的斜率${\displaystyle f(t_{j},y_{j})}$和隱式歐拉法的斜率平均，這是一個步階的二個端點^[14]，可以得到整個區間內較佳的近似值。事實上，可以證明隱式梯形法可以用下式得到: ${\displaystyle y_{j+1}=y_{j}+{\frac {h}{2}}{\Big (}f(t_{j},y_{j})+f(t_{j+1},y_{j+1}){\Big )}}$ 其收斂階數為2。此方法的穩定特性很好，不過是隱式法，因此在每一步階需要求解𝑦_𝑗+1的方程。若這個變數是用方程右邊，用顯式歐拉法來近似，其結果就是顯式Heun法^[15]

${\displaystyle y_{j+1}=y_{j}+{\frac {h}{2}}{\Big (}f(t_{j},y_{j})+f{\big (}t_{j+1},y_{j}+hf(t_{j},y_{j}){\big )}{\Big )}}$,

其收斂階數也是2。另一個二階的顯式方法，改良的歐拉法，可以用以下的考量得到，此方法中步階的「平均」斜率是在步階中間解𝑦的斜率。不過，其解未知，可以用顯式歐拉法，以及一半的步階來近似。其結果是用以下的方式

${\displaystyle y_{j+1}=y_{j}+hf{\big (}t_{j}+{\tfrac {h}{2}},y_{j}+{\tfrac {h}{2}}f(t_{j},y_{j}){\big )}}$.

1895年時，德國數學家Carl Runge發布了這些二階的單步法，當做歐拉法的改良版^[16]。

龍格-庫塔法

經典的四階龍格-庫塔法，在每一步中對四個輔助梯度（紅色）進行平均

上述單步法的概念，若在進一步推廣，就可以得到龍格-庫塔法。例如，Heun法可以更清楚的表示如下：

首先，先計算輔助斜率${\displaystyle k_{1}=f(t_{j},y_{j})}$，這是顯式歐拉法的斜率。

上述斜率會用在計算下一個輔助斜率${\displaystyle k_{2}=f(t_{j}+h,y_{j}+hk_{1})}$。

實際用的梯度${\displaystyle \Phi }$會用上述兩個梯度平均而成，也就是Heun法裡的${\displaystyle {\tfrac {1}{2}}k_{1}+{\tfrac {1}{2}}k_{2}}$。

上述程式可以再擴展，不只使用兩個輔助斜率。s階龍格-庫塔法先在適當的點計算𝑓，以求得輔助斜率${\displaystyle k_{1},\dotsc ,k_{s}}$，接著加權平均，產生${\displaystyle \Phi }$。在顯式龍格-庫塔法中，輔助斜率${\displaystyle k_{1},k_{2},k_{3},\dotsc }$是直接一個接一個計算，而在隱式方法中，用求解方程來求得輔助斜率。一個典型的例子是4階的顯式龍格-庫塔法，有時會直接用龍格-庫塔法來指稱此特定方法：先計算以下四個輔助斜率^[17]

${\displaystyle {\begin{aligned}k_{1}&=f(t_{j},y_{j})\\k_{2}&=f(t_{j}+{\tfrac {h}{2}},y_{j}+{\tfrac {h}{2}}k_{1})\\k_{3}&=f(t_{j}+{\tfrac {h}{2}},y_{j}+{\tfrac {h}{2}}k_{2})\\k_{4}&=f(t_{j}+h,y_{j}+hk_{3})\end{aligned}}}$

接著算以下的加權平均。

${\displaystyle {\tfrac {1}{6}}k_{1}+{\tfrac {1}{3}}k_{2}+{\tfrac {1}{3}}k_{3}+{\tfrac {1}{6}}k_{4}}$

此方法是由德國數學家馬丁·威廉·庫塔在1901年提出，而三年前Karl Heun有找到一個三階單步法，其收斂階數也是三^[18]。

外推法

用階數${\textstyle p=2}$的方法外推${\displaystyle h=0}$的值

外推的概念不限制在用單步法算初始值問題的解中，也可以應用在用步階${\displaystyle h}$將問題離散化，進而求解的問題。廣為人知的外推法是數值積分的Romberg積分（英語：Romberg integration）。一般來說，令${\displaystyle v}$是要知道其數值的值，在此條目中，是初始值問題解在某一點的值。數值方法（例如單步法），可以計算其近似值${\displaystyle {\tilde {v}}(h)}$，數值依步階大小${\displaystyle h>0}$而定。假設此方法收斂，也就是${\displaystyle {\tilde {v}}(h)}$收斂到${\displaystyle v}$，當${\displaystyle h}$收斂到零時。不過，此收斂只是純理論的說法，因為在有限個數的步階${\displaystyle h_{1}>h_{2}>\ldots >h_{m}}$下，可以計算近似值${\displaystyle {\tilde {v}}(h_{1}),{\tilde {v}}(h_{2}),\dotsc ,{\tilde {v}}(h_{m})}$，步階大小當然不能「收斂為0」。不過，不同步階的計算近似可以解讀為對未知函數的資訊${\displaystyle {\tilde {v}}}$：在外推法中，${\displaystyle {\tilde {v}}}$會用內插多項法來近似，也就是滿足以下性質的多項式${\displaystyle P}$^[19]

${\displaystyle P(h_{k})={\tilde {v}}(h_{k})}$

其中${\displaystyle k=1,2,\dotsc ,m}$。多項式在${\displaystyle h=0}$的值${\displaystyle P(0)}$可以用作在${\displaystyle h}$趨近於0時，${\displaystyle {\tilde {v}}(h)}$不可計算極限值的可計算近似。早期針對初始值問題的成功外推法是由Roland Bulirsch（英語：Roland Bulirsch）和Josef Stoer（英語：Josef Stoer）在1966年發表^[20]。

用一個${\displaystyle p}$階的單步法來說明外推法的步驟。在此種方法下，小步階ℎ下計算的近似值可以用以下的多項式表示

${\displaystyle P(h)=a+bh^{p}}$

其參數${\displaystyle a}$和${\displaystyle b}$未知。若分別用步階${\displaystyle h_{1}}$和${\displaystyle h_{2}={\tfrac {1}{2}}h_{1}}$，計算兩個近似值${\displaystyle y_{h_{1}}}$和${\displaystyle y_{h_{2}}}$，可以從內插條件${\displaystyle P(h_{1})=y_{h_{1}}}$ and ${\displaystyle P(h_{2})=y_{h_{2}}}$求得方程，得到${\displaystyle a}$和${\displaystyle b}$。

外推的數值為

${\displaystyle P(0)=a=y_{h_{2}}+{\frac {y_{h_{2}}-y_{h_{1}}}{2^{p}-1}}}$

會是比前兩個近似值更好的近似值。可以證明這個方法的一步法，其階數至少是${\displaystyle p+1}$，也就是比原來的階數至少再高一階^[21]。

步長控制方法

單步法的優點是每一步的步長可以獨立控制。因此就自然會有每一步的步長ℎ𝑗如何選擇的問題。實務的應用上，初始值問題的解一定會有允許的誤差範圍。若某個數值方法的精度遠比其問題的初始值和參數的精度（初始值和參數的精度會受測量誤差限制）要高，在實務上其實沒什麼用。因此其目標是選擇適當的步長大小，一方面符合誤差範圍，另一方面也讓運算量可以降到最低。配合初始條件的常微分方程，在自然科學和工程科學上相當重要，在經濟學和社會科學上也越來越重要^[22]。

針對良好條件數的初始值問題，可以證明其全域誤差大致等於每一步局部截尾誤差的和${\displaystyle \eta _{j}:=\|\eta (t_{j},h_{j})\|}$。因此，應該選押最大可能的${\displaystyle h_{j}}$作步長，而選擇${\displaystyle \eta _{j}}$低於選擇的容許誤差範圍。不過有一個問題：${\displaystyle \eta _{j}}$無法直接計算，${\displaystyle \eta _{j}}$和初始值問題在${\displaystyle t_{j}}$點的正確解（其解未知）有關。步長控制的基本概念是用一個比基本方法更精確的方法來近似${\displaystyle y(t_{j})}$^[23]。

步長控制的兩個基本概念是步長減半（step width halving）和embedded process。若是用步長減半的方法，會用步長減半後，兩步所得的值，和原始步長的單步法結果比較。透過二個值的外推法可以得到${\displaystyle y(t_{j})}$更準確的近似值，並且估計局部誤差${\displaystyle \eta _{j}}$。若是太大，就要再選擇更小的步長。若局部誤差明顯大於要求的允許誤差，下一步就可以將步長加大^[24]。步長減半法會增加的運算量較大，因此目前的實現方式會用embedded process來進行步長控制。其基本概念是每一步用二個不同的單步法近似${\displaystyle y(t_{j})}$，兩者有不同的收斂階數，因此可以估測局部誤差。為了減小運算量，這兩種方法會選擇有較多相近可共用的計算：讓二個近似「彼此嵌合」。Embedded Runge-Kutta法，就用了相同的輔助斜率，只是平均的係數不同。著名的embedded方法包括Runge-Kutta-Fehlberg方法（Erwin Fehlberg​（德語）, 1969）和Dormand-Prince法（J. R. Dormand和P. J. Prince, 1980）^[25]。

相關條目

• 常微分方程數值方法