数学建模：微分方程模型—常微分方程数值解算法及 Python 实现

考虑一阶常微分方程的初值问题
$$\{\begin{array}{l}\frac{dy}{dx}=f(x,y)x\in [a,b]\\ y(a)=y_0\end{array}$$

只要 $f(x,y)$ 在 $[a,b]\times {\mathbb{R}}^1$ 上连续，且关于 $y$ 满足 Lipschitz 条件，即存在与 $x,y$ 无关的常数 $L$, $s.t.$ $\left|f\left(x,y_1\right)-f\left(x,y_2\right)\right|\le L\left|y_1-y_2\right|$, 对任意定义在 $[a,b]$ 上的 $y_1(x)$ 和 $y_2(x)$ 成立，则上述初值问题的连续可微解 $y(x)$ 在 $[a,b]$ 上存在且唯一。

常微分方程数值解要计算出解函数 $y(x)$ 在一系列节点 $a=x_0<x_1<\dots <x_n=b$ 处的近似值 $y_i\approx y\left(x_i\right)(i=1,\dots ,n)$，节点间距 $h_i=x_{i+1}-x_i(i=0,\dots ,n-1)$ 称为步长，通常采用等距节点，即取 $h_i=h$ (常数)。

一、显式欧拉 (Euler) 法

在 $xy$ 平面上, 微分方程的解 $y=y(x)$ 称为积分曲线. 积分曲线上一点 $(x,y)$ 的切线斜率等于 函数 $f(x,y)$ 的值.

从初始点 $P_0\left(x_0,y_0\right)$ 出发, 沿切线方向推进到 $x=x_1$ 上一点 $P_1$, 再从 $P_1$ 沿切线方向推进到 $x=x_2$ 上一点 $P_2$, 循此前进作出一条折线 $P_0{P}_1{P}_2\cdots$, 作为积分曲线的一个近似. 即向前取差商近似导数
$$\{\begin{array}{l}{y}_0=y(x_0),\\ y_{i+1}=y_i+hf(x_i,y_i)(i=0,1,\cdots ,n-1)\end{array}$$

称为 (显式) 欧拉法

定义 在假设 $y_i=y\left(x_i\right)$ ，即第 $i$ 步计算是精确的前提下，$R_i=y\left(x_{i+1}\right)-y_{i+1}$ 称为局部截断误差。

定义 若某算法的局赔截断误差为 $O\left(h^{p+1}\right)$ ，则称该算法有 $p$ 阶精度。

显式欧拉格式具有 1 阶精度

Python 实现:

import numpy as np
#显式欧拉法的Python实现
def Euler(f,x0,y0,h,l):
    #f函数, x0,y0初值, h步长, l数量
    xn, yn = x0, y0
    n = 0
    ns, xs, ys = [n], [x0], [y0]
    while n <= l:
        n += 1
        xn += h
        yn = yn + f(xn,yn) * h
        ns.append(n)
        xs.append(xn)
        ys.append(yn)
    return (ns,xs,ys)

二、显式欧拉法的改进

隐式欧拉法 (后退欧拉法)

将常微分方程化为积分方程
$$y\left(x_{i+1}\right)=y\left(x_i\right)+{\int }_{x_i}^{x_{i+1}}f(x,y(x))dx.$$

采用左矩形公式计算积分
$$y_{i+1}=y_i+hf\left(x_i,y_i\right),$$

即是显式欧拉法的递推公式.

用右矩形公式计算积分, 即是向后取差商近似导数
$$y_{i+1}=y_i+hf\left(x_{i+1},y_{i+1}\right)(i=0,\dots ,n-1)$$

称为隐式欧拉法或后退欧拉法.

隐式欧拉格式具有 1 阶精度

梯形法

即显、隐式两种算法的平均

$$y_{i+1}=y_i+\frac{h}{2}[f(x_i,y_i)+f(x_{i+1},y_{i+1})](i=0,...,n-1)$$

称为梯形法.

梯形格式具有 2 阶精度

两步欧拉法 (中点法)

取中心差商近似导数
$$y_{i+1}=y_{i-1}+2hf(x_i,y_i)(i=1,...,n-1)$$

称为两步欧拉法或中点法.

中点格式具有 2 阶精度

预报 - 校正法 (改进欧拉法)

1. 先用显式欧拉法作预报, 算出 ${\bar{y}}_{i+1}=y_i+hf\left(x_i,y_i\right)$

2. 再将 ${\bar{y}}_{i+1}$ 代入隐式梯形法的右边作校正，得到
   $$y_{i+1}=y_i+\frac{h}{2}\left[f\left(x_i,y_i\right)+f\left(x_{i+1},{\bar{y}}_{i+1}\right)\right]$$

即
$$y_{i+1}=y_i+\frac{h}{2}\left[f\left(x_i,y_i\right)+f\left(x_{i+1},y_i+hf\left(x_i,y_i\right)\right)\right](i=0,\dots ,n-1)$$

预报 - 校正法通常只迭代一两次就转人下一步的计算, 避免了梯形法和中点法每次迭代, 都要重新计算函数 $f(x,y)$, 减少了计算量.

预报 - 校正法具有 2 阶精度

三、龙格 - 库塔 (Runge-Kutta) 法

基本思路

积分方程右端积分的求积公式可以表示为
$${\int }_{x_i}^{x_{i+1}}f(x,y(x))\mathrm{d}x\approx h\sum _{i=1}^r{c}_{i}f\left(x_i+{\lambda }_{i}h,y\left(x_i+{\lambda }_{i}h\right)\right)$$

将上式表示为
$$y_{n+1}=y_n+h\varphi \left(x_n,y_n,h\right)$$

其中
$$\begin{array}{cc}& \varphi \left(x_i,y_i,h\right)=\sum _{i=1}^r{c}_i{K}_i\\ & K_1=f\left(x_i,y_i\right)\\ & K_i=f\left(x_i+{\lambda }_{i}h,y_i+h\sum _{j=1}^{i-1}{\mu }_{ij}{K}_j\right)(i=2,\cdots ,r)\end{array}$$

称为 $r$ 阶显式龙格 - 库塔法. 当 $r=1$ 时, $\varphi \left(x_i,y_i,h\right)=f\left(x_i,y_i\right)$, 就是欧拉法.

2 阶龙格 - 库塔法

取 $r=2$ 时
$$\{\begin{array}{l}{y}_{i+1}=y_i+h\left[{\lambda }_1{K}_1+{\lambda }_2{K}_2\right]\\ K_1=f\left(x_i,y_i\right)\\ K_2=f\left(x_i+ph,y_i+ph{K}_1\right)\end{array}$$

其中 ${\lambda }_1+{\lambda }_2=1,{\lambda }_{2}p=\frac{1}{2}$, 算法有 2 阶精度

若取 ${\lambda }_1=0,{\lambda }_2=1$, 则 $p=\frac{1}{2}$, 则
$$\{\begin{array}{l}{y}_{i+1}=y_i+h{K}_2\\ K_1=f\left(x_i,y_i\right)\\ K_2=f\left(x_i+\frac{h}{2},y_i+\frac{h}{2}{K}_1\right)\end{array}$$

即是中点法

4 阶经典龙格 - 库塔法 (最常用)

4 阶龙格 - 库塔法是最常用的龙格 - 库塔法

$$\{\begin{array}{l}{y}_{i+1}=y_i+\frac{h}{6}\left(K_1+2K_2+2K_3+K_4\right)\\ K_1=f\left(x_i,y_i\right)\\ K_2=f\left(x_i+\frac{h}{2},y_i+\frac{h}{2}{K}_1\right)\\ K_3=f\left(x_i+\frac{h}{2},y_i+\frac{h}{2}{K}_2\right)\\ K_4=f\left(x_i+h,y_i+h{K}_3\right)\end{array}$$

Python 实现:

import numpy as np
#四阶龙格库塔法的Python实现
def RK4(f,x0,y0,h,maxx):
    #f函数, x0,y0初值, h步长, maxx: x最大值
    xn, yn = x0, y0
    n = 0
    ns, xs, ys = [n], [xn], [yn]
    while xn < maxx:
        n += 1
        xn += h
        k1 = f(xn,yn)
        k2 = f(xn + (h / 2),yn + (h * k1) / 2)
        k3 = f(xn + (h / 2),yn + (h * k2) / 2)
        k4 = f(xn + h,yn + h * k3)
        yn = yn + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
        ns.append(n)
        xs.append(xn)
        ys.append(yn)
    return (ns,xs,ys)

一般的龙格 - 库塔法

$$\{\begin{array}{rl}{y}_{i+1}& =y_i+h\left[{\lambda }_1{K}_1+{\lambda }_2{K}_2+\dots +{\lambda }_m{K}_m\right]\\ K_1& =f\left(x_i,y_i\right)\\ K_2& =f\left(x_i+{\alpha }_{2}h,y_i+{\beta }_{21}h{K}_1\right)\\ K_3& =f\left(x_i+{\alpha }_{3}h,y_i+{\beta }_{31}h{K}_1+{\beta }_{32}h{K}_2\right)\\ & \cdots \cdots \\ K_m& =f\left(x_i+{\alpha }_{m}h,y+{\beta }_{m1}h{K}_1+{\beta }_{m2}h{K}_2+\dots +{\beta }_{mm-1}h{K}_{m-1}\right)\end{array}$$

由于龙格-库塔法的导出基于泰勒展开, 故精度主要受解函数的光滑性影响. 对于光滑性不太好的解, 最好采用低阶算法而将步长 $h$ 取小.