python - Scipy odeint 非负解

Question

显然，getting a non-negative solution from an ODE solver is non-trivial .在 Matlab 中，有 NonNegative option对于某些求解器，以获得非负解。 scipy 中是否有类似的选项？

如果不是，那么施加非负约束的“最佳”方式是什么？目前，我有以下内容:

def f(x, t, params):
     ... ... ...
     ... ... ...
     x_dot[(x <= 0) * (x_dot <= 0)] = 0.0
     return x_dot
... ... ...
x = odeint(f, x0, t, args=params)

然而，这会导致数值不稳定。我需要将 mxstep 设置为 1e8 和 hmin=1e-15。

Answer 1

问题不仅在于您必须避免对负 x 取平方。问题是施加约束的“最佳”方式仍然取决于您的系统应用程序是什么以及您认为什么行为是“合理的”。如果您的系统在 0 处没有平衡，那么您的问题可能是不适定的。它以非零速度移动到负 x 域的意义是什么？如果解释是解应该保持为零，那么您实际上不再有 ODE 系统作为您的预期模型:您有一个具有非光滑分量的混合动力系统，即当轨迹 x(t) 在t = t_1，对于所有 t > t_1，它必须保持在 x(t)。这可以通过适当的动态系统包(例如 PyDSTool)轻松实现。

或者，x=0 是一个稳定的平衡，您只需要防止对 x<0 的 f 进行评估。这也可以通过事件检测来破解。

在任何一种情况下，当您的 f 未定义为 x<0 时，x=0 处的事件检测都很棘手。几乎没有标准的 ODE 求解器可以在所有情况下强制避免在子域中进行评估，并且大多数事件检测将涉及对边界任一侧的评估。一个实用的解决方案是为 x 选择一个较小的数字，低于该数字是安全的(在您的应用程序上下文中)声明 x = 0。然后让事件检测 x 何时达到该值(假设您可以控制步长保持足够小)应该防止 x 被评估为负值。然后，如果这是您想要的行为，那么您将在该点之后创建一个条件使 x = 0。同样，这在 scipy/python 中有点大惊小怪，但你可以做到。在 PyDSTool 中设置所需的行为也相当容易，如果您在其帮助论坛中发帖，我愿意为您提供建议。