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

Question

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

9

我目前正在使用 deSolve 解决一组普通微分方程组，想知道是否有任何方法可以防止微分变量值低于零。我看到了一些其他帖子关于在向量、数据框等中将负值设置为零，但由于这是生物模型（T细胞计数不可能为负），我需要阻止它发生，以便这些值不会影响结果，而不仅仅是在最终输出中替换负数。

- Dorian

你正在寻找参数约束吗？ - Roman Luštrik

这可能有所帮助——我的主要问题不在于参数，而在于独立变量。是否有任何类似于MATLAB的非负函数的等效物，或者我可以实现的代码，以使x [1]永远不会低于0？ - Dorian

当x[1]达到零时，方程会发生什么？零是一个稳定的固定点吗？如果是这样，您可以将第一个负值和比它大的所有值都设置为零。此外，在定义ODE的函数中，您可以设置x[1]=max(x[1],0)，然后即使需要后处理，动态也应该是正确的。 - goryh

2个回答

0

如果一个值在现实中不会变成负数，但在您的模型中变成了负数，那么您应该更改您的模型或者等效地修改您的微分方程，使其不可能出现这种情况。换句话说：不要试图限制您的动态变量，而是限制它们的导数。其他任何方法都只会导致求解器出现问题，而它不应该关心微分方程的变化。

举个简单的例子，假设：

您有一个一维微分方程ẏ = f(y)，
y不应该变成负数，
您的初始y为正数。

在这种情况下，y只有在f(0) < 0时才会变成负数。因此，您只需要修改f，使得f(0) ≥ 0（并且仍然平滑）即可。

为了证明这个原理，您可以将f与适当修改的sigmoid函数相乘（这允许您使用平滑函数组合每个逻辑操作）。这样，在大多数值y的情况下，没有任何改变，并且只有在y接近0时才更改微分方程，即当您要进行操作时。

但是，我不建议您在不考虑模型的情况下使用sigmoid。如果您的模型在y=0附近完全错误，那么它很可能已经对附近的值无用。如果您的模拟涉及到此地形，并且希望结果具有意义，则应修复此问题。

- Wrzlprmft

很遗憾，我认为这是不可能的，因为变量彼此相互作用。例如，从某种细胞类型推导出病毒载量的方程式为dVB <- (k_b * x [7]) - (d_vb * x [1])，其中k_b = 病毒释放速率，x [7] = 感染细胞，d_vb = 病毒死亡率，而x [1] = 病毒载量。也许我太无知了，但我看不到重写这些方程式以避免在一定点之后下降到0以下的方法（在这种情况下，当裂解感染的细胞用尽时），这就是为什么我希望可能有一些约束条件可以防止这种情况发生的原因。 - Dorian

@Dorian：你绝对可以实现这样的约束，即使是针对交互变量：Sigmoid函数可以让你实现所有逻辑操作。然而，如果你改进我的模型，很可能会有更合理的方法。我强烈建议你考虑这一点，因为如果你的模型在0处是错误的，那么它几乎肯定在接近0处也是错误的，因此你的结果很可能是无用的。另请参见我的编辑。 - Wrzlprmft

@Dorian：此外，如果我正确理解了您的模型，dVB是x [1]（即ẋ [1]）的时间导数，而k_b、x [7]和d_vb始终为正。因此，您的方程已经确保x [1]永远不会变为负数，因为对于x [1] = 0，您自动具有非负的dVB = ẋ [1]。因此，您的模型已经通过成为一个良好的模型来避免了所讨论的问题。 - Wrzlprmft

@Wrzlprmft：在数学上是正确的，但在数字环境中不一定如此。 - Ben Bolker

1

@BenBolker：是的和不是。例如，在提问者的示例中，稍微负值将被“推回”到0，因此这样的数值误差不是问题。对于大多数模型也应该如此（如果不是，您可以再次修改微分方程来解决这个问题）。此外，通过良好的积分器和合理选择的参数，您可以完全避免这种情况。使用自适应步长的Runge-Kutta积分器，相对误差受控制，提问者的示例永远不会低于0。 - Wrzlprmft

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ben Bolker · Accepted Answer

我的标准方法是将状态变量转换为不受限制的尺度。对于正变量来说，最明显/标准的做法是编写关于 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) 一次并将结果存储/重复使用比两次计算更为高效...