傳染病擴散的 Logistic 微分方程

封存文章

這是一個取材自 A First Course in Mathematical Modeling 的 ODE 教學例子。單一 Logistic 狀態不能表示易感、感染、康復與再感染，因此本文把 $N$ 稱為「模型中的累積規模」，不把它當成可用於公共衞生預測的流行病模型。

考慮自治方程

\frac{dN}{dt}=0.25N(10-N),\qquad N(0)=2,

其中 $N$ 以百人為單位。這條方程是 carrying capacity 為 10 的 Logistic model。它很適合示範定性分析、separation of variables 及數值積分，但它沒有接觸網絡、潛伏期、康復、人口流動或 intervention。

本文依次回答：

平衡點在哪裏，是否穩定？
何時增長最快？
解析解如何取得？
Euler 與 RK4 的誤差如何隨步長改變？

一、相線與平衡點

令

f(N)=0.25N(10-N).

平衡點滿足 $f(N)=0$ ，所以

N^\star=0,\qquad N^\star=10.

由符號判斷：

$0<N<10$ 時， $f(N)>0$ ，解向上增長；
$N>10$ 時， $f(N)<0$ ，解向下移動；
若容許 $N<0$ ， $f(N)<0$ ，但負人口沒有實際意義。

線性化導數為

f'(N)=2.5-0.5N.

因此

f'(0)=2.5>0,\qquad f'(10)=-2.5<0.

$N=0$ 在數學上是不穩定平衡點， $N=10$ 是局部漸近穩定平衡點。在非負初始條件下，除 $N_0=0$ 外的解都趨近 10。

import numpy as np
import matplotlib.pyplot as plt

def rhs(t, n):
    return 0.25 * n * (10 - n)

n = np.linspace(0, 15, 500)
plt.plot(n, rhs(0, n), color="black")
plt.axhline(0, color="0.55", linewidth=0.8)
plt.axvline(0, color="0.55", linestyle="--")
plt.axvline(10, color="0.55", linestyle="--")
plt.xlabel(r"$N$")
plt.ylabel(r"$dN/dt$")

由於 $f(N)$ 是開口向下的拋物線，其最大值位於頂點：

N_{\rm fastest} =-\frac{2.5}{2(-0.25)} =5.

亦可由 $f'(N)=0$ 得到同一結果。要注意，這是狀態為 5 時的最大正增長率；若考慮 $N>10$ 的負增長，「變化最快」要先說明是最大導數還是最大 $|dN/dt|$ 。

二、斜率場

對自治方程，斜率只依賴 $N$ ，因此同一水平線上的小線段具有相同斜率：

t_grid = np.linspace(0, 5, 24)
n_grid = np.linspace(0, 15, 24)
T, N = np.meshgrid(t_grid, n_grid)
S = rhs(T, N)

dt = np.ones_like(S)
length = np.sqrt(dt**2 + S**2)

plt.quiver(
    T, N, dt / length, S / length,
    angles="xy", pivot="mid", color="black",
)
plt.xlabel(r"$t$")
plt.ylabel(r"$N$")

圖形應顯示：

$N=0$ 和 $N=10$ 是水平解；
$0<N<10$ 的線段向上；
$N>10$ 的線段向下；
$N=5$ 附近的正斜率最大。

這與相線分析一致。斜率場提供所有初值的局部方向，但不代替解析解或誤差控制。

三、解析解

對 $0<N<10$ 分離變數：

\frac{dN}{N(10-N)}=0.25\,dt.

利用 partial fractions，

\frac{1}{N(10-N)} =\frac1{10} \left(\frac1N+\frac1{10-N}\right).

積分後：

\frac1{10} \ln\left|\frac{N}{10-N}\right| =0.25t+C.

整理得一般 Logistic 解

N(t) =\frac{10} {1+\left(\frac{10-N_0}{N_0}\right)e^{-2.5t}}.

對 $N_0=2$ ，

N(t) =\frac{10}{1+4e^{-2.5t}} =\frac{10e^{2.5t}}{4+e^{2.5t}}.

def exact(t, n0=2.0):
    t = np.asarray(t, dtype=float)
    return 10 / (1 + ((10 - n0) / n0) * np.exp(-2.5 * t))

這個公式亦說明：

$N_0=2$ 或 $7$ 時，解由下方單調接近 10；
$N_0=14$ 時，解由上方單調下降至 10；
不同初值要使用各自的常數，不能拿 $N_0=2$ 的解析曲線去比較 $N_0=7,14$ 並稱為「同一精確解」。

四、最大增長發生時間

當 $N_0=2$ ，最大增長在 $N(t)=5$ 。由解析解：

5=\frac{10}{1+4e^{-2.5t}},

所以

e^{-2.5t}=\frac14, \qquad t_{\rm fastest} =\frac{\ln4}{2.5} \approx0.5545.

無需用 generic root finder；解析式更精確，也展示初始條件如何影響轉折時間。一般 $0<N_0<5$ 時，

t_{\rm fastest} =\frac1{2.5} \ln\left(\frac{10-N_0}{N_0}\right).

若 $N_0\ge5$ ，往後時間區間內的最大正增長率已出現在 $t=0$ ，或解根本由 10 上方向下移動。

五、Euler 法

顯式 Euler 更新為

N_{k+1}=N_k+h\,0.25N_k(10-N_k).

def euler(f, y0, t0, tf, h):
    steps = int(round((tf - t0) / h))
    if not np.isclose(t0 + steps * h, tf):
        raise ValueError("tf - t0 必須是 h 的整數倍")

    t = t0 + h * np.arange(steps + 1)
    y = np.empty(steps + 1, dtype=float)
    y[0] = y0

    for k in range(steps):
        y[k + 1] = y[k] + h * f(t[k], y[k])
    return t, y

舊程式以 $h=1$ 計算後，用 linear interpolation 補出 $N(0.5)$ 。那個值不是 Euler 法在 $t=0.5$ 的一步結果，因為 0.5 根本不在網格上。若要比較 $N(0.5)$ ，應選可整除 0.5 的步長，或只把插值結果清楚標成 post-processing。

for h in (0.5, 0.1, 0.05):
    t, n_num = euler(rhs, 2.0, 0.0, 5.0, h)
    n_true = exact(t)
    max_error = np.max(np.abs(n_num - n_true))
    print(h, max_error)

Euler 法是 global first-order：

\max_k|N(t_k)-N_k|=O(h)

在解足夠光滑且方法位於穩定範圍時成立。把步長減半後，誤差大致減半才是較可靠的 order check。

平衡點附近的穩定性

在穩定平衡點 $N=10$ 附近令 $e=N-10$ 。線性化為

e'=-2.5e.

Euler 更新則是

e_{k+1}=(1-2.5h)e_k.

要有線性穩定性需

|1-2.5h|<1 \quad\Longrightarrow\quad 0<h<0.8.

因此 $h=1$ 不只「比較粗」；它在 carrying capacity 附近已超出 Euler 對線性化系統的穩定區，會產生振盪或不合理行為。舊文報告 $h=1$ 在 $t=5$ 得到 6，主要反映這個數值問題。

六、四階 Runge–Kutta

經典 RK4 以四個 stage 近似一步內的平均斜率：

\begin{aligned} k_1&=f(t_n,N_n),\\ k_2&=f(t_n+h/2,N_n+hk_1/2),\\ k_3&=f(t_n+h/2,N_n+hk_2/2),\\ k_4&=f(t_n+h,N_n+hk_3),\\ N_{n+1}&=N_n+\frac h6(k_1+2k_2+2k_3+k_4). \end{aligned}

def rk4(f, y0, t0, tf, h):
    steps = int(round((tf - t0) / h))
    if not np.isclose(t0 + steps * h, tf):
        raise ValueError("tf - t0 必須是 h 的整數倍")

    t = t0 + h * np.arange(steps + 1)
    y = np.empty(steps + 1, dtype=float)
    y[0] = y0

    for k in range(steps):
        tk, yk = t[k], y[k]
        k1 = f(tk, yk)
        k2 = f(tk + h/2, yk + h*k1/2)
        k3 = f(tk + h/2, yk + h*k2/2)
        k4 = f(tk + h, yk + h*k3)
        y[k + 1] = yk + h*(k1 + 2*k2 + 2*k3 + k4)/6
    return t, y

RK4 在固定終點、足夠光滑和穩定條件下有

\max_k|N(t_k)-N_k|=O(h^4).

同樣地，若 $h=1$ 的網格沒有 $t=0.5$ ，不能把兩個 RK4 節點線性插值後，再把該數字當作 RK4 的四階精度證據。應使用 $h=0.5,0.25,0.125$ 等嵌套網格，並計算 observed order：

p_{\rm obs} =\log_2\left( \frac{E_h}{E_{h/2}} \right).

Euler 應趨近 $p_{\rm obs}\approx1$ ，RK4 應在 round-off 主導前趨近 4。

七、一個可重現比較

def convergence_table(method, steps):
    rows = []
    for h in steps:
        t, y = method(rhs, 2.0, 0.0, 5.0, h)
        error = np.max(np.abs(y - exact(t)))
        rows.append({"h": h, "max_abs_error": error})
    return rows

print(convergence_table(euler, [0.4, 0.2, 0.1, 0.05]))
print(convergence_table(rk4, [0.4, 0.2, 0.1, 0.05]))

圖表應使用相同節點或以精確解在各自節點取值，並保留黑色文字：

t_exact = np.linspace(0, 5, 500)
plt.plot(t_exact, exact(t_exact), color="black", label="解析解")

for method, style, label in [
    (euler, "o--", "Euler"),
    (rk4, "s-.", "RK4"),
]:
    t, y = method(rhs, 2.0, 0.0, 5.0, 0.2)
    plt.plot(t, y, style, label=label)

plt.xlabel(r"$t$")
plt.ylabel(r"$N(t)$")
plt.legend()

八、模型能與不能說甚麼

這條 Logistic 方程可解釋為：

初期規模近似指數增長；
隨 $N$ 接近 10，增長率下降；
半 carrying capacity 時瞬時增長最快。

但在真實感染模型中，感染人數下降通常來自 susceptible depletion、康復、隔離或行為改變。單一 $N$ 方程沒有區分這些機制，也沒有直接產生 incidence、 prevalence 或 $R_0$ 。若研究問題涉及感染期與免疫，至少應考慮 SIR：

\begin{aligned} \dot S&=-\beta SI/M,\\ \dot I&=\beta SI/M-\gamma I,\\ \dot R&=\gamma I. \end{aligned}

即使如此，參數亦需由資料估計並做 identifiability、uncertainty 和 validation。

結論

這個小例子把 ODE 分析的四個層次連在一起：

相線先給出平衡點、方向與穩定性；
解析解確認長期極限與最大增長時間；
Euler 顯示步長與穩定區的重要性；
RK4 以較高階局部資訊改善固定步長精度。

最重要的數值教訓不是「RK4 永遠準確」，而是比較方法時要讓評估時間落在網格上、做步長收斂測試，並把離散誤差與模型誤差分開。

參考資料

Giordano, F. R., Fox, W. P., Horton, S. B., & Weir, M. D. A First Course in Mathematical Modeling. Brooks/Cole.
Hairer, E., Nørsett, S. P., & Wanner, G. (1993). Solving Ordinary Differential Equations I: Nonstiff Problems (2nd ed.). Springer.
Brauer, F., Castillo-Chavez, C., & Feng, Z. (2019). Mathematical Models in Epidemiology. Springer.