2

我发帖是因为我在数值上求解常微分方程并想约束状态变量。总而言之，我是一名模拟人口的生物学家，我希望当人口变得非常小时，他们就会灭绝。

当我运行我的模型时，种群可能会变得非常非常小（例如，10^{-166}），但永远不会变为 0。它们这样做对我来说很重要，因此我可以计算灭绝。但更现实的是，当人口的密度远小于地球上的原子密度时，这不太现实；）

即使在像下面这样简单的 2 种敌人-受害者模型中10^{-9}，我希望将其解释为 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”对消光进行建模，但是任何帮助指导我解决此类问题的一般答案/名称的帮助也将不胜感激。

PS 这类似于想用 0 （链接）替换负值的帖子，但不同，因为我希望人口不仅仅是非负数，而是无限地保持在 0。但是，在这篇文章中也没有关于如何做到这一点的好答案。

4

1 回答 1

2

这里举个例子。如前所述，非常务实。

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])
于 2019-12-03T21:42:30.923 回答