R 中的 deSolve 给出了不正确的中间值

Question

我在使用 deSolve 的 R 中的 ode 求解器时遇到问题。

我发现求解器为应该始终等于零的变量给出了非零中间值。在我开始尝试使用 if(R>0){browser()} 之类的测试进行调试之前，它看起来不像是一个问题，它会被触发。

我的代码如下。提前致谢！

艾伦

library(deSolve)

simpleSIR <- function(t,states,par){
    with(as.list(c(states,par)),{

    S=states[1]
    I=states[2]
    R=states[3] 

    newinfections = beta*I*S

    dS <- -newinfections
    dI <- newinfections - gamma*I
    dR <- gamma*I

    if(R>0)
    {
        print(paste("why is this not zero?",R))
    }

    return(list(c(dS,dI,dR)))   
})}


par=list(beta=0.3,gamma=0.0)

init=c(0.99,0.01,0)

times <- seq(0,500,by = 1)


out <- as.data.frame(ode(y = init, times = times, func = simpleSIR, parms = par,maxsteps=2000))


plot(out[,2],type="l")
lines(out[,3],type="l",col="red")
lines(out[,4],type="l",col="blue")

Answer 1

首先，不知道大家是否了解state-variables和parameters的区别

R （在你的情况下）是一个状态变量。这意味着，在模拟开始时（t = 0）R = 0。但随着时间的发展，这可能会改变。查看您的状态变量 S！

您的参数 beta 和 gamma，它们在计算过程中不会改变。

然后你可以简化你的代码。使用此行而不是您对 par 和 init:

的定义

par <- c(beta = 0.3, gamma = 0.0)

init <- c(S = 0.99, I = 0.01, R = 0)

现在您可以删除以下行，因为 with(as.list(c(states, par)) 参数和状态参数的名称可以使用它们的名称 S、I 和 R。

S = states[1]
I = states[2]
R = states[3]

你也可以简化你的情节：

matplot(out[,1], out[,2:4], type = "l", col = c("black", "red", "blue"))

此致，

J_F

Answer 2

我认为你的问题在于 ode() 中数值积分的默认方法是 lsoda。此方法可以在刚性和非刚性问题的求解器之间切换。如果您切换到刚性求解器，雅可比矩阵将以数值方式计算，这可能会导致您看到的数值错误。在下面的代码中可以看到可能是这个原因：

out <- deSolve::ode(y = init, times = times, func = simpleSIR, parms = par,maxsteps=2000)
deSolve::diagnostics.deSolve(out)

"[...]14 到目前为止雅可比计算和 LU 分解的数量：23 [...]"

这对应于您的原始代码生成的打印消息数 (23)。您将通过使用像 RK4 这样的非刚性求解器来解决这个问题：

out.rk4 <- deSolve::ode(y = init, times = times, func = simpleSIR,method = "rk4", parms = par,maxsteps=2000)

如果您坚持使用 lsoda，您可能想尝试为 lsoda 提供您分析计算的雅可比矩阵。

R 中的 deSolve 给出了不正确的中间值

deSolve in R gives incorrect intermediate values

r

ode