傳染病的數學模型是數學建模中的典型問題，常見的傳染病模型有 SI、SIR、SIRS、SEIR 模型，
SI 模型是最簡單的傳染病模型，適用于只有易感者和患病者兩類人群，
我們就從 SI 模型開始吧，從模型、例程、運行結果到模型分析，全都在這個系列中，
『Python小白的數學建模課 @ Youcans』帶你從數模小白成為國賽達人，
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 模型的假設
- 考察地區的總人數 N 不變，即不考慮生死或遷移；
- 人群分為易感者（S類）和患病者（I類）兩類；
- 易感者（S類）與患病者（I類）有效接觸即被感染，變為患病者，無潛伏期、無治愈情況、無免疫力；
- 每個患病者每天有效接觸的易感者的平均人數（日接觸數）是 \(\lambda\)，稱為日接觸率；
- 將第 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() 的編程步驟：
- 匯入 scipy、numpy、matplotlib 包；
- 定義導數函式 \(f(i,t)=\lambda i (1-i)\) ；
- 定義初值 \(i_0\) 和 \(i\) 的定義區間 \([t_0,\ t]\)；
- 呼叫 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，表明所有人都被感染成為患者，
- 不論日接觸率多大，都具有緩慢發展、爆發、增長放緩 3 個階段，進入爆發階段后患病者的比例急劇增長，疫情就很難控制了，
4.2 患病者比例的初值對 SI 模型的影響
對患病者比例初值的比較表明，患病者初值的人數或比例只影響疫情爆發期到來的快慢，對疫情傳播的程序和結果幾乎沒有影響，
這與我們直觀的經驗不太一致，一個原因是 SI 模型本身存在不足，另一方面也說明如果對傳染病不加控制，即使開始患病人數很少，經過一段時間的傳播后也終將會引起爆發，
4.3 SI 模型結果討論
- 在 \(i(t)=0.5,\ I(t) = N/2\) 時 $ di/dt$ 達到最大值，病人數目 \(I(t)\) 增加最快，由此可以預報傳染病高潮的到來，即為醫院的門診量最大的一天，衛生部門要重點關注，
- \(t_m\) 與 \(\lambda\) 成反比，日接觸率 \(\lambda\) 反映衛生水平、防控手段，提高衛生水平、強化防控手段，降低病人的日接觸率，可以推遲傳染病高潮的到來，
- 當 $ t \to \infty$ 時 \(i \to 1\) ，表明所有人最終都會被傳染而變成病人，這完全不符合實際情況，表明該模型太不講 politics 了，只能適用于美帝國家建模，
- 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數模筆記-模擬退火演算法
轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/288609.html
標籤：Python