Python小白的數學建模課-09 微分方程模型

youcans發表於2021-06-21

Python模型

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

本文介紹微分方程模型的建模與求解，通過常微分方程、常微分方程組、高階常微分方程 3個案例手把手教你搞定微分方程。

通過二階 RLC 電路問題，學習微分方程模型的建模、求解和討論。

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

1. 微分方程

1.1 基本概念

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

具體來說，微分方程是指含有未知函式及其導數的關係式。

微分方程按自變數個數分為：只有一個自變數的常微分方程（Ordinary Differential Equations）和包含兩個或兩個以上獨立變數的偏微分方程（Partial Differential Equations）。
微分方程按階數分為：一階、二階、高階，微分方程的階數取決於方程中最高次導數的階數。
微分方程還可以分為：（非）齊次，常（變）係數，（非）線性，初值問題/邊界問題...

以上內容看看就算了，看多了就嚇跑了。

1.2 微分方程的數學建模

微分方程的數學建模其實並不複雜，基本過程就是分析題目屬於哪一類問題、可以選擇什麼微分方程模型，然後如何使用現有的微分方程模型建模。

在數學、力學、物理、化學等各個學科領域的課程中，針對該學科的各種問題都會建立適當的數學模型。在中學課程中，各學科的數學模型主要是線性或非線性方程，而在大學物理和各專業的課程中，越來越多地出現用微分方程描述的數學模型。

數學建模中的微分方程問題，通常還是這些專業課程中相對簡單的模型，專業課程的教材在介紹一個模型時，往往都做了非常詳細的講解。只要搞清楚問題的型別、選擇好數學模型，建模和求解並不是很難，而且在撰寫論文時對問題背景、使用範圍、假設條件、求解過程有大量現成的內容可以複製參考。

小白之所以害怕，一是看到微分方程就心裡發怵，二是缺乏專業背景，不知道從哪裡查資料、不能判斷問題的型別、不知道選擇什麼模型、不善於從題目內容得出模型引數，也不知道如何程式設計求解。所以，老師說，一看這就是××問題，顯然就可以用××模型。小白說，我們還是換 B題吧。

本系列將會從簡單的微分方程模型入手，重點介紹微分方程數值解法的程式設計實現，並通過分析問題、建立模型的案例幫助小白樹立信心和動力。

希望你在學習本系列之後，會發現微分方程模型是數學建模中最容易的題型：模型找教材，建模找例題，求解有例程，討論有套路，論文夠檔次。

1.3 微分方程的數值解法

在學習專業課程時，經常會推導和求解微分方程的解析解，小白對微分方程模型的恐懼就是從高等數學“微分方程”開始，經過專業課的不斷強化而形成的。實際上，只有很少的微分方程可以解析求解，大多數的微分方程只能採用數值方法進行求解。

微分方程的數值求解是先把時間和空間離散化，然後將微分化為差分，建立遞推關係，然後反覆進行迭代計算，得到任意時間和空間的值。

如果你還是覺得頭暈目眩，我們可以說的更簡單一些。建模就是把專業課教材上的公式抄下來，求解就是把公式的引數輸入到 Python 函式中。

我們先說求解。求解常微分方程的基本方法，有尤拉法、龍格庫塔法等，可以詳見各種教材，撰寫數模競賽論文時還是可以抄幾段的。本文沿用“程式設計方案”的概念，不涉及這些演算法的具體內容，只探討如何使用 Python 的工具包、庫函式，零基礎求解微分方程模型。

我們的選擇是 Python 常用工具包三劍客：Scipy、Numpy 和 Matplotlib：

Scipy 是 Python 演算法庫和數學工具包，包括最優化、線性代數、積分、插值、特殊函式、傅立葉變換、訊號和影像處理、常微分方程求解等模組。有人介紹 Scipy 就是 Python 語言的 Matlab，所以大部分數學建模問題都可以用它搞定。
Numpy 提供了高維陣列的實現與計算的功能，如線性代數運算、傅立葉變換及隨機數生成，另外還提供了與 C/C++ 等語言的整合工具。
Matplotlib 是視覺化工具包，可以方便地繪製各種資料視覺化圖表，如折線圖、散點圖、直方圖、條形圖、箱形圖、餅圖、三維圖，等等。

順便說一句，還有一個 Python 符號運算工具包 SymPy，以解析方式求解積分、微分方程，也就是說給出的結果是微分方程的解析解表示式。很牛，但只能求解有解析解的微分方程，所以，你知道就可以了。

2. SciPy 求解常微分方程（組）

2.1 一階常微分方程（組）模型

給定初始條件的一階常微分方程（組）的標準形式是：

\[\begin{cases} \begin{aligned} &\frac{dy}{dt} = f(y,t)\\ &y(t_0) = y_0 \end{aligned} \end{cases} \]

式中的 y 在常微分方程中是標量，在常微分方程組中是陣列向量。

2.2 scipy.integrate.odeint() 函式

SciPy 提供了兩種方式求解常微分方程：基於 odeint 函式的 API 比較簡單易學，基於 ode 類的物件導向的 API 更加靈活。

**scipy.integrate.odeint() **是求解微分方程的具體方法，通過數值積分來求解常微分方程組。在 odeint 函式內部使用 FORTRAN 庫 odepack 中的 lsoda，可以求解一階剛性系統和非剛性系統的初值問題。官網介紹詳見： scipy.integrate.odeint — SciPy v1.6.3 Reference Guide 。

scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0, tfirst=False)

odeint 的主要引數：

求解標準形式的微分方程（組）主要使用前三個引數：

func: callable(y, t, …) 　　導數函式 $f(y,t)$ ，即 y 在 t 處的導數，以函式的形式表示
y0: array：　　初始條件 $y_0$，對於常微分方程組 $y_0$ 則為陣列向量
t: array：　　求解函式值對應的時間點的序列。序列的第一個元素是與初始條件 $y_0$ 對應的初始時間 $t_0$；時間序列必須是單調遞增或單調遞減的，允許重複值。

其它引數簡介如下：

args: 嚮導數函式 func 傳遞引數。當導數函式 $f(y,t,p1,p2,..)$ 包括可變引數 p1,p2.. 時，通過 args =(p1,p2,..) 可以將引數p1,p2.. 傳遞給導數函式 func。argus 的用法參見 2.4 中的例項2。
Dfun: func 的雅可比矩陣，行優先。如果 Dfun 未給出，則演算法自動推導。
col_deriv: 自動推導 Dfun的方式。
printmessg: 布林值。控制是否列印收斂資訊。
其它引數用於控制求解演算法的引數，一般情況可以忽略。

odeint 的主要返回值：

y: array 　　陣列，形狀為 (len(t),len(y0)，給出時間序列 t 中每個時刻的 y 值。

3. 例項1：Scipy 求解一階常微分方程

3.1 例題 1：求微分方程的數值解

\[\begin{cases} \begin{aligned} &\frac{dy}{dt} = sin(t^2)\\ &y(-10) = 1 \end{aligned} \end{cases} \]

3.2 常微分方程的程式設計步驟

以該題為例講解 scipy.integrate.odeint() 求解常微分方程初值問題的步驟：

匯入 scipy、numpy、matplotlib 包；
定義導數函式 $f(y,t)=sin(t^2)$ ；
定義初值 $y_0$ 和 $y$ 的定義區間 $[t_0,\ t]$；
呼叫 odeint() 求 $y$ 在定義區間 $[t_0,\ t]$ 的數值解。

3.3 Python 例程

# 1. 求解微分方程初值問題(scipy.integrate.odeint)
from scipy.integrate import odeint  # 匯入 scipy.integrate 模組
import numpy as np
import matplotlib.pyplot as plt

def dy_dt(y, t):  # 定義函式 f(y,t)
    return np.sin(t**2)

y0 = [1]  # y0 = 1 也可以
t = np.arange(-10,10,0.01)  # (start,stop,step)
y = odeint(dy_dt, y0, t)  # 求解微分方程初值問題

# 繪圖
plt.plot(t, y)
plt.title("scipy.integrate.odeint")
plt.show()

3.4 Python 例程執行結果

4. 例項2：Scipy 求解一階常微分方程組

4.1 例題 2：求洛倫茲（Lorenz）方程的數值解

洛倫茲（Lorenz）混沌吸引子的軌跡可以由如下的 3個微分方程描述：

\[\begin{cases} \begin{aligned} &\frac{dx}{dt} = \sigma (y-x)\\ &\frac{dy}{dt} = x (\rho-z) - y\\ &\frac{dz}{dt} = xy - \beta z\\ \end{aligned} \end{cases} \]

洛倫茲方程將大氣流體運動的強度 x 與水平和垂直方向的溫度變化 y 和 z 聯絡起來，進行大氣對流系統的模擬，現已廣泛應用於天氣預報、空氣汙染和全球氣候變化的研究。引數 $\sigma$ 稱為普蘭特數，$\rho$ 是規範化的瑞利數，$\beta$ 和幾何形狀相關。洛倫茲方程是非線性微分方程組，無法求出解析解，只能使用數值方法求解。

4.2 洛倫茲（Lorenz）方程問題的程式設計步驟

以該題為例講解 scipy.integrate.odeint() 求解常微分方程初值問題的步驟：

匯入 scipy、numpy、matplotlib 包；
定義導數函式 lorenz(W, t, p, r, b)

注意 odeint() 函式中定義導數函式的標準形式是 $f(y,t)$ ，對於微分方程組 y 表示向量。

為避免混淆，我們記為 $W=[x,y,z]$，函式 lorenz(W,t) 定義導數函式 $f(W,t)$ 。

用 p,r,b 分別表示方程中的引數 $\sigma、\rho、\beta$，則對導數定義函式程式設計如下：

# 導數函式，求 W=[x,y,z] 點的導數 dW/dt
def lorenz(W,t,p,r,b):
    x, y, z = W  # W=[x,y,z]
    dx_dt = p*(y-x)  # dx/dt = p*(y-x), p: sigma
    dy_dt = x*(r-z) - y  # dy/dt = x*(r-z)-y, r:rho
    dz_dt = x*y - b*z  # dz/dt = x*y - b*z, b;beta
    return np.array([dx_dt,dy_dt,dz_dt])

定義初值 $W_0$ 和 $W$ 的定義區間 $[t_0,\ t]$；
呼叫 odeint() 求 $W$ 在定義區間 $[t_0,\ t]$ 的數值解。

注意例程中通過 args=paras 或 args = (10.0,28.0,3.0) 將引數 (p,r,b) 傳遞給導數函式 lorenz(W,t,p,r,b)。引數 (p,r,b) 當然也可以不作為函式引數傳遞，而是在導數函式 lorenz() 中直接設定。但例程的引數傳遞方法，使導數函式結構清晰、更為通用。另外，對於可變引數問題，使用這種引數傳遞方式就非常方便。

4.3 洛倫茲（Lorenz）方程問題 Python 例程

# 2. 求解微分方程組初值問題(scipy.integrate.odeint)
from scipy.integrate import odeint    # 匯入 scipy.integrate 模組
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 導數函式, 求 W=[x,y,z] 點的導數 dW/dt
def lorenz(W,t,p,r,b):  # by youcans
    x, y, z = W  # W=[x,y,z]
    dx_dt = p*(y-x)  # dx/dt = p*(y-x), p: sigma
    dy_dt = x*(r-z) - y  # dy/dt = x*(r-z)-y, r:rho
    dz_dt = x*y - b*z  # dz/dt = x*y - b*z, b;beta
    return np.array([dx_dt,dy_dt,dz_dt])

t = np.arange(0, 30, 0.01)  # 建立時間點 (start,stop,step)
paras = (10.0, 28.0, 3.0)  # 設定 Lorenz 方程中的引數 (p,r,b)

# 呼叫ode對lorenz進行求解, 用兩個不同的初始值 W1、W2 分別求解
W1 = (0.0, 1.00, 0.0)  # 定義初值為 W1
track1 = odeint(lorenz, W1, t, args=(10.0, 28.0, 3.0))  # args 設定導數函式的引數
W2 = (0.0, 1.01, 0.0)  # 定義初值為 W2
track2 = odeint(lorenz, W2, t, args=paras)  # 通過 paras 傳遞導數函式的引數

# 繪圖
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(track1[:,0], track1[:,1], track1[:,2], color='magenta') # 繪製軌跡 1
ax.plot(track2[:,0], track2[:,1], track2[:,2], color='deepskyblue') # 繪製軌跡 2
ax.set_title("Lorenz attractor by scipy.integrate.odeint")
plt.show()

4.4 洛倫茲（Lorenz）方程問題 Python 例程執行結果

5. 例項3：Scipy 求解高階常微分方程

高階常微分方程，必須做變數替換，化為一階微分方程組，再用 odeint 求數值解。

5.1 例題 3：求二階 RLC 振盪電路的數值解

零輸入響應的 RLC 振盪電路可以由如下的二階微分方程描述：

\[\begin{cases} \begin{aligned} &\frac{d^2 u}{dt^2} + \frac{R}{L} * \frac{du}{dt} + \frac{1}{LC}*u = 0\\ &u(0) = U_0\\ &u'(0)= 0 \end{aligned} \end{cases} \]

令 $ \alpha = R/2L$、$\omega_0^2=1/LC$，在零輸入響應 $u_s=0$ 時上式可以寫成：

\[\begin{cases} \begin{aligned} &\frac{d^2 u}{dt^2} + 2 \alpha \frac{du}{dt} + \omega_0^2 u = 0\\ &u(0) = U_0\\ &u'(0)= 0 \end{aligned} \end{cases} \]

對二階微分方程問題，引入變數 $v = {du}/{dt}$，通過變數替換就把原方程化為如下的微分方程組：

\[\begin{cases} \begin{aligned} &\frac{du}{dt} = v \\ &\frac{dv}{dt} = -2\alpha v - \omega_0^2 u\\ &u(0)=U_0\\ &v(0)=0 \end{aligned} \end{cases} \]

這樣就可以用上節求解微分方程組的方法來求解高階微分方程問題。

5.2 二階微分方程問題的程式設計步驟

以RLC 振盪電路為例講解 scipy.integrate.odeint() 求解高階常微分方程初值問題的步驟：

匯入 scipy、numpy、matplotlib 包；
定義導數函式 deriv(Y, t, a, w)

注意 odeint() 函式中定義導數函式的標準形式是 $f(y,t)$ ，本問題中 y 表示向量，記為 $Y=[u,v]$

導數定義函式 deriv(Y, t, a, w) 程式設計如下，其中 a, w 分別表示方程中的引數 $\alpha、\omega$：

# 導數函式，求 Y=[u,v] 點的導數 dY/dt
def deriv(Y, t, a, w):
    u, v = Y  # Y=[u,v]
    dY_dt = [v, -2*a*v-w*w*u]
    return dY_dt

定義初值 $Y_0=[u_0,v_0]$ 和 $Y$ 的定義區間 $[t_0,\ t]$；
呼叫 odeint() 求 $Y=[u,v]$ 在定義區間 $[t_0,\ t]$ 的數值解。

例程中通過 args=paras 將引數 (a,w) 傳遞給導數函式 deriv(Y, t, a, w) 。本例要考察不同引數對結果的影響，這種引數傳遞方法使用非常方便。

5.3 二階微分方程問題 Python 例程

# 3. 求解二階微分方程初值問題(scipy.integrate.odeint)
# Second ODE by scipy.integrate.odeint
from scipy.integrate import odeint  # 匯入 scipy.integrate 模組
import numpy as np
import matplotlib.pyplot as plt

# 導數函式，求 Y=[u,v] 點的導數 dY/dt
def deriv(Y, t, a, w):
    u, v = Y  # Y=[u,v]
    dY_dt = [v, -2*a*v-w*w*u]
    return dY_dt

t = np.arange(0, 20, 0.01)  # 建立時間點 (start,stop,step)
# 設定導數函式中的引數 (a, w)
paras1 = (1, 0.6)  # 過阻尼：a^2 - w^2 > 0
paras2 = (1, 1)  # 臨界阻尼：a^2 - w^2 = 0
paras3 = (0.3, 1)  # 欠阻尼：a^2 - w^2 < 0

# 呼叫ode對進行求解, 用兩個不同的初始值 W1、W2 分別求解
Y0 = (1.0, 0.0)  # 定義初值為 Y0=[u0,v0]
Y1 = odeint(deriv, Y0, t, args=paras1)  # args 設定導數函式的引數
Y2 = odeint(deriv, Y0, t, args=paras2)  # args 設定導數函式的引數
Y3 = odeint(deriv, Y0, t, args=paras3)  # args 設定導數函式的引數
# W2 = (0.0, 1.01, 0.0)  # 定義初值為 W2
# track2 = odeint(lorenz, W2, t, args=paras)  # 通過 paras 傳遞導數函式的引數

# 繪圖
plt.plot(t, Y1[:, 0], 'r-', label='u1(t)')
plt.plot(t, Y2[:, 0], 'b-', label='u2(t)')
plt.plot(t, Y3[:, 0], 'g-', label='u3(t)')
plt.plot(t, Y1[:, 1], 'r:', label='v1(t)')
plt.plot(t, Y2[:, 1], 'b:', label='v2(t)')
plt.plot(t, Y3[:, 1], 'g:', label='v3(t)')
plt.axis([0, 20, -0.8, 1.2])
plt.legend(loc='best')
plt.title("Second ODE by scipy.integrate.odeint")
plt.show()

5.4 二階方程問題 Python 例程執行結果

結果討論：

RLC串聯電路是典型的二階系統，在零輸入條件下根據 $\alpha$ 與 $\omega$ 的關係，電路的輸出響應存在四種情況：

過阻尼： $\alpha^2 - \omega^2>0$ ，有 2 個不相等的負實數根；
臨界阻尼： $\alpha^2 - \omega^2 = 0$，有 2 個相等的負實數根；
欠阻尼： $\alpha^2 - \omega^2 <0$，有一對共軛複數根；
無阻尼：$R=0$，有一對純虛根。

例程中所選擇的 3 組引數分別對應過阻尼、臨界阻尼和欠阻尼的條件，微分方程的數值結果很好地體現了不同情況的相應曲線。

6. 小結

小白首先要有信心，用 Scipy 工具包求解標準形式的微分方程（組），程式設計實現是非常簡單、容易掌握的。
其次要認識到，由於微分方程的解的特性是與模型引數密切相關的，不同引數的解可能具有完全不同的形態。這就涉及模型穩定性、靈敏度的分析，很容易使論文寫的非常豐富和精彩。
不會從問題建立微分方程模型怎麼辦，不會展開引數對穩定性、靈敏度的影響進行討論怎麼辦？誰讓你自己做呢，當然是先去找相關專業的教材、論文，從中選擇比較接近、比較簡單的理論和模型，然後通過各種假設強行將題目簡化為模型中的條件，這就可以照貓畫虎了。

【本節完】

版權宣告：

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

原創作品，轉載必須標註原文連結：https://www.cnblogs.com/youcans/p/14912966.html。

Crated：2021-06-08

歡迎關注『Python小白的數學建模課 @ Youcans』，每週更新數模筆記
Python小白的數學建模課-01.新手必讀
 Python小白的數學建模課-02.資料匯入
 Python小白的數學建模課-03.線性規劃
 Python小白的數學建模課-04.整數規劃
 Python小白的數學建模課-05.0-1規劃
 Python小白的數學建模課-06.固定費用問題
 Python小白的數學建模課-07.選址問題
 Python小白的數學建模課-09.微分方程模型
 Python數模筆記-PuLP庫
 Python數模筆記-StatsModels統計迴歸
 Python數模筆記-Sklearn
Python數模筆記-NetworkX
Python數模筆記-模擬退火演算法

Python小白的數學建模課-B5. 新冠疫情 SEIR模型
2021-07-10
Python模型
Python小白的數學建模課-04.整數規劃
2021-06-03
Python
Python小白的數學建模課-05.0-1規劃
2021-06-06
Python
Python小白的數學建模課-07 選址問題
2021-06-13
Python
Python小白的數學建模課-06 固定費用問題
2021-06-10
Python
Python小白的數學建模課-03.線性規劃
2021-06-01
Python
Python小白的數學建模課-16.最短路徑演算法
2021-08-06
Python演算法
Python小白的數學建模課-19.網路流優化問題
2021-08-24
Python優化
微分方程在遊戲經濟建模中的應用
2023-05-15
遊戲
【小白學PyTorch】4 構建模型三要素與權重初始化
2020-09-03
PyTorch模型
Python數學建模-01.新手必讀
2021-05-29
Python
Python數學建模-02.資料匯入
2021-05-30
Python
微分方程數值解法的matlab程式
2020-09-24
Matlab
小白學python系列-(3)基礎數量型別
2018-04-05
Python型別
引數匹配模型——Python學習之引數（二）
2020-04-04
模型Python
Arduino和Python螺旋槳角度數學模型
2020-12-03
UIPython模型
遊戲數學建模工程手冊
2020-05-18
遊戲
基於大資料模型的數字孿生建模方法
2022-10-17
大資料模型
用 Python 進行貝葉斯模型建模（4）
2016-08-11
Python模型
小白學 Python（1）：開篇
2019-10-18
Python
python適合小白學嗎？
2021-09-11
Python
小白學python系列-(4)list
2018-04-06
Python
小白學python系列-(5) tuple
2018-04-06
Python
小白學python系列-(8)dict
2018-04-06
Python
漢羅塔問題java數學模型（從數學原型到數學模型）
2017-10-28
Java模型原型
Python入門_給小白的學習指南
2020-02-06
Python
華中杯數學建模 A題簡單覆盤（附Python原始碼）
2022-08-09
Python原始碼
實驗四常微分方程的數值解（android）
2016-05-11
Android
小白學python遇到問題了.......
2019-09-28
Python
新手小白學Python有前途嗎?
2019-08-15
Python
python 裝飾器小白學習
2018-05-28
Python
小白學python系列-(2)hello word
2018-04-05
Python
小白學python系列-(7)迴圈
2018-04-06
Python
小白學python系列-(9)函式
2018-04-14
Python函式
小白的演算法課-方法論
2021-04-26
演算法
數學模型——數學與人類文明的橋樑
2021-01-02
模型
python小白之路：第十九章 Boosting模型（二）
2020-10-14
Python模型
數學建模 # 1 論文撰寫技巧
2020-12-12