Python小白的数学建模课-10 微分方程边值问题
Posted youcans
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了Python小白的数学建模课-10 微分方程边值问题相关的知识，希望对你有一定的参考价值。
小白往往听到微分方程就觉得害怕，其实数学建模中的微分方程模型不仅没那么复杂，而且很容易写出高水平的数模论文。
本文介绍微分方程模型边值问题的建模与求解，不涉及算法推导和编程，只探讨如何使用 Python 的工具包，零基础求解微分方程模型边值问题。
通过 3个 BVP 案例层层深入，手把手教你搞定微分方程边值问题。
欢迎关注『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.微分方程边值问题
1.2 常微分方程边值问题的数学模型
只含边界条件作为定解条件的常微分方程求解问题，称为常微分方程的边值问题（boundary value problem）。
一般形式的二阶常微分方程边值问题：
y
′
′
=
f
(
x
,
y
,
y
′
)
，
a
<
x
<
b
y{\\ ''} = f(x,y,y{\\ '})，\\; a<x<b
y ′′=f(x,y,y ′)，a<x<b
有三种情况的边界条件：
（1）第一类边界条件（两点边值问题）：
y ( a ) = y a ， y ( b ) = y b y(a)=ya，y(b)=yb y(a)=ya，y(b)=yb
（2）第二类边界条件：
y ′ ( a ) = y a ， y ′ ( b ) = y b y\\ '(a)=ya，y\\ '(b)=yb y ′(a)=ya，y ′(b)=yb
（3）第三类边界条件：
{
y
′
(
a
)
−
a
0
y
(
a
)
=
a
1
y
′
(
b
)
−
b
0
y
(
b
)
=
b
1
\\begin{cases} y\\ '(a)-a_0\\ y(a) = a_1\\\\ y\\ '(b)-b_0\\ y(b) = b_1 \\end{cases}
{y ′(a)−a0 y(a)=a1y ′(b)−b0 y(b)=b1
其中： a 0 ≥ 0 ， b 0 ≥ 0 ， a 0 + b 0 > 0 a_0 \\geq 0，b_0 \\geq 0，a_0+b_0>0 a0≥0，b0≥0，a0+b0>0
1.3 常微分方程边值问题的数值解法
简单介绍求解常微分方程边值问题的数值解法，常用方法有：打靶算法、有限差分法和有限元法。打靶算法把边值问题转化为初值问题求解，是根据边界条件反复迭代调整初始点的斜率，使初值问题的数值解在边界上“命中”问题的边值条件。有限差分法把空间离散为网格节点，用差商代替微商，将微分方程离散化为线性或非线性方程组来求解。 有限元法将微分方程离散化，有限元就是指近似连续域的离散单元，对每一单元假定一个近似解，然后推导求解域满足条件，从而得到问题的解。
按照本系列“编程方案”的概念，不涉及这些算法的具体内容，只探讨如何使用 Python 的工具包、库函数，零基础求解微分方程模型边值问题。我们的选择还是 Python 常用工具包三剑客：Scipy、Numpy 和 Matplotlib。
2. SciPy 求解常微分方程边值问题
2.1 BVP 问题的标准形式
Scipy 用 solve_bvp() 函数求解常微分方程的边值问题，定义微分方程的标准形式为：
{
y
′
=
f
(
x
,
y
)
，
a
<
x
<
b
g
(
y
(
a
)
,
y
(
b
)
=
0
)
\\begin{cases} y{\\ '} = f(x,y)，\\; a<x<b\\\\ g(y(a),y(b)=0) \\end{cases}
{y ′=f(x,y)，a<x<bg(y(a),y(b)=0)
因此要将第一类边界条件
y
(
a
)
=
y
a
，
y
(
b
)
=
y
b
y(a)=ya，y(b)=yb
y(a)=ya，y(b)=yb 改写为：
{
y
(
a
)
−
y
a
=
0
y
(
b
)
−
y
b
=
0
\\begin{cases} y(a)-ya=0\\\\ y(b)-yb=0 \\end{cases}
{y(a)−ya=0y(b)−yb=0
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 ) f(y,x) 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 ∣ = 0 y ( x = 0 ) = 0.5 y ( 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 ′′+∣y∣=0y(x=0)=0.5y(x=4)=−1.5
引入变量 y 0 = y ， y 1 = y ′ y0 = y，y1 = y\\ ' y0=y，y1=y ′，通过变量替换就把原方程化为如下的标准形式的微分方程组：
{ y 0 ′ = y 1 y 1 ′ = − ∣ y 0 ∣ y ( x = 0 ) − 0.5 = 0 y ( 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} ⎩⎪⎪⎪⎨⎪⎪⎪⎧y0′=y1y1′=−∣y0∣y(x=0)−0.5=0y(x=4)+1.5=0
这样就可以用 solve_bvp() 求解该常微分方程的边值问题。
3.2 常微分方程的编程步骤
以该题为例讲解scipy.integrate.solve_bvp 求解常微分方程边值问题的步骤：
-
导入 scipy、numpy、matplotlib 包；
-
定义导数函数 dydx(x,y)
注意本问题中 y 表示向量，记为 y = [ y 0 , y 1 ] y=[y_0,y_1] y=[y0,以上是关于Python小白的数学建模课-10 微分方程边值问题的主要内容，如果未能解决你的问题，请参考以下文章
Python小白的数学建模课-B5. 新冠疫情 SEIR模型