SIR相平面r

yeotifhr 于 2023-03-27 发布在其他

关注(0)|答案(1)|浏览(139)

1.[~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（函数2，时间，y，rho）中出错：模型函数必须返回列表

r

来源：https://stackoverflow.com/questions/75842745/sir-phase-plane-r

1条答案

按热度按时间

uwopmtnx1#

原题中的代码不完整，附加注解中的代码有一些错误，所以我试着做了一个综合。
有几点：

变量N是冗余的
参数向量包含N作为第一个元素，因此参数被混淆。通常，通过命名向量而不是通过它们的位置来传递参数和状态变量更健壮。
有些乘法符号不见了
增加了steps的数量，以获得平滑的图形。

除了2D状态图，还可以使用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)

赞(0）回复(0）举报 2023-03-27

我来回答

SIR相平面r

1条答案

相关问题

热门标签

最新问答