Python小白的数学建模课-B5. 新冠疫情 SEIR模型
傳染病的數(shù)學(xué)模型是數(shù)學(xué)建模中的典型問題,常見的傳染病模型有 SI、SIR、SIRS、SEIR 模型。
考慮存在易感者、暴露者、患病者和康復(fù)者四類人群,適用于具有潛伏期、治愈后獲得終身免疫的傳染病。
本文詳細(xì)給出了 SEIR 模型微分方程的建模、例程、結(jié)果和分析,讓小白都能懂。
『Python小白的數(shù)學(xué)建模課 @ Youcans』帶你從數(shù)模小白成為國賽達人。
歡迎關(guān)注『Python小白的數(shù)學(xué)建模課 @ Youcans』系列,每周持續(xù)更新
Python小白的數(shù)學(xué)建模課-B2. 新冠疫情 SI模型
Python小白的數(shù)學(xué)建模課-B3. 新冠疫情 SIS模型
Python小白的數(shù)學(xué)建模課-B4. 新冠疫情 SIR模型
Python小白的數(shù)學(xué)建模課-B5. 新冠疫情 SEIR模型
Python小白的數(shù)學(xué)建模課-B6. 新冠疫情 SEIR改進模型
Python小白的數(shù)學(xué)建模課-09.微分方程模型
Python小白的數(shù)學(xué)建模課-10.微分方程邊值問題
1. SEIR 模型
1.1 SEIR 模型的提出
建立傳染病的數(shù)學(xué)模型來描述傳染病的傳播過程,要根據(jù)傳染病的發(fā)病機理和傳播規(guī)律, 結(jié)合疫情數(shù)據(jù)進行擬合分析,可以認(rèn)識傳染病的發(fā)展趨勢,預(yù)測疫情持續(xù)時間和規(guī)模,分析和模擬各種防控措施對疫情發(fā)展的影響程度, 為傳染病防控工作提供決策指導(dǎo),具有重要的理論意義和現(xiàn)實意義。
SI 模型是最簡單的傳染病傳播模型,把人群分為易感者(S 類)和患病者(I 類)兩類,通過 SI 模型可以預(yù)測傳染病高潮的到來;提高衛(wèi)生水平、強化防控手段,降低病人的日接觸率,可以推遲傳染病高潮的到來。在 SI 模型基礎(chǔ)上發(fā)展的 SIS 模型考慮患病者可以治愈而變成易感者,SIS 模型表面?zhèn)魅酒诮佑|數(shù) σ\sigmaσ 是傳染病傳播和防控的關(guān)鍵指標(biāo),決定了疫情終將清零或演變?yōu)榈胤讲¢L期存在。在 SI 模型基礎(chǔ)上考慮病愈免疫的康復(fù)者(R 類)就得到 SIR 模型,通過 SIR 模型也揭示傳染期接觸數(shù) σ\sigmaσ 是傳染病傳播的閾值,滿足 s0>1/σs_0>1/\sigmas0?>1/σ 才會發(fā)生傳染病蔓延,由此可以分析各種防控措施,如:提高衛(wèi)生水平來降低日接觸率λ\lambdaλ、提高醫(yī)療水平來提高日治愈率 μ\muμ,通過預(yù)防接種達到群體免疫來降低 s0s_0s0? 等。
傳染病大多具有潛伏期(incubation period),也叫隱蔽期,是指從被病原體侵入肌體到最早臨床癥狀出現(xiàn)的一段時間。在潛伏期的后期一般具有傳染性。不同的傳染病的潛伏期長短不同,從短至數(shù)小時到長達數(shù)年,但同一種傳染病有固定的(平均)潛伏期。例如,流感的潛伏期為 1~3天,冠狀病毒感染的潛伏期為4~7天,新型冠狀病毒肺炎傳染病(COVID-19)的潛伏期為1-14天(* 來自:新型冠狀病毒肺炎診療方案試行第八版,潛伏時間 1~14天,多為3~7天,在潛伏期具有傳染性),肺結(jié)核的潛伏期從數(shù)周到數(shù)十年。
SEIR 模型考慮存在易感者(Susceptible)、暴露者(Exposed)、患病者(Infectious)和康復(fù)者(Recovered)四類人群,適用于具有潛伏期、治愈后獲得終身免疫的傳染病。易感者(S 類)被感染后成為潛伏者(E類),隨后發(fā)病成為患病者(I 類),治愈后成為康復(fù)者(R類)。這種情況更為復(fù)雜,也更為接近實際情況。
SEIR 模型的倉室結(jié)構(gòu)示意圖如下:
1.2 SEIR 模型假設(shè)
考察地區(qū)的總?cè)藬?shù) N 不變,即不考慮生死或遷移;
人群分為易感者(S 類)、暴露者(E 類)、患病者(I 類)和康復(fù)者(R 類)四類;
易感者(S 類)與患病者(I 類)有效接觸即變?yōu)楸┞墩?#xff08;E 類),暴露者(E 類)經(jīng)過平均潛伏期后成為患病者(I 類);患病者(I 類)可被治愈,治愈后變?yōu)榭祻?fù)者(R 類);康復(fù)者(R類)獲得終身免疫不再易感;
將第 t 天時 S 類、E 類、I 類、R 類人群的占比記為 s(t)s(t)s(t)、e(t)e(t)e(t)、i(t)i(t)i(t)、r(t)r(t)r(t),數(shù)量分別為 S(t)S(t)S(t)、E(t)E(t)E(t)、I(t)I(t)I(t)、R(t)R(t)R(t);初始日期 t=0t=0t=0 時,各類人群占比的初值為 s0s_0s0?、e0e_0e0?、i0i_0i0?、r0r_0r0?;
日接觸數(shù) λ\lambdaλ,每個患病者每天有效接觸的易感者的平均人數(shù);
日發(fā)病率 δ\deltaδ,每天發(fā)病成為患病者的暴露者占暴露者總數(shù)的比例;
日治愈率 μ\muμ,每天被治愈的患病者人數(shù)占患病者總數(shù)的比例,即平均治愈天數(shù)為 1/μ1/\mu1/μ;
傳染期接觸數(shù) σ=λ/μ\sigma = \lambda / \muσ=λ/μ,即每個患病者在整個傳染期內(nèi)有效接觸的易感者人數(shù)。
1.3 SEIR 模型的微分方程
由
{Ndsdt=?NλsiNdedt=Nλsi?NδeNdidt=Nδe?NμiNdrdt=Nμi\begin{cases} N \frac{ds}{dt} = - N \lambda s i\\ N \frac{de}{dt} = N \lambda s i - N \delta e\\ N \frac{di}{dt} = N \delta e - N \mu i\\ N \frac{dr}{dt} = N \mu i\\ \end{cases} ??????????Ndtds?=?NλsiNdtde?=Nλsi?NδeNdtdi?=Nδe?NμiNdtdr?=Nμi?
得:
{dsdt=?λsi,s(0)=s0dedt=λsi?δe,e(0)=e0didt=δe?μi,i(0)=i0\begin{cases} \frac{ds}{dt} = -\lambda s i, &s(0)=s_0\\ \frac{de}{dt} = \lambda s i - \delta e, &e(0)=e_0\\ \frac{di}{dt} = \delta e - \mu i, &i(0)=i_0 \end{cases} ??????dtds?=?λsi,dtde?=λsi?δe,dtdi?=δe?μi,?s(0)=s0?e(0)=e0?i(0)=i0??
SEIR 模型不能求出解析解,可以通過數(shù)值計算方法求解。
2. SEIR 模型的 Python 編程
2.1 Scipy 工具包求解微分方程組
SIS 模型是常微分方程初值問題,可以使用 Scipy 工具包的 scipy.integrate.odeint() 函數(shù)求數(shù)值解。
scipy.integrate.odeint(func, y0, t, args=())**scipy.integrate.odeint() **是求解微分方程的具體方法,通過數(shù)值積分來求解常微分方程組。
odeint() 的主要參數(shù):
- func: callable(y, t, …) 導(dǎo)數(shù)函數(shù) f(y,t)f(y,t)f(y,t) ,即 y 在 t 處的導(dǎo)數(shù),以函數(shù)的形式表示
- y0: array: 初始條件 y0y_0y0?,注意 SEIR模型是二元常微分方程組, 初始條件為數(shù)組向量 y0=[i0,s0]y_0=[i_0, s_0]y0?=[i0?,s0?]
- t: array: 求解函數(shù)值對應(yīng)的時間點的序列。序列的第一個元素是與初始條件 y0y_0y0? 對應(yīng)的初始時間 t0t_0t0?;時間序列必須是單調(diào)遞增或單調(diào)遞減的,允許重復(fù)值。
- args: 向?qū)?shù)函數(shù) func 傳遞參數(shù)。當(dāng)導(dǎo)數(shù)函數(shù) f(y,t,p1,p2,..)f(y,t,p1,p2,..)f(y,t,p1,p2,..) 包括可變參數(shù) p1,p2… 時,通過 args =(p1,p2,…) 可以將參數(shù)p1,p2… 傳遞給導(dǎo)數(shù)函數(shù) func。
odeint() 的返回值:
- y: array 數(shù)組,形狀為 (len(t),len(y0),給出時間序列 t 中每個時刻的 y 值。
2.2 odeint() 求解 SEIR 模型的編程步驟
常微分方程的導(dǎo)數(shù)定義(SIS模型)
def dySIS(y, t, lamda, mu): # SIS 模型,導(dǎo)數(shù)函數(shù)dy_dt = lamda*y*(1-y) - mu*y # di/dt = lamda*i*(1-i)-mu*ireturn dy_dt常微分方程組的導(dǎo)數(shù)定義(SEIR模型)
def dySEIR(y, t, lamda, delta, mu): # SEIR 模型,導(dǎo)數(shù)函數(shù)s, e, i = yds_dt = - lamda*s*i # ds/dt = -lamda*s*ide_dt = lamda*s*i - delta*e # de/dt = lamda*s*i - delta*edi_dt = delta*e - mu*i # di/dt = delta*e - mu*ireturn np.array([ds_dt,de_dt,di_dt])Python 可以直接對向量、向量函數(shù)進行定義和賦值,使程序更為簡潔。但考慮讀者主要是 Python 小白,又涉及到看著就心煩的微分方程組,所以我們寧愿把程序?qū)懙美圪樢恍?#xff0c;便于讀者將程序與前面的微分方程組逐項對應(yīng)。
SEIR 模型是三元常微分方程,返回值 y 是 len(t)*3 的二維數(shù)組。
2.3 Python例程:SEIR 模型
# modelCovid4_v1.py # Demo01 of mathematical modeling for Covid2019 # SIR model for epidemic diseases. # Copyright 2021 Youcans, XUPT # Crated:2021-06-13 # Python小白的數(shù)學(xué)建模課 @ Youcans# 1. SEIR 模型,常微分方程組 from scipy.integrate import odeint # 導(dǎo)入 scipy.integrate 模塊 import numpy as np # 導(dǎo)入 numpy包 import matplotlib.pyplot as plt # 導(dǎo)入 matplotlib包def dySIS(y, t, lamda, mu): # SI/SIS 模型,導(dǎo)數(shù)函數(shù)dy_dt = lamda*y*(1-y) - mu*y # di/dt = lamda*i*(1-i)-mu*ireturn dy_dtdef dySIR(y, t, lamda, mu): # SIR 模型,導(dǎo)數(shù)函數(shù)s, i = y # youcansds_dt = -lamda*s*i # ds/dt = -lamda*s*idi_dt = lamda*s*i - mu*i # di/dt = lamda*s*i-mu*ireturn np.array([ds_dt,di_dt])def dySEIR(y, t, lamda, delta, mu): # SEIR 模型,導(dǎo)數(shù)函數(shù)s, e, i = y # youcansds_dt = -lamda*s*i # ds/dt = -lamda*s*ide_dt = lamda*s*i - delta*e # de/dt = lamda*s*i - delta*edi_dt = delta*e - mu*i # di/dt = delta*e - mu*ireturn np.array([ds_dt,de_dt,di_dt])# 設(shè)置模型參數(shù) number = 1e5 # 總?cè)藬?shù) lamda = 0.3 # 日接觸率, 患病者每天有效接觸的易感者的平均人數(shù) delta = 0.03 # 日發(fā)病率,每天發(fā)病成為患病者的潛伏者占潛伏者總數(shù)的比例 mu = 0.06 # 日治愈率, 每天治愈的患病者人數(shù)占患病者總數(shù)的比例 sigma = lamda / mu # 傳染期接觸數(shù) fsig = 1-1/sigma tEnd = 300 # 預(yù)測日期長度 t = np.arange(0.0,tEnd,1) # (start,stop,step) i0 = 1e-3 # 患病者比例的初值 e0 = 1e-3 # 潛伏者比例的初值 s0 = 1-i0 # 易感者比例的初值 Y0 = (s0, e0, i0) # 微分方程組的初值# odeint 數(shù)值解,求解微分方程初值問題 ySI = odeint(dySIS, i0, t, args=(lamda,0)) # SI 模型 ySIS = odeint(dySIS, i0, t, args=(lamda,mu)) # SIS 模型 ySIR = odeint(dySIR, (s0,i0), t, args=(lamda,mu)) # SIR 模型 ySEIR = odeint(dySEIR, Y0, t, args=(lamda,delta,mu)) # SEIR 模型# 輸出繪圖 print("lamda={}\tmu={}\tsigma={}\t(1-1/sig)={}".format(lamda,mu,sigma,fsig)) plt.title("Comparison among SI, SIS, SIR and SEIR models") plt.xlabel('t-youcans') plt.axis([0, tEnd, -0.1, 1.1]) plt.plot(t, ySI, 'cadetblue', label='i(t)-SI') plt.plot(t, ySIS, 'steelblue', label='i(t)-SIS') plt.plot(t, ySIR[:,1], 'cornflowerblue', label='i(t)-SIR') # plt.plot(t, 1-ySIR[:,0]-ySIR[:,1], 'cornflowerblue', label='r(t)-SIR') plt.plot(t, ySEIR[:,0], '--', color='darkviolet', label='s(t)-SIR') plt.plot(t, ySEIR[:,1], '-.', color='orchid', label='e(t)-SIR') plt.plot(t, ySEIR[:,2], '-', color='m', label='i(t)-SIR') plt.plot(t, 1-ySEIR[:,0]-ySEIR[:,1]-ySEIR[:,2], ':', color='palevioletred', label='r(t)-SIR') plt.legend(loc='right') # youcans plt.show()2.4 SI /SIS/SIR 模型與SEIR 模型的比較
例程 2.3 的參數(shù)和初值為:λ=0.3,δ=0.03,μ=0.06,(s0,e0,i0)=(0.001,0.001,0.998)\lambda=0.3,\delta=0.03,\mu=0.06,(s_0,e_0,i_0)=(0.001,0.001,0.998)λ=0.3,δ=0.03,μ=0.06,(s0?,e0?,i0?)=(0.001,0.001,0.998),上圖為例程的運行結(jié)果。
曲線 i(t)-SI 是 SI 模型的結(jié)果,患病者比例急劇增長到 1.0,所有人都被傳染而變成患病者。
曲線 i(t)-SIS 是 SIS 模型的結(jié)果,患病者比例快速增長并收斂到某個常數(shù),即穩(wěn)態(tài)特征值 i∞=1?μ/λ=0.8i_\infty=1-\mu/\lambda = 0.8i∞?=1?μ/λ=0.8,表明疫情穩(wěn)定,并將長期保持一定的患病率,稱為地方病平衡點。
曲線 i(t)-SIR 是 SIR 模型的結(jié)果,患病者比例 i(t) 先上升達到峰值,然后再逐漸減小趨近于常數(shù)。
曲線 s(t)-SEIR、e(t)-SEIR、i(t)-SEIR、r(t)-SEIR 分別表示 SEIR模型中易感者(S類)、潛伏者(E類)、患病者(I類)和康復(fù)者(R 類)人群的占比。
圖中易感者比例 s(t) 單調(diào)遞減并收斂到接近于 0 的穩(wěn)定值。潛伏者比例 e(t) 曲線存在波峰,先逐漸上升而達到峰值,然后再逐漸減小,最終趨于 0?;疾≌弑壤?i(t) 曲線與潛伏者比例曲線類似,上升達到峰值后逐漸減小,最終趨于 0;但患病者比例曲線發(fā)展、達峰的時間比潛伏者曲線要晚一些,峰值強度也較低??祻?fù)者比例 r(i) 單調(diào)遞增并收斂到非零的穩(wěn)態(tài)值。以上分析只是對本圖進行的討論,并非普遍結(jié)論,取決于具體參數(shù)條件。
比較相同參數(shù)條件下 SIR 和 SEIR 模型的結(jié)果,SIR 模型中患病者比例 i(t) 的波形起點、峰值和終點到來的時間都顯著早于 SEIR 模型,峰值強度也高于 SEIR 模型。這表明具有潛伏期的傳染病,疫情發(fā)生和峰值的到來要晚于沒有潛伏期的傳染病,而且持續(xù)時間更長。
3. SEIR 模型參數(shù)的影響
SEIR 模型中有日接觸率 λ\lambdaλ 、日發(fā)病率 δ\deltaδ 和日治愈率 μ\muμ 三個參數(shù),還有 i0、e0、s0i_0、e_0、s_0i0?、e0?、s0? 等初始條件,我們先用單因素分析的方法來觀察參數(shù)條件對于疫情傳播的影響。
3.1 初值條件 i0、e0、s0i_0、e_0、s_0i0?、e0?、s0? 初始條件的影響
SEIR 模型中有 i0、e0、s0i_0、e_0、s_0i0?、e0?、s0? 等 3個初始條件,組合眾多無法窮盡。考慮實際情況中,疫情初始階段尚無康復(fù)者,而潛伏者比例往往高于確診的發(fā)病者,我們假定 e0/i0=2,r0=0e_0/i_0=2,r_0=0e0?/i0?=2,r0?=0,考察不同 i0i_0i0? 時的疫情傳播情況。
通過對于該參數(shù)下不同的患病者、潛伏者比例初值條件的模擬,可以看到患病者、潛伏者比例的初值條件對疫情發(fā)生、達峰、結(jié)束的時間早晚具有直接影響,但對疫情曲線的形態(tài)和特征影響不大。不同初值條件下的疫情曲線,幾乎是沿著時間指標(biāo)平移的。
這說明如果不進行治療防控等人為干預(yù),疫情傳播過程與初始患病者、潛伏者比例關(guān)系并不大,該來的總會來。
圖中患病率達到高峰后逐步降低,直至趨近于 0;易感率在疫情爆發(fā)后迅速下降,直至趨近于 0。但這一現(xiàn)象是基于具體的參數(shù)條件的觀察,僅由該圖并不能確定其是否普遍規(guī)律。
3.2 日接觸率 λ\lambdaλ 的影響
首先考察日接觸率 λ\lambdaλ 的影響。
保持參數(shù) δ=0.1,μ=0.06,(i0=0.001,e0=0.002,s0=0.997)\delta =0.1,\mu=0.06, (i_0=0.001,e_0=0.002,s_0=0.997)δ=0.1,μ=0.06,(i0?=0.001,e0?=0.002,s0?=0.997) 不變,$\lambda = [0.12, 0.25, 0.5, 1.0, 2.0] $ 時 i(t),s(t)i(t), s(t)i(t),s(t) 的變化曲線如下圖所示。
通過對于該條件下日接觸率的單因素分析,可以看到隨著日接觸率 λ\lambdaλ 的增大,患病率比例 i(t)i(t)i(t) 出現(xiàn)的峰值更早、更強,而易感者比例 s(t)s(t)s(t) 逐漸降低,但最終都趨于穩(wěn)定。
3.3 日發(fā)病率 δ\deltaδ 的影響
下面考察日發(fā)病率 δ\deltaδ 的影響。保持參數(shù) λ=0.25,μ=0.06,(i0=0.001,e0=0.002,s0=0.997)\lambda =0.25,\mu=0.06, (i_0=0.001,e_0=0.002,s_0=0.997)λ=0.25,μ=0.06,(i0?=0.001,e0?=0.002,s0?=0.997) 不變,$\delta = [0.02, 0.05, 0.1, 0.2, 0.4] $ 時 i(t),s(t)i(t), s(t)i(t),s(t) 的變化曲線如下圖所示。
通過對于該條件下日接觸率的單因素分析,可以看到隨著日接觸率 λ\lambdaλ 的增大,患病率比例 i(t)i(t)i(t) 出現(xiàn)的峰值更早、更強,而易感者比例 s(t)s(t)s(t) 逐漸降低,但最終都趨于穩(wěn)定。
3.4 日治愈率 μ\muμ 的影響
下面考察日治愈率 μ\muμ 的影響。保持參數(shù) λ=0.25,δ=0.1,(i0=0.001,e0=0.002,s0=0.997)\lambda =0.25,\delta=0.1, (i_0=0.001,e_0=0.002,s_0=0.997)λ=0.25,δ=0.1,(i0?=0.001,e0?=0.002,s0?=0.997) 不變,$\mu = [0.02, 0.05, 0.1, 0.2, 0.4] $ 時 i(t),s(t)i(t), s(t)i(t),s(t) 的變化曲線如下圖所示。
通過對該條件下日治愈率的單因素分析,可以看到在日治愈率 μ=0.4\mu=0.4μ=0.4 時,患病者比例始終非常低并趨于 0,易感者比例幾乎不變,表明疫情不會傳播,這是因為患病者治愈的速度快于易感者被感染的速度。
在日治愈率 μ=0.2\mu=0.2μ=0.2 時,患病者比例也始終非常低(接近 0),易感者比例緩慢降低并趨于穩(wěn)定值,表明疫情的傳播十分緩慢微弱,這是因為患病者治愈的速度較快,在易感者比例逐漸降低到某一水平后治愈人數(shù)與感染人數(shù)將達到平衡。
在日治愈率較低時 (μ<0.2\mu<0.2μ<0.2 ),患病者比例曲線存在波峰,然后再逐漸減小,最終趨于 0。隨著日治愈率 μ\muμ 的減小,患病率比例 i(t)i(t)i(t) 曲線的峰值更強、也稍早。易感者比例 s(t)s(t)s(t) 隨著患病者比例上升而迅速降低,最終趨于某個穩(wěn)定值,也達到治愈與感染的平衡。
通過對不同參數(shù)的單因素實驗和結(jié)果分析,可以發(fā)現(xiàn)雖然各個模型參數(shù)和初始條件都會影響疫情曲線,但日治愈率 與日接觸率的關(guān)系更為重要。進一步的分析和模擬可以發(fā)現(xiàn),傳染期接觸數(shù) σ=λ/μ\sigma = \lambda / \muσ=λ/μ是傳染病蔓延的閾值,滿足 s0>1/σs_0>1/\sigmas0?>1/σ 才會發(fā)生傳染病蔓延。
這說明具有決定性影響的特征參數(shù),往往不是模型中的某個參數(shù),而是多個參數(shù)特定關(guān)系的組合,僅從單因素實驗很難充分反映模型中的內(nèi)在特征。
4. SEIR 模型的相空間分析
4.1 Python例程:SEIR 模型的相軌跡
# modelCovid4_v1.py # Demo01 of mathematical modeling for Covid2019 # SIR model for epidemic diseases. # Copyright 2021 Youcans, XUPT # Crated:2021-06-15 # Python小白的數(shù)學(xué)建模課 @ Youcans# 7. SEIR 模型,常微分方程組 相空間分析: e(t)~i(t) from scipy.integrate import odeint # 導(dǎo)入 scipy.integrate 模塊 import numpy as np # 導(dǎo)入 numpy包 import matplotlib.pyplot as plt # 導(dǎo)入 matplotlib包def dySEIR(y, t, lamda, delta, mu): # SEIR 模型,導(dǎo)數(shù)函數(shù)s, e, i = yds_dt = -lamda*s*i # ds/dt = -lamda*s*ide_dt = lamda*s*i - delta*e # de/dt = lamda*s*i - delta*edi_dt = delta*e - mu*i # di/dt = delta*e - mu*ireturn np.array([ds_dt,de_dt,di_dt])# 設(shè)置模型參數(shù) number = 1e5 # 總?cè)藬?shù) lamda = 0.3 # 日接觸率, 患病者每天有效接觸的易感者的平均人數(shù) delta = 0.1 # 日發(fā)病率,每天發(fā)病成為患病者的潛伏者占潛伏者總數(shù)的比例 mu = 0.1 # 日治愈率, 每天治愈的患病者人數(shù)占患病者總數(shù)的比例 sigma = lamda / mu # 傳染期接觸數(shù) tEnd = 500 # 預(yù)測日期長度 t = np.arange(0.0, tEnd, 1) # (start,stop,step)# e0List = np.arange(0.01,0.4,0.05) # (start,stop,step)e0List = np.arange(0.01,0.4,0.05) # (start,stop,step) for e0 in e0List:# odeint 數(shù)值解,求解微分方程初值問題i0 = 0 # 潛伏者比例的初值s0 = 1 - i0 - e0 # 易感者比例的初值ySEIR = odeint(dySEIR, (s0,e0,i0), t, args=(lamda,delta,mu)) # SEIR 模型plt.plot(ySEIR[:,1], ySEIR[:,2]) # (e(t),i(t))print("lamda={}\tdelta={}\mu={}\tsigma={}\ti0={}\te0={}".format(lamda,delta,mu,lamda/mu,i0,e0))# 輸出繪圖 plt.title("Phase trajectory of SEIR models: e(t)~i(t)") plt.axis([0, 0.4, 0, 0.4]) plt.plot([0,0.4],[0,0.35],'y--') #[x1,x2][y1,y2] plt.plot([0,0.4],[0,0.18],'y--') #[x1,x2][y1,y2] plt.text(0.02,0.36,r"$\lambda=0.3, \delta=0.1, \mu=0.1$",color='black') plt.xlabel('e(t)') plt.ylabel('i(t)') plt.show()4.2 SEIR 模型的相軌跡分析
上圖為例程 4.2 的運行結(jié)果,是 SEIR 模型的相軌跡圖。
圖中每一條 e-s 曲線,從直線 i(t)+s(t)=1 上的某一初值點出發(fā),最終收斂于 s軸上的某一點,對應(yīng)著某一個初值條件下的患病者與易感者比例隨時間的變化關(guān)系。
SEIR 模型的相軌跡圖比較復(fù)雜,難以在本文展開進行討論。但是,相軌跡圖中的虛線還是反映出了時間變化曲線圖中難以表達的一些重要特征,以此為線索可以進行更深入的研究。
5. SEIR 模型結(jié)果討論
最后,我們簡單總結(jié)一下 SEIR 模型的特點:
SEIR 模型是一個單向模型,易感者(S)不斷轉(zhuǎn)變?yōu)闈摲?#xff08;E),潛伏者(E)發(fā)病后成為患病者(I),患病者(I)不斷轉(zhuǎn)變?yōu)榭祻?fù)者(R),因此易感者比例 s(t) 單調(diào)遞減,康復(fù)者比例 r(i) 單調(diào)遞增。
SEIR 模型假設(shè)潛伏期無傳染性,因此潛伏期延遲了傳染期的開始,疫情發(fā)生和峰值的到來要晚于沒有潛伏期的 SIR 模型。但潛伏期不會改變感染人群的累計數(shù)量,而且持續(xù)時間更長。
1/σ1/\sigma1/σ 是傳染病蔓延的閾值,滿足 s0>1/σs_0>1/\sigmas0?>1/σ 才會發(fā)生傳染病蔓延。因此,為了控制傳染病的蔓延:一方面要提高閾值 1/σ1/\sigma1/σ,這可以通過提高衛(wèi)生水平來降低日接觸率λ\lambdaλ、提高醫(yī)療水平來提高日治愈率 μ\muμ;另一方面要降低 s0s_0s0?,這可以通過預(yù)防接種達到群體免疫來實現(xiàn)。
在 SEIR 模型的基礎(chǔ)上,可以根據(jù)不同傳染病病理特征及疫情傳播特點,對模型進行進一步的改進,使模型與實際情況更加吻合,以便更準(zhǔn)確地預(yù)測疫情發(fā)展趨勢。
在 SEIR 模型的基礎(chǔ)上,可以結(jié)合實際的疫情數(shù)據(jù)來擬合和估計模型參數(shù),進而用來模擬和分析不同治療方案和防控措施對疫情發(fā)展的影響,為新冠疫情的防控工作提供決策指導(dǎo)。
【本節(jié)完】
版權(quán)聲明:
歡迎關(guān)注『Python小白的數(shù)學(xué)建模課 @ Youcans』 原創(chuàng)作品
原創(chuàng)作品,轉(zhuǎn)載必須標(biāo)注原文鏈接:(https://blog.csdn.net/youcans/article/details/117932162)。
Copyright 2021 Youcans, XUPT
Crated:2021-06-15
歡迎關(guān)注 『Python小白的數(shù)學(xué)建模課 @ Youcans』 系列,持續(xù)更新
Python小白的數(shù)學(xué)建模課-01.新手必讀
Python小白的數(shù)學(xué)建模課-02.數(shù)據(jù)導(dǎo)入
Python小白的數(shù)學(xué)建模課-03.線性規(guī)劃
Python小白的數(shù)學(xué)建模課-04.整數(shù)規(guī)劃
Python小白的數(shù)學(xué)建模課-05.0-1規(guī)劃
Python小白的數(shù)學(xué)建模課-06.固定費用問題
Python小白的數(shù)學(xué)建模課-07.選址問題
Python小白的數(shù)學(xué)建模課-09.微分方程模型
Python小白的數(shù)學(xué)建模課-10.微分方程邊值問題
Python小白的數(shù)學(xué)建模課-A1.國賽賽題類型分析
Python小白的數(shù)學(xué)建模課-A2.2021年數(shù)維杯C題探討
Python小白的數(shù)學(xué)建模課-A3.12個新冠疫情數(shù)模競賽賽題及短評
Python小白的數(shù)學(xué)建模課-B2. 新冠疫情 SI模型
Python小白的數(shù)學(xué)建模課-B3. 新冠疫情 SIS模型
Python小白的數(shù)學(xué)建模課-B4. 新冠疫情 SIR模型
Python小白的數(shù)學(xué)建模課-B5. 新冠疫情 SEIR模型
Python小白的數(shù)學(xué)建模課-B6. 新冠疫情 SEIR改進模型
Python數(shù)模筆記-PuLP庫
Python數(shù)模筆記-StatsModels統(tǒng)計回歸
Python數(shù)模筆記-Sklearn
Python數(shù)模筆記-NetworkX
Python數(shù)模筆記-模擬退火算法
總結(jié)
以上是生活随笔為你收集整理的Python小白的数学建模课-B5. 新冠疫情 SEIR模型的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: mysql如何查看表拥有的键_如何查看表
- 下一篇: mac python运行按哪个键_#ma