在 R 中使用函工具的非负 ODE 解决方案?

Question

我正在尝试实现一个处理非负常微分方程系统的R算法。我需要像 MATLAB 中的 ode45 这样的东西来定义必须为非负的状态。

我在 3 年前就已经讨论过这个问题，但没有真正的解决方案。 deSolve仍然不是要走的路。我发现了一些看起来很有前途的python代码。也许这在 R 中也是可能的。最后，我必须定义一个函数包装器，作为python中的functools。它的作用非常简单。这是python包装器的代码：

def wrap(f):
@wraps(f)
def wrapper(t, y, *args, **kwargs):
low = y < 0
y = np.maximum(y, np.ones(np.shape(y))*0)
result = f(t, y, *args, **kwargs)
result[too_low]  = np.maximum(result[low], np.ones(low.sum())*0)

return result
return wrapper
return wrap

我的意思是在python中这是直截了当的。包装器将在 调用的集成的每个步骤中使用

solver = scipy.integrate.odeint(f, y0)
solution = solver.solve()

在 R 中可能相同吗？我知道还有一个functools包和functools函数。但我不知道这是否真的有效。我可以为此使用 deSolve 中的事件吗？

我现在在这个项目上工作了 5 年，我没有想法。我使用了MATLAB，C++和Python接口，但所有这些都是慢下来的，我在R中需要它。 非常感谢您的帮助！

Answer 1

deSolve 不支持自动非负性约束，这是有充分理由的。我们过去曾多次遇到过这样的问题，但事实证明，在所有这些情况下，负值的原因是模型规格不完整。典型情况是从空池中导出某些内容。由于不需要的负值通常是模型规范不足的指标，因此我们(目前)不考虑将来添加"非负"约束。

示例：在以下等式中，X 可以通过模型设计变为负数：

dX/dt = -k

而以下不能：

dX/dt = -k * X

如果您需要"大部分时间"的线性递减，在 X 变为零之前不久减少到零，您可以使用 Monod 类型的保护措施(或类似的东西)：

dX/dt = -k*X/(k2 + X)

k2的选择相对不重要。与求解器的数值精度相比，它应该足够小，不会影响整体行为，也不要太小。

避免负值的另一种方法是在对数转换的空间中工作。以下是一些相关线程：

https://stat.ethz.ch/pipermail/r-sig-dynamic-models/2010q2/000028.html

https://stat.ethz.ch/pipermail/r-sig-dynamic-models/2013q3/000222.html

https://stat.ethz.ch/pipermail/r-sig-dynamic-models/2016/000437.html

此外，当然也可以在 R 中编写自己的包装器。

希望对你有帮助