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

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

本文介紹微分方程模型的建模與求解，通過常微分方程、常微分方程組、高階常微分方程 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數模筆記-模擬退火演算法

轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/288244.html

標籤：Python

上一篇：Python bytearray/bytes/string區別 - Python零基礎入門教程

下一篇：Python bytes 和 string 相互轉換 - Python零基礎入門教程