数计算数学

用Scipy解微分方程（组）

字数统计: 3.4k字   |   阅读时长: 15分
, , , ,

摘要: 用 Scipy 解常微分方程（组）

【对数据分析、人工智能、金融科技、风控服务感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：潮汐朝夕
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

当我们对业务进行建模的时候，有时会考虑微分方程。实践中当我们拍了一个微分方程的模型之后，可以通过 Scipy 快速求解微分方程的解并画出 $y = f(x), y = f’(x)$ 的图像，看看是不是符合预期，验证我们拍的模型是不是对的。

本文我们大致看一下 Scipy 解微分方程并画出图像的流程。

Scipy 常微分方程（组）求解组件

在 Scipy 中，scipy.integrate 子模块可以对函数进行定积分、二重积分、三重积分、n重积分，以及对离散数据进行积分计算，以及对微分方程组进行求解。

其中求解微分方程的组件主要是 solve_ivp 以及 solve_bvp。

solve_ivp

在 Scipy 中，solve_ivp 用于求解一阶常微分方程组。

一阶常微分方程组的形式如下：

其中 $y(0)$ 为初值，$\mathbf{y}$ 是是 n 维方程组，$t_{1}, \cdots, t_{m}$ 是额外参数。

solve_ivp 解常微分方程（组）的函数接口如下：

1
2
3
4
5
6
7
8
9
10
11
scipy.integrate.solve_ivp(fun
                         ,t_span
                         ,y0
                         ,method='RK45'
                         ,t_eval=None
                         ,dense_output=False
                         ,events=None
                         ,vectorized=False
                         ,args=None
                         ,**options
                         )

其中 method='RK45' 表示用龙格-库塔（Runge-Kutta）法进行数值求解。其它参数的含义如下：

  • fun: 可调用函数 f(t, y)，如果 y 的形状是 (n,) 或 (n, k) 的数组，则返回值的形状也与 y 的形状相同。
  • t_span: 由两个浮点数构成的元组 (t1, t2)设置 t 的积分下限和上限。
  • y0: 设置 y 的初始值，其形状是 (n, )。
  • method: 积分方法，可以取 RK45, Radau 等。
  • t_eval: 设置 t 的离散点，必须在 t_span 范围内。
  • dense_output: 是否进行连续求解。
  • vectorized: 是否把 y 的列当成向量处理，此时 y 的形状是 (n, k)。
  • args: 设置额外参数 t0, …, tm。
  • fiest_step: 步长。
  • max_step: 最大步长。
  • rtol, atol: 相对误差、绝对误差。
  • jac: 设置雅可比矩阵，形状是 (n, n)，None 时通过有限差分法获取。
  • jac_sparsity: 如果雅可比矩阵通过有限差分法获取，定义雅可比矩阵的稀疏结构，形状 (n, n)，元素是 0 的位置其值始终是 0。与 method 设置的积分方法有关。

该方法返回的的对象中。比较关键的属性如下：

  • t: 变量 t 的离散数据 ，形状 (n_points, )
  • y: 与 t 对应的 y 离散数据，形状 (n, n_points)

solve_bvp

solve_bvp 是求解一阶微分方程（组）的两点边值问题的方法（第一类边界条件）。

1
2
3
4
5
6
7
8
9
10
11
12
13
scipy.integrate.solve_bvp(fun
                         ,bc
                         ,x
                         ,y
                         ,p=None
                         ,S=None
                         ,fun_jac=None
                         ,bc_jac=None
                         ,tol=0.001
                         ,max_nodes=1000
                         ,verbose=0
                         ,bc_tol=None
                         )

参数含义如下：

  • func: 导数函数 f(t, y)，可以有参数。
  • bc: 边界条件，$y$ 在两点边界的函数，以函数形式表示，可以有参数。
  • x: 初始网格序列，必须是单调增额实数序列，起止于两点边界值。
  • y: 网格节点处函数值的初值。
  • p: 向导数函数func、边界条件函数 bc 传递参数。

其中 p 这个参数需要注意，如果原方程有待定参数，那么 func 和 bc 中都要带着 p，比较复杂，一般还是尽量不带着参数。

该方法返回的对象中，主要的属性如下，其中 n 的含义是方程组的方程个数：

  • sol: 通过插值求出网格节点处的 $y$ 值。
  • x: 数组，形状为 (m,)，最终输出的网格节点。
  • y: 二维数组，形状为 (n, m)，输出的网格节点处的 $y$ 值。
  • yp: 二维数组，形状为 (n ,m)，输出的网格节点处的 $y’$ 值。

solve_ivp 例子

一阶方程组 初值问题（混沌蝴蝶）

洛伦兹方程组是基于流体力学中的Navier-Stokes方程、热传导方程和连续性方程构建的，属于耗散系统。相空间中，耗散系统的终态都将收缩到吸引子的状态上。

洛伦茨吸引子由以下微分方程组定义：

$(y_{0}, y_{1}, y_{2})$ 定义了三维空间中各个坐标点上的速度矢量。从某个坐标开始，沿着速度矢量进行积分，可以得出无质量点在该空间中的运动轨迹。

其中 $\sigma, \rho, \beta$ 为三个参数，当参数为某些值时，轨迹出现混沌现象，即初值的微小差别也会显著影响运动轨迹。

这已经是一个一阶微分方程组，可以直接用 solve_ivp 求解。我们可以先求解该一阶微分方程组，画出 $\mathbf{y}$，以及其导数 $\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t}$ 的曲线。然后画出洛伦茨吸引子三维空间中的轨迹。

对比 solve_ivp 的接口，代码流程大致是先定义方程组的右端。返回 $\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t}$，然后设置画图的区域（积分上下限），生成自变量的离散点，求解，画图。完整代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp

# 返回 yi 对 t 的的一阶导数
def f(t, y, sigma, rho, beta):
    dy0_dt = sigma * (y[1] - y[0])
    dy1_dt = y[0] * (rho - y[2]) - y[1]
    dy2_dt = y[0] * y[1] - beta * y[2]
    return [dy0_dt, dy1_dt, dy2_dt]

# 画图区域
a = 0
b = 30

# 独立变量的离散点
t_series = np.linspace(a, b, 1501)

# 初始值
initial = [0.0, 1.0, 0.0]

# 其它参数
c = (10.0, 28.0, 3.0)

# 求解
sol = solve_ivp(fun=f
               ,t_span=[a,b]
               ,y0=initial
               ,method="RK45"
               ,args=c
               ,t_eval=t_series
               )

# 独立变量处的一阶导数
dy_dt = f(t_series, sol.y, *c)

if sol.success:
    fig = plt.figure()
    ax1 = fig.add_subplot(2, 1, 1)
    ax1.plot(sol.t, sol.y[0], "-", label="y0(t)")
    ax1.plot(sol.t, sol.y[1], "--", label="y1(t)")
    ax1.plot(sol.t, sol.y[2], "-.", label="y2(t)")
    ax1.legend()
    ax2 = fig.add_subplot(2, 1, 2)
    ax2.plot(sol.t, dy_dt[0], "-", label="y0'(t)")
    ax2.plot(sol.t, dy_dt[1], "--", label="y1'(t)")
    ax2.plot(sol.t, dy_dt[2], "-.", label="y2'(t)")
    ax2.legend()
    ax2.set_xlabel("t")
    plt.show()

下面画出 (0.0, 1.0, 0.0) 和 (0.0, 1.01, 0.0) 这两组初值下，吸引子的轨迹，观察混沌现象。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp

# 返回 yi 对 t 的的一阶导数
def f(t, y, sigma, rho, beta):
    dy0_dt = sigma * (y[1] - y[0])
    dy1_dt = y[0] * (rho - y[2]) - y[1]
    dy2_dt = y[0] * y[1] - beta * y[2]
    return [dy0_dt, dy1_dt, dy2_dt]

# 画图区域
a = 0
b = 30

# 独立变量的离散点
t_series = np.linspace(a, b, 1501)

# 初始值
initials = [[0.0, 1.0, 0.0]
           ,[0.0, 1.01, 0.0]
           ]

# 其它参数
c = (10.0, 28.0, 3.0)

fig = plt.figure()
ax3d = fig.add_subplot(1, 1, 1, projection="3d")

for initial in initials:
    # 求解
    sol = solve_ivp(fun=f
                   ,t_span=[a,b]
                   ,y0=initial
                   ,method="RK45"
                   ,args=c
                   ,t_eval=t_series
                   )

    if sol.success:
        ax3d.plot3D(sol.y[0], sol.y[1], sol.y[2])

ax3d.set_xlabel("x")
ax3d.set_ylabel("y")
ax3d.set_zlabel("z")
plt.show()

二阶方程组 初值问题（3自由度质量块-弹簧系统）

对于二阶或更高阶的微分方程，需要想办法将其转化成由一阶微分方程组成的方程组组的形式。在文章 n维线性空间中的微分方程 中我们知道将高阶微分方程转换为由一阶微分方程组成的方程组的一般处理方法。

下面以3自由度质量块-弹簧系统为例，看一下由高阶微分方程组成的方程组怎么处理。

如下图，$x_1, x_2, x_3$ 为质量块的位移，$m_1, m_2, m_3$ 时质量块的质量，$k_1, k_2, k_3$ 为弹簧的刚度，$c_1, c_2, c_3$ 为弹簧的阻尼，$F_1(t), F_2(t), F_3(t)$ 为作用在质量块上的外力。

由胡克定律、阻尼的定义、牛顿第二定律，列出动力学方程：

简写为：

进一步写为：

用标准方法（参考 入门级常微分方程小结）将二阶方程变为一阶方程组，令 $\mathbf{Y} = \mathbf{\dot{X}}$，则 $\mathbf{\dot{Y}} = \mathbf{\ddot{X}}$，于是以上方程可以写为：

进一步写为：

令 $\mathbf{Z} = \begin{bmatrix}\mathbf{X} \\ \mathbf{Y}\end{bmatrix}$，上式写为：

也就是：

其中 $\mathbf{M}, \mathbf{K}, \mathbf{C}$ 是额外参数。$\mathbf{F}(t)$ 如下：

该方程组可以用 solve_ivp 求解。代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp

# 返回 Zi 对 t 的的一阶导数
def f(t, Z, M, K, C):
    F1 = 10 * np.sin(10 * t)
    F2 = np.zeros_like(F1)
    F3 = 10 * np.sin(20 * t)
    F = np.array([[F1], [F2], [F3]]) # F 向量
    M_inv = np.linalg.inv(M) # 质量矩阵 M 的逆矩阵

    D = np.array([[Z[0]], [Z[1]], [Z[2]]]) # 位移
    V = np.array([[Z[3]], [Z[4]], [Z[5]]]) # 速度
    A = M_inv @ (F - C @ V - K @ D) # 加速度

    return [V[0, 0], V[1, 0], V[2, 0], A[0, 0], A[1, 0], A[2, 0]]

# 其它参数
m1 = 2; m2 = 2.5; m3 = 1
k1 = 10; k2 = 20; k3 = 30;
c1 = 0.4; c2 = 0.8; c3 = 1
M = np.diag([m1, m2, m3])
K = np.array([[k1 + k2, -k2, 0]
             ,[-k2, k2 + k3, -k3]
             ,[0, -k3, k3]]
             )
C = np.array([[c1 + c2, -c2, 0]
             ,[-c2, c2 + c3, -c3]
             ,[0, -c3, c3]]
             )

# 画图区域
a = 0
b = 15

# 独立变量的离散点
t_series = np.linspace(a, b, 15001)

# 初始值
initial = [10, -10, 0, 0, 0, 0]

# 求解
sol = solve_ivp(fun=f
               ,t_span=[a,b]
               ,y0=initial
               ,method="RK45"
               ,args=(M, K, C)
               ,t_eval=t_series
               )

if sol.success:
    fig = plt.figure()
    ax1 = fig.add_subplot(1, 2, 1)
    ax1.plot(sol.t, sol.y[0], "-", label="x1(t)")
    ax1.plot(sol.t, sol.y[1], "--", label="x2(t)")
    ax1.plot(sol.t, sol.y[2], "-.", label="x3(t)")
    ax1.legend()
    ax2 = fig.add_subplot(1, 2, 2)
    ax2.plot(sol.t, sol.y[3], "-", label="v1(t) = x1'(t)")
    ax2.plot(sol.t, sol.y[4], "--", label="v2(t) = x2'(t)")
    ax2.plot(sol.t, sol.y[5], "-.", label="v3(t) = x3'(t)")
    ax2.legend()
    plt.show()

solve_bvp

二阶微分方程 边值问题 (悬链线)

考虑悬链线方程的边值问题如下：

该方程的推导过程可以参考 入门级常微分方程小结。其中 $a$ 是参数，与线的长度有关。

参考文章 n维线性空间中的微分方程 中将二阶或更高阶的微分方程转化成由一阶微分方程组成的方程组的方法。

令 $y_{1} = y, y_{2} = \frac{\mathrm{d}y}{\mathrm{d}x}$，于是 $\frac{\mathrm{d}y_{2}}{\mathrm{d}x} = \frac{\mathrm{d}^{2}y}{\mathrm{d}x^{2}}$，原二阶微分方程的边值问题转化为微分方程组的边值问题：

如果我们已知参数 a，则方程中就没有未知参数了，可以直接求解。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
import numpy as np

# 导数函数
def f(x, y):
    dy1_dx = y[1]
    dy2_dx = 2 * np.sqrt(1 + y[1] ** 2)
    return [dy1_dx, dy2_dx]

# 边界条件
def bound_cond(ya, yb):
    fa = 8 # y(0) = 8
    fb = 5 # y(4) = 5
    return np.array([ya[0] - fa, yb[0] - fb])

# 边界点 (xa,xb)
xa, xb = 0, 4

# 用于迭代的离散点
x_series = np.linspace(xa, xb, 11)
y = np.zeros((2, x_series.size))  # 网格点上 vec(y) 的迭代初值

# 求解
res = solve_bvp(fun=f
               ,bc=bound_cond
               ,x=x_series
               ,y=y
               )

# 画图
x_plot = np.linspace(0, 4, 101)
y_plot = res.sol(x_plot)[0]
plt.plot(x_plot, y_plot, label="a=2.0")
plt.ylabel("y")
plt.xlabel("x")
plt.legend()
plt.show()

下面我们看一下参数 a 取不同值时的结果：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
import numpy as np

# 导数函数
def get_f(a):
    def f(x, y):
        dy1_dx = y[1]
        dy2_dx = a * np.sqrt(1 + y[1] ** 2)
        return [dy1_dx, dy2_dx]
    return f

# 边界条件
def bound_cond(ya, yb):
    fa = 8 # y(0) = 8
    fb = 5 # y(4) = 5
    return np.array([ya[0] - fa, yb[0] - fb])

# 边界点 (xa,xb)
xa, xb = 0, 4

# 用于迭代的离散点
x_series = np.linspace(xa, xb, 11)

for a in [0.2, 0.5, 1.0, 1.5, 2.0]:
    # 参数已知
    f = get_f(a)
    y = np.zeros((2, x_series.size))  # 网格点上 vec(y) 的迭代初值

    # 求解
    res = solve_bvp(fun=f
                   ,bc=bound_cond
                   ,x=x_series
                   ,y=y
                   )

    # 画图
    x_plot = np.linspace(0, 4, 101)
    y_plot = res.sol(x_plot)[0]
    plt.plot(x_plot, y_plot, label="a = {:.3f}".format(a))

plt.ylabel("y")
plt.xlabel("x")
plt.legend()
plt.show()

Share