r - 微分方程解中非常小的总体/状态变量为零

Question

我发帖是因为我正在数值求解常微分方程并且想要约束状态变量。总之，我是一名对种群进行建模的生物学家，希望种群在变得非常小时就灭绝。

当我运行模型时，种群数量可能会变得非常非常小(例如 )，但永远不会变成 0。它们这样做对我来说很重要，这样我就可以计算灭绝。但更现实的是，当人口密度远小于地球上原子的密度时，这不太现实;)

即使在像下面这样的简单的 2 种敌人-受害者模型中，密度也达到 我希望被解释为 0。

library(package = "deSolve")

lv <- function(times, state, parms) {
    with(as.list(c(state, parms)), {
    dR <- 2*R - 0.5*R*C
    dC <- 0.2*R*C - 0.6*C
    return(list(c(dR, dC)))
    })
}

time_vec <- seq(from = 0, to = 100, length.out = 1e4)
y_0 <- c(R = 50, C = 20)

out <- ode(y = y_0, times = time_vec, func = lv, parms = NULL, method = "lsoda")
min(out[,-1])
plot(x = out[,2], y = out[,3], type = "l")

我非常希望了解如何使用 R 中的解算器“deSolve”对灭绝进行建模，但如果有任何帮助指导我找到此类问题的一般答案/名称，我也将非常感激。

附注这类似于想要用 0 替换负值 ( link ) 的帖子，但有所不同，因为我希望总体不只是非负数，而是无限地保持在 0。而且，这篇文章中没有关于如何做到这一点的良好答案。

Answer 1

这是一个例子。正如所说，非常务实。

library(package = "deSolve")

lv <- function(times, state, parms) {
  with(as.list(c(state, parms)), {
    dR <- 2*R - 0.5*R*C
    dC <- 0.2*R*C - 0.6*C
    return(list(c(dR, dC)))
  })
}

time_vec <- seq(from = 0, to = 100, length.out = 1e4)
y_0 <- c(R = 50, C = 20)
out <- ode(y = y_0, times = time_vec, func = lv, parms = NULL, method = "lsoda")
min(out[,-1])

eps <- 1e-2
## event triggered if state variable <= eps
rootfun <- function (t, y, pars) {
  return(y - eps)
}

## sets state variable = 0
eventfun <- function(t, y, pars) {
  if (y[1] <= eps) y[1] <- 0
  if (y[2] <= eps) y[2] <- 0
  return(y)
}

out1 <- ode(y = y_0, times = time_vec, func = lv, parms = NULL, method = "lsoda",
           rootfun = rootfun,
           events = list(func = eventfun, root = TRUE))

plot(out, out1)
min(out1[,-1])