用零替换模型(ODE系统)中的负值 [英] Replacing negative values in a model (system of ODEs) with zero
问题描述
我目前正在使用deSolve
求解常微分方程组,并且想知道是否有任何方法可以防止微分变量值降至零以下.我还看过其他一些有关在向量,数据帧等中将负值设置为零的文章,但是由于这是一个生物学模型(T细胞计数变为负数没有意义),需要从一开始就阻止它发生,以便这些值不会使结果产生偏差,而不仅仅是替换最终输出中的负数.
I'm currently working on solving a system of ordinary differential equations using deSolve
, and was wondering if there's any way of preventing differential variable values from going below zero. I've seen a few other posts about setting negative values to zero in a vector, data frame, etc., but since this is a biological model (and it doesn't make sense for a T cell count to go negative), I need to stop it from happening to begin with so these values don't skew the results, not just replace the negatives in the final output.
推荐答案
我的标准方法是将状态变量转换为不受限制的比例.对正变量执行此操作的最明显/标准方法是为log(x)
而不是x
的动力学写下方程.
My standard approach is to transform the state variables to an unconstrained scale. The most obvious/standard way to do this for positive variables is to write down equations for the dynamics of log(x)
rather than of x
.
例如,对于传染病流行的敏感传染病恢复(SIR)模型,当方程为dS/dt = -beta*S*I; dI/dt = beta*S*I-gamma*I; dR/dt = gamma*I
时,我们会天真地将梯度函数写为
For example, with the Susceptible-Infected-Recovered (SIR) model for infectious disease epidemics, where the equations are dS/dt = -beta*S*I; dI/dt = beta*S*I-gamma*I; dR/dt = gamma*I
we would naively write the gradient function as
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-gamma
;其余等式需要使用exp(logI)
来引用I
.所以:
If we make log(I) rather than I be the state variable (in principle we could do this with S as well, but in practice S is much less likely to approach the boundary), then we have d(log(I))/dt = (dI/dt)/I = beta-gamma
; the rest of the equations need to use exp(logI)
to refer to I
. So:
gfun_log <- function(time, y, params) {
g <- with(as.list(c(y,params)),
c(-beta*S*exp(logI),
beta-gamma,
gamma*exp(logI))
)
return(list(g))
}
(一次计算exp(logI)
并存储/重用它比两次计算的效率更高……)
(it would be slightly more efficient to compute exp(logI)
once and store/re-use it rather than computing it twice ...)
这篇关于用零替换模型(ODE系统)中的负值的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!