Scipy odeint 非负解

Question

显然，从ODE求解器中得到非负解并不容易。在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。

- carmichael561

---

正如您提供的链接所指出的那样，解决这个问题并不一定容易。您的非负解的性质是什么，为什么会引起问题？至少有三种可能性：（1）解收敛于0的稳定平衡状态，但由于常规数值误差，它会稍微低于0（并可能在0周围振荡）（例如dx/dt = -2*x）；（2）0是“半稳定”的平衡状态，因此如果解变为负数，它会爆炸（例如dx/dt = -x ** 2）；（3）微分方程在负x处未定义（例如dx/dt = -sqrt（x）。 - Warren Weckesser

（3）由于存在平方根，ODE在负x处未定义。 - carmichael561

我对你施加了类似的限制，并没有遇到任何问题。你能发布你的ODE系统吗？ - Daniel Farrell

Answer 1

问题不仅仅在于要避免平方根的负x，而是“最佳”实施约束的方式依赖于您系统的应用和您认为“合理”的行为。如果您的系统在0处没有平衡点，则您的问题可能存在问题。将其移动到负x域中非零速度的意义是什么？如果解释是解决方案应保持在零，则实际上您的预期模型不再具有ODE系统：您拥有具有非平滑组件的混合动力系统，即当轨迹x（t）在t = t_1时命中0时，它必须在所有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中设置所需的行为也相当容易，如果您在其帮助论坛中发布，我愿意为您提供建议。

- RHC