Runge-Kutta方法

常微分方程数值格式

张瑞
中国科学技术大学数学科学学院

Runge-Kutta方法

  • 在高阶多步格式中，起步计算和预估计算都需要格式具有相应的高阶精度。
  • 但已知的单步格式中，精度最高只有2阶（梯形公式）。
  • 如何构造高阶的单步格式？

函数$y(x+h)$在$x$处Taylor展开，则有

\[\begin{aligned} y(x+h)=y(x)+hy'(x)+\frac{h^2}{2}y''(x)+\cdots+\frac{h^k}{k!}f^{(k)}(x)\\ +\frac{h^{k+1}}{(k+1)!}y^{(k+1)}(\xi) \end{aligned} \]

由微分方程，有

\[\begin{aligned} y'(x)=&f(x,y(x)) \\ y''(x)=&f'_x(x,y)+f'_y(x,y)y'(x)\\ =&f'_x(x,y)+f'_y(x,y)f(x,y) \\ y'''(x)=&\cdots \end{aligned} \]

得到

\[\begin{aligned} y(x_{n+1})=y(x_n)+hf(x_n,y_n)+\frac{h^2}2(f'_x(x_n,y_n)+f'_y(x_n,y_n)f(x_n,y_n))\\ +\cdots+\frac{h^{k+1}}{(k+1)!}y^{(k+1)}(\xi) \end{aligned} \]

这样，可以得到单步的高阶格式：

\[y_{n+1}=y_n+hf(x_n,y_n)+\frac{h^2}{2}(f'_x(x_n,y_n)+f'_y(x_n,y_n)f(x_n,y_n))+\cdots \]

局部截断误差为$\displaystyle\frac{h^{k+1}}{(k+1)!}y^{(k+1)}(\xi)$。

这种格式可以达到任意阶的精度，但用到了$f(x,y)$的各阶偏导数，使用不便。

Runge-Kutta方法

前面改写为

\[y(x_{n+1})=y(x_n)+hF(x_n,y_n,f,h)+hT_{n+1} \]

其中

\[T_{n+1}=\frac{h^{k}}{(k+1)!}y^{(k+1)}(\xi)=O(h^{k}) \]

上式中，$F$的计算因为涉及到$f(x,y)$的各阶偏导数，使用极为不便。

  • 取$f(x,y)$在$(x_n,y_n)$及其附近点的线性组合$Q$来近似表达$F$，使得误差也为$O(h^k)$。
  • 这样，不会影响到整体的局部截断误差，又避免了对$f(x,y)$的高阶偏导数的计算。

以2阶为例，

\[F(x_n,y_n,f,h)=f(x_n,y_n)+\frac{h}{2}(f'_x(x_n,y_n)+f'_y(x_n,y_n)f(x_n,y_n)) \]

则取合适的$c_1$, $c_2$, $a_2$, $b_{21}$，使得

\[Q=c_1f(x_n,y_n)+c_2f(x_n+a_2h,y_n+b_{21}hf(x_n,y_n)) \]

满足

\[Q=F+O(h^2) \]

利用Taylor展开

\[\begin{aligned} Q=&c_1f(x_n,y_n)+c_2f(x_n,y_n)\\ &+c_2a_2hf'_x(x_n,y_n)+c_2b_{21}hf(x_n,y_n)f'_y(x_n,y_n)+O(h^2) \end{aligned} \]

比较系数，得到如下非线性方程组

\[\begin{cases} c_1+c_2=1 \\ c_2a_2=\frac12 \\ c_1b_{21}=\frac12 \end{cases} \]

方程有无穷多个解

改进的Euler格式（Heun格式）

\[\begin{cases} c_1=c_2=\frac12 \\ a_2=1 , b_{21}=1 \end{cases} \Rightarrow \begin{cases} y_{n+1}=y_n+h\frac12(k_1+k_2) \\ k_1=f(x_n,y_n) \\ k_2=f(x_n+h,y_n+hk_1) \end{cases} \]

例 1. 用改进的Euler格式解方程

\[\begin{cases} y'(x)=-\sin(x) , \\ y(0)=1 \end{cases}, x\in[0,\pi] \]

解. 如图

num-ode-euler-improved

一个N级的Runge-Kutta格式的一般形式

\[\begin{cases} y_{n+1}=\displaystyle y_n+h\sum_{i=1}^Nc_iK_i \\ K_i=\displaystyle f(x_n+a_ih,y_n+h\sum_{j=1}^Nb_{ij}K_j) , i=1,\cdots,N \end{cases} \]

通常用如下的Butcher表来记录这些系数

\[\begin{array}{c|cccc} a_1 & b_{11} & b_{12} & \cdots & b_{1N} \\ a_2 & b_{21} & b_{22} & \cdots & b_{2N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_N & b_{N1} & b_{N2} & \cdots & b_{1N} \\ \hline & c_1 & c_2 & \cdots & c_{1N} \end{array} =\begin{array}{c|c} a & B \\ \hline & c^T \end{array} \]

显式Runge-Kutta格式

一个显式Runge-Kutta法具有如下形式

\[\begin{cases} y_{n+1}=\displaystyle y_n+h\sum_{i=1}^N c_iK_i \\ K_1=f(x_n,y_n) \\ K_2=f(x_n+a_2h,y_n+h(b_{21}K_1))\\ K_3=f(x_n+a_3h,y_n+h(b_{31}K_1+b_{32}K_2))\\ \cdots \\ K_N=\displaystyle f(x_n+a_Nh,y_n+h\sum_{j=1}^{N-1}b_{ij}K_j) \end{cases} \]

Butcher表

\[\begin{array}{c|ccccc} 0 & 0 & 0 & 0 & \cdots & 0 \\ a_2 & b_{21} & 0 & 0 & \cdots & 0 \\ a_3 & b_{31} & b_{32} & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ a_N & b_{N1} & b_{N2} & b_{N3} & \cdots & 0 \\ \hline & c_1 & c_2 & c_3 & \cdots & c_{1N} \end{array} \]

可以看到，此时矩阵$B$是一个对角元素全为0的下三角阵。

例 2. 格式

\[\begin{matrix} 0 & | & 0 \\ \hline & | & 1 \end{matrix} \]

对应格式为

\[\begin{cases} y_{n+1} = y_n+h K_1 \\ K_1 = f(x_n, y_n) \end{cases} \]

即为向前Euler格式

\[y_{n+1}=y_n+hf(x_n, y_n) \]

显式中点格式

\[\begin{array}{c|cc} 0 & 0 & 0 \\ 1/2 & 1/2 & 0 \\ \hline & 0 & 1 \\ \end{array} \]
\[\begin{cases} y_{n+1}=y_n+h K_2 \\ K_1=f(x_n,y_n) \\ K_2=f(x_n+\frac{h}2,y_n+h\frac12K_1) \end{cases} \]

Heun格式(改进的Euler格式)

\[\begin{cases} y_{n+1}=y_n+h(\frac12K_1+\frac12 K_2) \\ K_1=f(x_n,y_n) \\ K_2=f(x_n+{h},y_n+hK_1) \end{cases} \]
\[\begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1 & 0 \\ \hline & 1/2 & 1/2 \\ \end{array} \]

Ralston格式(在二阶格式中，具有最小的局部截断误差)

\[\begin{array}{c|cc} 0 & 0 & 0 \\ 2/3 & 3/4 & 0 \\ \hline & 1/4 & 3/4 \\ \end{array} \]
\[\begin{cases} y_{n+1}=y_n+h(\frac14K_1+\frac34 K_2) \\ K_1=f(x_n,y_n) \\ K_2=f(x_n+\frac23{h},y_n+\frac34hK_1) \end{cases} \]

3阶Kutta格式

\[\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 \\ 1 & -1 & 2 & 0 \\ \hline & 1/6 & 2/3 & 1/6 \\ \end{array} \]
\[\begin{cases} y_{n+1}=y_n+h(\frac16K_1+\frac23 K_2+\frac16K_3) \\ K_1=f(x_n,y_n) \\ K_2=f(x_n+\frac12{h},y_n+h\frac12K_1) \\ K_3=f(x_n+h,y_n+h(-K_1+2K_2)) \end{cases} \]

4阶经典Runge-Kutta格式

\[\begin{cases} y_{n+1}=y_n+h(\frac16K_1+\frac26 K_2+\frac26K_3+\frac16K_4) \\ K_1=f(x_n,y_n) \\ K_2=f(x_n+\frac12{h},y_n+h\frac12K_1) \\ K_3=f(x_n+\frac12{h},y_n+h\frac12K_2) \\ K_4=f(x_n+h,y_n+hK_3) \end{cases} \]
\[\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\ 1/2 & 1/2 & 0 & 0 & 0\\ 1/2 & 0 & 1/2 & 0 & 0\\ 1 & 0 & 0 & 1 & 0\\ \hline & 1/6 & 1/3 & 1/3 & 1/6\\ \end{array} \]

例 3. 用4阶经典Runge-Kutta格式

\[\begin{cases} y'(x)=-\sin(x) , \\ y(0)=1 \end{cases}, x\in[0,\pi] \]

解. 如图

num-ode-rk4

隐式Runge-Kutta格式

当Runge-Kutta的Butcher表中矩阵$B$不再是下三角阵时， 就是隐式的Runge-Kutta格式

例 4. N=1时，

\[\begin{array}{c|c} \frac{1}{2} & \frac{1}{2} \\ \hline & 1\\ \end{array} \]

得到的隐式中点格式

\[\begin{cases} y_{n+1}=y_n+hK_1 \\ K_1=f(x_n+\frac{h}2,y_n+\frac{h}2K_1) \end{cases} \]

可以看到，解$K_1$时，可能需要解非线性方程。

对应的显式中点格式

\[\begin{cases} y_{n+1}=y_n+h K_2 \\ K_1=f(x_n,y_n) \\ K_2=f(x_n+\frac{h}2,y_n+h\frac12K_1) \end{cases} \]

例 5. N=2时，得到2级4阶的隐式Runge-Kutta格式, Gauss–Legendre：

\[\begin{array}{c|cc} \frac{1}{2}-\frac{\sqrt3}{6} & \frac{1}{4} & \frac{1}{4}-\frac{\sqrt3}{6} \\ \frac{1}{2}+\frac{\sqrt3}{6} & \frac{1}{4}+\frac{\sqrt3}{6} &\frac{1}{4} \\ \hline & \frac{1}{2} & \frac{1}{2}\\ & \frac12+\frac12 \sqrt3 & \frac12-\frac12 \sqrt3 \\ \end{array} \]
\[\begin{cases} y_{n+1}=y_n+\frac{h}2(K_1+K_2) \\ K_1=f(x_n+(\frac12-\frac{\sqrt3}6)h,y_n+\frac{h}4K_1+(\frac14-\frac{\sqrt3}6)hK_2) \\ K_2=f(x_n+(\frac12+\frac{\sqrt3}6)h,y_n+(\frac14+\frac{\sqrt3}6)hK_1+\frac{h}4K_2) \end{cases} \]

在解$K_1$, $K_2$时，需要解一个非线性方程组。

程序作业

常微分方程

  1. 经典4阶 Runge-Kutta 方法的程序
  2. Adams 隐式3阶方法的程序，预估步分别用向前Euler格式，改进的Euler格式，4阶Runge-Kutta格式
  3. 求解如下的问题
    \[\begin{cases} y'(x)&=-x^2y^2 \\ y(0)&=3 \end{cases}, \quad 0\leq x\leq 1.5 \]
    分别取步长$h$为$0.1$, $\frac{0.1}2$, $\frac{0.1}4$, $\frac{0.1}8$计算$y(1.5)$
  4. 给出$x=1.5$处的误差，并计算误差阶。（精确解为$y=\frac{3}{1+x^3}$）

输出 

Runge-Kutta法，误差和误差阶为
  k=0, 0.244934066848e00
  k=1, 0.534607244904e-01,  3.90
  ...
Adams隐格式，误差和误差阶为
向前Euler做预估计算
  k=0, 0.244934066848e00
  k=1, 0.534607244904e-01 , 3.01
  ...
改进Euler做预估计算
  k=0, 0.244934066848e00
  k=1, 0.534607244904e-01 , 3.01
  ...
4阶Runge-Kutta做预估计算
  k=0, 0.244934066848e00
  k=1, 0.534607244904e-01 , 3.01
  ...

谢谢

目录

本节读完

参考：

  1. https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods

例 6.

6.