音乐网站前台模板,网站 chat now怎么做,陕西网站建设报价,企业网站备案域名信息小白往往听到微分方程就觉得害怕#xff0c;其实数学建模中的微分方程模型不仅没那么复杂#xff0c;而且很容易写出高水平的数模论文。
本文介绍微分方程模型边值问题的建模与求解#xff0c;不涉及算法推导和编程#xff0c;只探讨如何使用 Python 的工具包#xff0c;…
小白往往听到微分方程就觉得害怕其实数学建模中的微分方程模型不仅没那么复杂而且很容易写出高水平的数模论文。
本文介绍微分方程模型边值问题的建模与求解不涉及算法推导和编程只探讨如何使用 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′)axby{\ } f(x,y,y{\ })\; axb y ′′f(x,y,y ′)axb
有三种情况的边界条件
1第一类边界条件两点边值问题
y(a)yay(b)yby(a)yay(b)yb y(a)yay(b)yb
2第二类边界条件
y′(a)yay′(b)yby\ (a)yay\ (b)yb y ′(a)yay ′(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} {y ′(a)−a0 y(a)a1y ′(b)−b0 y(b)b1
其中a0≥0b0≥0a0b00a_0 \geq 0b_0 \geq 0a_0b_00a0≥0b0≥0a0b00 1.3 常微分方程边值问题的数值解法
简单介绍求解常微分方程边值问题的数值解法常用方法有打靶算法、有限差分法和有限元法。打靶算法把边值问题转化为初值问题求解是根据边界条件反复迭代调整初始点的斜率使初值问题的数值解在边界上“命中”问题的边值条件。有限差分法把空间离散为网格节点用差商代替微商将微分方程离散化为线性或非线性方程组来求解。 有限元法将微分方程离散化有限元就是指近似连续域的离散单元对每一单元假定一个近似解然后推导求解域满足条件从而得到问题的解。
按照本系列“编程方案”的概念不涉及这些算法的具体内容只探讨如何使用 Python 的工具包、库函数零基础求解微分方程模型边值问题。我们的选择还是 Python 常用工具包三剑客Scipy、Numpy 和 Matplotlib。 2. SciPy 求解常微分方程边值问题
2.1 BVP 问题的标准形式
Scipy 用 solve_bvp() 函数求解常微分方程的边值问题定义微分方程的标准形式为 {y′f(x,y)axbg(y(a),y(b)0)\begin{cases} y{\ } f(x,y)\; axb\\ g(y(a),y(b)0) \end{cases} {y ′f(x,y)axbg(y(a),y(b)0)
因此要将第一类边界条件 y(a)yay(b)yby(a)yay(b)yby(a)yay(b)yb 改写为 {y(a)−ya0y(b)−yb0\begin{cases} y(a)-ya0\\ y(b)-yb0 \end{cases} {y(a)−ya0y(b)−yb0
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, pNone, SNone, fun_jacNone, bc_jacNone, tol0.001, max_nodes1000, verbose0, bc_tolNone)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,)。必须是单调递增的实数序列起止于两点边界值 xaxb。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(x0)0.5y(x4)−1.5\begin{cases} \begin{aligned} y\ \left| y \right| 0\\ y(x0) 0.5\\ y(x4) -1.5 \end{aligned} \end{cases} ⎩⎪⎨⎪⎧y ′′∣y∣0y(x0)0.5y(x4)−1.5
引入变量 y0yy1y′y0 yy1 y\ y0yy1y ′通过变量替换就把原方程化为如下的标准形式的微分方程组
{y0′y1y1′−∣y0∣y(x0)−0.50y(x4)1.50\begin{cases} y_0^{} y_1\\ y_1^{} -\left| y_0 \right| \\ y(x0) - 0.5 0\\ y(x4) 1.5 0 \end{cases} ⎩⎪⎪⎪⎨⎪⎪⎪⎧y0′y1y1′−∣y0∣y(x0)−0.50y(x4)1.50
这样就可以用 solve_bvp() 求解该常微分方程的边值问题。 3.2 常微分方程的编程步骤
以该题为例讲解scipy.integrate.solve_bvp 求解常微分方程边值问题的步骤 导入 scipy、numpy、matplotlib 包 定义导数函数 dydx(x,y) 注意本问题中 y 表示向量记为 y[y0,y1]y[y_0,y_1]y[y0,y1]导数定义函数 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(xa0) 0.5fb -1.5 # 边界条件 y(xb4) -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(xa0) 0.5fb -1.5 # 边界条件 y(xb4) -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, labely)
# 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/20h(x−1)h(x1)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(x1) 0 \end{aligned} \end{cases} ⎩⎨⎧dx2d2h[1−h]∗[1(dxdh)2]3/20h(x−1)h(x1)0
引入变量 h0hh1h′h0 hh1 h\ h0hh1h ′通过变量替换就把原方程化为如下的标准形式的微分方程组
{h0′h1h1′(h0−1)∗[1h12]3/2h0(x−1)h0(x1)0\begin{cases} h_0^{} h_1\\ h_1^{} (h_0-1) * [1h_1^2]^{3/2}\\ h_0(x-1) h_0(x1) 0 \end{cases} ⎩⎪⎨⎪⎧h0′h1h1′(h0−1)∗[1h12]3/2h0(x−1)h0(x1)0
这样就可以用 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)*(1h[1]*h[1])**1.5 # 计算 dh1/dxreturn np.vstack((dh0, dh1))# 计算 边界条件
def boundCond(ha,hb):# ha 0 # 边界条件h0(x-1) 0# hb 0 # 边界条件h0(x1) 0return np.array([ha[0],hb[0]])xa, xb -1, 1 # 边界点 (xa0, xb1)
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)初始网格的序列起止于两点边界值 xaxb。
# y: shape (n,m)网格节点处函数值的初值第 i 列对应于 x[i]。xSol np.linspace(xa, xb, 100) # 输出的网格节点
hSol res.sol(xSol)[0] # 网格节点处的 h 值
plt.plot(xSol, hSol, labelh(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 例题 3Mathieu 方程的特征函数
Mathieu 在研究椭圆形膜的边界值问题时导出了一个二阶常微分方程其形式为 d2ydx2[λ−2qcos(2x)]y0\frac{d^2 y}{dx^2} [\lambda-2q\ cos(2x)]\ y 0 dx2d2y[λ−2q cos(2x)] y0
用这种形式的数学方程可以描述自然中的物理现象包括振动椭圆鼓、四极质谱仪和四极离子阱、周期介质中的波动、强制振荡器参数共振现象、广义相对论中的平面波解决方案、量子摆哈密顿函数的本征函数、旋转电偶极子的斯塔克效应。
式中 λ、q\lambda、qλ、q 是两个实参数方程的系数是以 π\piπ 或 2π2\pi2π 为周期的但只有在 λ、q\lambda、qλ、q 满足一定关系时 Mathieu 方程才有周期解。
引入变量 y0yy1y′y0 yy1 y\ y0yy1y ′通过变量替换就把原方程化为如下的标准形式的微分方程组
{y0′y1y1′−[λ−2qcos(2x)]y0y0(x0)1y1(x0)0y1(xπ)0\begin{cases} y_0^{} y_1\\ y_1^{} -[\lambda-2q\ cos(2x)]\ y_0 \\ y_0(x0) 1\\ y_1(x0) 0\\ y_1(x\pi)0 \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧y0′y1y1′−[λ−2q cos(2x)] y0y0(x0)1y1(x0)0y1(xπ)0
这样就可以用 solve_bvp() 求解该常微分方程的边值问题。
5.2 常微分方程的编程步骤
以该题为例讲解scipy.integrate.solve_bvp 求解常微分方程边值问题的步骤。
需要注意的是1本案例涉及一个待定参数 λ\lambdaλ 需要通过 solve_bvp(fun, bc, x, y, pNone) 中的可选项 p 传递到导数函数和边界条件函数2本案例涉及 3 个边界条件要注意边界条件函数的定义。 导入 scipy、numpy、matplotlib 包 定义导数函数 dydx(x,y,p) 本问题中 y 表示向量记为 y[y0,y1]y[y_0,y_1]y[y0,y1]定义函数 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,-,labelA {:.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(locbest)
plt.text(2,-4,youcans-xupt,colorwhitesmoke)
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)。
Copyright 2021 Youcans, XUPT
Crated2021-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数模笔记-模拟退火算法
