常微分方程（ODE）的数值计算方法

常微分方程组的求解比较麻烦，通常在计算机上使用数值计算的方式去进行。

假设一阶常微分方程组（ODEs）由下式给出

$$ \frac{dx}{dt}=f_i(x),~i=1,2,\dots,n $$

其中 $x$ 是一个变量构成的向量 $x=[x_1,x_2,\dots,x_n]$。任意高阶微分方程都可以通过变量替换变成一阶微分方程组。

求解 ODE 的数值方法有很多种：欧拉法（Euler）、龙格-库塔法（Runge-Kutta）、Adams 线性多步法、Gear 法等 ¹。

1/ 欧拉法（Euler Method）²

在数学和计算科学中，欧拉方法（也叫前向欧拉方法）是一种用于求解给定初值的 ODE 的一阶数值过程，是对常微分方程进行数值积分的最基本的显式方法，也是最简单的 Runge-Kutta 方法。

欧拉方法以列昂哈德·欧拉（Leonhard Euler）的名字命名，在他的《计算积分学》（1768-1870年出版）一书中首次提出。

欧拉方法是一阶方法，它的局部误差（每步误差）与步长的平方成正比，全局误差（给定时间的误差）与步长成正比。 欧拉方法通常作为构造更复杂方法的基础，例如预测器-校正器（Predictor-Corrector）方法。

假设有一个 ODE：

$$ y^{\prime}(t)=f(t,y(t)),\quad y(t_0)=y_0, $$

尽管曲线是未知的，但是起点是已知的。根据微分方程，可以计算出 $A_0$ 点的斜率，因此可以计算出切线，我们根据该切线走一小步到 $A_1$ 点。然后再根据 $A_1$ 点得到该位置的斜率，根据这个斜率走一小步去下一个点 $A_2$......

如此反复，如果步长取足够的小，那么得到的结果将会很接近确切解。

Euler 方法的迭代公式为：

$$ \begin{aligned} y_{n+1}&=y_n+hf(t_n,y_n)\\ t_{n+1}&=t_n+h \end{aligned} $$

2/ 龙格-库塔法（Runge-Kutta Method）

2.1/ 四阶 Runge-Kutta 方法

在数值分析中，Runge-Kutta 方法是一组隐式和显式迭代方法，其中包括欧拉方法，用于时间离散化中的常微分方程的近似解。由国数学家卡尔 • 朗格（Carl Runge）和威廉 • 库塔（Wilhelm Kutta）于 1900 年左右开发。

设一个初值问题如下定义：

$$ \frac{dy}{dt}=y^{\prime}(t)=f(t,y),\quad y(t_0)=y_0 $$

其中 $y$ 是一个关于时间 $t$ 的未知函数（标量或向量），也是需要近似的函数。已知 $y^{\prime}(t)$ 是关于 $t$ 和 $y$ 自身的函数，在时间 $t_0$ 对应的 $y$ 值为 $y_0$。 函数 $f$ 和初始条件 $t_0,y_0$ 已经给出。选择一个步长 $h>0$ 并且对于 $n=0,1,2,3,\dots,$ 定义：

$$ \begin{aligned} y_{n+1}&=y_{n}+\frac{1}{6} h\left(k_{1}+2 k_{2}+2 k_{3}+k_{4}\right)\\ t_{n+1}&=t_{n}+h \end{aligned} $$

其中，

$$ \begin{aligned} k_{1}&=f\left(t_{n}, y_{n}\right),\\ k_{2}&=f\left(t_{n}+\frac{h}{2}, y_{n}+h \frac{k_{1}}{2}\right),\\ k_{3}&=f\left(t_{n}+\frac{h}{2}, y_{n}+h \frac{k_{2}}{2}\right),\\ k_{4}&=f\left(t_{n}+h, y_{n}+h k_{3}\right). \end{aligned} $$

在不同的资料中可能会有不同的但是等价的公式表示。

其中 $y_{n+1}$ 为 $y(t_{n+1})$ 的四阶 Runge-Kutta （RK4）近似，下一个值 $(y_{n+1})$ 由当前值 $(y_n)$ 加上四个增量的加权平均确定，每个增量为间隔大小 $h$ 和微分方程右侧的函数 $f$ 指定的估计斜率的乘积。

• $k_1$ 是间隔开始处的斜率，依赖于 $y$（这第一步是欧拉法）；
• $k_2$ 是间隔中点处的斜率，依赖于 $y$ 和 $k_1$；
• $k_3$ 依旧是间隔中点处的斜率，只是依赖于 $y$ 和 $k_2$ 了；
• $k_4$ 是间隔终点处的斜率，依赖于 $y$ 和 $k_3$.

RK4 是四阶的方法，意味着局部截断误差的阶数为 $\mathcal{O}(h^5)$，总积累误差的阶数为 $\mathcal{O}(h^4)$.

在实际应用的时候，特别是在物理上，函数 $f$ 与 $t$ 无关，即所谓的自治系统（Autonomous system）或者时不变系统，它们的增量不会被计算，也不会传递给函数 $f$，而仅仅使用 $t_{n+1}$ 的最终公式。

2.2/ Runge-Kutta 的一般形式

上面介绍的是 RK4，它只是 Runge-Kutta 的其中一种。实际上，Runge-Kutta 有很多种，如果把它写成一般的形式，有

$$ y_{n+1}=y_{n}+h \sum_{j=1}^{s} b_{\mathbf{z}} k_{i} $$

其中

$$ \begin{aligned} k_{1} &=f\left(t_{n}, y_{n}\right) \\ k_{2} &=f\left(t_{n}+c_{2} h, y_{n}+h\left(a_{21} k_{1}\right)\right), \\ k_{3} &=f\left(t_{n}+c_{3} h, y_{n}+h\left(a_{31} k_{1}+a_{32} k_{2}\right)\right), \\ & \vdots \\ k_{s} &=f\left(t_{n}+c_{s} h, y_{n}+h\left(a_{s 1} k_{1}+a_{s 2} k_{2}+\cdots+a_{s, s-1} k_{s-1}\right)\right) . \end{aligned} $$

从公式中可以知道，如果要指定某个特定的方法，就需要提供一个整数 $s$（代表级数的项数），还有系数 $a_{ij},~\text{for }1\le j<i\le s$, $b_i,~\text{for }i=1,2,\dots,s$ 和 $c_i~\text{for }i=2,3,\dots,s$.

其中矩阵 $[a_{ij}]$ 称为 Runge-Kutta 矩阵，而 $b_i$ 和 $c_i$ 是权重和节点。

Taylor 级数展开式表明，Runge-Kutta 方法是一致的，当且仅当

$$ \sum_{i=1}^s b_i=1 $$

参考