常微分方程的数值解法

欧拉方法

回到 Euler 方法的基本思想——用差商代替导数——上来。实际上，按照微分中值定理应有

$$\frac{y(x_{n + 1}) - y(x_{n})}{h} = y'(x_{n}+\theta h),0 < \theta < 1$$

注意到方程 $y’ = f(x,y)$ 就有

$$y(x_{n + 1}) = y(x_{n})+hf(x_{n}+\theta h,y(x_{n}+\theta h)) \tag{13}$$

不妨记 $\overline{K}=f(x_{n}+\theta h,y(x_{n}+\theta h))$，称为区间 $[x_{n},x_{n + 1}]$ 上的平均斜率。可见给出一种斜率 $\overline{K}$，（13）式就对应地导出一种算法。

向前 Euler 公式简单地取 $f(x_{n},y_{n})$ 为 $\overline{K}$，精度自然很低。改进的 Euler 公式可理解为 $\overline{K}$ 取 $f(x_{n},y_{n})$，$f(x_{n + 1},\overline{y}_{n + 1})$ 的平均值，其中 $\overline{y}_{n + 1}=y_{n}+hf(x_{n},y_{n})$，这种处理提高了精度。

龙格库塔方法

基本思想

基于欧拉方法的分析启示我们，在区间 $[x_{n},x_{n + 1}]$ 内多取几个点，将它们的斜率加权平均作为 $\overline{K}$，就有可能构造出精度更高的计算公式。这就是龙格 - 库塔方法的基本思想。

确定系数以提高精度

首先不妨在区间$[x_{n},x_{n + 1}]$内仍取 2 个点，仿照（13）式用以下形式试一下

$$\begin{cases} y_{n + 1}=y_{n}+h(\lambda_{1}k_{1}+\lambda_{2}k_{2}) \\ k_{1}=f(x_{n},y_{n}) \\ k_{2}=f(x_{n}+\alpha h,y_{n}+\beta hk_{1}), 0 < \alpha,\beta < 1 \end{cases} \tag{14}$$

其中$\lambda_{1},\lambda_{2},\alpha,\beta$为待定系数，看看如何确定它们使（14）式的精度尽量高。为此我们分析局部截断误差$y(x_{n + 1}) - y_{n + 1}$，因为$y_{n}=y(x_{n})$，所以（14）可以化为

$$\begin{cases} y_{n + 1}=y(x_{n})+h(\lambda_{1}k_{1}+\lambda_{2}k_{2}) \\ k_{1}=f(x_{n},y(x_{n})) = y'(x_{n}) \\ k_{2}=f(x_{n}+\alpha h,y(x_{n})+\beta hk_{1}) \\ \quad=f(x_{n},y(x_{n}))+\alpha hf_{x}(x_{n},y(x_{n})) \\ \quad\quad+\beta hk_{1}f_{y}(x_{n},y(x_{n})) + O(h^{2}) \end{cases} \tag{15}$$

其中$k_{2}$在点$(x_{n},y(x_{n}))$作了 Taylor 展开。（15）式又可表为

$$y_{n + 1}=y(x_{n})+(\lambda_{1}+\lambda_{2})hy'(x_{n})+\lambda_{2}\alpha h^{2}(f_{x}+\frac{\beta}{\alpha}ff_{y})+O(h^{3})$$

注意到

$$y(x_{n + 1})=y(x_{n})+hy'(x_{n})+\frac{h^{2}}{2}y''(x_{n})+O(h^{3})$$

中$y’ = f$，$y’’ = f_{x}+ff_{y}$（多元复合函数求偏导），可见为使误差$y(x_{n + 1}) - y_{n + 1}=O(h^{3})$，只须令

$$\lambda_{1}+\lambda_{2}=1,\lambda_{2}\alpha=\frac{1}{2},\frac{\beta}{\alpha}=1 \tag{16}$$

待定系数满足（16）的（15）式称为 2 阶龙格—库塔公式。由于（16）式有 4 个未知数而只有 3 个方程，所以解不唯一。不难发现，若令$\lambda_{1}=\lambda_{2}=\frac{1}{2}$，$\alpha = \beta = 1$，即为改进的 Euler 公式。可以证明，在$[x_{n},x_{n + 1}]$内只取 2 点的龙格—库塔公式精度最高为 2 阶。

RK方法 : 4 阶龙格—库塔公式

$$\begin{cases} y_{n + 1}=y_{n}+h(\lambda_{1}k_{1}+\lambda_{2}k_{2}+\lambda_{3}k_{3}+\lambda_{4}k_{4}) \\ k_{1}=f(x_{n},y_{n}) \\ k_{2}=f(x_{n}+\alpha_{1}h,y_{n}+\beta_{1}hk_{1}) \\ k_{3}=f(x_{n}+\alpha_{2}h,y_{n}+\beta_{2}hk_{1}+\beta_{3}hk_{2}) \\ k_{4}=f(x_{n}+\alpha_{3}h,y_{n}+\beta_{4}hk_{1}+\beta_{5}hk_{2}+\beta_{6}hk_{3}) \end{cases} \tag{17}$$

其中待定系数$\lambda_{i},\alpha_{i},\beta_{i}$共 13 个，经过与推导 2 阶龙格—库塔公式类似、但更复杂的计算，得到使局部误差$y(x_{n + 1}) - y_{n + 1}=O(h^{5})$的 11 个方程。取既满足这些方程、又较简单的一组$\lambda_{i},\alpha_{i},\beta_{i}$，可得：

$$\begin{cases} y_{n + 1}=\frac{h}{6}(k_{1}+2k_{2}+2k_{3}+k_{4}) \\ k_{1}=f(x_{n},y_{n}) \\ k_{2}=f(x_{n}+\frac{h}{2},y_{n}+\frac{hk_{1}}{2}) \\ k_{3}=f(x_{n}+\frac{h}{2},y_{n}+\frac{hk_{2}}{2}) \\ k_{4}=f(x_{n}+h,y_{n}+hk_{3}) \end{cases} \tag{18}$$

这就是常用的 4 阶龙格—库塔方法（简称 RK 方法）.

线性多步法

多步法的基本思想 、增量函数

以上所介绍的各种数值解法都是单步法，这是因为它们在计算$y_{n + 1}$时，都只用到前一步的值$y_{n}$，单步法的一般形式是

$$y_{n + 1}=y_{n}+h\varphi(x_{n},y_{n},h)\quad (n = 0,1,\cdots,N - 1) \tag{19}$$

其中$\varphi(x,y,h)$称为增量函数，例如 Euler 方法的增量函数为$f(x,y)$，改进 Euler 法的增量函数为

$$\varphi(x,y,h)=\frac{1}{2}[f(x,y)+f(x + h,y + hf(x,y))]$$

如何通过较多地利用前面的已知信息，如$y_{n},y_{n - 1},\cdots,y_{n - r}$，来构造高精度的算法计算$y_{n + 1}$，这就是多步法的基本思想。经常使用的是线性多步法。
让我们试验一下$r = 1$，即利用$y_{n},y_{n - 1}$计算$y_{n + 1}$的情况。

从用数值积分方法离散化方程的（4）式

$$y(x_{n + 1}) - y(x_{n})=\int_{x_{n}}^{x_{n + 1}}f(x,y(x))dx$$

出发，记$f(x_{n},y_{n}) = f_{n}$，$f(x_{n - 1},y_{n - 1}) = f_{n - 1}$，式中被积函数$f(x,y(x))$用二节点$(x_{n - 1},f_{n - 1})$，$(x_{n},f_{n})$的插值公式得到（因$x\geq x_{n}$），所以是外插。

$$\begin{align*} f(x,y(x))&=f_{n}\frac{x - x_{n - 1}}{x_{n}-x_{n - 1}}+f_{n - 1}\frac{x - x_{n}}{x_{n - 1}-x_{n}}\\ &=\frac{1}{h}[(x - x_{n - 1})f_{n}-(x - x_{n})f_{n - 1}] \end{align*} \tag{20}$$

此式在区间$[x_{n},x_{n + 1}]$上积分可得

$$\int_{x_{n}}^{x_{n + 1}}f(x,y(x))dx=\frac{3h}{2}f_{n}-\frac{h}{2}f_{n - 1}$$

于是得到

$$y_{n + 1}=y_{n}+\frac{h}{2}(3f_{n}-f_{n - 1}) \tag{21}$$

注意到插值公式（20）的误差项含因子$(x - x_{n - 1})(x - x_{n})$，在区间$[x_{n},x_{n + 1}]$上积分后得出$h^{3}$，故公式（21）的局部截断误差为$O(h^{3})$，精度比向前 Euler 公式提高 1 阶。

若取$r = 2,3,\cdots$可以用类似的方法推导公式，如对于$r = 3$有

$$y_{n + 1}=y_{n}+\frac{h}{24}(55f_{n}-59f_{n - 1}+37f_{n - 2}-9f_{n - 3}) \tag{22}$$

其局部截断误差为$O(h^{5})$。

如果将上面代替被积函数$f(x,y(x))$用的插值公式由外插改为内插，可进一步减小误差。内插法用的是$y_{n + 1},y_{n},\cdots,y_{n - r + 1}$，取$r = 1$时得到的是梯形公式，取$r = 3$时可得

$$y_{n + 1}=y_{n}+\frac{h}{24}(9f_{n + 1}+19f_{n}-5f_{n - 1}+f_{n - 2}) \tag{23}$$

与（22）式相比，虽然其局部截断误差仍为$O(h^{5})$，但因它的各项系数（绝对值）大为减小，误差还是小了。当然，（23）式右端的$f_{n + 1}$未知，需要如同向后 Euler 公式一样，用迭代或校正的办法处理。