tiancanrong的个人博客分享 http://blog.sciencenet.cn/u/tiancanrong

博文

一种SEIR模型的R代码

已有 7619 次阅读 2020-2-21 19:11 |系统分类:科研笔记

SEIR是很经典的模型,现在随着城市之间交通的日益频繁,考虑人口城市交通工具扩散的SEIR模型似乎更加符合实际。微分方程模型如下:

\begin{eqnarray*}

\frac{\mathrm{d}\vec{S}_t}{\mathrm{d}t} & = & -\vec{S}_t \odot \vec{\lambda}_t \\

\frac{\mathrm{d}\vec{E}_t}{\mathrm{d}t} & = & \vec{S}_t \odot \vec{\lambda}_t - \alpha \vec{E}_t \\

\frac{\mathrm{d}\vec{I}_t}{\mathrm{d}t} & = & \alpha \vec{E}_t - \gamma \vec{I}_t \\

\frac{\mathrm{d}\vec{R}_t}{\mathrm{d}t} & = & \gamma \vec{I}_t

\end{eqnarray*}

where $\vec{S}_t = (S_{1t},\dots,S_{nt})^T$ and likewise for $\vec{E}_t$, $\vec{I}_t$, and $\vec{R}_t$.  Furthermore,

$$ \vec{\lambda}_t = \beta \left(\vec{I}/\vec{N} + (K \cdot (\vec{I_t}/\vec{N}))/\vec{N}\right) $$

解模型的R语言如下:

test <- load("china_population.rda")

K=diag(187)  #%% 城市之间交通流矩阵

N<<-china_population  #%% 城市人口向量

K<<-K

alpha=1/4

beta=2

gammaa=1/6

I0W = 1

param=c(beta, gammaa)

func = function(t, y, parms) {

  n <- length(N) # number of cities

  alpha=1/4

  ymat = matrix(y, nrow=n, byrow=FALSE)

  S = ymat[,1]

  E = ymat[,2]

  I = ymat[,3]

  R = ymat[,4]

  beta = parms[1]

  gamma = parms[2]

  

  lambda = beta * (I/N + (t(K/N)%*%I/N))

  dS = -lambda * S

  dE = lambda * S - alpha * E

  dI = alpha * E - gamma * I

  dR = gamma * I

  

  list(c(dS, dE, dI, dR))

}


max_t=22

n = length(N) # number of cities

t = 0:max_t

 # Initial infectives in Wuhan

  I0 = rep(0, length(N))

  I0[151] = I0W

  y = c(N-I0, rep(0, length(N)), I0, rep(0, length(N)))

  sim = deSolve::ode(y=y, times=t, func=func, parms=param)

  t=sim[,1]

  S=sim[,seq(2,len=n)]

  E=sim[,seq(2+n, len=n)]

  I=sim[,seq(2+2*n, len=n)]

  R=sim[,seq(2+3*n, len=n)]

  simulate=list(t=t, S=S, E=E, I=I, R=R)

  write.csv(I,"mydata_I.csv")

  write.csv(S,"mydata_S.csv")

  write.csv(E,"mydata_E.csv")

  write.csv(R,"mydata_R.csv")

    plot(1:23, I[,151],type="o",

         main="I")

    grid()

 box()




https://blog.sciencenet.cn/blog-306734-1219685.html

上一篇:2013-2017年Journal of differential equations高被引文章
下一篇:动力系统的解能不能用动态规划来做

2 张鹰 宁利中

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2022-5-27 07:28

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部