常微分方程初值问题: Runge-Kutta 方法

Taylor 方法

预备知识: 二维函数的 Taylor 展开.

假设对于初值问题\[ y’=f(t,y),\quad 0\leq t\leq L,\quad y(0)=y_0,\tag{1}\]足够光滑, 那么它的 Taylor 展开具有如下形式:\[ y(t_{i+1})=y(t_i)+hy'(t_i)+\frac{h^2}{2}y^{\prime\prime}(t_i) +\cdots+\frac{h^n}{n!}y^{(n)}(t_i)+\frac{h^{n+1}}{(n+1)!}y^{(n+1)}(\xi_i). \tag{2}\]通过微分方程 (1),\[\begin{align*} y’&=f,\\ y^{\prime\prime}&=\frac{d}{dt} f(t,y)=f_t+f_yy’=f_t+f_yf,\\ y^{\prime\prime\prime}&=\frac{d^2}{dt^2} f(t,y)=f_{tt}+f_{ty}f+(f_{ty}+f_{yy}f)f+f_y(f_t+f_yf),\\ & \cdots \\ y^{(k)} &= \frac{d^{k-1}}{dt^{k-1}}f(t,y).\end{align*}\]将上面的等式代入 (2) 并去掉余项就得到了高阶 Taylor 方法:\[ \begin{equation} \begin{split} y_{i+1}&=y_i+hf(t_i,y_i)+\frac{h^2}{2}\left[ \frac{d}{dt} f(t,y) \right]_{|t=t_i, y=y_i} +\cdots+\frac{h^n}{n!}\left[ \frac{d^{n-1}}{dt^{n-1}} f(t,y) \right]_{|t=t_i, y=y_i} \\ &=y_i+h \left[ f(t_i,y_i)+\frac{h}{2} \frac{d}{dt} f(t,y) +\cdots+\frac{h^{n-1}}{n!}\frac{d^{n-1}}{dt^{n-1}} f(t,y) \right]_{|t=t_i, y=y_i}, \end{split} \end{equation} \tag{3}\]并且容易验证高阶 Taylor 方法 (3) 具有局部截断误差 (高阶导数一致有界)\[ \tau_{i+1} = \frac{y(t_{i+1})-y_{i+1}}{h} = O(h^{n}).\]这里的局部误差表明了具体每一步计算产生的新的误差, 通常依赖于微分方程, 步长大小以及数值格式. 当 \(n=1\) 时就是最基础的 Euler 方法.

二阶 Runge-Kutta 方法

由上一节可知 Taylor 方法具有高阶局部误差, 因此理论上可以达到高阶收敛. 但是, 它的实现需要计算右端函数的高阶导数, 这样的过程通常是极其复杂的, 从而很少在实际应用中使用. 我们以二阶方法为例阐述 Runge-Kutta 方法的基本思想, 它可以避免计算高阶导数但是仍可以达到高阶收敛. Runge-Kutta 方法最早是由 Carl Runge (卡尔龙格) 在 1895 提出的, 后来由 Martin Wilhelm Kutta 在 1901 年推广到微分方程组的求解上. 他们提出的方法最后发展成现在成熟的 Runge-Kutta 方法. 关于 Runge 的生平, 读者可以参考另外一篇文章: 数值分析中的数学家 – 龙格(C. Runge).

修正的 Euler 法

首先, 我们希望找到合适的参数 \(a_1,\alpha_1, \beta_1\) 使得\[ a_1 f(t+\alpha_1 h, y+\beta_1 f(t,y)) – \left[ f(t,y)+\frac{h}{2} \frac{d}{dt} f(t,y) \right]=O(h^2). \tag{4}\]从而, 当 \(n=2\) 时, (3) 可以写为\[ y_{i+1} = y_i+ h \left[ a_1 f(t_i+\alpha_1 h, y_i+\beta_1f(t_i,y_i)) \right]. \tag{5}\]从而得到一个特殊二阶 Runge-Kutta 方法–中点方法 (Midpoint method).

下面, 我们来具体证明存在合适的参数使得 (4) 成立. 直接微分可得\[ f+\frac{h}{2} \frac{d}{dt} f = f+\frac{h}{2} f_t + \frac{h}{2} f_y f. \tag{6}\]另一方面, 由二维函数的 Taylor 展开公式可知 (一阶展开)\[ a_1 f(t+\alpha_1 h, y+\beta_1f) = a_1 f + a_1 \alpha_1 hf_t + a_1 \beta_1f_y f + a_1 R_1. \tag{7}\]其中, 如果我们假设函数足够光滑 (任意导数一致有界), 那么 \(R_1=O(\alpha_1^2h^2+\beta_1^2)\). 对比 (6)(7) 两个等式中的系数, 得到三个方程:\[ a_1=1;\quad a_1 \alpha_1=\frac{1}{2};\quad a_1 \beta_1=\frac{h}{2}.\]求解可知 \(a_1=1, \alpha_1=\frac{h}{2}, \beta_1=\frac{h}{2}\). 由 (7) 中的残差可知, 这些参数满足 (4), 证明结束. 从而得到中点方法或修正的 Euler 法:\[ y_{i+1} = y_i+ h f \left( t_i+\frac{h}{2}, y_i+\frac{h}{2} f(t_i,y_i) \right).\]

改进的 Euler 法

我们希望找到合适的参数 \(a_1,a_2, c_2,e_2\) 使得\[ a_1 f(t,y)+a_2 f(t+c_2h, y+e_{2}hf(t,y)) – \left[ f(t,y)+\frac{h}{2} \frac{d}{dt} f(t,y) \right]=O(h^2).\]由二维函数的展开公式, 我们得到\[ a_1 f+a_2 f(t+c_2h, y+e_{2}hf) =(a_1+a_2) f+ a_2c_2h f_t +a_2e_2h f_yf + ….\tag{8}\]对比 (6)(8) 两个等式中的系数, 得到三个方程:\[ a_1+a_2=1,\quad a_2c_2=\frac{1}{2},\quad a_2e_2=\frac{1}{2}.\]三个方程中包含四个未知数, 因此我们对于参数有多种选择, 通常由如下的参数值\[ a_1=a_2=\frac{1}{2},\quad c_2= e_2=1,\]得到的方法称为改进的 Euler 方法:\[ y_{i+1} = y_i+ \frac{h}{2} \left[ f(t_i,y_i)+ f(t_i+h, y_i+hf(t_i,y_i)) \right]. \tag{9}\]

高阶 Runge-Kutta 方法

由二阶方法的求解过程可以知道, Runge-Kutta 方法的基本思想是利用合适的线性组合来逼近对应的 Taylor 方法. 我们考虑经典的四阶方法, 寻找参数 \(a_i, c_i,e_{ij}\) 使得\[\begin{align*} k_1 &= f(t,y)\\ k_2 &= f(t+c_2h, y+he_{21}k_1)\\ k_3 &= f(t+c_3h, y+h(e_{31}k_1+e_{32}k_2)),\\ k_4 &= f(t+c_4h, y+h(e_{41}k_1+e_{42}k_2+e_{43}k_3)),\\ \sum_{i=1}^4 a_ik_i &- \left[ f(t_i,y_i)+\frac{h}{2} \frac{d}{dt} f(t,y) +\cdots+\frac{h^{3}}{4!}\frac{d^{3}}{dt^{3}} f(t,y) \right] = O(h^4)\end{align*}\]最后一个等式的左边第二项就是高阶 Taylor 方法的前四项, 参考 (3). 在通过系数对比并确定参数之后就得到如下的四阶方法:\[ y_{i+1}=y_i+h\sum_{i=1}^4 a_ik_i.\]上述推导过程相当繁琐, 有兴趣的读者可以作为练习自己推导一下. 最经典的四阶方法 (RK4) 为:\[ y_{i+1}=y_i+\frac{h}{6}(k_1+2 k_2+2 k_3+k_4), \tag{10}\]其中,\[\begin{align*} k_1 &= f(t_i,y_i),\\ k_2 &= f \left( t_i+\frac{h}{2}, y_i+\frac{h}{2} k_1 \right),\\ k_3 &= f \left( t_i+\frac{h}{2}, y_i+ \frac{h}{2} k_2 \right),\\ k_4 &= f(t_i+h, y_i+hk_3),\end{align*}\]

RK4 的几何意义

我们考虑区间 \([t_i,t_{i+1}]\) 上的函数 \(y(t)\), 由积分基本定理可知\[ y(t_{i+1}) = y(t_i) + \int_{t_i}^{t_{i+1}} f(t,y)dt.\tag{11}\]利用数值积分方法–Simpson 公式–可以得到积分的一个近似值如下\[ \int_{t_i}^{t_{i+1}} f(t,y)dt \simeq \frac{h}{6} \left( f(t_i,y(t_i))+4f(t_i+\frac{h}{2},y(t_i+\frac{h}{2}))+f(t_{i+1},y(t_{i+1})) \right). \tag{12}\]我们知道在 RK4 方法中, \(k_1,k_2,k_3,k_4\) 分别是曲线在不同位置斜率的近似值, 具体来说, \(k_1\) 是左端点处的近似, \(k_2,k_3\) 为区间中点处的近似, \(k_4\) 是右端点处的近似. 如果我们在 (12) 中令 \(k_1=f(t_i,y(t_i))\), \(k_4=f(t_{i+1},y(t_{i+1}))\), 最后令 \(k_2+k_3=2f(t_i+\frac{h}{2},y(t_i+\frac{h}{2}))\), 那么方程 (11) 可以写成\[ y(t_{i+1}) \simeq y(t_i) + \frac{h}{6} \left( k_1+2k_2+2k_3+k_4 \right).\]它具有与 RK4 方法相同的形式.

总结

由前面的介绍我们了解到 Runge-Kutta 方法的基本思想是利用合适的 “导数” 的线性组合来逼近高阶 Taylor 方法得到的改进的高阶方法, 它相比于 Taylor 方法的优点是不需要实际计算 \(f(t,y)\) 的各阶导数, 而只需要计算函数 \(f\) 在特定点处的值. 当利用 Euler 方法 (一阶) 时, 需要计算的次数为 1, 这里的截断误差为 \(O(h)\). 当利用四阶方法时, 需要的计算次数为 4, 对应的截断误差为 \(O(h^4)\). 但是, 如何比较它们的优劣呢? 一个直接的方法是计算如下情况下的误差:

1. 利用步长 \(h\), Euler 方法计算 4n 次;
2. 利用步长 \(2h\), 改进或修正的 Euler 方法计算 2n 次;
3. 利用步长 \(4h\), 四阶方法 RK4 计算 n 次.

上述情况下的计算次数都 4n 次, 如果 RK4 得到的误差较小, 那么就说明四阶方法是优于 Euler 和其它方法的. 我们利用下面的数值实验来回答这个问题.

例 1. 考虑如下问题\[ y’=y-t^2+1,\quad 0\leq t\leq 1, \quad y(0)=0.5.\]我们在下面的表格中列出了四种方法在不同步长情况下的计算结果, 从中可以看出四阶 Runge-Kutta 方法在相同量级的计算复杂度的情况下得到最精确的数值结果, 因此是优于其它方法的.

参考

[1] Burden, R. L., & Faires, J. D. (2010). Numerical Analysis. In Cengage Learning (9th ed.).
[2] 常微分方程初值问题: 若干数学模型.
[3] Python 科学计算: NumPy.
[4] Python 科学计算: SymPy.

附录

二维函数的 Taylor 展开

假设 \(f(t,y)\) 为区域 \(\Omega=[a,b]\times[c,d]\) 上 的\(n+1\) 阶连续可微函数. 令 \((t_0,y_0)\in \Omega\), 那么对任意 \((t,y)\in \Omega\), 都还在 \(\xi, \mu\) 使得如下等式成立\[ f(t,y)=T_n(t,y)+\frac{1}{(n+1)!}\sum_{j=0}^{n+1} \begin{pmatrix} n+1\\j \end{pmatrix} (t-t_0)^{n+1}(y-y_0)^j \frac{\partial^{n+1} f}{\partial t^{n+1-j} \partial y^j} (\xi,\mu),\]其中,\[\begin{align*} T_n(t,y) =& f(t_0,y_0)+ \left[ (t-t_0)\frac{\partial f}{\partial t} (t_0,y_0)+(y-y_0)\frac{\partial f}{\partial y} (t_0,y_0) \right] +\cdots \\ &+\frac{1}{n!}\sum_{j=0}^n \begin{pmatrix} n\\j \end{pmatrix} (t-t_0)^{n-j}(y-y_0)^j \frac{\partial^n f}{\partial t^{n-j}\partial y^j}(t_0,y_0).\end{align*}\]

Runge-Kutta 方法的 Python 实现

# Runge-Kutta method
import numpy as np
import matplotlib.pyplot as plt

class ode_sol:
    def __init__(self):
        self.y = None
        self.err = None

class my_ode_ivp:
    def __init__(self, fun, t_span, y0, t_eval=None, method='Euler'):
        self.N = 100    # if t_eval is none
        self.fun = fun
        self.t_span = t_span
        self.y0 = y0
        self.sol = ode_sol()

        if t_eval is None:
            self.t_eval = np.linspace(t_span[0], t_span[1], self.N)
        else:
            self.t_eval = t_eval

        self.sol.y = np.zeros(self.t_eval.shape[0])
        self.sol.y[0] = y0

        if method == 'Euler':
            self.euler()
        elif method == 'ModifiedEuler':
            self.modified_euler()
        elif method == 'ImprovedEuler':
            self.improved_euler()
        elif method == 'RK4':
            self.rk4()
        else:
            pass

    def euler(self):
        for i in range(1, self.t_eval.shape[0]):
            self.sol.y[i] = self.sol.y[i-1] + (self.t_eval[i]-self.t_eval[i-1])*self.fun(self.t_eval[i-1], self.sol.y[i-1])
    
    def modified_euler(self):
        for i in range(1, self.t_eval.shape[0]):
            h = self.t_eval[i]-self.t_eval[i-1]
            k1 = self.fun(self.t_eval[i-1], self.sol.y[i-1])
            self.sol.y[i] = self.sol.y[i-1] + h*self.fun(self.t_eval[i-1]+.5*h, self.sol.y[i-1]+.5*h*k1)
            
    def improved_euler(self):
        for i in range(1, self.t_eval.shape[0]):
            h = self.t_eval[i]-self.t_eval[i-1]
            k1 = self.fun(self.t_eval[i-1], self.sol.y[i-1])
            k2 = self.fun(self.t_eval[i-1]+h, self.sol.y[i-1]+h*k1)
            self.sol.y[i] = self.sol.y[i-1]+.5*h*(k1+k2)

    def rk4(self):
        for i in range(1, self.t_eval.shape[0]):
            h = self.t_eval[i]-self.t_eval[i-1]
            k1 = self.fun(self.t_eval[i-1], self.sol.y[i-1])
            k2 = self.fun(self.t_eval[i-1]+.5*h, self.sol.y[i-1]+.5*h*k1)
            k3 = self.fun(self.t_eval[i-1]+.5*h, self.sol.y[i-1]+.5*h*k2)
            k4 = self.fun(self.t_eval[i-1]+h, self.sol.y[i-1]+h*k3)
            self.sol.y[i] = self.sol.y[i-1]+h*(k1+2*k2+2*k3+k4)/6


# setting for ode: f(t,y), y0, t_span, t_eval
# two test examples: ex1, ex2
def f(t, y):
    return y-t**2+1 #ex1
    # return -t*y/(1-t*t) #ex2

# analytic solution
def yfun(t):
    return (t+1)**2 - .5*np.exp(t) #ex1
    # return np.sqrt(1-t*t) #ex2

y0 = .5 #ex1
# y0 = 1. #ex2
t_span = (0, 1.)


N = 11
t_eval = np.linspace(t_span[0], t_span[1], N)

# solve, compute error and plot
ode = my_ode_ivp(fun=f, t_span=t_span, y0=y0, t_eval=t_eval, method='RK4')
# 'RK4', 'Euler', 'ModifiedEuler', 'ImprovedEuler'

y_exact = yfun(ode.t_eval)
print(y_exact[-1], ode.sol.y[-1])
ode.sol.err = abs(ode.sol.y[-1]-y_exact[-1])
print(ode.sol.err)

plt.plot(ode.t_eval, y_exact,'.-', label='Exact sol')
plt.plot(ode.t_eval, ode.sol.y, 'o-', label='Euler sol')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.show()