原先写了篇论文，现在把代码分享给大家。

代码

import scipy.integrate as spi

import numpy as np

import matplotlib.pyplot as plt

import math

import xlrd

# https://zhuanlan.zhihu.com/p/104091330

# I_0为感染者的初始人数 524*1.5

I_0 = 109

# E_0为潜伏者的初始人数

E_0 = 223

# R_0为治愈者的初始人数

R_0 = 4

# S_0为易感者的初始人数

# S_0 = N - I_0 - E_0 - R_0

S_0=21536000

# Sq 隔离易感染者

S_q=1000

# Eq隔离潜伏者

E_q=200

# H住院患者

H=I_0+E_q

# T为传播时间

T = 35

#θ潜伏者相对感染者传播能力的比值-假设潜伏期患者与已经表现出症状的患者传染能力相同

cita=1

#γ隔离接触速率-隔离14天

gamma=1/14

# σ潜伏者向感染者转化速率 潜伏期7天

sigama=1/7

# α病死率----需基于实际死亡人数进行调整

alpha=0.013

# δI感染者隔离速率

delta_I=0.13

# γI感染者恢复速率

gamma_I=0.043

# δq 隔离潜伏者向隔离感染者的转化速率

delta_q=0.013

# γH隔离感染者的恢复速率

gamma_H=0.043

# β 传染概率

beta = 3.3*math.pow(10,-9)

# q 隔离比例

q=59*math.pow(10,-6)

# c 接触率

c=2.9

# ρ有效接触系数 ρc有效接触率

rou=1

#λ隔离解除速率

lamada=1/14

# INI为初始状态下的数组

INI = (S_0,E_0,I_0,R_0,S_q,E_q,H)

def funcSEIR(inivalue,_):

Y = np.zeros(7)

X = inivalue

Y[0]=-(rou*c*beta+rou*c*q*(1-beta))*X[0]*(X[2]+cita*X[1])+lamada*X[4]

Y[1]=rou*c*beta*(1-q)*X[0]*(X[2]+cita*X[1])-sigama*X[1]

Y[2]=sigama*X[1]-(delta_I+alpha+gamma_I)*X[2]

Y[3]=rou*c*q*(1-beta)*X[0]*(X[2]+cita*X[1])-lamada*X[4]

Y[4]=rou*c*beta*q*X[0]*(X[2]+cita*X[1])-delta_I*X[4]

Y[5]=delta_I*X[2]+delta_q*X[5]-(alpha+gamma_H)*X[6]

Y[6]=gamma_I*X[2]+gamma_H*X[6]

return Y

T_range = np.arange(0,T + 1)

RES = spi.odeint(funcSEIR,INI,T_range)

# plt.plot(RES[:,0],color = 'darkblue',label = 'Susceptible',marker = '.')

# plt.plot(RES[:,1],color = 'orange',label = 'Exposed',marker = '.')

plt.plot(RES[:,2],color = 'red',label = '模拟感染人数',marker = '.')

# plt.plot(RES[:,3],color = 'green',label = 'S_q',marker = '.')

# plt.plot(RES[:,4],color = 'blue',label = 'E_q',marker = '.')

# plt.plot(RES[:,5],color = 'black',label = 'H',marker = '.')

# plt.plot(RES[:,6],color = 'darkgreen',label = 'Recovery',marker = '.')

filename=r'paint1.xlsx'

file= xlrd.open_workbook(filename) # 打开文件

table=file.sheets()[3]#通过索引顺序获取 从0开始

x=table.row_values(2)

plt.rcParams['font.sans-serif']=['SimHei']#用来正常显示中文标签

plt.plot(x,label='现实感染人数')

plt.title('北京')

plt.legend()

plt.xlabel('天数')

plt.ylabel('人数')

plt.show()

效果图

参考文献

2. 理论参考-曹盛力,冯沛华,时朋朋.修正SEIR传染病动力学模型应用于湖北省2019冠状病毒病(COVID-19)疫情预测和评估[J/OL].浙江大学学报(医学版):1-13[2020-03-29].http://kns.cnki.net/kcms/detail/33.1248.R.20200303.1722.004.html.