小白往往聽到微分方程就覺得害怕，其實數學建模中的微分方程模型不僅沒那么復雜，而且很容易寫出高水平的數模論文，
本文介紹微分方程模型的建模與求解，通過常微分方程、常微分方程組、高階常微分方程 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，
Copyright 2021 Youcans, XUPT
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/288257.html
標籤：其他
上一篇：作業流專案一般都有哪些模塊？