Runge-Kutta方法

在工程算法上Runge-Kutta(龙格库塔)数值积分方法应用十分广泛。Runge-Kutta方法最早出现在1895年,Runge-Kutta的经典算法是4阶Runge-Kutta法,简称为RK4。RK4是为了解决没有解析解或是难以求出解析解的积分问题。RK4相对于欧拉积分(或称为三角积分)精度更高;相对于泰勒展开,不需要再去求导。

Runge-Kutta方法的主要思想来源于泰勒展开,但是泰勒展开往往要求高阶导数,很多时候求高阶导数并不容易。Runge-Kutta方法用一阶导数的组合形式对泰勒展开中的高阶项进行逼近,使得该方法不需要求任何导数就可以达到不错的精度,而且计算量还挺小的——因为是迭代式的计算,所有计算结果重复使用,非常精简。

关于Runge-Kutta还有一些问题没有搞清楚,请带有自己的思考和判断进行阅读

推导

符号与说明

现有如下微分方程

$$ y'(t)=f(t,y) $$

我们已知$t$时刻的值$y(t)$,现在求$y(t+h)$的值,其中$h$为某小量。

根据泰勒展开,我们可以知道

$$ \begin{align*} y(t+h)&=y(t)+h\cdot y'(t)+\frac{1}{2}h^2\cdot y''(t)+\cdots+\frac{1}{n!}h^n\cdot y^{(n)}(t)+O(h^{n+1}) \\ &=y(t)+h\cdot f(t,y)+\frac{1}{2}h^2\cdot f'(t,y)+\cdots+\frac{1}{n!}h^n\cdot f^{(n-1)}(t,y)+O(h^{n+1}) \tag{1} \end{align*} $$

现在有Runge-Kutta的一般形式

$$ y(t+h)=y(t)+\sum_{i=1}^{q}\omega_i K_i \tag{2} $$

其中

$$ K_i=h\cdot f(t+\alpha_i,y+\sum_{j=1}^{i-1}\beta_{i,j}K_j) $$

且$q$为阶数,通常

$$ \alpha_1=0 $$

前置知识

  • $$ f(x+\delta x,y+\delta y)\approx f(x,y)+\delta x\cdot f'_x(x,y)+\delta\cdot f'_y(x,y) \tag{3}$ $$

  • $$ f'(x,y(x))=f'_x(x,y)+f'_y(x,y)\cdot f(x,y) \tag{4}$ $$

1阶Runge-Kutta推导

1阶Runge-Kutta方法(简称RK1)就是欧拉积分方法,这里为了给读者表达Runge-Kutta的思想,就先用这个简单的例子让读者快速明白其思想。

当$q=1$时,式$(2)$可以写为

$$ \begin{align*} y(t+h)=y(t)+w_1\cdot h\cdot f(t,y) \tag{5} \end{align*} $$

对比式$(4)$和式$(1)$,我们可以发现RK1就是泰勒展开的一阶导数的截断项,其中$\omega_1=1$,则RK1方法的表达式即为

$$ y(t+h)=y(t)+h\cdot f(t,y) $$

在这个推导中,我们可以发现Runge-Kutta方法总是往泰勒展开上靠拢,下面我们继续推导2阶Runge-Kutta方法。

2阶Runge-Kutta推导

当$q=2$时,式$(2)$可以写为

$$ y(t+h)=y(t)+\omega_1K_1+\omega_2K_2 \tag{6} $$

其中

$$ \begin{align*} K_1=&h\cdot f(t,y) \tag{7} \\ K_2=&h\cdot f(t+\alpha_2h,y+\beta_{2,1}K_1) \tag{8} \end{align*} $$

将式$(3)$的近似方法带入式$(8)$中,可得

$$ K_2\approx h\cdot\left[f(t,y)+\alpha_2h\cdot f'_t(t,y)+\beta_{2,1}K_1\cdot f'_y(t,y)\right] \tag{9} $$

再将式$(7)$和式$(9)$代回式$(6)$中

$$ \begin{align*} y(t+h)&=y(t)+\omega_1h\cdot f(t,y)+\omega_2h\cdot\left[f(t,y)+\alpha_2\cdot f'_t(t,y)+\beta_{2,1}K_1\cdot f'_y(t,y)\right] \\ &=y(t)+(\omega_1+\omega_2)h\cdot f(t,y)+\omega_2\alpha_2h^2\cdot f'_t(t,y)+\omega_2\beta_{2,1}hK_1\cdot f'_y(t,y) \\ &=y(t)+(\omega_1+\omega_2)h\cdot f(t,y)+\omega_2\alpha_2h^2\cdot f'_t(t,y)+\omega_2\beta_{2,1}h^2\cdot f(t,y)\cdot f'_y(t,y) \tag{10} \end{align*} $$

下面我们再回过头来看式$(1)$的泰勒展开,我们这时候只关注泰勒展开二阶截断项,并且利用式$(4)$进行化简展开,可得

$$ \begin{align*} y(t+h)&\approx y(t)+h\cdot f(t,y)+\frac{1}{2}h^2\cdot f'(t,y) \\ &=y(t)+h\cdot f(t,y)+\frac{1}{2}h^2\cdot f'(t,y(t)) \\ &=y(t)+h\cdot f(t,y)+\frac{1}{2}h^2\left[f'_t(t,y)+f'_y(t,y)\cdot f(t,y)\right] \\ &=y(t)+h\cdot f(t,y)+\frac{1}{2}h^2\cdot f'_t(t,y)+\frac{1}{2}h^2\cdot f(t,y)\cdot f'_y(t,y) \tag{11} \end{align*} $$

对比式$(10)$的Runge-Kutta方法和式$(11)$的泰勒展开,我们可以发现

$$ \left\{ \begin{align*} y(t+h)&=&y(t)&+&(\omega_1+\omega_2)h\cdot f(t,y)&+&\omega_2\alpha_2h^2\cdot f'_t(t,y)&+&\omega_2\beta_{2,1}h^2\cdot f(t,y)\cdot f'_y(t,y) \\ y(t+h)&=&y(t)&+&h\cdot f(t,y)&+&\frac{1}{2}h^2\cdot f'_x(t,y)&+&\frac{1}{2}h^2\cdot f(t,y)\cdot f'_y(t,y) \end{align*} \right. $$

因此我们可以得到如下方程

$$ \left\{ \begin{align*} &\omega_1+\omega_2=1\\ &\omega_2\alpha_2=\frac{1}{2}\\ &\omega_2\beta_{2,1}=\frac{1}{2} \end{align*} \right. \tag{12} $$

在这个方程组中,有四个未知数、三个方程,所以有无穷组解。

我们可以取$\omega_1=\omega_2=\frac{1}{2}$,也就是Runge Kutta中取平均数,可得

$$ \left\{ \begin{align*} &\omega_1=\frac{1}{2}\\ &\omega_2=\frac{1}{2}\\ &\alpha_2=1 \\ &\beta_{2,1}=1 \end{align*} \right. $$

将$\omega_1$、$\omega_2$、$\alpha_2$和$\beta_{2,1}$代入到式$(6)$中可得RK2的表达式

$$ \left\{ \begin{aligned} &y(t+h)=y(t)+\frac{1}{2}K_1+\frac{1}{2}K_2 \\ &K_1=h\cdot f(t,y)\\ &K_2=h\cdot f(t+h,y+K_1) \end{aligned} \right. $$

在2阶Runge-Kutta的推导中,我们可以清晰的看到RK方法和泰勒展开的关系和区别。当然在最后求解式$(12)$的方程组时,我们也可以取$\alpha_2=\frac{1}{2}$,那么最后的RK2的表达式则为

$$ \left\{ \begin{aligned} &y(t+h)=y(t)+K_2 \\ &K_1=h\cdot f(t,y)\\ &K_2=h\cdot f(t+\frac{h}{2},y+\frac{K_1}{2}) \end{aligned} \right. $$

4阶Runge-Kutta推导

有了2阶Runge-Kutta的推导,我们可以参照2阶的推导方法进而推导,但由于涉及到$f(t,y)$的二阶导和三阶导,推导和计算量比较大。但跟2阶的原理是一致的。这里给出经典的RK4表达式

$$ \left\{ \begin{aligned} &y(t+h)=y(t)+\frac{1}{6}(K_1+2K_2+2K_3+K_4) \\ &K_1=h\cdot f(t,y)\\ &K_2=h\cdot f(t+\frac{h}{2},y+\frac{K_1}{2}) \\ &K_3=h\cdot f(t+\frac{h}{2},y+\frac{K_2}{2}) \\ &K_4=h\cdot f(t+h,y+K_3) \end{aligned} \right. $$

还有一种4阶Runge-Kutta的变体,是Kutta在1901年提出的名为3/8法则的表达式

$$ \left\{ \begin{aligned} &y(t+h)=y(t)+\frac{1}{8}(K_1+3K_2+3K_3+K_4) \\ &K_1=h\cdot f(t,y)\\ &K_2=h\cdot f(t+\frac{h}{3},y+\frac{K_1}{3}) \\ &K_3=h\cdot f(t+\frac{2h}{3},y-\frac{K_1}{3}+K_2) \\ &K_4=h\cdot f(t+h,y+K_1-K_2+K_3) \end{aligned} \right. $$

这种4阶Runge-Kutta方法的精度比经典的RK4精度要高一些,但是相对的,计算量会大一些。

Runge-Kutta一些特定情况的简化

一些情况下,$f(t,y)$只和$t$相关,那么这个时候经典RK4就可以化简为

$$ \left\{ \begin{aligned} &y(t+h)=y(t)+\frac{1}{6}(K_1+2K_2+2K_3+K_4) \\ &K_1=h\cdot f(t)\\ &K_2=h\cdot f(t+\frac{h}{2}) \\ &K_3=h\cdot f(t+\frac{h}{2}) \\ &K_4=h\cdot f(t+h) \end{aligned} \right. \tag{13} $$

再将方程组化简

$$ y(t+h)=y(t)+h\cdot\frac{1}{6}\left[f(t)+2f(t+\frac{h}{2})+2f(t+\frac{h}{2})+f(t+h)\right] $$

这时就可以理解为$y$在不同$t$时的一阶导数的加权平均值来作为整个积分的斜率,如下所示

$$ \left\{ \begin{align*} &\bar{k}=\frac{1}{6}f(t)+\frac{2}{6}f(t+\frac{h}{2})+\frac{2}{6}f(t+\frac{h}{2})+\frac{1}{6}f(t+h) \\ &y(t+h)=y(t)+h\bar{k} \end{align*} \right. $$

式$(13)$这个形式在SLAM或是机器人学中经常会用到。

Reference

  1. https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods
  2. http://fourier.eng.hmc.edu/e176/lectures/NM/node42.html
  3. https://arxiv.org/abs/1711.02508
Last modification:June 2, 2021
If you think my article is useful to you, please feel free to appreciate