SIR 相平面 r

Question

4. [~15%] 使用 deSolve 包中的函数 ode()，计算 SIR 模型的解 N = 100，β = 0.02，γ = 1，初始条件 S(0) = 99.9，I(0) = 0.1，时间点向量 范围从 0 到 20。在 (S,I) 相平面中绘制解。

对这个问题有帮助吗？

library(deSolve)

beta<-0.02
gamma<-1
N<-100
initial<-c(S=99.9, I=0.1, R=1000)
parms<-c(beta=beta, gamma=gamma, N=N)

times<-seq(0, 20 ,by=0.1)

output<-ode(y=initial, times=times, func=SIR_dydt, parms=parms)

plot(output[, "S", "I"], type="l", xlab="t", ylab="No. individuals, main="SIR solutions")
legend("topright", legend = c("S", "I"), lty = 1, col = c("blue", "red"))

试过这个代码但得到这个错误

checkFunc(Func2, times, y, rho) 错误： 模型函数必须返回一个列表

Answer 1

原题中的代码不完整，补充注释中的代码有一些错误，所以尝试综合一下

几点：

• 变量

  N

  是多余的
• 参数向量包含

  N

  作为第一个元素，因此参数被混淆了。一般来说，通过命名向量而不是它们的位置来传递参数和状态变量更健壮。
• 缺少一些乘号
• 增加了

  steps

  的数量以获得平滑的图形。

除了二维状态图外，包phaseR可用于添加可选的流场。

library("deSolve")

SIR_dydt <- function(t, y, parms) { 
   beta  <- parms[1] 
   gamma <- parms[2] 
   S <- y[1] 
   I <- y[2] 
   dSdt <- -beta * S * I 
   dIdt <-  beta * S * I - gamma * I 
   return(list(c(dSdt, dIdt))) 
} 
 
beta <- 0.02 
gamma <- 1 
## inital conditions: 
Szero <- 99.9 
Izero <- 0.1 

## time points 
steps <- seq(from = 0, to = 20, by = .1) 
SIRsolution <- ode(c(S = Szero, I = Izero), steps, SIR_dydt, c(beta, gamma)) 
 
## time series
plot(SIRsolution)

## state diagram
par(mfrow = c(1, 1)) # reset plotting parameters
plot(SIRsolution[,2], SIRsolution[,3], 
  type = "l", xlab = "S", ylab="I", 
  col = "magenta", main = "SI phase plane for SIR model")

### === additional goody, add flow field ========================
library("phaseR")
 
init <- c(s = Szero, I = Izero)
parms <- c(beta, gamma)

plot(SIRsolution[, c(2, 3)], type = "l", col = "magenta", 
     xlim = c(0, 100), ylim=c(0, 20),
     main="SI phase plane for SIR model")
flow <- flowField(SIR_dydt, xlim = c(0, 100), ylim = c(0, 20),
          parameters = parms)