用零替换模型（ODE 系统）中的负值

Question

我目前正在使用 deSolve 求解常微分方程组，想知道是否有任何方法可以防止微分变量值低于零。我看过其他一些关于在向量、数据框等中将负值设置为零的帖子，但由于这是一个生物模型（T 细胞计数变为负数没有意义），我需要从一开始就阻止它发生，这样这些值就不会扭曲结果，而不仅仅是替换最终输出中的底片。

Answer 1

如果一个值在现实中没有变成负数，但在您的模型中变成了负数，您应该更改您的模型，或者，等同地，修改您的微分方程，这样就不可能发生这种情况。换句话说：不要试图约束你的动力学变量，而是它们的导数。 其他一切只会导致你的求解器出现问题，而它不应该关心微分方程的变化。

举个简单的例子，假设：

你有一个 one-dimensional 微分方程 ẏ = f(y),
y不得变为负数,
你的初始 y 是正值。

在这种情况下，y 只能在 f(0) < 0 时变为负数。因此，您所要做的就是修改 f 使得 f(0) ≥ 0（并且仍然平滑）。

为了证明原理，您可以将 f 与适当修改的 sigmoid function 相乘（这样您就可以用平滑的函数组合每个逻辑运算）。这样，y, 的大多数值都不会发生任何变化，并且只有在 y 接近 0 时才更改微分方程，即当您无论如何我们都会操纵事情。

但是，如果不考虑您的模型，我真的不会推荐使用 sigmoid。如果您的模型在 y = 0 附近完全错误，它很可能已经对附近的值无用。如果您的模拟在这个领域冒险，并且您希望结果有意义，您应该解决这个问题。

Answer 2

我的标准方法是将状态变量转换为无约束尺度。对于正变量，最 obvious/standard 的方法是写下 log(x) 的动力学方程，而不是 x.

的动力学方程

例如，对于传染病流行的 Susceptible-Infected-Recovered (SIR) 模型，方程为 dS/dt = -beta*S*I; dI/dt = beta*S*I-gamma*I; dR/dt = gamma*I 我们会天真地将梯度函数写为

gfun <- function(time, y, params) {
   g <- with(as.list(c(y,params)),
       c(-beta*S*I,
          beta*S*I-gamma*I,
          gamma*I)
       )
   return(list(g))
}

如果我们让 log(I) 而不是 I 成为状态变量（原则上我们也可以用 S 做到这一点，但实际上 S 很多不太可能接近边界），那么我们有 d(log(I))/dt = (dI/dt)/I = beta*S-gamma；其余等式需要用exp(logI)代指I。所以：

gfun_log <- function(time, y, params) {
   g <- with(as.list(c(y,params)),
       c(-beta*S*exp(logI),
          beta*S-gamma,
          gamma*exp(logI))
       )
   return(list(g))
}

（计算一次 exp(logI) 和计算一次 store/re-use 比计算两次更有效 ...）

用零替换模型（ODE 系统）中的负值

Replacing negative values in a model (system of ODEs) with zero

r

ode

negative-number