我目前正在使用deSolve 求解常微分方程组，并且想知道是否有任何方法可以防止微分变量值低于零。我看过其他一些关于在向量、数据框等中将负值设置为零的帖子，但由于这是一个生物模型（T 细胞计数变为负数没有意义），我需要从一开始就阻止它发生，这样这些值就不会扭曲结果，而不仅仅是替换最终输出中的负数。
【问题讨论】：
x[1]=max(x[1],0) 然后动态应该是正确的，即使它需要后处理。
标签： r ode negative-number
我的标准方法是将状态变量转换为不受约束的规模。对正变量执行此操作的最明显/标准方法是写下 log(x) 而不是 x 的动力学方程。
例如，对于传染病流行的易感感染恢复 (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) 并存储/重用它会比计算两次更有效...）
【讨论】：
如果某个值在现实中没有变为负值，但在您的模型中变为负值，则您应该更改模型，或者等效地修改微分方程，这样就不可能了。换句话说：不要试图约束你的动态变量，而是它们的导数。 其他一切只会导致你的求解器出现问题，而它不应该关心微分方程的变化。
举个简单的例子，假设：
在这种情况下，y 只能在 f(0) f em> 使得 f(0) ≥ 0（并且它仍然是平滑的）。
为了证明原理，您可以将 f 与经过适当修改的 sigmoid function 相乘（它允许您用平滑函数组合每个逻辑运算）。这样，y, 的大多数值都不会发生任何变化，并且只有在 y 接近 0 时（即，当您无论如何都要操作事物时）才更改微分方程.
但是，如果不考虑您的模型，我不会真的推荐使用 sigmoid。如果您的模型在 y = 0 附近完全错误，那么它很可能对于附近的值已经无用。如果您的模拟在这种情况下冒险，并且您希望结果有意义，您应该解决这个问题。
【讨论】：