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

1/ 欧拉法（Euler Method）²

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

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

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

^{假设有一个 ODE：

y′(t)=f(t,y(t)),y(t0)=y0,y^{\prime}(t)=f(t,y(t)),\quad y(t_0)=y_0, y′(t)=f(t,y(t)),y(t0)=y0,

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

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

^{Euler 方法的迭代公式为：

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

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

^{2.1/ 四阶 Runge-Kutta 方法}

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

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

dydt=y′(t)=f(t,y),y(t0)=y0\frac{dy}{dt}=y^{\prime}(t)=f(t,y),\quad y(t_0)=y_0 dtdy=y′(t)=f(t,y),y(t0)=y0

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

yn+1=yn+16h(k1+2k2+2k3+k4)tn+1=tn+h\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} yn+1tn+1=yn+61h(k1+2k2+2k3+k4)=tn+h

其中，

k1=f(tn,yn),k2=f(tn+h2,yn+hk12),k3=f(tn+h2,yn+hk22),k4=f(tn+h,yn+hk3).\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} k1k2k3k4=f(tn,yn),=f(tn+2h,yn+h2k1),=f(tn+2h,yn+h2k2),=f(tn+h,yn+hk3).}

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

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

^{k1k_1k1 是间隔开始处的斜率，依赖于 yyy（这第一步是欧拉法）；}
^{k2k_2k2 是间隔中点处的斜率，依赖于 yyy 和 k1k_1k1；}
^{k3k_3k3 依旧是间隔中点处的斜率，只是依赖于 yyy 和 k2k_2k2 了；}
^{k4k_4k4 是间隔终点处的斜率，依赖于 yyy 和 k3k_3k3.}

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

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

^{2.2/ Runge-Kutta 的一般形式}

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

yn+1=yn+h∑j=1sbzkiy_{n+1}=y_{n}+h \sum_{j=1}^{s} b_{\mathbf{z}} k_{i} yn+1=yn+hj=1∑sbzki

其中

k1=f(tn,yn)k2=f(tn+c2h,yn+h(a21k1)),k3=f(tn+c3h,yn+h(a31k1+a32k2)),⋮ks=f(tn+csh,yn+h(as1k1+as2k2+⋯+as,s−1ks−1)).\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} k1k2k3ks=f(tn,yn)=f(tn+c2h,yn+h(a21k1)),=f(tn+c3h,yn+h(a31k1+a32k2)),⋮=f(tn+csh,yn+h(as1k1+as2k2+⋯+as,s−1ks−1)).

从公式中可以知道，如果要指定某个特定的方法，就需要提供一个整数 sss（代表级数的项数），还有系数 aij,for 1≤j<i≤sa_{ij},~\text{for }1\le j<i\le saij, for 1≤j<i≤s, bi,for i=1,2,…,sb_i,~\text{for }i=1,2,\dots,sbi, for i=1,2,…,s 和 cifor i=2,3,…,sc_i~\text{for }i=2,3,\dots,sci for i=2,3,…,s.}

^{其中矩阵 [aij][a_{ij}][aij] 称为 Runge-Kutta 矩阵，而 bib_ibi 和 cic_ici 是权重和节点。}

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

∑i=1sbi=1\sum_{i=1}^s b_i=1 i=1∑sbi=1}

^参考

^{知乎：常微分法数值计算方法 ↩︎}
^{Wikipedia: Euler method ↩︎}