微分方程解中非常小的 populations/state 变量为零

Question

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

当我运行我的模型数量可能变得非常非常小（例如，），但永远不会变为 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”对消光进行建模，但是任何帮助指导我解决此类问题的一般 answers/names 也将不胜感激。

P.S。这类似于 post 想要将负值替换为 0 ()，但有所不同，因为我希望人口不仅是非负的，而且无限地保持在 0。但是，在 post.

中也没有关于如何做到这一点的好答案

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])

微分方程解中非常小的 populations/state 变量为零

Very small populations/state variables to zero in differential equation solutions

r

ode

numerical-integration

population