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

Runge-Kutta 方法是求解常微分方程初值问题最经典的方法之一, 自适应 Runge-Kutta 通常称为 Runge-Kutta-Fehlberg 方法, 它是由 Erwin Fehlberg 在十九世纪 60 年代为 NASA 工作时提出的一系列误差控制方法的统称. 本文主要介绍这种通过调整步长大小将初值问题的数值方法的局部截断误差控制在一定范围内的重要思想和方法.

问题表述

考虑初值问题\[ y’=f(t,y),\quad 0\leq t\leq L,\quad y(0)=y_0,\tag{1}\]的单步方法 (如典型的 Runge-Kutta) 具有如下形式:\[ y_{i+1}=y_i+h_i \phi(t_i, y_i, h_i),\quad i=0,1,\cdots,N-1, \tag{2}\]其中, 步长 \(h_i\) 不是固定不变的. 如果网格大小是一致的, RK4 方法的局部截断误差为 \(O(h^4)\), 也就是说存在一个正常数 \(C>0\) 满足\[ \tau_{i+1}=\frac{|y(t_{i+1})-y(t_i)-h \phi(t_i, y(t_i), h)|}{h}\approx Ch^4.\]如果我们计算的 \(y_i\) 足够精确, 那么可以认为 \(y_i=y(t_i)\), 从而局部误差可以表述如下\[ \tau_{i+1}(h)=\frac{|y(t_{i+1})-y_{i+1}|}{h}\approx Ch^4. \tag{3}\]

现在, 一个有趣且重要的问题是:

问题 1. 给定一个误差上限 \(\epsilon>0\), 怎样找到最少的网格点使得全局误差 \(|y(t_i)-y_i|\leq \epsilon\)?

如何理解上面问题的 “有趣和重要” 呢? 我们有必要考虑一下背后的动机.

• 在实际应用, 如热传输问题 (见附录) 中, 初值问题 (1) 通常表示一个向量值问题, 即 \(y\in \mathbb{R}^N\), 其中整数 \(N\) 依赖于问题的规模.
• 当问题规模比较大时, 每计算一次函数 \(f\) 的值就变得相当耗时, 进而影响整个数值离散方法的计算效率.
• 在问题规模比较大时, 如何在保证数值解满足一定精度的前提下减少 \(f\) 的计算次数从而提高算法效率呢? 答案就是采用自适应的 Runge-Kutta-Fehlberg 方法.

对问题 1的回答对数值算法本身和实际应用都具有重要的意义.

RKF 方法概述

现假设 RK 方法 (2) 具有精度 \(O(h^n)\), 另外我们还有一个高一阶的 RK 方法\[ w_{i+1} = w_i + h \psi(t_i, w_i,h).\tag{4}\]首先, 我们需要假设 \(y_i \approx w_i \approx y(t_i)\), 然后考虑两种方法下一步时间迭代产生的误差 (时间步长依旧为 \(h\)). 容易知道\[ \tau_{i+1}(h) = \frac{y(t_{i+1})-y_{i+1}}{h} = O(h^n),\quad \tilde{\tau}_{i+1}(h) = \frac{y(t_{i+1})-w_{i+1}}{h} = O(h^{n+1}),\]从而,\[ \tau_{i+1}(h) = \frac{y(t_{i+1})-w_{i+1}}{h} + \frac{w_{i+1}-y_{i+1}}{h}=O(h^{n+1})+\frac{w_{i+1}-y_{i+1}}{h}\approx \frac{w_{i+1}-y_{i+1}}{h}.\]因此, RK 方法的局部截断误差可以很容易地近似计算出来. 我们计算这个误差有什么用呢, 下面可以看到它对新的步长选取非常关键. 具体地, 如果用 \(qh\) 代替 (3) 中的 \(h\) 可以得到如下关系式\[ \tau_{i+1}(qh) \approx C(qh)^n =q^n(Ch^n) \approx q^n \tau_{i+1}(h) \approx q^n\frac{w_{i+1}-y_{i+1}}{h}.\]给定精度 \(\epsilon>0\), 就得到新的步长 \(qh\) 满足的不等式 \(q^n\frac{|w_{i+1}-y_{i+1}|}{h} \leq \epsilon\), 求解可得\[ q\leq \left( \frac{\epsilon h}{|w_{i+1}-y_{i+1}|} \right)^{1/n}.\tag{5}\]

想要实现新步长的选取, 我们必须使用新的高一阶的 RK 方法计算对应的数值解. Fehlberg 给出的方法 [1] RK 方法的非常精妙, 因为两种 RK 方法共享部分 \(k_i\), 减少了近 40% 的计算量. 这里我们列出经典的 RKF45:\[ \begin{align} y_{i+1}&=y_i+\frac{25}{216}k_1+\frac{1408}{2565}k_3+\frac{2197}{4104}k_4-\frac{1}{5}k_5, && (O(h^4))\\ w_{i+1}&=y_i+\frac{16}{135}k_1+\frac{6656}{12825}k_3+\frac{28561}{56430}k_4-\frac{9}{50}k_5+\frac{2}{55}k_6 && (O(h^5)) \end{align}\tag{6}\]其中\[\begin{align*} k_1 &= hf(t_i,y_i),\\ k_2 &= hf \left( t_i+\frac{h}{4}, y_i+\frac{1}{4}k_1 \right),\\ k_3 &= hf \left( t_i+\frac{3h}{8}, y_i+\frac{3}{32}k_1 +\frac{9}{32}k_2 \right),\\ k_4 &= hf \left( t_i+\frac{12h}{13}, y_i+\frac{1932}{2197}k_1-\frac{7200}{2197}k_2 + \frac{7296}{2197}k_3 \right),\\ k_5 &= hf \left( t_i+h, y_i+\frac{439}{216}k_1-8k_2 + \frac{3680}{513}k_3-\frac{845}{4104}k_4 \right),\\ k_6 &= hf \left( t_i+\frac{h}{2}, y_i-\frac{8}{27}k_1+2k_2 – \frac{3544}{2565}k_3+\frac{1859}{4104}k_4-\frac{11}{40}k_5 \right).\end{align*}\]

算法实现

现在我们总结 RKF45 方法的算法实现如下. 给定条件: [0,L]; 初值 y0; 精度 \(\epsilon\); 最大最小步长 hmax, hmin.

1. 当前节点为 \(t_i,y_i\), 步长为 \(h\);
2. 计算 \(k_i,i=1,2,\cdots,6\), 并利用 (6) 中两个方法的差值 \[ D=w_{i+1}-y_{i+1}=\frac{1}{360}k_1-\frac{128}{4275}k_3-\frac{2197}{75240}k_4+\frac{1}{50}k_5+\frac{2}{55}k_6. \] 这里我们需要注意, Fehlberg 在其论文 [1, (30)] 中给出的这个差值是有个书写错误的, 2197 被写成了 2187.
3. 如果 \(D/h<\epsilon\), 说明当前步长是合适的, 利用 RK4 计算出 \(y_{i+1}\) 并更新时间 \(t_i \leftarrow t_i+h\).
4. 利用 \(q=(\epsilon/D)^{1/4}\) 得到新的步长为 \(qh\); 如果 \(t_i+qh>L\), \(h=L-t_i\).
5. 重复步骤 1-4 直到达到区间的右端点 \(L\).

改进的 RKF45

我们知道在经典的 RKF45 方法中, 当步长合适的时候, 数值解是由 RK4 计算得到. 现在的问题是为什么不利用 RK5 来计算数值解呢? 理由如下:

1. RK5 是 RK4 的高阶方法, 得到的解理论上更加精确;
2. 计算 RK5 对应的解 \(w_{i+1}\) 非常简单, 只需要稍微对经典的 RKF45 算法进行修改. 具体地, 因为我们已经计算出了 \(D=w_{i+1}-y_{i+1}\), 并且也已经知道 \(w_{i+1} = y_{i+1}+D\), 所以只需在算法第 3 步中 \(y_{i+1}\) 的基础上加上 \(D\) 的值.

数值算例

例 1. 考虑初值问题\[ y’=y-t^2+1,\quad 0\leq t\leq 1, \quad y(0)=0.5.\]

在数值计算中, 我们给定初始条件为: 精度 \(\epsilon=10^{-5}, hmin=0.001, hmax=0.25\).下面两个表格中分别列出了经典的 RKF45 和改进的 RKF45 的数值结果. 可以看出改进的方法得到了更加精确的结果.

参考

[1] Erwin Fehlberg (1969) Low-order classical Runge-Kutta formulas with stepsize control and their application to some heat transfer problems. Vol. 315. National aeronautics and space administration.
[2] Erwin Fehlberg (1968) Classical fifth-, sixth-, seventh-, and eighth-order Runge-Kutta formulas with step size control. NASA Technical Report 287.
[3] Burden, R. L., & Faires, J. D. (2010). Numerical Analysis. In Cengage Learning (9th ed.).
[4] 常微分方程初值问题: Runge-Kutta 方法.

附录

热传导方程

为简化表述, 考虑区间 \((0,1)\) 上的一维热传导方程\[ \begin{align*} \frac{\partial u}{\partial t} &=\frac{\partial^2 u}{\partial x^2} , && (\text{方程})\\ u(x,t) &= u_0(x), && (\text{初值})\\ u(0,t) &= u(1,t) = 0. && (\text{边值}) \end{align*}\]上面这个方程是最简单的一种抛物偏微分方程. 前向 Euler 方法 (最简单的显式方法) 离散格式如下:\[ \frac{u_{i,j+1}-u_{i,j}}{k} = \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2}. \tag{A.1}\]其中, \(k,h\) 分别是时间和空间步长. 可以证明, 为保证方法 (A.1) 稳定, 时间和空间步长必须满足如下条件:\[ \frac{k}{h^2}\leq \frac{1}{2}.\tag{A.2}\]稳定性条件 (A.2) 告诉我们, 如果想要利用 Euler 方法 (A.1) 得到一个令人满意的数值解, 通常需要足够小的空间步长 \(h\); 这也意味着时间步长 \(k\) 必须足够小, 关于时间的迭代就会更加耗时. 一个可行的方法就是采用比 Euler 法更高阶的 Runge-Kutta 方法. 具体来说就是求解向量 \(\textbf{u}(t)\in \mathbb{R}^n\):\[ \frac{d u_i}{dt} = \frac{u_{i+1}-2u_{i}+u_{i-1}}{h^2} \tag{A.3}\]为了更容易地利用稳定性条件 (A.2), 令 \(u(\tau)=u(t/h^2)\), 上面的数值格式等价为\[ \frac{d u_i}{d\tau} = u_{i+1}-2u_{i}+u_{i-1}.\]

热方程的 RK4 方法讨论

Fehlberg [1] 在其论文中指出热方程的 RK4 方法相比于 Euler 方法 (A.1) 的优点并不明显. 原因如下: 1. RK4 方法的稳定性条件是 \(k/h^2< 0.7\), 优越性并不明显; 2. RK4 每一次时间迭代需要更多的计算量. 因此反观 Euler 方法, 如果我们要求的精度不太高, 其优点是明显的.

但是, 自适应地 RKF 方法并不是单纯的显式 Runge-Kutta 方法, 它具有时间步长自动调节的能力. Fehlberg 做了这样一个实验: 在 8 位精度的计算机上采用尽可能高的局部截断误差要求, 利用低价 RKF 方法求解热方程. 数值结果表明, 在相当多次时间迭代后, 数值解的有效数字依旧可以保持在 5-6 位, 而且不会出现稳定性问题.

在文献 [2] 中 Fehlberg 引入了高阶 Runge-Kutta 方法, 但是相对流行的还是低阶方法 (<5), 这是因为高阶方法带来的计算复杂度的问题.

RKF34

我们列出一种经典的 RKF34:\[ \begin{align} y_{i+1}&=y_i+\frac{1}{6}k_1+\frac{27}{52}k_3+\frac{49}{156}k_4, && (O(h^3))\\ w_{i+1}&=w_i+\frac{43}{288}k_1+\frac{243}{416}k_3+\frac{343}{1872}k_4+\frac{1}{12}k_5, && (O(h^4)) \end{align}\tag{6}\]其中\[\begin{align*} k_1 &= hf(t_i,y_i),\\ k_2 &= hf \left( t_i+\frac{h}{4}, y_i+\frac{1}{4}k_1 \right),\\ k_3 &= hf \left( t_i+\frac{4h}{9}, y_i+\frac{4}{81}k_1 +\frac{32}{81}k_2 \right),\\ k_4 &= hf \left( t_i+\frac{6h}{7}, y_i+\frac{57}{98}k_1-\frac{432}{343}k_2 + \frac{1053}{686}k_3 \right),\\ k_5 &= hf \left( t_i+h, y_i+\frac{1}{6}k_1+ \frac{27}{52}k_3-\frac{49}{156}k_4 \right).\end{align*}\]

RKF45 的 Python 实现

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

class my_ode_ivp:
    def __init__(self, fun, t_span, y0, yfun=None):
        self.fun = fun
        self.t_span = t_span
        self.y0 = y0

    def rkf45(self, TOL=1e-6, hmin=1e-7, hmax=0.2):
        flag = True
        h = hmax
        yi = self.y0
        ti = self.t_span[0]
        # print(ti, yi, h)

        while flag:
            # print(h)
            # 
            k1=h*self.fun(ti, yi)
            k2=h*self.fun(ti+.25*h, yi+.25*k1)
            k3=h*self.fun(ti+3*h/8, yi+3*k1/32+9*k2/32)
            k4=h*self.fun(ti+12*h/13, yi+1932*k1/2197-7200*k2/2197+7296*k3/2197)
            k5=h*self.fun(ti+h, yi+439*k1/216-8*k2+3680*k3/513-845*k4/4104)
            k6=h*self.fun(ti+.5*h, yi-8*k1/27+2*k2-3544*k3/2565+1859*k4/4104-11*k5/40)
            
            # compute |y_{i+1}-w_{i+1}|/h
            D = k1/360-128*k3/4275-2197*k4/75240+k5/50+2*k6/55
            R = abs(D)/h
            
            if R<TOL: # accurate is acceptable
                ti += h
                # w_{i+1} = y_{i+1} + (w_{i+1}-y_{i+1}) = y_{i+1} + D
                yi += 25*k1/216+1408*k3/2565+2197*k4/4104-k5/5 + D
                if yfun is None:
                    print(h, ti, yi)
                else:
                    print("{:11.7f} {:11.7f} {:11.5e} {:11.5e} {:11.4e}".format(h, ti,yi, yfun(ti), abs(yfun(ti)-yi)))
            
            # compute q
            q = .84*(TOL/R)**(.25)
            # print( q, h, q*h)
            # new size for step i-th or (i+1)-th
            if q<=.1:
                h=.1*h
            elif q>=4:
                h=4*h
            else:
                h *= q

            if h>hmax:
                h=hmax

            if ti >= self.t_span[1]:
                flag=False
            elif ti+h > self.t_span[1]:
                h = self.t_span[1]-ti
                if h<hmin:
                    flag=False



# setting for ode: f(t,y), y0, t_span
def f(t, y):
    return y-t**2+1 

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

y0 = .5 
t_span = (0, 2.)

# solve, compute error
ode = my_ode_ivp(fun=f, t_span=t_span, y0=y0, yfun=yfun)
ode.rkf45(TOL=1.0e-5, hmin=.001, hmax=.25)