小白往往聽到微分方程就覺得害怕，其實數學建模中的微分方程模型不僅沒那么復雜，而且很容易寫出高水平的數模論文。

本文介紹微分方程模型邊值問題的建模與求解，不涉及算法推導和編程，只探討如何使用 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{}^{\prime \prime }=f(x,y,y{}^{\prime })，a<x<b$$

有三種情況的邊界條件：

（1）第一類邊界條件（兩點邊值問題）：

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

（2）第二類邊界條件：

$$y{}^{\prime }(a)=ya，y{}^{\prime }(b)=yb$$

（3）第三類邊界條件：
$$\{\begin{array}{l}y{}^{\prime }(a)?a_{0}y(a)=a_1\\ y{}^{\prime }(b)?b_{0}y(b)=b_1\end{array}$$

其中：$a_0\ge 0，b_0\ge 0，a_0+b_0>0$

1.3 常微分方程邊值問題的數值解法

簡單介紹求解常微分方程邊值問題的數值解法，常用方法有：打靶算法、有限差分法和有限元法。打靶算法把邊值問題轉化為初值問題求解，是根據邊界條件反復迭代調整初始點的斜率，使初值問題的數值解在邊界上“命中”問題的邊值條件。有限差分法把空間離散為網格節點，用差商代替微商，將微分方程離散化為線性或非線性方程組來求解。 有限元法將微分方程離散化，有限元就是指近似連續域的離散單元，對每一單元假定一個近似解，然后推導求解域滿足條件，從而得到問題的解。

按照本系列“編程方案”的概念，不涉及這些算法的具體內容，只探討如何使用 Python 的工具包、庫函數，零基礎求解微分方程模型邊值問題。我們的選擇還是 Python 常用工具包三劍客：Scipy、Numpy 和 Matplotlib。

---

2. SciPy 求解常微分方程邊值問題

2.1 BVP 問題的標準形式

Scipy 用 solve_bvp() 函數求解常微分方程的邊值問題，定義微分方程的標準形式為：
$$\{\begin{array}{l}y{}^{\prime }=f(x,y)，a<x<b\\ g(y(a),y(b)=0)\end{array}$$

因此要將第一類邊界條件 $y(a)=ya，y(b)=yb$ 改寫為：
$$\{\begin{array}{l}y(a)?ya=0\\ y(b)?yb=0\end{array}$$

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：一階常微分方程邊值問題

求常微分方程邊值問題的數值解。

$$?\begin{array}{l}\begin{array}{rl}& y{}^{\prime \prime }+\left|y\right|=0\\ & y(x=0)=0.5\\ & y(x=4)=?1.5\end{array}\end{array}$$

引入變量 $y0=y，y1=y{}^{\prime }$，通過變量替換就把原方程化為如下的標準形式的微分方程組：

$$?\begin{array}{l}{y}_0^{{}^{\prime }}=y_1\\ y_1^{{}^{\prime }}=?\left|y_0\right|\\ y(x=0)?0.5=0\\ y(x=4)+1.5=0\end{array}$$

這樣就可以用 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：水滴橫截面形狀問題

水平面上的水滴橫截面形狀，可以用如下的微分方程描述：
$$?\begin{array}{l}\begin{array}{rl}& \frac{d^{2}h}{d{x}^2}+[1?h]?[1+(\frac{dh}{dx}{)}^2{]}^{3/2}=0\\ & h(x=?1)=h(x=1)=0\end{array}\end{array}$$

引入變量 $h0=h，h1=h{}^{\prime }$，通過變量替換就把原方程化為如下的標準形式的微分方程組：

$$?\begin{array}{l}{h}_0^{{}^{\prime }}=h_1\\ h_1^{{}^{\prime }}=(h_0?1)?[1+h_1^2{]}^{3/2}\\ h_0(x=?1)=h_0(x=1)=0\end{array}$$

這樣就可以用 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 在研究橢圓形膜的邊界值問題時，導出了一個二階常微分方程，其形式為：
$$\frac{d^{2}y}{d{x}^2}+[\lambda ?2qcos(2x)]y=0$$

用這種形式的數學方程可以描述自然中的物理現象，包括振動橢圓鼓、四極質譜儀和四極離子阱、周期介質中的波動、強制振蕩器參數共振現象、廣義相對論中的平面波解決方案、量子擺哈密頓函數的本征函數、旋轉電偶極子的斯塔克效應。

式中 $\lambda 、q$ 是兩個實參數，方程的系數是以 $\pi$ 或 $2\pi$ 為周期的，但只有在 滿足一定關系時 Mathieu 方程才有周期解。

引入變量 $y0=y，y1=y{}^{\prime }$，通過變量替換就把原方程化為如下的標準形式的微分方程組：

$$?\begin{array}{l}{y}_0^{{}^{\prime }}=y_1\\ y_1^{{}^{\prime }}=?[\lambda ?2qcos(2x)]y_0\\ y_0(x=0)=1\\ y_1(x=0)=0\\ y_1(x=\pi )=0\end{array}$$

這樣就可以用 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)。

Copyright 2021 Youcans, XUPT

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數模筆記-模擬退火算法

總結

以上是生活随笔為你收集整理的Python小白的数学建模课-10.微分方程边值问题的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。