微分方程的数值解法

本文将研究两种常见的常微分方程的数值解法：欧拉法和龙格-库塔法。

问题描述

本文我们主要关心寻找以下柯西问题的数值解： \(\left\{ \begin{aligned} y^\prime = f(x,y) \\ y(x_0) = y_0 \end{aligned} \right.\)

所谓寻找数值解，就是在一系列离散的点上给出微分方程解$y(x_n)$的近似值。 两个相邻的点之间的距离称为步长$h$，若步长为常数，则每个点都可写为 \(x_n = x_0 + h \cdot n\)

最常用的解法为步进法，即沿离散点从小到大的顺序依次计算，这种算法的核心在于从之前的点处的近似值计算当前点的近似值的递推公式。 如果递推公式中只需要之前一个点的值，那么这种方法称为单步法，否则称为多步法。

单步法容易实现且易于分析，因此我们本文中主要关心的方法都是单步法。

解法又可分为显式法和隐式法，其区别在于求解$y_n$的方程中是否需要$y_n$的值。 举例来说，对显示的单步法，递推公式总是可以写为： \(y_{k+1} = y_k + h \phi (x_k, y_k, h)\) 而对隐式的单步法，则有： \(y_{k+1} = y_k + h \phi (x_k, \textcolor{red}{y_{k+1}}, y_k, h)\) 隐式法通常数值稳定性好于显式法，但实现难度更高。

微分方程解的唯一性

（柯西-利普希茨定理） 对柯西初值问题： \(\left\{ \begin{aligned} y^\prime = f(x,y) \\ y(x_0) = y_0 \end{aligned} \right.\) 若函数关于$y$局部利普西茨连续，即存在常数$K$，满足 \(\forall y_1, y_2, \; \left| f(x,y_1) - f(x,y_2) \right| \le K | y_1 - y_2 |\) 那么这个初值问题的局部解以至最大解是唯一的。

在之后的研究中，我们总是假设这个条件成立。

局部截断误差

设$y(x)$是初值问题的精确解，而$y_n$为其数值解，那么： \(T_{n} = y(x_{n}) - y_{n}\) 称为其局部截断误差。

对显示单步法，这个误差可以化为： \(T_{n+1} = y(x_{n+1}) - y(x_n) - h \phi (x_n, y(x_n), h)\)

若局部截断误差是步长的$p$阶无穷小，那么一个数值方法可以说具有$p$阶精度，更形式化地说：

若存在最大整数$p$，使得局部截断误差可以写为： \(T_n = \psi h^{p+1} + \mathcal{O}(h^{p+2})\) 那么该方法称为具有$p$阶精度，其中$\psi h^{p+1}$称为其局部截断误差主项。

欧拉法解微分方程

我们首先关注最简单的数值解方法，即欧拉法。

显式、隐式与梯形欧拉法

欧拉法基于一个非常简单的变换： \(y^\prime = f(x, y) \implies y(x_{n+1}) - y(x_n) = \int_{x_{n}}^{x_{n+1}} f(x,y(x)) \mathrm d x\) 对于右侧的积分，我们可以使用三种方法进行近似。

如果使用左矩形积分，则有： \(\begin{multline} \int_{x_{n}}^{x_{n+1}} f(x,y(x)) \mathrm d x = h f(x_n, y(x_n)) \\ \implies y_{n+1} - y_n = h f(x_n, y_n) \end{multline}\) 这种方法称为显式欧拉法。

如果使用右矩形积分，则有： \(\begin{multline} \int_{x_{n}}^{x_{n+1}} f(x,y(x)) \mathrm d x = h f(x_{n+1}, y(x_{n+1})) \\ \implies y_{n+1} - y_n = h f(x_{n+1}, y_{n+1}) \end{multline}\) 这种方法称为隐式欧拉法，或称后退欧拉法。

如果使用梯形公式积分，则有： \(\begin{multline} \int_{x_{n}}^{x_{n+1}} f(x,y(x)) \mathrm d x = h \frac{f(x_n, y(x_n)) + f(x_{n+1}, y(x_{n+1}))}{2} \\ \implies y_{n+1} - y_n = \frac{h}{2} (f(x_n, y_n) + f(x_{n+1}, y_{n+1})) \end{multline}\) 这种方法称为梯形欧拉法，也是一种隐式方法。

隐式法的迭代计算

使用隐式法的一大问题在于如何计算$y_{n+1}$，这里给出一种使用迭代方法进行计算的思路。

以后退欧拉法为例子，我们首先利用显式欧拉法给出一个估计值： \(y_{n+1,0} = y_n + h f(x_n, y_n)\) 然后代入隐式欧拉法的公式，进行迭代： \(y_{n+1, k+1} = y_n + h f(x_{n+1}, y_{n+1, k})\) 经过数次迭代之后，就能达到需要的精度。

在满足李普希茨条件的情况下，对于足够小的步长，迭代方法计算出的数值解收敛，且收敛至 \(y_{n+1} = y_n + h f(x_{n+1}, y_{n+1})\)

我们有： \(\begin{aligned} |y_{n+1, k+1} - y_{n+1}| &= | [y_n + h f(x_{n+1}, y_{n+1, k})] - [y_n + h f(x_{n+1}, y_{n+1})] | \\ &= h | f(x_{n+1}, y_{n+1, k}) - f(x_{n+1}, y_{n+1}) | \\ &\le hK |y_{n+1, k} - y_{n+1}| \end{aligned}\) 从而，只要$hK < 1$，迭代法求出的值就能收敛至$y_{n+1}$

注意到这种方法实际上是求出函数$g(y_{n+1}) = y_n + h f(x_{n+1}, y_{n+1})$的不动点，因此这种方法也称不动点迭代法。 实际上，这个函数的李普希茨常数为$hK$，若这个常数小于一，则巴拿赫不动点定理指出这个不动点唯一存在。 也可以使用其他的迭代方法，比如牛顿迭代法，来求解这个方程。

对梯形法，我们使用以下迭代公式： \(\begin{cases} y_{n+1, 0} = y_n + h f(x_n, y_n) \\ y_{n+1, k+1} = y_n + \frac{h}{2} \left[ f(x_n, y_n) + f(x_{n+1}, y_{n+1,k}) \right] \end{cases}\) 不难用类似的方法证明其收敛。

截断误差分析

对于显式法，其局部截断误差为： \(\begin{aligned} T_{n+1} &= y(x_{n+1}) - y(x_n) - h \phi(x_n, y(x_n)) \\ &= y(x_n + h) - y(x_n) - h y^\prime(x_n) \\ &= y(x_n) + h y^\prime(x_n) + \frac{h^2}{2} y^{\prime\prime}(x_n) + \mathcal O(h^3) - y(x_n) - h y^\prime(x_n) \\ &= \frac{h^2}{2} y^{\prime\prime}(x_n) + \mathcal O(h^3) \end{aligned}\) 从而其具有一阶精度，局部截断误差主项为$\frac{h^2}{2} y^{\prime\prime}(x_n)$。

对隐式法和梯形法，相同的分析给出：

• 隐式法的主项为$-\frac{h^2}{2} y^{\prime\prime}(x_n)$，具有一阶精度。
• 梯形法的主项为$-\frac{h^3}{12} y^{\prime\prime\prime}(x_n)$，具有二阶精度。

改进欧拉法

隐式法相较于显式法通常数值稳定性更好，而梯形法相较于后退法的精度更高。 然而，这种方法却要求迭代地求解一个代数方程，这极大的提高了每一次求解的计算量。

改进的欧拉法只进行一次迭代，以此在速度和精度之间取得平衡。

改进的欧拉法使用以下公式： \(\begin{cases} y_p = y_n + hf(x_n, y_n) \\ y_c = y_n + hf(x_{n+1}, y_p) \\ y_{n+1} = \frac{1}{2}(y_p + y_c) \end{cases}\)

其中，$p$表示预测，$c$表示校正。 改进欧拉法亦称休恩方法，以德国数学家卡尔·休恩命名。

龙格-库塔法解微分方程

欧拉法的根本思路和局限性在于数值积分公式的使用。 这为我们提供了提高其精度的一种思路，即使用更精确的公式进行数值积分。 为了达成这种目标，我们需要提高在区间$[x_n, x_{n+1}]$上进行数值积分的选点数，此时我们继续进行近似： \(\int_{x_n}^{x_{n+1}} f(x, y(x)) \mathrm d x \approx \sum_{i=1}^r c_i f(x_{n,i}, y(x_{n,i}))\) 更进一步地，我们可以把$x_{n,i}$写成到$x_n$的距离的形式： \(x_{n,i} = x_n + \lambda_i \cdot h, \; \lambda_i \le 1\)

把上述公式带回微分方程中，并提出一个步长$h$，可得： \(y_{n+1} = y_n + h \sum_{i=1}^r c_i f(x_{n,i}, y(x_{n,i}))\) 最后把$f(x_{n,i}, y(x_{n,i}))$写成$K_i$的形式，即可得： \(\left\{ \begin{aligned} y_{n+1} &= y_n + h \sum_{i=1}^r c_i K_i \\ K_1 &= f(x_n, y_n) \\ K_i &= f(x_n + \lambda_i h, y_n + \sum_{j=1}^{i-1} \mu_{i,j} K_j) \end{aligned} \right.\) 其中$c_i, \lambda_i, \mu_{i,j}$都是待确定的常数。 这就是龙格库塔法的一般形式。

二阶显式R-K方法

二阶的龙格库塔方法使用以下通式： \(\begin{cases} y_{n+1} &= y_n + h(c_1 K_1 + c_2 K_2) \\ K_1 &= f(x_n, y_n) \\ K_2 &= f(x_n+\lambda_2 h, y_n + \mu_{21} h K_1) \end{cases}\) 我们需要指定这些常数的值。

利用局部截断误差的定义，有： \(T_{n+1} = y(x_{n+1}) - y(x_n) - h[c_1 f(x_n, y_n) + c_2 f(x_n + \lambda_2 h, y_n + \mu_{21} h f_n)]\) 接下来使用泰勒展开，将$y(x_{n+1})$和$f$展开到三阶，可得： \(\begin{aligned} T_{n+1} &= (1-c_1-c_2)f_n h \\ &+ (\frac{1}{2} - c_2 \lambda_2) f_x^\prime(x_n, y_n) h^2 \\ &+ (\frac{1}{2} - c_2 \mu_{21}) f_y^\prime(x_n, y_n) h^2 \\ &+ \mathcal O (h^3) \end{aligned}\) 其中$f_n = f(x_n, y_n)$，$f_x^\prime, f_y^\prime$分别代表两个偏导数。

欲使该方法具有二阶精度，那么$h$的次数小于等于二的项必须都为零，这给出以下方程： \(\begin{cases} c_2 \lambda_2 &= \frac{1}{2} \\ c_2 \mu_{21} &= \frac{1}{2} \\ c_1 + c_2 &= 1 \end{cases}\) 所有常数满足以上方程的方法都是二阶龙格库塔法。

值得注意的是，若取 \(c_1 = c_2 = \dfrac{1}{2}, \; \lambda_2 = \mu_{21} = 1\) 得出的公式就是改进欧拉法。

如果进一步提高泰勒展开的阶数，不难注意到如果要求更高阶的截断误差，则需要满足的方程数目会增加，从而出现无解的情况。 因此，该方法不能提供更小的局部截断误差。 如果需要更高阶的精度，则需要提高每一步中估计的次数。

MATLAB中常用的数值解方法ode45使用的曾经就是四阶龙格库塔方法，其截断误差可达$\mathcal O(h^5)$¹。 其阶数正是其名称的由来。

1. https://blogs.mathworks.com/cleve/2014/05/26/ordinary-differential-equation-solvers-ode23-and-ode45 ↩