Python小白的數學建模課-B2. 新冠疫情 SI模型

傳染病的數學模型是數學建模中的典型問題，常見的傳染病模型有 SI、SIR、SIRS、SEIR 模型。

SI 模型是最簡單的傳染病模型，適用於只有易感者和患病者兩類人羣。

咱們就從 SI 模型開始吧，從模型、例程、運行結果到模型分析，全都在這個系列中。

『Python小白的數學建模課 @ Youcans』帶你從數模小白成爲國賽達人。

Python小白的數學建模課-09.微分方程模型
Python小白的數學建模課-B2.新冠疫情 SI模型
Python小白的數學建模課-B3.新冠疫情 SIS模型
Python小白的數學建模課-B4.新冠疫情 SIR模型
Python小白的數學建模課-B5.新冠疫情 SEIR模型
Python小白的數學建模課-B6.改進 SEIR疫情模型

---

1. 前言

新冠疫情不只嚴重影響到全球的政治和經濟，深入和全面地影響着社會和生活的方方面面，也已經成爲數學建模競賽的背景帝。

傳染病的數學模型是數學建模中的典型問題，標準名稱是流行病的數學模型（Mathematical models of epidemic diseases）。創建傳染病的數學模型來描述傳染病的傳播過程，研究傳染病的傳播速度、空間範圍、傳播途徑、動力學機理等問題，以指導對傳染病的有效地預防和控制，具備重要的現實意義。

不一樣類型傳染病的傳播具備不一樣的特色，傳染病的傳播模型不是從醫學角度分析傳染病的傳播過程，而是按照傳播機理創建不一樣的數學模型。

首先，把傳染病流行範圍內的人羣分爲 S、E、I、R 四類，具體含義以下：

• S 類（Susceptible），易感者，指缺少免疫能力的健康人，與感染者接觸後容易受到感染；

• E 類（Exposed），暴露者，指接觸過感染者但暫無傳染性的人，適用於存在潛伏期的傳染病；

• I 類（Infectious），患病者，指具備傳染性的患病者，能夠傳播給 S 類成員將其變爲 E 類或 I 類成員；

• R 類（Recovered），康復者，指病癒後具備免疫力的人。若是免疫期有限，仍能夠從新變爲 S 類成員，進而被感染；若是是終身免疫，則不能再變爲 S類、E類或 I 類成員。

常見的傳染病模型按照傳染病類型分爲 SI、SIR、SIRS、SEIR 模型等，就是由以上四類人羣根據不一樣傳染病的特徵進行組合而產生的不一樣模型。

2. 疫情傳播 SI 模型

2.1 SI 模型的適用範圍

SI 模型適用於只有易感者和患病者兩類人羣，且沒法治癒的疾病，例如 T型病、殭屍。

2.2 SI 模型的假設

1. 考察地區的總人數 N 不變，即不考慮生死或遷移；
2. 人羣分爲易感者（S類）和患病者（I類）兩類；
3. 易感者（S類）與患病者（I類）有效接觸即被感染，變爲患病者，無潛伏期、無治癒狀況、無免疫力；
4. 每一個患病者天天有效接觸的易感者的平均人數（日接觸數）是 \(\lambda\)，稱爲日接觸率；
5. 將第 t 天時 S類、I 類人羣的佔比記爲 \(s(t)\)、\(i(t)\)，數量爲 \(S(t)\)、\(I(t)\)；初始日期 \(t=0\) 時， S類、I 類人羣佔比的初值爲 \(s_0\)、\(i_0\)。

2.3 SI 模型的微分方程

由

\[\begin{align*} N\frac{di}{dt} = N\lambda s i \end{align*} \]

得：

\[\begin{align*} \frac{di}{dt} = \lambda i (1-i)，\ i(0) = i_0 \end{align*} \]

這是 Logistic 模型，用分離變量法能夠求出其解析解爲：

\[\begin{align*} & i(t)=\frac{1}{1+(1/i_0 - 1)\ e^{-\lambda t}}\\ & I(t)= N\ i(t) \end{align*} \]

---

3. SI 模型的 Python 編程

3.1 SI 模型的解析解

上文已經獲得 SI 模型的解析解，對此很容易經過 Python 編程實現，詳見本文例程。

雖然 SI 模型的解析解並不複雜，並且解的精度固然是最好的，但咱們仍然不鼓勵用解析解的方法。緣由在於，一是對於小白求解析解的過程相對複雜困難，並且可能出錯，二是對於更復雜的模型是沒有解析解的，即使大神也只能用數值方法求解。既然如此，不如從一開始就學習、掌握數值求解方法，熟悉數值解法的編程實現。

3.2 SI 模型的數值解

SI 模型是常微分方程初值問題，可使用 Scipy 工具包的 scipy.integrate.odeint() 函數求數值解，具體方法能夠參考前文《Python小白的數學建模課-09 微分方程模型》。

scipy.integrate.odeint(func, y0, t, args=())

**scipy.integrate.odeint() **是求解微分方程的具體方法，經過數值積分來求解常微分方程組。

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。

odeint() 的返回值：

• y: array 　　數組，形狀爲 (len(t),len(y0)，給出時間序列 t 中每一個時刻的 y 值。

odeint() 的編程步驟：

1. 導入 scipy、numpy、matplotlib 包；
2. 定義導數函數 \(f(i,t)=\lambda i (1-i)\) ；
3. 定義初值 \(i_0\) 和 \(i\) 的定義區間 \([t_0,\ t]\)；
4. 調用 odeint() 求 \(i\) 在定義區間 \([t_0,\ t]\) 的數值解。

3.3 Python例程：SI 模型的解析解與數值解

# 1. SI 模型，常微分很是，解析解與數值解的比較
from scipy.integrate import odeint  # 導入 scipy.integrate 模塊
import numpy as np  # 導入 numpy包
import matplotlib.pyplot as plt  # 導入 matplotlib包

def dy_dt(y, t, lamda, mu):  # 定義導數函數 f(y,t)
    dy_dt = lamda*y*(1-y)  # di/dt = lamda*i*(1-i)
    return dy_dt

# 設置模型參數
number = 1e7  # 總人數
lamda = 1.0  # 日接觸率, 患病者天天有效接觸的易感者的平均人數
mu1 = 0.5  # 日治癒率, 天天被治癒的患病者人數佔患病者總數的比例
y0 = i0 = 1e-6  # 患病者比例的初值
tEnd = 50  # 預測日期長度
t = np.arange(0.0,tEnd,1)  # (start,stop,step)

yAnaly = 1/(1+(1/i0-1)*np.exp(-lamda*t))  # 微分方程的解析解
yInteg = odeint(dy_dt, y0, t, args=(lamda,mu1))  # 求解微分方程初值問題
yDeriv = lamda * yInteg *(1-yInteg)

# 繪圖
plt.plot(t, yAnaly, '-ob', label='analytic')
plt.plot(t, yInteg, ':.r', label='numerical')
plt.plot(t, yDeriv, '-g', label='dy_dt')
plt.title("Comparison between analytic and numerical solutions")
plt.legend(loc='right')
plt.axis([0, 50, -0.1, 1.1])
plt.show()

3.4 解析解與數值解的比較

本圖爲例程 2.3 的運行結果，圖中對解析解（藍色）與使用 odeint() 獲得的數值解（紅色）進行比較。在該例中，沒法觀察到解析解與數值解的差別，代表數值解的偏差很小。

圖中 \(di/dt\) 具備最大值，最大值表示疫情增加的高潮，達到最大值後 \(di/dt\) 逐漸減少，但患病者比例很快增加到 100%，代表全部人都被感染成爲患者。

這是特定參數的結果，仍是模型的必然趨勢，須要對參數的影響進行更詳細的研究。

---

4. SI 模型參數的影響

對於 SI 模型，只有日接觸率 \(\lambda\) 和患病者比例的初值 \(i_0\) 會影響模型的結果，其它參數如總人數 N 並無影響。

4.1 日接觸率對 SI 模型的影響

對不一樣日接觸率的比較代表：

1. 日接觸率越大，疫情從發生到爆發的時間越短，爆發過程的增加速度也越快。
2. 不論日接觸率多大，患病者的比例最終都會增加到 1，代表全部人都被感染成爲患者。
3. 不論日接觸率多大，都具備緩慢發展、爆發、增加放緩 3 個階段，進入爆發階段後患病者的比例急劇增加，疫情就很難控制了。

4.2 患病者比例的初值對 SI 模型的影響

對患病者比例初值的比較代表，患病者初值的人數或比例隻影響疫情爆發期到來的快慢，對疫情傳播的過程和結果幾乎沒有影響。

這與咱們直觀的經驗不太一致，一個緣由是 SI 模型自己存在不足，另外一方面也說明若是對傳染病不加控制，即便開始患病人數不多，通過一段時間的傳播後也終將會引發爆發。

4.3 SI 模型結果討論

1. 在 \(i(t)=0.5,\ I(t) = N/2\) 時 $ di/dt$ 達到最大值，病人數目 \(I(t)\) 增長最快。由此能夠預報傳染病高潮的到來，即爲醫院的門診量最大的一天，衛生部門要重點關注。
2. \(t_m\) 與 \(\lambda\) 成反比。日接觸率 \(\lambda\) 反映衛生水平、防控手段，提升衛生水平、強化防控手段，下降病人的日接觸率，能夠推遲傳染病高潮的到來。
3. 當 $ t \to \infty$ 時 \(i \to 1\) ，代表全部人最終都會被傳染而變成病人。這徹底不符合實際狀況，代表該模型太不講 politics 了，只能適用於美帝國家建模。
4. SI 模型很是明顯而嚴重的缺陷，是該模型沒有考慮患病者能夠治癒，所以只能是健康人患病，而患病者不能恢復健康（甚至也不會死亡，而是不斷傳播疫情），因此終將所有被傳染。

【本節完】

版權聲明：

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

原創做品，轉載必須標註原文連接：(https://www.cnblogs.com/youcans/p/14944259.html)。

Copyright 2021 Youcans, XUPT

Crated：2021-06-09

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