r - 常微分方程 (ODE) - 有什么方法可以防止出现负值吗？

Question

我正在尝试在每个空间网格单元上应用常微分方程组 (ODE)。因此，每个景观单元都有一个关联的 ODE 模型。易感蚊子 ( Sv )、暴露蚊子 ( Se ) 和受感染蚊子 ( St ) 的数量在每个时间步都会更新，并且 ODE 模型与基于离散代理的动物运动模型相结合。以下是运行 ODE 模型的示例:

library(deSolve)

mod1 <- function(out_tab, time_step, var){
 Sv <- out_tab[time_step,c("Sv")]
  Ev <- out_tab[time_step,c("Ev")]
  Iv <- out_tab[time_step,c("Iv")]

  Nh <- out_tab[time_step,c("Nh")]
  Ih <- out_tab[time_step,c("Ih")]
  bv <- 100
  dv <- 0.07
  betav <- 0.33
  av <- 0.5
  muv <- 0.1

  mod <- function(times, states, parameters) {
    with(as.list(c(states, parameters)), {

      dSv <- bv - dv*Sv - betav*av*(Ih/Nh)*Sv
      dEv <- betav*av*(Ih/Nh)*Sv - dv*Ev - muv*Ev
      dIv <- muv*Ev - dv*Iv

      return(list(c(dSv, dEv, dIv)))

    })
  }

  states <- c(Sv = Sv, Ev = Ev, Iv = Iv)
  input_parameters <- c(bv = bv, dv = dv, betav = betav, av = av, muv = muv)

  ## Solve ODEs
  out <- ode(func=mod, y=states, times=seq(time_step, time_step + 1, by=1), parms=input_parameters, method="iteration")
  out <- as.data.frame(out)

  out_tab[time_step + 1, c("Sv", "Ev", "Iv")] <- out[dim(out)[1], c("Sv", "Ev", "Iv")]
  # out_tab[time_step + 1, c("Nh")] <- var
  # out_tab[time_step + 1, c("Ih")] <- var

  return(out_tab)

}
out_tab <- data.frame(Sv = 10, Ev = 3, Iv = 2, Nh = 2, Ih = 1)
Nh <- c(5, 1, 8, 0, 5)
Ih <- c(2, 0, 4, 0, 1)

for(time in 1:2){
  out_tab <- mod1(out_tab = out_tab, time_step = time)
  ## print(out_tab)
  out_tab[time + 1, c("Nh")] <- Nh[time + 1]
  out_tab[time + 1, c("Ih")] <- Ih[time + 1]
}

在time = 2 (以天为单位)，Ih (单元格中的动物数量)等于 0，因此 Ev 的值是负数。有什么办法可以防止 ODE 出现负值吗？我用的方法"iteration"因为 ODE 模型被合并到基于离散主体的模型中。

Answer 1

这里有几个主要问题；长话短说，您错误地定义了模型系统以与 method = "iteration" 一起使用。 ，因此无论进行多少轻微的修改都不会得到合理的结果。我将在答案的第二部分中讨论这一点，但首先，我将回答您原来的问题。

使用事件获取非负状态值

您可以强制deSolve中的非负人口规模使用事件。我建议您阅读deSolve进一步文档，因为触发事件的方法有很多，但我们会将一种方法集成到在每个时间步骤期间触发的代码中。因为它是基于时间的，所以需要某种最大时间来引用，所以我创建了一个 maxtime您也可以在 for 循环中使用该值。您应该在其他函数可以访问的地方定义这个值；我只是把它放在你的mod1之前函数声明。

# Here we declare the maximum time to which your system will evaluate
maxtime <- 2

# This is where we define your event function
# Add this directly above your call to ode()
  posfun <- function(t, y, parms){
    with(as.list(y), {
      y[which(y<0)] <- 0  
      return(y)
    })
  }

# Here's your original call to ode(), with a small addition
# Notice that we added events; iteration is missing, more on that later
  out <- ode(func=mod, 
             y=states, 
             times=seq(time_step, time_step + 1, by=1), 
             parms=input_parameters, 
             events=list(func = posfun, time = c(0:maxtime)))
  
  out <- as.data.frame(out)

# Don't forget to add maxtime to your for-loop
for(time in 1:maxtime){
  out_tab <- mod1(out_tab = out_tab, time_step = time)
  ## print(out_tab)
  out_tab[time + 1, c("Nh")] <- Nh[time + 1]
  out_tab[time + 1, c("Ih")] <- Ih[time + 1]

查看posfun ，我们看到它只是在每个时间步检查每个状态变量，并将任何负值设置为零。如果我们检查输出，我们会发现它给出了非负人口密度:

 out_tab
        Sv       Ev       Iv Nh Ih
1  10.0000  3.00000 2.000000  2  1
3 179.7493 16.67784 3.244288  1  0
4 416.2958 10.01499 6.133576  8  4

太棒了，对吧？嗯，不完全是。不幸的是，目前还没有办法在 method = iteration 时使用事件。 。鉴于到目前为止您所分享的有关模型的内容，对于连续时间建模肯定存在争议。仅仅因为传播事件是离散的，并不一定意味着出生、死亡和感染事件也必须是离散的。重要的是要概念化这些现象在现实生活中发生的不同时间尺度，并问问自己是否用模型准确地捕捉到了它们。然而，这已经超出了 StackOverflow 的范围，所以进入第二部分:

使用 deSolve 将离散时间模型编码到 R 中

查看 iteration 的文档deSolve中的方法，我们看到:“方法“迭代”的特殊之处在于，这里函数 func 应该返回状态变量的新值而不是变化率。”

您已经清楚地在返回导数值的连续时间模型中进行了编码，而不是返回状态变量的值的离散时间模型。您的出生成分在 dSv是一个纯粹的常数，因此您使用的是内在出生率；在离散时间模型中，您的出生分量将是某个常数数量的后代(几乎总是整数)乘以上一个时间步的“ parent ”数量。

看看你的方程组，我们可以看到这是如何产生问题的:而不是让 Sv 的每个个体都存在。在每个时间步产生 100 个个体，您设置 Sv到 100。 不可避免的是 Sv当您快速累积净亏损时，该值将暴跌至零以上。

离散时间模型示例如下:

# discrete-time host-parasite model
Parasite <- function(t, y, ks) {
P <- y[1]
H <- y[2]
f <- A * P / (ks + H)
Pnew <- H * (1 - exp(-f))
Hnew <- H * exp(rH * (1 - H) - f)
list (c(Pnew, Hnew))
}

请注意我们如何不断引用 P 和 H 的“当前”值，以便计算下一个时间步的种群动态。我希望这对您有所帮助，并祝您在建模冒险中一切顺利!