Python小白的数学建模课-B4. 新冠疫情 SIR模型

Posted 2021-06-27 youcans

传染病的数学模型是数学建模中的典型问题，常见的传染病模型有 SI、SIR、SIRS、SEIR 模型。

SIR 模型将人群分为易感者（S类）、患病者（I类）和康复者（R 类），考虑了患病者治愈后的免疫能力。

本文详细给出了 SIR 模型微分方程、相空间分析的建模、例程、结果和分析，让小白都能懂。

『Python小白的数学建模课 @ Youcans』带你从数模小白成为国赛达人。

---

1. 疫情传播 SIR 模型

传染病的传播特性不可能通过真实的试验开展研究，因此需要针对不同的传染病传播方式和流行特点建立相应的数学模型，并对模型进行理论研究和数值模拟。通过研究发现传染病传播的特征阈值，就可以为预防和控制传染病提供数据支撑和防控策略。

1927年，W. Kermack 在论文 “Contributions to the mathcmatical theory of epidemics” 中研究了伦敦黑死病和孟买瘟疫的流行过程，创造性地提出了 SIR 模型。

SI 模型和 SIS 模型将人群分为感染者（S类）和患病者（I类）两类人群，但大多数传染病的患者在治愈后就有很强的免疫力，终身或在一段时期内不再会被感染而变成病人。这类人群称为病愈免疫的康复者（R 类）。康复者已经退出传染系统，对于致死性疾病的死亡者也可以用该类别描述其传播特性。

SIR 模型适用于具有易感者、患病者和康复者三类人群，可以治愈，且治愈后终身免疫不再复发的疾病，例如天花、肝炎、麻疹等免疫力很强的传染病。

SIR 模型假设：

1. 考察地区的总人数 N 不变，即不考虑生死或迁移；
2. 人群分为易感者（S类）、患病者（I类）和康复者（R 类）三类；
3. 易感者（S类）与患病者（I类）有效接触即被感染，变为患病者（I类）；患病者（I类）可被治愈，治愈后变为康复者；康复者（R类）获得终身免疫不再易感；无潜伏期；
4. 将第 t 天时 S类、I 类、R 类人群的占比记为 $s(t)$、$i(t)$、$r(t)$，数量为 $S(t)$、$I(t)$、$R(t)$；初始日期 $t=0$ 时， S类、I 类、R 类人群占比的初值为 $s_0$、$i_0$、$r_0$。
5. 日接触数 $\lambda$，每个患病者每天有效接触的易感者的平均人数；
6. 日治愈率 $\mu$，每天被治愈的患病者人数占患病者总数的比例，即平均治愈天数为 $1/\mu$；
7. 传染期接触数 $\sigma =\lambda /\mu$，即每个患病者在整个传染期内有效接触的易感者人数。

SIR 模型的微分方程：
由
$$\{\begin{array}{ll}& N\frac{di}{dt}=N\lambda si\\ & N\frac{dr}{dt}=N\mu i\\ & s(t)+i(t)+r(t)=1\end{array}$$
得：
$$\{\begin{array}{llc}& \frac{di}{dt}=\lambda si-\mu i,& i(0)=i_0\\ & \frac{ds}{dt}=-\lambda si,& s(0)=s_0\end{array}$$

SIR 模型不能求出解析解，只能通过数值计算方法求解。

---

2. SIR 模型的 Python 编程

2.1 Scipy 工具包求解微分方程组

SIS 模型是常微分方程初值问题，可以使用 Scipy 工具包的 scipy.integrate.odeint() 函数求数值解。

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

**scipy.integrate.odeint() **是求解微分方程的具体方法，通过数值积分来求解常微分方程组。

odeint() 的主要参数：

• func: callable(y, t, …) 　　导数函数 $f(y,t)$ ，即 y 在 t 处的导数，以函数的形式表示
• y0: array：　　初始条件 $y_0$，注意 SIR模型是二元常微分方程组， 初始条件为数组向量 $y_0=[i_0,s_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 值。

2.2 odeint() 求解 SIR 模型的编程步骤

1. 导入 scipy、numpy、matplotlib 包。
2. 定义导数函数 $f(y,t)$。注意对于常微分方程（例如 SIS模型）和常微分方程组（SIR模型），y 分别表示标量和向量，函数定义略有不同，以下给出两种情况的例程以供对比。

常微分方程的导数定义（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 小白，又涉及到看着就心烦的微分方程组，所以我们宁愿把程序写得累赘一些，便于读者将程序与前面的微分方程组逐项对应。

3. 定义初值 $y_0$ 和 $y$ 的定义区间 $[t_0,t]$，注意初值为数组向量 $y_0=[i_0,s_0]$。
4. 调用 odeint() 求 $y$ 在定义区间 $[t_0,t]$ 的数值解。

SIR 模型是二元常微分方程，返回值 y 是 len(t)*2 的二维数组。

2.3 Python例程：SI 模型、SIS 模型与SIR 模型的比较

# 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=Python小白的数学建模课-A3. 12个新冠疫情数模竞赛赛题与点评

 Python小白的数学建模课-A3. 10个新冠疫情数模竞赛赛题与点评
 Python小白的数学建模课-B2. 新冠疫情 SI模型
 Python小白的数学建模课-B5. 新冠疫情 SEIR模型
 Python小白的数学建模课-B5. 新冠疫情 SEIR模型
 Python小白的数学建模课-B6. 新冠疫情 SEIR 改进模型