10.2 專案二：ODE 建模#

本節目標#

用三個實際建模專案學會微分方程工作流： 數學建模 → SymPy 求解 → 畫圖 → 解讀
人口成長：logistic 方程的 dsolve 解與畫圖
彈簧振動： \(my'' + cy' + ky = 0\) 的三種阻尼行為
SIR 疫情模型：方程組求解與參數實驗
學會「改參數、看結果」的實驗方法

專案 1：人口成長（Logistic 方程）#

數學建模#

馬爾薩斯模型 \(\dfrac{dP}{dt} = rP\) 預測指數爆炸，不符現實。 Logistic 模型加入環境承載量 \(K\) ：

\[ \frac{dP}{dt} = rP\left(1 - \frac{P}{K}\right) \]

\(P\) 很小時近似指數成長； \(P\) 接近 \(K\) 時成長放緩。

程式求解#

import sympy as sp

t = sp.Symbol('t', positive=True)
P = sp.Function('P')
r, K, P0 = sp.symbols('r K P0', positive=True)

ode = sp.Eq(P(t).diff(t), r*P(t)*(1 - P(t)/K))
sol = sp.dsolve(ode, P(t), ics={P(0): P0})
print(sp.simplify(sol.rhs))
# P0*K*exp(r*t) / (K + P0*(exp(r*t) - 1))

解讀：通解是 S 形曲線—— \(t \to \infty\) 時 \(P \to K\) 。驗證：對解取極限。

print(sp.limit(sol.rhs, t, sp.oo))   # K

畫圖#

import matplotlib.pyplot as plt
import numpy as np

sol_num = sp.lambdify(t, sol.rhs.subs({r: 0.5, K: 1000, P0: 10}))
tt = np.linspace(0, 25, 300)
plt.plot(tt, sol_num(tt))
plt.axhline(1000, ls='--', color='gray', label='K = 1000')
plt.xlabel('t'); plt.ylabel('P(t)'); plt.legend()
plt.title('Logistic 成長曲線')
plt.show()

曲線先緩升、再陡增、最後貼平到 \(K\) ——標準 S 形。

專案 2：彈簧振動（三種阻尼）#

數學建模#

質量 \(m\) 的物體掛在彈簧（彈性係數 \(k\) ）上，阻尼係數 \(c\) ：

\[ my'' + cy' + ky = 0 \]

特徵方程 \(m\lambda^2 + c\lambda + k = 0\) 的判別式決定三種行為：

條件	名稱	行為
\(c^2 < 4mk\)	欠阻尼	來回振動、振幅遞減
\(c^2 = 4mk\)	臨界阻尼	最快回到平衡、不振盪
\(c^2 > 4mk\)	過阻尼	緩慢爬回平衡、不振盪

程式求解#

import sympy as sp

t = sp.Symbol('t', positive=True)
y = sp.Function('y')
m, c, k = sp.symbols('m c k', positive=True)

ode = sp.Eq(m*y(t).diff(t, 2) + c*y(t).diff(t) + k*y(t), 0)
sol = sp.dsolve(ode, y(t), ics={y(0): 1, y(t).diff(t).subs(t, 0): 0})
print(sol.rhs)

SymPy 會以指數與餘弦/正弦的形式給出通解，形式隨判別式而變。

三種阻尼各解一次（數值代入）#

def solve_case(m, c, k):
    ode = sp.Eq(m*y(t).diff(t, 2) + c*y(t).diff(t) + k*y(t), 0)
    s = sp.dsolve(ode, y(t), ics={y(0): 1, y(t).diff(t).subs(t, 0): 0})
    return sp.lambdify(t, s.rhs, 'numpy')

under  = solve_case(1, 0.5, 4)    # c² = 0.25 < 16 → 欠阻尼
crit   = solve_case(1, 4,   4)    # c² = 16 = 16 → 臨界阻尼
over   = solve_case(1, 10,  4)    # c² = 100 > 16 → 過阻尼

畫圖比較#

import numpy as np
import matplotlib.pyplot as plt

tt = np.linspace(0, 10, 400)
plt.plot(tt, under(tt),  label='欠阻尼 c=0.5')
plt.plot(tt, crit(tt),   label='臨界 c=4')
plt.plot(tt, over(tt),   label='過阻尼 c=10')
plt.axhline(0, color='gray', lw=0.5)
plt.xlabel('t'); plt.ylabel('y(t)'); plt.legend()
plt.title('m=1, k=4，三種阻尼')
plt.show()

解讀：欠阻尼有明顯振盪；臨界阻尼最快且不越過平衡點；過阻尼反而更慢——這就是為什麼車門緩衝器設計在臨界附近。

專案 3：SIR 疫情模型#

數學建模#

人口分成三類：易感 \(S\) 、感染 \(I\) 、康復 \(R\) ：

\[ \frac{dS}{dt} = -\beta SI, \quad \frac{dI}{dt} = \beta SI - \gamma I, \quad \frac{dR}{dt} = \gamma I \]

\(\beta\) 是傳染率、 \(\gamma\) 是康復率。 \(R_0 = \beta S_0 / \gamma\) 決定疫情是否爆發。

程式求解#

非線性方程組通常沒有解析解，SymPy 解不出閉式時，改用 dsolve 嘗試後退回數值解（sp.dsolve 搭配 ics 失敗時，實務上直接用 lambdify + 自行積分，或改用 SciPy，見 10.4）。這裡用符號推導基本再生數，再以數值方法畫圖：

import sympy as sp

S, I, R = sp.symbols('S I R', positive=True)
beta, gamma = sp.symbols('beta gamma', positive=True)

# 符號分析：疫情初期 I 增長的條件
dI = beta*S*I - gamma*I
# 爆發條件：dI/dt > 0 ⟺ beta*S > gamma
R0_expr = beta*S/gamma
print(sp.simplify(R0_expr))   # beta*S/gamma，即 R₀

數值畫圖（SymPy 符號化 + 手寫 Euler 步進）#

import numpy as np
import matplotlib.pyplot as plt

def simulate(beta_v, gamma_v, S0=990, I0=10, days=120):
    Sv, Iv, Rv = S0, I0, 0.0
    hist = [(Sv, Iv, Rv)]
    for _ in range(days):
        new = beta_v*Sv*Iv/1000          # 標準化人口 N=1000
        rec  = gamma_v*Iv
        Sv, Iv, Rv = Sv-new, Iv+new-rec, Rv+rec
        hist.append((Sv, Iv, Rv))
    return np.array(hist)

h = simulate(beta_v=0.3, gamma_v=0.1)    # R₀ = 3 > 1 → 爆發
plt.plot(h[:, 0], label='S')
plt.plot(h[:, 1], label='I')
plt.plot(h[:, 2], label='R')
plt.xlabel('天'); plt.ylabel('人數'); plt.legend()
plt.title('SIR 模型（R₀=3）')
plt.show()

實驗：改參數看疫情峰值#

for b in [0.1, 0.2, 0.3, 0.5]:
    h = simulate(beta_v=b, gamma_v=0.1)
    peak_day = h[:, 1].argmax()
    print(f"beta={b}: 峰值感染 {h[:, 1].max():.0f} 人，第 {peak_day} 天")

預期結果與解讀：

beta=0.1： \(R_0 = 1\) ，感染人數緩慢下降，幾乎沒有峰值
beta=0.2： \(R_0 = 2\) ，峰值數百人
beta=0.3： \(R_0 = 3\) ，峰值更高、更早
beta=0.5： \(R_0 = 5\) ，峰值逼近全人口，快速燒完

結論：降低 \(\beta\) （戴口罩、保持距離）能把峰值壓平拉後—— 這正是防疫政策「壓平曲線」背後的數學。

建模紀律：每次只改一個參數，記下結論再改下一個。同時改兩個參數，你會分不清是誰造成結果。

三專案共通的工作流回顧#

數學：寫下微分方程與初始條件
符號：dsolve 求解析解；解不出就轉數值
畫圖：lambdify + Matplotlib 看行為
解讀：對照理論預測（S 形、阻尼分類、 \(R_0\) 條件）
實驗：掃參數，找出「什麼改變會造成什麼」

小結與練習#

本節完成三個 ODE 建模專案：logistic 人口成長（解析解 + 極限驗證）、彈簧振動（三種阻尼的系統性比較）、SIR 疫情模型（符號分析 + 數值模擬 + 參數實驗）。核心技能是把數學模型翻譯成 SymPy 程式，並用圖形與參數實驗建立直覺。

練習

修改 logistic 模型加入收穫項： \(\dfrac{dP}{dt} = rP(1 - P/K) - h\) ，用 dsolve 求解並討論 \(h\) 多大時人口會崩潰。
彈簧加上外力 \(F_0\cos(\omega t)\) ，求受迫振動的解，掃 \(\omega\) 觀察共振（提示： \(\omega \to \sqrt{k/m}\) ）。
在 SIR 模型中加入疫苗：每天把固定比例的 \(S\) 直接移到 \(R\) ，模擬並找出「多少接種率能防止爆發」（對照 \(R_0\) 的理論門檻 \(1 - 1/R_0\) ）。
用「一次只改一個參數」的紀律，對彈簧模型掃 \(k = 1, 4, 16\) ，畫三條欠阻尼曲線在同一張圖，驗證勁度越大振動越快。