传染病的数学模型是数学建模中的典型问题,常见的传染病模型有 SI、SIR、SIRS、SEIR 模型。
SIR 模型将人群分为易感者(S类)、患病者(I类)和康复者(R 类),考虑了患病者治愈后的免疫能力。
本文详细给出了 SIR 模型微分方程、相空间分析的建模、例程、结果和分析,让小白都能懂。
『Python小白的数学建模课 @ Youcans』带你从数模小白成为国赛达人。
传染病的传播特性不可能通过真实的试验开展研究,因此需要针对不同的传染病传播方式和流行特点建立相应的数学模型,并对模型进行理论研究和数值模拟。通过研究发现传染病传播的特征阈值,就可以为预防和控制传染病提供数据支撑和防控策略。
1927年,W. Kermack 在论文 “Contributions to the mathcmatical theory of epidemics” 中研究了伦敦黑死病和孟买瘟疫的流行过程,创造性地提出了 SIR 模型。
SI 模型和 SIS 模型将人群分为感染者(S类)和患病者(I类)两类人群,但大多数传染病的患者在治愈后就有很强的免疫力,终身或在一段时期内不再会被感染而变成病人。这类人群称为病愈免疫的康复者(R 类)。康复者已经退出传染系统,对于致死性疾病的死亡者也可以用该类别描述其传播特性。
SIR 模型适用于具有易感者、患病者和康复者三类人群,可以治愈,且治愈后终身免疫不再复发的疾病,例如天花、肝炎、麻疹等免疫力很强的传染病。
SIR 模型假设:
SIR 模型的微分方程:
由
{
N
d
i
d
t
=
N
λ
s
i
N
d
r
d
t
=
N
μ
i
s
(
t
)
+
i
(
t
)
+
r
(
t
)
=
1
\begin{cases} & N\frac{di}{dt} = N\lambda s i\\ & N \frac{dr}{dt} = N \mu i\\ & s(t) + i(t) + r(t) = 1 \end{cases}
⎩⎪⎨⎪⎧Ndtdi=NλsiNdtdr=Nμis(t)+i(t)+r(t)=1
得:
{
d
i
d
t
=
λ
s
i
−
μ
i
,
i
(
0
)
=
i
0
d
s
d
t
=
−
λ
s
i
,
s
(
0
)
=
s
0
\begin{cases} & \frac{di}{dt} = \lambda s i - \mu i , &i(0)=i_0\\ & \frac{ds}{dt} = -\lambda s i , &s(0)=s_0 \end{cases}
{dtdi=λsi−μi,dtds=−λsi,i(0)=i0s(0)=s0
SIR 模型不能求出解析解,只能通过数值计算方法求解。
SIS 模型是常微分方程初值问题,可以使用 Scipy 工具包的 scipy.integrate.odeint() 函数求数值解。
scipy.integrate.odeint(func, y0, t, args=())
**scipy.integrate.odeint() **是求解微分方程的具体方法,通过数值积分来求解常微分方程组。
odeint() 的主要参数:
odeint() 的返回值:
常微分方程的导数定义(SIS模型)
def dySIS(y, t, lamda, mu): # SIS 模型,导数函数 dy_dt = lamda*y*(1-y) - mu*y # di/dt = lamda*i*(1-i)-mu*i return dy_dt
常微分方程组的导数定义(SIR模型)
def dySIR(y, t, lamda, mu): # SIR 模型,Y=[i,s] 点的导数dy/dt i, s = y di_dt = lamda*s*i - mu*i # di/dt = lamda*s*i-mu*i ds_dt = -lamda*s*i # ds/dt = -lamda*s*i return np.array([di_dt,ds_dt])
Python 可以直接对向量、向量函数进行定义和赋值,使程序更为简洁。但考虑读者主要是 Python 小白,又涉及到看着就心烦的微分方程组,所以我们宁愿把程序写得累赘一些,便于读者将程序与前面的微分方程组逐项对应。
SIR 模型是二元常微分方程,返回值 y 是 len(t)*2 的二维数组。
# modelCovid3_v1.py # Demo01 of mathematical modeling for Covid2019 # SIR model for epidemic diseases. # Copyright 2021 Youcans, XUPT # Crated:2021-06-12 # Python小白的数学建模课 @ Youcans # 1. SIR 模型,常微分方程组 from scipy.integrate import odeint # 导入 scipy.integrate 模块 import numpy as np # 导入 numpy包 import matplotlib.pyplot as plt # 导入 matplotlib包 def dySIS(y, t, lamda, mu): # SI/SIS 模型,导数函数 dy_dt = lamda*y*(1-y) - mu*y # di/dt = lamda*i*(1-i)-mu*i return dy_dt def dySIR(y, t, lamda, mu): # SIR 模型,导数函数 i, s = y di_dt = lamda*s*i - mu*i # di/dt = lamda*s*i-mu*i ds_dt = -lamda*s*i # ds/dt = -lamda*s*i return np.array([di_dt,ds_dt]) # 设置模型参数 number = 1e5 # 总人数 lamda = 0.2 # 日接触率, 患病者每天有效接触的易感者的平均人数 sigma = 2.5 # 传染期接触数 mu = lamda/sigma # 日治愈率, 每天被治愈的患病者人数占患病者总数的比例 fsig = 1-1/sigma tEnd = 200 # 预测日期长度 t = np.arange(0.0,tEnd,1) # (start,stop,step) i0 = 1e-4 # 患病者比例的初值 s0 = 1-i0 # 易感者比例的初值 Y0 = (i0, s0) # 微分方程组的初值 print("lamda={}\tmu={}\tsigma={}\t(1-1/sig)={}".format(lamda,mu,sigma,fsig)) # odeint 数值解,求解微分方程初值问题 ySI = odeint(dySIS, i0, t, args=(lamda,0)) # SI 模型 ySIS = odeint(dySIS, i0, t, args=(lamda,mu)) # SIS 模型 ySIR = odeint(dySIR, Y0, t, args=(lamda,mu)) # SIR 模型 # 绘图 plt.title("Comparison among SI, SIS and SIR models") plt.xlabel('t-youcans') plt.axis([0, tEnd, -0.1, 1.1]) plt.axhline(y=0,ls="--",c='c') # 添加水平直线 plt.plot(t, ySI, ':g', label='i(t)-SI') plt.plot(t, ySIS, '--g', label='i(t)-SIS') plt.plot(t, ySIR[:,0], '-r', label='i(t)-SIR') plt.plot(t, ySIR[:,1], '-b', label='s(t)-SIR') plt.plot(t, 1-ySIR[:,0]-ySIR[:,1], '-m', label='r(t)-SIR') plt.legend(loc='best') # youcans plt.show()
本图为例程 2.3 的运行结果,参数和初值为: λ = 0.2 , μ = 0.08 , ( i 0 , s 0 , r 0 ) = ( 0.0001 , 0.9999 , 0 ) \lambda =0.2,\mu=0.08,(i_0,s_0,r_0)=(0.0001,0.9999,0) λ=0.2,μ=0.08,(i0,s0,r0)=(0.0001,0.9999,0)。
曲线 i(t)-SI 是 SI 模型的结果,患病者比例急剧增长到 1.0,所有人都被传染而变成患病者。
曲线 i(t)-SIS 是 SIS 模型的结果,患病者比例快速增长并收敛到某个常数,即稳态特征值 i ∞ = 1 − μ / λ = 0.6 i_\infty=1-\mu/\lambda = 0.6 i∞=1−μ/λ=0.6,表明疫情稳定,并将长期保持一定的患病率,称为地方病平衡点。
曲线 i(t)-SIS、s(t)-SIS、r(t)-SIS 分别是 SIS 模型的易感者(S类)、患病者(I 类)、康复者(R 类)人群的占比。图中易感者比例 s(t) 单调递减并收敛到非零的稳态值 s ∞ s_\infty s∞,康复者比例 r(i) 单调递增并收敛到非零的稳态值 r ∞ r_\infty r∞,患病者比例 i(t) 先上升达到峰值,然后再逐渐减小趋近于 0 。
SIR 模型中有日接触率 λ \lambda λ 与日治愈率 μ \mu μ 两个参数,还有 i 0 、 s 0 i_0、s_0 i0、s0 两个初始条件,共有 4 个可以调整的参数条件都会影响微分方程的解,也就是会影响患病者、易感者比例的时间变化曲线。其中的各种组合无穷无尽,如果没有恰当的研究方法、不能把握内在的规律,即使在几十、几百组参数条件下进行模拟,仍然只是盲人摸象、管中窥豹。
下面考察初值条件
i
0
,
s
0
i_0, s_0
i0,s0 的影响。固定参数
λ
=
0.2
,
μ
=
0.02
\lambda=0.2, \mu=0.02
λ=0.2,μ=0.02不变,不同初值条件下
i
(
t
)
,
s
(
t
)
i(t), s(t)
i(t),s(t) 的变化曲线如下图所示。
通过对于该参数条件下不同初值条件的单因素分析,可以看到患病者比例、易感者比例的初值条件对疫情发生、达峰、结束的时间早晚具有直接影响,但对疫情曲线的形态和特征影响不大。不同初值条件下的疫情曲线,几乎是沿着时间指标平移的。
这说明如果不进行治疗防控等人为干预,疫情传播过程与初始患病率无关,该来的总会来。
图中患病率达到高峰后逐步降低,直至趋近于 0;易感率在疫情爆发后迅速下降,直至趋近于 0。但这一现象是基于具体的参数条件 λ = 0.2 , μ = 0.02 \lambda=0.2, \mu=0.02 λ=0.2,μ=0.02 的观察,无法确定其是否普遍规律。
首先考察日接触率
λ
\lambda
λ 的影响。固定参数
μ
=
0.25
,
i
0
=
0.002
,
s
0
=
1
−
i
0
\mu=0.25, i_0=0.002,s_0=1-i_0
μ=0.25,i0=0.002,s0=1−i0不变,$\lambda = [0.2, 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) 出现的峰值更早、更强,而易感者比例 s ( t ) s(t) s(t) 从几乎不变到迅速降低,但最终都趋于稳定。
对本图我们好像感觉到存在一些规律,但又似乎说不清,很难总结出来。即便总结出某些特征,也只能说是在该固定参数条件下的特征,不能说是 SIR 模型的共有特征。即便反复地改变固定参数的取值或日接触率的范围,情况也差不多。
下面考察日治愈率
μ
\mu
μ 的影响。固定参数
λ
=
0.2
,
i
0
=
0.002
,
s
0
=
1
−
i
0
\lambda=0.2, i_0=0.002,s_0=1-i_0
λ=0.2,i0=0.002,s0=1−i0不变,
μ
=
[
0.4
,
0.2
,
0.1
,
0.05
,
0.025
]
\mu = [0.4, 0.2, 0.1, 0.05, 0.025]
μ=[0.4,0.2,0.1,0.05,0.025] 时
i
(
t
)
,
s
(
t
)
i(t), s(t)
i(t),s(t) 的变化曲线如下图所示。
通过对于该条件下日治愈率的单因素分析,可以看到随着日治愈率 μ \mu μ 的减小,患病率比例 i ( t ) i(t) i(t) 出现的峰值更强、也稍早,而易感者比例 s ( t ) s(t) s(t) 从几乎不变到迅速降低,但最终都趋于稳定。
对于本图的观察和分析情况与上图是差不多的,看起来内容更丰富,似乎也有规律可循,但很难说的清,只能做一些简单的描述。即便进行更多的模拟,情况也差不多。
这是因为,对于SIR 模型这类微分方程,对结果具有决定性影响的特征参数,往往不是模型中的某个参数,而是多个参数特定关系的组合,因此仅从单因素实验很难充分反映模型中的内在特征。
不过,对于数学建模,通过几组单因素实验获得一系列曲线、图表,再从各个角度对结果进行一些描述和解读,就已经足够了。
SIR 模型不能求出解析解,可以通过相空间方法来研究解的周期性、稳定性。
由于患病者比例 i ( t ) i(t) i(t) 和易感者比例 s ( t ) s(t) s(t) 都是时间 t 的函数,因此当 t 取任意值时都对应着 i − s i-s i−s 平面上的一个点,当 t 连续变化时对应着 i − s i-s i−s 平面上的一条轨迹,称为相轨迹。通过相轨迹图可以分析微分方程的性质。
对于 SIR 模型,消去 dt 可以得到:
d
i
d
s
=
1
σ
s
−
1
,
i
(
s
=
s
0
)
=
i
0
\frac{di}{ds} = \frac{1}{\sigma s} - 1 ,\; i(s=s_0)=i_0
dsdi=σs1−1,i(s=s0)=i0
该微分方程的解为:
i
=
(
s
0
+
i
0
)
−
s
+
1
σ
l
n
s
s
0
i = (s_0 + i_0) - s +\frac{1}{\sigma} ln \frac{s}{s_0}
i=(s0+i0)−s+σ1lns0s
# modelCovid3_v1.py # Demo01 of mathematical modeling for Covid2019 # SIR model for epidemic diseases. # Copyright 2021 Youcans, XUPT # Crated:2021-06-12 # Python小白的数学建模课 @ Youcans # 2. SIR 模型,常微分方程组 相空间分析 from scipy.integrate import odeint # 导入 scipy.integrate 模块 import numpy as np # 导入 numpy包 import matplotlib.pyplot as plt # 导入 matplotlib包 def dySIR(y, t, lamda, mu): # SIR 模型,导数函数 i, s = y di_dt = lamda*s*i - mu*i # di/dt = lamda*s*i-mu*i ds_dt = -lamda*s*i # ds/dt = -lamda*s*i return np.array([di_dt,ds_dt]) # 设置模型参数 number = 1e5 # 总人数 lamda = 0.2 # 日接触率, 患病者每天有效接触的易感者的平均人数 sigma = 2.5 # 传染期接触数 mu = lamda/sigma # 日治愈率, 每天被治愈的患病者人数占患病者总数的比例 fsig = 1-1/sigma print("lamda={}\tmu={}\tsigma={}\t(1-1/sig)={}".format(lamda, mu, sigma, fsig)) # odeint 数值解,求解微分方程初值问题 tEnd = 200 # 预测日期长度 t = np.arange(0.0,tEnd,1) # (start,stop,step) s0List = np.arange(0.01,0.91,0.1) # (start,stop,step) for s0 in s0List: # s0, 易感者比例的初值 i0 = 1 - s0 # i0, 患病者比例的初值 Y0 = (i0, s0) # 微分方程组的初值 ySIR = odeint(dySIR, Y0, t, args=(lamda,mu)) # SIR 模型 plt.plot(ySIR[:,1], ySIR[:,0]) # 绘图 plt.title("Phase trajectory of SIR models") plt.axis([0, 1, 0, 1]) plt.plot([0,1],[1,0],'b-') plt.plot([1/sigma,1/sigma],[0,1-1/sigma],'b--') plt.xlabel('s(t)-youcans') plt.ylabel('i(t)-xupt') plt.text(0.8,0.9,r"$1/\sigma$ = {}".format(1/sigma),color='b') plt.show()
上图为例程 3.2 的运行结果( λ = 0.2 , μ = 0.08 , 1 / σ = 0.4 \lambda =0.2,\mu=0.08,1/\sigma=0.4 λ=0.2,μ=0.08,1/σ=0.4),是 SIR 模型的相轨迹图。
图中每一条 i-s 曲线,从直线 i(t)+s(t)=1 上的某一初值点出发,最终收敛于 s轴上的某一点,对应着某一个初值条件下的患病者与易感者比例随时间的变化关系。
利用相轨迹图可以分析和讨论 SIR 模型的性质:
对于小白来说,比较容易理解 2.4 节图中变量随时间的变化曲线,而对于本节相轨迹方法的思想、方法和图形都会觉得不容易理解,甚至感到困惑。虽然相轨迹的每一条线也对应着 t 从 t0 到 tend 的过程,但为什么要这么画,为什么轨迹这么怪怪的呢?
相轨迹图看上去比较怪,也不容易理解,是因为忽略时间轴而着重关注两个变量之间的关系,这种视角与我们日常观察问题和思考问题的习惯完全不同。也正是因为这个原因,相轨迹图能反映出时间变化曲线图中难以表达的一些重要特征。
例如,患病者比例在 s = 1 / σ s=1/\sigma s=1/σ 时达到峰值,即使把不同 σ \sigma σ 值下的患病者比例的时间变化曲线放在一张图中也无法观察到这一特征。进一步地,既然在 s = 1 / σ s=1/\sigma s=1/σ 时达到峰值,那么 s 0 s_0 s0 与 1 / σ 1/\sigma 1/σ 的关系自然就成为重要的分界线,并在图中可以观察到分界线两侧具有明显不同的特征。
有了对这些特征的认识和把握,才能选择不同的参数条件,在时间变化曲线图上进行比较系统的比较。要知道 SIR 模型中有 λ 、 μ \lambda、\mu λ、μ 两个参数,还有 i 0 、 s 0 i_0、s_0 i0、s0 两个初始条件,共有 4 个可以设置的参数都会影响微分方程的解,也就是会影响患病者、易感者比例的时间变化曲线。其中的各种组合无穷无尽,如果没有恰当的研究方法、不能把握内在的规律,即使在几十、几百组参数条件下进行模拟,仍然只是盲人摸象、管中窥豹。
看到这里,小白同学可能会对相轨迹研究的意义有所认识,但还是对这种分析方法难以理解望而却步。
没关系,还记得我们在”09 微分方程模型“中说的吗:
不会从问题建立微分方程模型怎么办,不会展开参数对稳定性、灵敏度的影响进行讨论怎么办?谁让你自己做呢,当然是先去找相关专业的教材、论文,从中选择比较接近、比较简单的理论和模型,然后通过各种假设强行将题目简化为模型中的条件,这就可以照猫画虎了。
最后,我们简单总结一下 SIR 模型的特点:
【本节完】
版权声明:
欢迎关注『Python小白的数学建模课 @ Youcans』 原创作品
原创作品,转载必须标注原文链接:https://blog.csdn.net/youcans/article/details/117843875。
Copyright 2021 Youcans, XUPT
Crated:2021-06-12
欢迎关注『Python小白的数学建模课 @ Youcans』系列,每周持续更新
Python小白的数学建模课-A1.国赛赛题类型分析
Python小白的数学建模课-A2.2021年数维杯C题探讨
Python小白的数学建模课-A3.12个新冠疫情数模竞赛赛题及短评
Python小白的数学建模课-B2. 新冠疫情 SI模型
Python小白的数学建模课-B3. 新冠疫情 SIS模型
Python小白的数学建模课-B3. 新冠疫情 SIR模型
Python小白的数学建模课-01.新手必读
Python小白的数学建模课-02.数据导入
Python小白的数学建模课-03.线性规划
Python小白的数学建模课-04.整数规划
Python小白的数学建模课-05.0-1规划
Python小白的数学建模课-06.固定费用问题
Python小白的数学建模课-07.选址问题
Python小白的数学建模课-09.微分方程模型
Python数模笔记-StatsModels 统计回归
Python数模笔记-Sklearn(5)支持向量机
Python数模笔记-NetworkX(5)关键路径法
Python数模笔记-模拟退火算法(4)旅行商问题