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()