常微分方程数值格式

Runge-Kutta方法

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

Runge-Kutta方法

函数$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 \]

局部截断误差为

\[\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} \]

显式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|cccc} a_1 & b_{11} & b_{12} & \cdots & b_{1N} \\ a_2 & b_{21} & b_{22} & \cdots & b_{2N} \\ a_3 & b_{31} & b_{32} & \cdots & b_{3N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_N & b_{N1} & b_{N2} & \cdots & b_{1N} \\ \hline & c_1 & c_2 & \cdots & c_{1N} \end{array} \]

例 1. 向前Euler格式可以记为

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

例 2. 显式中点格式

\[\begin{array}{c|cc} 0 & 0 & 0 \\ 1/2 & 1/2 & 0 \\ \hline & 0 & 1 \\ \end{array} \]

Heun格式

\[\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\times 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格式

\[\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} \]

Ralston格式

\[\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} \]

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阶经典格式

\[\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} \]

隐式Runge-Kutta格式

一个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} \]

例 3. N=1时，得到的隐式中点格式

\[\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} \]

例 4. N=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} \]

程序作业

常微分方程

  1. 经典4阶 Runge-Kutta 方法与 Adams 隐式3阶方法解微分方程的程序
  2. 求解如下的问题
    \[\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)$
  3. 给出$x=1.5$处的误差，并计算误差阶。（精确解为$y=\frac{3}{1+x^3}$）

输出 

Runge-Kutta法，误差和误差阶为
  k=0, 0.244934066848e00
  k=1, 0.534607244904e-01,  3.90
  ...
Adams隐格式，误差和误差阶为
  k=0, 0.244934066848e00
  k=1, 0.534607244904e-01 , 3.01
  ...

谢谢

参考：

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

vertical slide 2

目录

本节读完

例 5.

[#ex9-1-0].