r - 在R中的操纵()内部求解微分方程

我正在尝试分析一组微分方程(描述疾病在人群中的进展)中变化的起始条件和变量值如何影响系统的动力学(如通过图表所示)。我已经编写了代码，它可以很好地工作，但是迫使我更改值，然后重新运行整个代码。

因此，我试图将这段代码放到操纵()中，这样我就可以操纵变量并立即看到对产生的图的影响:

library(deSolve)
library(manipulate)

xyz <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    dt <- 1
    dv <- v0*v1*cos(2*pi*t/33)
    dX <- (v*N-B*X*Y-(mu+N/K)*X)
    dY <- (B*X*Y-(mu+m+g+N/K)*Y)
    dZ <- (g*Y-(mu+N/K)*Z)
    dN <- (v-mu-N/K)*N-m*Y

    return(list(c(dt,dv,dX, dY, dZ, dN)))
  })
}

times <- seq(0,365*3,by = 1)
init <- c(t=0,v=0.02,X=995,Y=5,Z=0,N=1000)
parameters <- c(B=0.14,mu=.01,m=.075,g=.025,K=10000,v0=.02,v1=.5)


manipulate(
  out <- as.data.frame(ode(y = init, times = seq(0,365*3,by=1), func = xyz, parms = parameters)),
  matplot(times,out[4:6],type="l",xlab="Time",ylab="Susceptibles and Recovereds",main="SIR Model",lwd=1,lty=1,bty="l",col=2:4),
  B=slider(0,1,initial=0.14,step=0.01)
)

无论我是否将全部或部分代码包含在操纵()中，在其外部或内部定义变量或其他任何内容，我都会不断收到错误消息。任何帮助将不胜感激!

最佳答案

当我第一次运行您的代码时，遇到的错误是:

Error in manipulate(out <- as.data.frame(ode(y = init, times = seq(0,  :
all controls passed to manipulate must be named

发生此错误的原因是，要操作的第二个参数是matplot()命令，而不是命名的控制参数(例如滑块)。因此，我将前两行放在花括号内，使它们成为单个表达式:

manipulate( {
    out <- as.data.frame(ode(y = init, times = seq(0,365*3,by=1), func = xyz, parms = parameters))
    matplot(times,out[4:6],type="l",xlab="Time",ylab="Susceptibles and Recovereds",main="SIR Model",lwd=1,lty=1,bty="l",col=2:4)
  }, B=slider(0,1,initial=0.14,step=0.01)
)

这样可以消除错误，但是移动滑块不会对绘图产生任何影响。为什么？因为名为B的滑块在传递给manipulate()的表达式中未引用任何内容。通过将parameters <- ...行移动到操纵表达式中，然后更改该行，以便有一个变量B(不仅仅是列表中的名称)，我解决了该问题。换句话说，我们需要B=B而不是B=0.14。现在，当您移动滑块时，图将发生变化，我相信这是您想要的:

manipulate( {
    parameters <- c(B=B,mu=.01,m=.075,g=.025,K=10000,v0=.02,v1=.5)
    out <- as.data.frame(ode(y = init, times = seq(0,365*3,by=1), func = xyz, parms = parameters))
    matplot(times,out[4:6],type="l",xlab="Time",ylab="Susceptibles and Recovereds",main="SIR Model",lwd=1,lty=1,bty="l",col=2:4)
  }, B=slider(0,1,initial=0.14,step=0.01)
)