单步法

数值分析里的单步法（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]。

相关条目

• 常微分方程数值方法