Python小白的數學建模課-10.微分方程邊值問題

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

本文介紹微分方程模型邊值問題的建模與求解，不涉及算法推導和編程，只探討如何使用 Python 的工具包，零基礎求解微分方程模型邊值問題。

通過 3個 BVP 案例層層深入，手把手教你搞定微分方程邊值問題。

歡迎關注『Python小白的數學建模課 @ Youcans』系列，每周持續更新

---

1. 常微分方程的邊值問題（BVP）

1.1 基本概念

微分方程是指含有未知函數及其導數的關系式。

微分方程是描述系統的狀態隨時間和空間演化的數學工具。物理中許多涉及變力的運動學、動力學問題，如空氣的阻力為速度函數的落體運動等問題，很多可以用微分方程求解。微分方程在化學、工程學、經濟學和人口統計等領域也有廣泛應用。

微分方程分為初值問題和邊值問題。初值問題是已知微分方程的初始條件，即自變量為零時的函數值，一般可以用歐拉法、龍哥庫塔法來求解。邊值問題則是已知微分方程的邊界條件，即自變量在邊界點時的函數值。

邊值問題的提出和發展，與流體力學、材料力學、波動力學以及核物理學等密切相關，並且在現代控制理論等學科中有重要應用。例如，力學問題中的懸鏈線問題、彈簧振動問題，熱學問題中的導熱細桿問題、細桿端點冷卻問題，流體力學問題、結構強度問題。

上節我們介紹的常微分方程，主要是微分方程的初值問題。本節介紹二階常微分方程邊值問題的建模與求解。

1.2 常微分方程邊值問題的數學模型

只含邊界條件作為定解條件的常微分方程求解問題，稱為常微分方程的邊值問題（boundary value problem）。

一般形式的二階常微分方程邊值問題：

\[y{\ ''} = f(x,y,y{\ '})，\; a<x<b \]

有三種情況的邊界條件：

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

\[y(a)=ya，y(b)=yb \]

（2）第二類邊界條件：

\[y\ '(a)=ya，y\ '(b)=yb \]

（3）第三類邊界條件：

\[\begin{cases} y\ '(a)-a_0\ y(a) = a_1\\ y\ '(b)-b_0\ y(b) = b_1 \end{cases} \]

其中：\(a_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() 函數求解常微分方程的邊值問題，定義微分方程的標准形式為：

\[\begin{cases} y{\ '} = f(x,y)，\; a<x<b\\ g(y(a),y(b)=0) \end{cases} \]

因此要將第一類邊界條件 \(y(a)=ya，y(b)=yb\) 改寫為：

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

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

\[\begin{cases} \begin{aligned} &y\ ''+ \left| y \right| = 0\\ &y(x=0) = 0.5\\ &y(x=4) = -1.5 \end{aligned} \end{cases} \]

引入變量 \(y0 = y，y1 = y\ '\)，通過變量替換就把原方程化為如下的標准形式的微分方程組：

\[\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 求解常微分方程邊值問題的步驟：

1. 導入 scipy、numpy、matplotlib 包；

2. 定義導數函數 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))

3. 定義邊界條件函數 boundCond(ya,yb)

   本問題中邊界條件為常數，具體編程如下。注意對照 3.1中的邊值條件，沒有為什么，函數就是這么定義的。

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

4. 設置 x、y 的初值

5. 調用 solve_bvp() 求解常微分方程在區間 [xa,xb] 的數值解

6. 由 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.5
    fb = -1.5  # 邊界條件 y(xb=4) = -1.5
    return 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)  # 求解 BVP

xSol = 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{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} \]

引入變量 \(h0 = h，h1 = h\ '\)，通過變量替換就把原方程化為如下的標准形式的微分方程組：

\[\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/dx
    dh1 = (h[0]-1)*(1+h[1]*h[1])**1.5  # 計算 dh1/dx
    return np.vstack((dh0, dh1))

# 計算 邊界條件
def boundCond(ha,hb):
    # ha = 0  # 邊界條件：h0(x=-1) = 0
    # hb = 0  # 邊界條件：h0(x=1) = 0
    return 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}{dx^2} + [\lambda-2q\ cos(2x)]\ y= 0 \]

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

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

引入變量 \(y0 = y，y1 = y\ '\)，通過變量替換就把原方程化為如下的標准形式的微分方程組：

\[\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 個邊界條件，要注意邊界條件函數的定義。

1. 導入 scipy、numpy、matplotlib 包；

2. 定義導數函數 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))

3. 定義邊界條件函數 boundCond(ya,yb,p)

   注意，雖然邊界條件定義函數並沒有用到參數 p，但也必須寫在輸入變量中，函數就是這么要求的。

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

4. 設置 x、y 的初值

5. 調用 solve_bvp() 求解常微分方程在區間 [xa,xb] 的數值解

6. 由 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 = 10
    dy0 = 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*k
    y0ini = 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])  # 求解 BVP
    y0 = 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. 小結

1. 微分方程的邊值問題相對初值問題來說更為復雜，但是用 Scipy 工具包求解標准形式的微分方程邊值問題，編程實現還是不難掌握的。
2. 關於邊值問題的模型穩定性、靈敏度的分析，是更為專業的問題。除非找到專業課程教材或范文中有相關內容可以參考套用，否則不建議小白自己摸索，這些問題不是調整參數試試就能試出來的。

【本節完】

版權聲明：

歡迎關注『Python小白的數學建模課 @ Youcans』 原創作品

原創作品，轉載必須標注原文鏈接：。
Copyright 2021 Youcans, XUPT
Crated：2021-06-23

歡迎關注 『Python小白的數學建模課 @ Youcans』，每周更新數模筆記
Python小白的數學建模課-01.新手必讀
Python小白的數學建模課-02.數據導入
Python小白的數學建模課-03.線性規划
Python小白的數學建模課-04.整數規划
Python小白的數學建模課-05.0-1規划
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數模筆記-模擬退火算法