将刚性 ODE 与 Python 集成

Question

我正在寻找一个能够在 Python 中集成严格 ODE 的优秀库。问题是，scipy 的 odeint 有时会给我很好的解决方案，但是初始条件的最轻微的变化都会导致它崩溃并放弃。 MATLAB 的刚性求解器（ode15s 和 ode23s）很好地解决了同样的问题，但我无法使用它（即使是在 Python 中，因为 MATLAB C API 的 Python 绑定都没有实现回调，并且我需要传递一个函数到 ODE 求解器）。我正在尝试 PyGSL，但它非常复杂。任何建议将不胜感激。

编辑：我在使用 PyGSL 时遇到的具体问题是选择正确的阶跃函数。其中有几个，但没有与 ode15s 或 ode23s 直接类似的东西（bdf 公式和修改后的 Rosenbrock，如果有意义的话）。那么对于刚性系统来说，选择什么好的阶跃函数呢？我必须对这个系统求解很长时间才能确保它达到稳态，而 GSL 求解器要么选择很小的时间步长，要么选择太大的时间步长。

Answer 1

如果您可以使用Matlab的ode15s解决您的问题，那么您应该也可以使用scipy的vode求解器解决它。为了模拟 ode15s，我使用了以下设置：

ode15s = scipy.integrate.ode(f)
ode15s.set_integrator('vode', method='bdf', order=15, nsteps=3000)
ode15s.set_initial_value(u0, t0)

然后您就可以愉快地使用 ode15s.integrate(t_final) 解决您的问题。它应该可以很好地解决棘手的问题。

（另请参阅链接）

Answer 2

Python可以调用C。行业标准是ODEPACK中的LSODE。它是公共领域的。您可以下载C 版本。这些求解器非常棘手，因此最好使用一些经过良好测试的代码。

补充：确保您确实拥有一个刚性系统，即速率（特征值）是否相差超过 2 或 3 个数量级。此外，如果系统是刚性的，但您只想寻找稳态解，这些求解器可以让您选择以代数方式求解某些方程。否则，像 DVERK 这样的优秀 Runge-Kutta 求解器将是一个很好且更简单的解决方案。

添加到此处是因为它不适合注释：这是来自 DLSODE 标题文档：

C     T     :INOUT  Value of the independent variable.  On return it
C                   will be the current value of t (normally TOUT).
C
C     TOUT  :IN     Next point where output is desired (.NE. T).

此外，是的 Michaelis-Menten 动力学是非线性的。不过，艾特肯加速可以与之配合。 （如果您想要简短的解释，请首先考虑 Y 为标量的简单情况。运行系统以获取 3 个 Y(T) 点。通过它们拟合指数曲线（简单代数）。然后将 Y 设置为渐近线并重复。现在假设 3 个点在一个平面上，也没关系。）此外，除非有强制函数（如恒定的 IV 滴注），否则 MM 消除将会衰减。远离，系统将接近线性。希望有帮助。

Answer 3

PyDSTool 封装了 Radau 求解器，它是一个优秀的隐式刚性积分器。它比 odeint 有更多的设置，但比 PyGSL 少很多。最大的好处是，您的 RHS 函数被指定为字符串（通常，尽管您可以使用符号操作构建系统）并转换为 C，因此没有缓慢的 python 回调，整个过程将非常快。

Answer 4

我目前正在研究一些 ODE 及其求解器，所以你的问题对我来说非常有趣......

根据我所听到和读到的，对于刚性问题，正确的方法是选择隐式方法作为阶跃函数（如果我错了请纠正我，我仍在学习 ODE 求解器的奥秘）。我无法引用您在哪里读到此内容，因为我不记得了，但这里有一个 线程 提出了类似的问题。

因此，简而言之，bsimp 方法似乎值得一试，尽管它需要雅可比函数。如果您无法计算雅可比行列式，我将尝试使用 rk2imp、rk4imp 或任何 gear 方法。