常微分方程 (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" 一起使用，因此再多的修修补补也不会给您带来明智的结果。我会在回答的第二部分谈到这一点，但首先，我会回答你原来的问题。

使用事件获取Non-Negative状态值

您可以使用事件在deSolve中强制non-negative人口规模。我建议您进一步阅读 deSolve 文档，因为触发事件的方法有很多种，但我们会将一种方法集成到您的代码中，在每个时间步长触发。因为它是 time-based，它需要某种最大时间来引用，所以我创建了一个 maxtime 值，您也可以在 for-loop 中使用它。您应该在其他功能可访问的地方定义此值；我只是把它放在你的 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，我们看到它只是在每个时间步检查我们的每个状态变量，并将所有负值设置为零。如果我们检查我们的输出，我们会看到它给了我们 non-negative 人口密度：

 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 时使用事件。鉴于您到目前为止所分享的关于您的模型的内容，肯定有关于在连续时间内对其建模的争论。仅仅因为传播事件是离散的并不一定意味着出生、死亡和感染事件也必须是离散的。重要的是要概念化这些现象在现实生活中发生的不同时间尺度，并问问自己是否用模型准确地捕捉到了它们。不过这已经超出了Whosebug的范畴，所以进入第二部分：

使用 deSolve

将 discrete-time 模型编码到 R 中

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

您已经清楚地在 continuous-time 模型中编码 return 是 衍生品 的值，而不是 discrete-time 模型 return 是您的 状态变量 的值。您在 dSv 中的出生成分是一个基本常数，因此您使用的是内在出生率 rate；在 discrete-time 模型中，您的出生成分将是后代的某个常数数量（几乎总是一个整数）乘以前一个时间步长的“parents”数量.

看看你的方程组，我们可以看出这是如何产生问题的：不是让 Sv 的每个个体在每个时间步长都生出 100 个个体，你是设置 Sv 到 100。 当您快速积累净亏损时，Sv 将不可避免地跌破零。

示例 discrete-time 模型如下所示：

# 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 的“当前”值，以便计算下一个时间步的种群动态。我希望这对你有所帮助，祝你在建模冒险中好运！

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

Ordinary differential equations (ODEs) - Is there any way to prevent negative values?

r

ode

differential-equations

使用事件获取Non-Negative状态值

使用 deSolve