Python小白的数学建模课-10.微分方程边值问题

小白往往听到微分方程就觉得害怕，其实数学建模中的微分方程模型不仅没那么复杂，而且很容易写出高水平的数模论文。

本文介绍微分方程模型边值问题的建模与求解，不涉及算法推导和编程，只探讨如何使用 Python 的工具包，零基础求解微分方程模型边值问题。

通过 3个 BVP 案例层层深入，手把手教你用 Python 搞定微分方程边值问题。

欢迎关注『Python小白的数学建模课 @ Youcans』系列，每周持续更新

1. 常微分方程的边值问题（BVP）

1.1 基本概念

微分方程是指含有未知函数及其导数的关系式。

微分方程是描述系统的状态随时间和空间演化的数学工具。物理中许多涉及变力的运动学、动力学问题，如空气的阻力为速度函数的落体运动等问题，很多可以用微分方程求解。微分方程在化学、工程学、经济学和人口统计等领域也有广泛应用。

微分方程分为初值问题和边值问题。初值问题是已知微分方程的初始条件，即自变量为零时的函数值，一般可以用欧拉法、龙哥库塔法来求解。边值问题则是已知微分方程的边界条件，即自变量在边界点时的函数值。

边值问题的提出和发展，与流体力学、材料力学、波动力学以及核物理学等密切相关，并且在现代控制理论等学科中有重要应用。例如，力学问题中的悬链线问题、弹簧振动问题，热学问题中的导热细杆问题、细杆端点冷却问题，流体力学问题、结构强度问题。

上节我们介绍的常微分方程，主要是微分方程的初值问题。本节介绍二阶常微分方程边值问题的建模与求解。

欢迎关注『Python小白的数学建模课 @ Youcans』系列，每周持续更新
Python小白的数学建模课-01.新手必读
Python小白的数学建模课-02.数据导入
Python小白的数学建模课-03.线性规划
Python小白的数学建模课-04.整数规划
Python小白的数学建模课-05.0-1规划
Python小白的数学建模课-06.固定费用问题
Python小白的数学建模课-07.选址问题
Python小白的数学建模课-09.微分方程模型
Python小白的数学建模课-10.微分方程边值问题
Python小白的数学建模课-12.非线性规划
Python小白的数学建模课-15.图论的基本概念
Python小白的数学建模课-16.最短路径算法
Python小白的数学建模课-17.条件最短路径算法
Python小白的数学建模课-18.最小生成树问题
Python小白的数学建模课-19.网络流优化问题
Python小白的数学建模课-20.网络流优化案例
Python小白的数学建模课-A3.12个新冠疫情数模竞赛赛题及短评
Python小白的数学建模课-B2. 新冠疫情 SI模型
Python小白的数学建模课-B3. 新冠疫情 SIS模型
Python小白的数学建模课-B4. 新冠疫情 SIR模型
Python小白的数学建模课-B5. 新冠疫情 SEIR模型
Python小白的数学建模课-B6. 新冠疫情 SEIR改进模型

1.2 常微分方程边值问题的数学模型

只含边界条件作为定解条件的常微分方程求解问题，称为常微分方程的边值问题（boundary value problem）。

一般形式的二阶常微分方程边值问题：
$y{\ ''} = f(x,y,y{\ '})，\; a<x<b$

有三种情况的边界条件：

（1）第一类边界条件（两点边值问题）：

$y (a) = y a ， y (b) = y b$

（2）第二类边界条件：

$y′(a)=ya，y′(b)=yby\ '(a)=ya，y\ '(b)=yb$

（3）第三类边界条件：
${y′(a)−a0y(a)=a1y′(b)−b0y(b)=b1\begin{cases} y\ '(a)-a_0\ y(a) = a_1\\ y\ '(b)-b_0\ y(b) = b_1 \end{cases}$

其中： $a0≥0，b0≥0，a0+b0>0a_0 \geq 0，b_0 \geq 0，a_0+b_0>0$

1.3 常微分方程边值问题的数值解法

简单介绍求解常微分方程边值问题的数值解法，常用方法有：打靶算法、有限差分法和有限元法。打靶算法把边值问题转化为初值问题求解，是根据边界条件反复迭代调整初始点的斜率，使初值问题的数值解在边界上“命中”问题的边值条件。有限差分法把空间离散为网格节点，用差商代替微商，将微分方程离散化为线性或非线性方程组来求解。有限元法将微分方程离散化，有限元就是指近似连续域的离散单元，对每一单元假定一个近似解，然后推导求解域满足条件，从而得到问题的解。

按照本系列“编程方案”的概念，不涉及这些算法的具体内容，只探讨如何使用 Python 的工具包、库函数，零基础求解微分方程模型边值问题。我们的选择还是 Python 常用工具包三剑客：Scipy、Numpy 和 Matplotlib。

2. SciPy 求解常微分方程边值问题

2.1 BVP 问题的标准形式

Scipy 用 solve_bvp() 函数求解常微分方程的边值问题，定义微分方程的标准形式为：
${y′=f(x,y)，a<x<bg(y(a),y(b)=0)\begin{cases} y{\ '} = f(x,y)，\; a<x<b\\ g(y(a),y(b)=0) \end{cases}$

因此要将第一类边界条件 $y (a) = y a ， y (b) = y b$ 改写为：
${y(a)−ya=0y(b)−yb=0\begin{cases} y(a)-ya=0\\ y(b)-yb=0 \end{cases}$

2.2 scipy.integrate.solve_bvp() 函数

**scipy.integrate.solve_bvp() **是求解微分方程边值问题的具体方法，可以求解一阶微分方程（组）的两点边值问题（第一类边界条件）。在 odeint 函数内部使用 FORTRAN 库 odepack 中的 lsoda，可以求解一阶刚性系统和非刚性系统的初值问题。官网介绍详见： scipy.integrate.solve_bvp — SciPy v1.7.0 Manual 。

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)

solve_bvp 的主要参数：

求解标准形式的微分方程（组）主要使用前 4 个参数：

func: callable fun(x, y, …) 　　导数函数 $f (y, x)$ ， y 在 x 处的导数，以函数的形式表示。可以带有参数 p。
bc: callable bc(ya, yb, …) 　　边界条件，y 在两点边界的函数，以函数的形式表示。可以带有参数 p。
x: array：　　初始网格的序列，shape (m,)。必须是单调递增的实数序列，起止于两点边界值 xa，xb。
y: array：　　网格节点处函数值的初值，shape (n,m)，第 i 列对应于 x[i]。
p: array：　　可选项，向导数函数 func、边界条件函数 bc 传递参数。

其它参数用于控制求解算法的参数，一般情况可以忽略。

solve_bvp 的主要返回值：

sol: PPoly 　　通过 PPoly （如三次连续样条函数）插值求出网格节点处的 y 值。
x: array 　　数组，形状为 (m,)，最终输出的网格节点。
y: array 　　二维数组，形状为 (n,m)，输出的网格节点处的 y 值。
yp: array 　　二维数组，形状为 (n,m)，输出的网格节点处的 y’ 值。

3. 实例 1：一阶常微分方程边值问题

3.1 例题 1：一阶常微分方程边值问题

求常微分方程边值问题的数值解。

${y′′+∣y∣=0y(x=0)=0.5y(x=4)=−1.5\begin{cases} \begin{aligned} &y\ ''+ \left| y \right| = 0\\ &y(x=0) = 0.5\\ &y(x=4) = -1.5 \end{aligned} \end{cases}$

引入变量 $y\ '$ ，通过变量替换就把原方程化为如下的标准形式的微分方程组：

${y0′=y1y1′=−∣y0∣y(x=0)−0.5=0y(x=4)+1.5=0\begin{cases} y_0^{'} = y_1\\ y_1^{'} = -\left| y_0 \right| \\ y(x=0) - 0.5 = 0\\ y(x=4) + 1.5 = 0 \end{cases}$

这样就可以用 solve_bvp() 求解该常微分方程的边值问题。

3.2 常微分方程的编程步骤

以该题为例讲解scipy.integrate.solve_bvp 求解常微分方程边值问题的步骤：

导入 scipy、numpy、matplotlib 包；
定义导数函数 dydx(x,y)

注意本问题中 y 表示向量，记为 $y=[y_0,y_1]$ ，导数定义函数 dydx(x,y) 编程如下：

# 导数函数，计算导数 dY/dx
def dydx(x, y):dy0 = y[1]dy1 = -abs(y[0])return np.vstack((dy0, dy1))

定义边界条件函数 boundCond(ya,yb)

本问题中边界条件为常数，具体编程如下。注意对照 3.1中的边值条件，没有为什么，函数就是这么定义的。

# 计算 边界条件
def boundCond(ya, yb):fa = 0.5  # 边界条件 y(xa=0) = 0.5fb = -1.5  # 边界条件 y(xb=4) = -1.5return np.array([ya[0]-fa,yb[0]-fb])

设置 x、y 的初值
调用 solve_bvp() 求解常微分方程在区间 [xa,xb] 的数值解
由 solve_bvp() 的返回值 sol，获得网格节点的处的 y值。

3.3 Python 例程

# mathmodel11_v1.py
# Demo10 of mathematical modeling algorithm
# Solving ordinary differential equations (boundary value problem) with scipy.from scipy.integrate import odeint, solve_bvp
import numpy as np
import matplotlib.pyplot as plt# 1. 求解微分方程组边值问题，DEMO
# y'' + abs(y) = 0, y(0)=0.5, y(4)=-1.5# 导数函数，计算导数 dY/dx
def dydx(x, y):dy0 = y[1]dy1 = -abs(y[0])return np.vstack((dy0, dy1))# 计算 边界条件
def boundCond(ya, yb):fa = 0.5  # 边界条件 y(xa=0) = 0.5fb = -1.5  # 边界条件 y(xb=4) = -1.5return np.array([ya[0]-fa,yb[0]-fb])xa, xb = 0, 4  # 边界点 (xa,xb)
# fa, fb = 0.5, -1.5  # 边界点的 y值
xini = np.linspace(xa, xb, 11)  # 确定 x 的初值
yini = np.zeros((2, xini.size))  # 确定 y 的初值
res = solve_bvp(dydx, boundCond, xini, yini)  # 求解 BVPxSol = np.linspace(xa, xb, 100)  # 输出的网格节点
ySol = res.sol(xSol)[0]  # 网格节点处的 y 值plt.plot(xSol, ySol, label='y')
# plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.title("scipy.integrate.solve_bvp")
plt.show()

3.4 Python 例程运行结果

在这里插入图片描述

4. 实例 2：水滴横截面的形状

4.1 例题 2：水滴横截面形状问题

水平面上的水滴横截面形状，可以用如下的微分方程描述：
${d2hdx2+[1−h]∗[1+(dhdx)2]3/2=0h(x=−1)=h(x=1)=0\begin{cases} \begin{aligned} &\frac{d^2 h}{dx^2} + [1-h]*[1+(\frac{dh}{dx})^2]^{3/2}= 0\\ &h(x=-1) = h(x=1) = 0 \end{aligned} \end{cases}$

引入变量 $h\ '$ ，通过变量替换就把原方程化为如下的标准形式的微分方程组：

${h0′=h1h1′=(h0−1)∗[1+h12]3/2h0(x=−1)=h0(x=1)=0\begin{cases} h_0^{'} = h_1\\ h_1^{'} = (h_0-1) * [1+h_1^2]^{3/2}\\ h_0(x=-1) = h_0(x=1) = 0 \end{cases}$

这样就可以用 solve_bvp() 求解该常微分方程的边值问题。

注：本问题来自司守奎等，数学建模算法与应用（第2版），国防工业出版社，2015

4.2 Python 例程：水滴横截面形状问题

# mathmodel11_v1.py
# Demo10 of mathematical modeling algorithm
# Solving ordinary differential equations (boundary value problem) with scipy.from scipy.integrate import odeint, solve_bvp
import numpy as np
import matplotlib.pyplot as plt# 3. 求解微分方程边值问题，水滴的横截面
# 导数函数，计算 h=[h0,h1] 点的导数 dh/dx
def dhdx(x,h):# 计算 dh0/dx, dh1/dx 的值dh0 = h[1]  # 计算 dh0/dxdh1 = (h[0]-1)*(1+h[1]*h[1])**1.5  # 计算 dh1/dxreturn np.vstack((dh0, dh1))# 计算 边界条件
def boundCond(ha,hb):# ha = 0  # 边界条件：h0(x=-1) = 0# hb = 0  # 边界条件：h0(x=1) = 0return np.array([ha[0],hb[0]])xa, xb = -1, 1  # 边界点 (xa=0, xb=1)
xini = np.linspace(xa, xb, 11)  # 设置 x 的初值
hini = np.zeros((2, xini.size))  # 设置 h 的初值res = solve_bvp(dhdx, boundCond, xini, hini)   # 求解 BVP
# scipy.integrate.solve_bvp(fun, bc, x, y,..)
#   fun(x, y, ..), 导数函数 f(y,x)，y在 x 处的导数。
#   bc(ya, yb, ..), 边界条件，y 在两点边界的函数。
#   x: shape (m)，初始网格的序列，起止于两点边界值 xa，xb。
#   y: shape (n,m)，网格节点处函数值的初值，第 i 列对应于 x[i]。xSol = np.linspace(xa, xb, 100)  # 输出的网格节点
hSol = res.sol(xSol)[0]  # 网格节点处的 h 值
plt.plot(xSol, hSol, label='h(x)')
plt.xlabel("x")
plt.ylabel("h(x)")
plt.axis([-1, 1, 0, 1])
plt.title("Cross section of water drop by BVP xupt")
plt.show()

4.3 Python 例程运行结果

在这里插入图片描述

5. 实例 3：带有未知参数的微分方程边值问题

5.1 例题 3：Mathieu 方程的特征函数

Mathieu 在研究椭圆形膜的边界值问题时，导出了一个二阶常微分方程，其形式为：
$d2ydx2+[λ−2qcos(2x)]y=0\frac{d^2 y}{dx^2} + [\lambda-2q\ cos(2x)]\ y= 0$

用这种形式的数学方程可以描述自然中的物理现象，包括振动椭圆鼓、四极质谱仪和四极离子阱、周期介质中的波动、强制振荡器参数共振现象、广义相对论中的平面波解决方案、量子摆哈密顿函数的本征函数、旋转电偶极子的斯塔克效应。

式中 $λ、q\lambda、q$ 是两个实参数，方程的系数是以 $π\pi$ 或 $2π2\pi$ 为周期的，但只有在 $λ、q\lambda、q$ 满足一定关系时 Mathieu 方程才有周期解。

引入变量 $y\ '$ ，通过变量替换就把原方程化为如下的标准形式的微分方程组：

${y0′=y1y1′=−[λ−2qcos(2x)]y0y0(x=0)=1y1(x=0)=0y1(x=π)=0\begin{cases} y_0^{'} = y_1\\ y_1^{'} = -[\lambda-2q\ cos(2x)]\ y_0 \\ y_0(x=0) = 1\\ y_1(x=0) = 0\\ y_1(x=\pi)=0 \end{cases}$

这样就可以用 solve_bvp() 求解该常微分方程的边值问题。

5.2 常微分方程的编程步骤

以该题为例讲解scipy.integrate.solve_bvp 求解常微分方程边值问题的步骤。

需要注意的是：（1）本案例涉及一个待定参数 $λ\lambda$ 需要通过 solve_bvp(fun, bc, x, y, p=None) 中的可选项 p 传递到导数函数和边界条件函数，（2）本案例涉及 3 个边界条件，要注意边界条件函数的定义。

导入 scipy、numpy、matplotlib 包；
定义导数函数 dydx(x,y,p)

本问题中 y 表示向量，记为 $y=[y_0,y_1]$ ，定义函数 dydx(x,y,p) 中的 p 是待定参数。

# 导数函数，计算导数 dY/dx
def dydx(x, y, p): # p 是待定参数lam = p[0]  # 待定参数，从 solve-bvp() 传递过来q = 10  # 设置参数dy0 = y[1]dy1 = -(lam-2*q*np.cos(2*x))*y[0]return np.vstack((dy0, dy1))

定义边界条件函数 boundCond(ya,yb,p)

注意，虽然边界条件定义函数并没有用到参数 p，但也必须写在输入变量中，函数就是这么要求的。

# 计算 边界条件
def boundCond(ya, yb, p):lam = p[0]return np.array([ya[0]-1,ya[0],yb[0]])

设置 x、y 的初值
调用 solve_bvp() 求解常微分方程在区间 [xa,xb] 的数值解
由 solve_bvp() 的返回值 sol，获得网格节点的处的 y值。

5.3 Python 例程

# mathmodel11_v1.py
# Demo10 of mathematical modeling algorithm
# Solving ordinary differential equations (boundary value problem) with scipy.from scipy.integrate import odeint, solve_bvp
import numpy as np
import matplotlib.pyplot as plt# 4. 求解微分方程组边值问题，Mathieu 方程
# y0' = y1, y1' = -(lam-2*q*cos(2x))y0)
# y0(0)=1, y1(0)=0, y1(pi)=0# 导数函数，计算导数 dY/dx
def dydx(x, y, p): # p 是待定参数lam = p[0]q = 10dy0 = y[1]dy1 = -(lam-2*q*np.cos(2*x))*y[0]return np.vstack((dy0, dy1))# 计算 边界条件
def boundCond(ya, yb, p):lam = p[0]return np.array([ya[0]-1,ya[0],yb[0]])xa, xb = 0, np.pi  # 边界点 (xa,xb)
xini = np.linspace(xa, xb, 11)  # 确定 x 的初值
xSol = np.linspace(xa, xb, 100)  # 输出的网格节点for k in range(5):A = 0.75*ky0ini = np.cos(8*xini)  # 设置 y0 的初值y1ini = -A*np.sin(8*xini)  # 设置 y1 的初值yini = np.vstack((y0ini, y1ini))  # 确定 y=[y0,y1] 的初值res = solve_bvp(dydx, boundCond, xini, yini, p=[10])  # 求解 BVPy0 = res.sol(xSol)[0]  # 网格节点处的 y 值y1 = res.sol(xSol)[1]  # 网格节点处的 y 值plt.plot(xSol, y0, '--')plt.plot(xSol, y1,'-',label='A = {:.2f}'.format(A))plt.xlabel("xupt")
plt.ylabel("y")
plt.title("Characteristic function of Mathieu equation")
plt.axis([0, np.pi, -5, 5])
plt.legend(loc='best')
plt.text(2,-4,"youcans-xupt",color='whitesmoke')
plt.show()

5.4 Python 例程运行结果

在这里插入图片描述

初值 A从 0~3.0 变化时，y-x 曲线（图中虚线）几乎不变，但 y’-x 的振幅增大；当 A 再稍微增大，系统就进入不稳定区， y-x 曲线振荡发散（图中未表示）。

关于 Mathieu 方程解的稳定性的讨论，已经不是数学建模课的内容，不再讨论。

6. 小结

微分方程的边值问题相对初值问题来说更为复杂，但是用 Scipy 工具包求解标准形式的微分方程边值问题，编程实现还是不难掌握的。
关于边值问题的模型稳定性、灵敏度的分析，是更为专业的问题。除非找到专业课程教材或范文中有相关内容可以参考套用，否则不建议小白自己摸索，这些问题不是调整参数试试就能试出来的。

更多微分方程数学模型案例，参见新冠疫情模型系列：

Python小白的数学建模课-B2. 新冠疫情 SI模型
Python小白的数学建模课-B3. 新冠疫情 SIS模型
Python小白的数学建模课-B4. 新冠疫情 SIR模型
Python小白的数学建模课-B5. 新冠疫情 SEIR模型
Python小白的数学建模课-B6. 新冠疫情 SEIR改进模型

【本节完】

欢迎关注『Python小白的数学建模课 @ Youcans』原创作品

原创作品，转载必须标注原文链接：(https://blog.csdn.net/youcans/article/details/118162990)。

Crated：2021-06-23

欢迎关注『Python小白的数学建模课 @ Youcans』系列，持续更新
Python小白的数学建模课-01.新手必读
Python小白的数学建模课-02.数据导入
Python小白的数学建模课-03.线性规划
Python小白的数学建模课-04.整数规划
Python小白的数学建模课-05.0-1规划
Python小白的数学建模课-06.固定费用问题
Python小白的数学建模课-07.选址问题
Python小白的数学建模课-09.微分方程模型
Python小白的数学建模课-10.微分方程边值问题
Python小白的数学建模课-12.非线性规划
Python小白的数学建模课-15.图论的基本概念
Python小白的数学建模课-16.最短路径算法
Python小白的数学建模课-17.条件最短路径算法
Python小白的数学建模课-18.最小生成树问题
Python小白的数学建模课-19.网络流优化问题
Python小白的数学建模课-20.网络流优化案例
Python小白的数学建模课-A1.国赛赛题类型分析
Python小白的数学建模课-A2.2021年数维杯C题探讨
Python小白的数学建模课-A3.12个新冠疫情数模竞赛赛题及短评
Python小白的数学建模课-B2. 新冠疫情 SI模型
Python小白的数学建模课-B3. 新冠疫情 SIS模型
Python小白的数学建模课-B4. 新冠疫情 SIR模型
Python小白的数学建模课-B5. 新冠疫情 SEIR模型
Python小白的数学建模课-B6. 新冠疫情 SEIR改进模型
Python数模笔记-PuLP库
Python数模笔记-StatsModels统计回归
Python数模笔记-Sklearn
Python数模笔记-NetworkX
Python数模笔记-模拟退火算法