数值计算方法 Chapter8. 常微分方程的数值解

• 数值计算方法 Chapter8. 常微分方程的数值解
  • 0. 问题描述
  • 1. Euler公式
    • 1. 向前Euler公式
    • 2. 向后Euler公式
    • 3. 梯形公式
  • 2. Runge-Kutta方法
    • 1. 二阶Runge-Kutta方法
    • 2. 高阶Runge-Kutta方法
      • 1. 三阶Runge-Kutta方法
      • 2. 四阶Runge-Kutta方法
    • 3. python伪代码实现
  • 3. 线性多步法
    • 1. 基本思路
    • 2. Adams公式
  • 4. 常微分方程组的数值解法
    • 1. 一阶常微分方程组的数值解法
    • 2. 高阶微分方程数值方法

0. 问题描述

这一章节考察的问题如标题所述，即常微分方程的数值求解：

$$\{\begin{array}{rl}\frac{dy}{dx}& =f(x,y)\\ y(x_0)& =y_0\end{array}$$

1. Euler公式

1. 向前Euler公式

Euler公式算是一个求解常微分方程数值解问题的一个比较直接的思路：

$$\frac{dy}{dx}=\frac{\delta y}{\delta x}=f(x,y)$$

从而有：

$$y+\delta y=f(x,y)\cdot \delta x$$

由此，我们不断地求解就能够近似的得到各个$x$取值下的$y$值。

$$y_{n+1}=y_n+h\cdot f(x_n,y_n)$$

给出python伪代码实现如下：

def fwd_euler_fn(f, x0, y0, step=1e-3):
    def fn(x):
        h = step if x >= x0 else -step
        n = math.ceil((x - x0) / h)
        x, y = x0, y0
        for _ in range(n):
            y += h * f(x, y)
            x += h
        return y
    return fn

2. 向后Euler公式

向后Euler公式和向前Euler公式并无本质区别，不过微调公式为：

$$y_{n+1}=y_n+h\cdot f(x_{n+1},y_{n+1})$$

但是显然的在计算$y_{n+1}$时事实上无法得知$f(x_{n+1},y_{n+1})$的值，因此这里需要叠加上迭代的思路：

$$y_{n+1}^{(k+1)}=y_n^{(k+1)}+h\cdot f(x_{n+1},y_{n+1}^{(k)})$$

上述迭代公式称为Picard迭代。

可以证明，当$h$足够小时，Picard迭代收敛。

给出python伪代码如下：

def bwd_euler_fn(f, x0, y0, step=1e-3):
    def fn(x, epsilon = 1e-6):
        h = step if x >= x0 else -step
        n = math.ceil((x - x0) / h)
        xlist = [x0 + i*h for i in range(n)]
        ylist = [y0 for i in range(n)]
        for i in range(n-1):
            ylist[i+1] = ylist[i] + h * f(xlist[i], ylist[i])
        while True:
            zlist = deepcopy(ylist)
            for i in range(n-1):
                zlist[i+1] = zlist[i] + h * f(xlist[i+1], ylist[i+1])
            delta = abs(zlist[-1] - ylist[-1])
            ylist = zlist
            if delta < epsilon:
                break
        return zlist[-1]
    return fn

3. 梯形公式

梯形公式本质上依然还是基于微分差商，不过不同于之前直接使用微分的形式，这里更加严格的使用了积分的表达，即：

$$y_{n+1}=y_n+{\int }_{x_n}^{x_{n+1}}f(x,y)dx$$

然后，这里使用梯形公式来近似掉其中的积分过程，有：

$$y_{n+1}=y_n+\frac{h}{2}(f(x_n,y_n)+f(x_{n+1},y_{n+1}))$$

这里同样会涉及到等号右侧的$y_{n+1}$的求解问题，不过，不同于向后Euler公式中的纯迭代思路，这里只使用一次迭代来近似，即：

$$y_{n+1}=y_n+\frac{h}{2}(f(x_n,y_n)+f(x_{n+1},y_n+h\cdot f(x_n,y_n)))$$

同样，我们给出python伪代码实现如下：

def tapezoid_fn(f, x0, y0, step=1e-3):
    def fn(x):
        h = step if x >= x0 else -step
        n = math.ceil((x - x0) / h)
        x, y = x0, y0
        for _ in range(n):
            y += h/2 * (f(x, y) + f(x+h, y+h*f(x, y)))
            x += h
        return y
    return fn

2. Runge-Kutta方法

1. 二阶Runge-Kutta方法

Runge-Kutta方法较之之前的Euler公式是一个相对而言精度更高的方法。

Euler公式来源于微分的定义，而Runge-Kutta方法则考察了Taylor展开，有：

$$\begin{array}{rl}{y}_{n+1}& =y_n+\sum _i\frac{h^i}{i!}{y}^{(i)}(x_n)\\ & =y_n+h\cdot y^{\prime }(x_n)+\frac{h^2}{2!}\cdot y^{\prime \prime }(x_n)+...\end{array}$$

我们保留二阶项即有：

$$\begin{array}{rl}{y}_{n+1}& =y_n+h\cdot y^{\prime }(x_n)+\frac{h^2}{2}{y}^{\prime \prime }(x_n)\\ & =y_n+h\cdot f(x_n,y_n)+\frac{h^2}{2}(\frac{\partial f}{\partial x}(x_n,y_n)+\frac{\partial f}{\partial y}(x_n,y_n)\cdot f(x_n,y_n))\end{array}$$

但是，这里的一个问题在于说两个偏导函数事实上也不一定好解，因此上述表达式没有办法直接调用。

而Runge-Kutta方法则是使用一个近似的位移公式来对其进行估计，使得二者在二阶导范围内没有误差。

具体而言：

$$\begin{array}{rl}\Delta & =y_{n+}-y_n\\ & =h\cdot y^{\prime }(x_n)+\frac{h^2}{2}{y}^{\prime \prime }(x_n)\\ & =h\cdot f(x_n,y_n)+\frac{h^2}{2}(\frac{\partial f}{\partial x}(x_n,y_n)+\frac{\partial f}{\partial y}(x_n,y_n)\cdot f(x_n,y_n))\\ & =c_1\cdot h\cdot f(x_n,y_n)+c_2\cdot h\cdot f(x_n+ah,y_n+bhf(x_n,y_n))\\ & =c_1\cdot h\cdot f(x_n,y_n)+c_2\cdot h\cdot [f(x_n,y_n)+\frac{\partial f}{\partial x}(x_n,y_n)\cdot ah+\frac{\partial f}{\partial y}(x_n,y_n)\cdot bhf(x_n,y_n)]\end{array}$$

或者说，给出参数$c_1,c_2,a,b$使得满足如下条件：

$$\{\begin{array}{rl}{y}_{n+1}& =y_n+h(c_1\cdot k_1+c_2\cdot k_2)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+ah,y_n+bh{k}_1)\end{array}$$

使得$y_{n+1}$在二阶导数范围内没有误差。

亦即：

$$\{\begin{array}{rl}(c_1+c_2)h\cdot f(x_n,y_n)& =h\\ c_{2}a{h}^2\cdot \frac{\partial f}{\partial x}(x_n,y_n)& =\frac{h^2}{2}\frac{\partial f}{\partial x}(x_n,y_n)\\ c_{2}b{h}^2\cdot \frac{\partial f}{\partial y}(x_n,y_n)& =\frac{h^2}{2}\frac{\partial f}{\partial y}(x_n,y_n)\end{array}$$

简化上式即可得到：

$$\{\begin{array}{rl}{c}_1+c_2& =1\\ 2c_{2}a& =1\\ 2c_{2}b& =1\end{array}$$

只要满足上述方程组，对应的参数$a,b,c_1,c_2$均可以令上两式在二阶导范围内没有误差。

给出两组常见的二阶Runge-Kutta公式如下：

1. $c_1=\frac{1}{2},c_2=\frac{1}{2},a=1,b=1$

   $$\{\begin{array}{rl}{k}_1& =f(x_n,y_n)\\ k_2& =f(x_n+h,y_n+h{k}_1)\\ y_{n+1}& =y_n+\frac{h}{2}(k_1+k_2)\end{array}$$

2. $c_1=0,c_2=1,a=\frac{1}{2},b=\frac{1}{2}$

   $$\{\begin{array}{rl}{k}_1& =f(x_n,y_n)\\ k_2& =f(x_n+\frac{h}{2},y_n+\frac{h}{2}{k}_1)\\ y_{n+1}& =y_n+h{k}_2\end{array}$$

2. 高阶Runge-Kutta方法

同样的，我们仿照上述的思路，给出一般情况下的高阶Runge-Kutta方法的表达式如下：

$$\{\begin{array}{rl}{y}_{n+1}& =y_n+h(c_1\cdot k_1+c_2\cdot k_2+...+c_m\cdot k_m)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+a_{1}h+y_n+b_{1,1}h\cdot k_1)\\ & ...\\ k_m& =f(x_n+a_{m-1}h+y_n+\sum _{i=1}^{m-1}{b}_{m-1,i}h\cdot k_i)\end{array}$$

参数$a,b,c$可以使得$y_n$具有$O(h^{m+1})$阶精度。

亦即，Taylor展开直到$h^m$都有两侧的展开系数相同。

1. 三阶Runge-Kutta方法

我们给出三阶Runge-Kutta方法的两组典型的系数如下：

1. 系数组合（一）

   $$\{\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{6}(k_1+4k_2+k_3)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{k}_1)\\ k_3& =f(x_n+h,y_n-h{k}_1+2h{k}_2)\end{array}$$

2. 系数组合（二）

   $$\{\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{4}(k_1+3k_3)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+\frac{1}{3}h,y_n+\frac{1}{3}h{k}_1)\\ k_3& =f(x_n+\frac{2}{3}h,y_n+\frac{2}{3}h{k}_2)\end{array}$$

3. 系数组合（三）

   $$\{\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{9}(2k_1+3k_2+4k_3)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{k}_1)\\ k_3& =f(x_n+\frac{3}{4}h,y_n+\frac{3}{4}h{k}_2)\end{array}$$

2. 四阶Runge-Kutta方法

同样的，我们可以给出两组典型的四阶Runge-Kutta公式如下：

1. 系数组合（一）

   $$\{\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{k}_1)\\ k_3& =f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{k}_2)\\ k_4& =f(x_n+h,y_n+h{k}_3)\end{array}$$

2. 系数组合（二）

   $$\{\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{8}(k_1+3k_2+3k_3+k_4)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+\frac{1}{3}h,y_n+\frac{1}{3}h{k}_1)\\ k_3& =f(x_n+\frac{2}{3}h,y_n+\frac{1}{3}h{k}_1+h{k}_2)\\ k_4& =f(x_n+h,y_n+h{k}_1-h{k}_2+h{k}_3)\end{array}$$

3. python伪代码实现

下面，我们来给出Runge-Kutta方法的python伪代码实现如下：

def runge_kutta_fn(f, x0, y0, a, b, c, step=1e-3):
    def fn(x):
        h = step if x >= x0 else -step
        n = math.ceil((x - x0) / h)
        m = len(c)
        assert(len(a) == m-1 and all(len(b[i]) == i+1 for i in range(m-1)))
        k = [0 for _ in range(m)]
        x, y = x0, y0
        for _ in range(n):
            k[0] = f(x, y)
            for i in range(m-1):
                xi = x + a[i] * h
                yi = y
                for j in range(i+1):
                    yi += h * b[i][j] * k[j]
                k[i+1] = f(xi, yi)
            x += h
            y += sum(ci * ki * h for ci, ki in zip(c, k))
        return y
    return fn

3. 线性多步法

1. 基本思路

线性多步法的思路来源依然还是积分公式：

$$y(x)=y(x_0)+{\int }_{x_0}^x{y}^{\prime }(t)dt$$

之前，Euler公式和Runge-Kutta方法都是直接对积分的值进行估计。

其中，Euler公式是直接令${\int }_x^{x+h}{y}^{\prime }(t)dt=y^{\prime }(x)h$，而Runge-Kutta方法则是通过引入一系列的偏移量使得$y(x)$在$n$阶Taylor展开当中没有误差。

而线性多步法的近似思路则是用采用之前的插值公式的思路，来对$y^{\prime }(x)$来进行拟合，然后用这个拟合函数来计算后面这个积分值。

给出书中的描述如下：

若用积分节点$x_n,x_{n-1},...,x_{n-q}$构造插值多项式近似$y^{\prime }(x)$，在区间$[x_{n-p},x_{n+1}]$上计算数值积分${\int }_{x_{n-p}}^{x_{n+1}}{y}^{\prime }(x)dx$，则称构造计算$y_{n+1}$的方法为线性多步法。

特别的：

• 如果使用$x_n,x_{n-1},...,x_{n-q}$构造插值多项式，则称拟合函数为显式Adams公式；
• 如果使用$x_{n+1},x_n,...,x_{n+1-q}$构造插值多项式，则称拟合函数为隐式Admas公式；

2. Adams公式

这里，我们舍去推导，直接给出一些常见的Adams公式如下：

1. 二阶显式Adams公式

   $$y_{n+1}=y_n+\frac{h}{2}(3f(x_n,y_n)-f(x_{n-1},y_{n-1}))$$

2. 二阶隐式Adams公式

   $$y_{n+1}=y_n+\frac{h}{2}(f(x_n,y_n)+f(x_{n+1},y_{n+1}))$$

3. 三阶显示Adams公式

   $$y_{n+1}=y_n+\frac{h}{12}(23f(x_n,y_n)-16f(x_{n-1},y_{n-1})+5f(x_{n-2},y_{n-2}))$$

4. 三阶隐示Adams公式

   $$y_{n+1}=y_n+\frac{h}{12}(5f(x_{n+1},y_{n+1}+8f(x_n,y_n)-f(x_{n-1},y_{n-1}))$$

5. 四阶显式Adams公式

   $$y_{n+1}=y_n+\frac{h}{24}(55f(x_n,y_n)-59f(x_{n-1},y_{n-1})+37f(x_{n-2},y_{n-2})-9f(x_{n-3},y_{n-3}))$$

6. 四阶隐式Adams公式

   $$y_{n+1}=y_n+\frac{h}{24}(9f(x_{n+1},y_{n+1}+19f(x_n,y_n)-5f(x_{n-1},y_{n-1})+f(x_{n-2},y_{n-2}))$$

4. 常微分方程组的数值解法

1. 一阶常微分方程组的数值解法

我们给出一阶常微分方程的初值问题表达如下：

$$\{\begin{array}{rl}\frac{d{y}_1}{dx}& =f_1(x,y_1,y_2,...,y_m)\\ \frac{d{y}_2}{dx}& =f_2(x,y_1,y_2,...,y_m)\\ & ...\\ \frac{d{y}_m}{dx}& =f_m(x,y_1,y_2,...,y_m)\\ y_1(x_0)& ={\eta }_1\\ y_2(x_0)& ={\eta }_2\\ & ...\\ y_m(x_0)& ={\eta }_m\end{array}$$

这类方程组的解法问题其实完全可以原模原样照搬常微分方程的数值解法，倒是也没啥必要详细展开就是了。

我们直接以二元方程组为例，给出一些常见的解：

$$\{\begin{array}{rl}\frac{dy}{dx}& =f(x,y,z)\\ \frac{dz}{dx}& =g(x,y,z)\\ y(x_0)& =y_0\\ z(x_0)& =z_0\end{array}$$

1. 前向Euler公式

   $$\{\begin{array}{rl}{y}_{n+1}& =y_n+hf(x_n,y_n,z_n)\\ z_{n+1}& =z_n+hg(x_n,y_n,z_n)\end{array}$$

2. 梯形公式

   $$\{\begin{array}{rl}\widehat{y_{n+1}}& =y_n+hf(x_n,y_n,z_n)\\ \widehat{z_{n+1}}& =z_n+hg(x_n,y_n,z_n)\\ y_{n+1}& =y_n+\frac{h}{2}(f(x_n,y_n,z_n)+f(x_{n+1},\widehat{y_{n+1}},\widehat{z_{n+1}}))\\ z_{n+1}& =z_n+\frac{h}{2}(g(x_n,y_n,z_n)+g(x_{n+1},\widehat{y_{n+1}},\widehat{z_{n+1}}))\end{array}$$

3. 四阶Runge-Kutta方法

   $$\{\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{6}(k_{11}+2k_{12}+2k_{13}+k_{14})\\ z_{n+1}& =z_n+\frac{h}{6}(k_{21}+2k_{22}+2k_{23}+k_{24})\\ k_{11}& =f(x_n,y_n,z_n)\\ k_{21}& =g(x_n,y_n,z_n)\\ k_{12}& =f(x_n+\frac{h}{2},y_n+\frac{h}{2}{k}_{11},z_n+\frac{h}{2}{k}_{21})\\ k_{22}& =g(x_n+\frac{h}{2},y_n+\frac{h}{2}{k}_{11},z_n+\frac{h}{2}{k}_{21})\\ k_{13}& =f(x_n+\frac{h}{2},y_n+\frac{h}{2}{k}_{12},z_n+\frac{h}{2}{k}_{22})\\ k_{23}& =g(x_n+\frac{h}{2},y_n+\frac{h}{2}{k}_{12},z_n+\frac{h}{2}{k}_{22})\\ k_{14}& =f(x_n+h,y_n+h{k}_{13},z_n+h{k}_{23})\\ k_{24}& =g(x_n+h,y_n+h{k}_{13},z_n+h{k}_{23})\end{array}$$

2. 高阶微分方程数值方法

这里，我们再来考察一下一元高阶微分方程的数值解法。

$$\{\begin{array}{rl}\frac{d^{m}y(x)}{d{x}^m}& =f(x,y,y^{\prime },...,y^{(m-1)})\\ y(x_0)& ={\eta }_0\\ y^{\prime }(x_0)& ={\eta }_1\\ & ...\\ y^{(m-1)}(x_0)& ={\eta }_{m-1}\end{array}$$

这一类问题事实上可以作为上述一阶常微分方程组的一个应用实例，我们只需要做如下变换就可以将问题完全转换为一个一阶常微分方程组，然后就可以运用之前的一阶常微分方程组的数值解法进行求解了。

$$\{\begin{array}{rl}\frac{d{y}_0}{dx}& =y_1(x)\\ \frac{d{y}_1}{dx}& =y_2(x)\\ & ...\\ \frac{d{y}_{m-2}}{dx}& =y_{m-1}(x)\\ \frac{d{y}_{m-1}}{dx}& =f(x,y_0,y_1,...,y_{m-1})\\ y_0(x_0)& ={\eta }_0\\ y_1(x_0)& ={\eta }_1\\ & ...\\ y_{m-1}(x_0)& ={\eta }_{m-1}\end{array}$$