我正在寻找一个可以在 Python 中集成僵硬 ODE 的好库。问题是，scipy 的 odeint有时会给我很好的解决方案，但初始条件的最轻微变化会导致它倒下并放弃。MATLAB 的刚性求解器（ode15s 和 ode23s）非常愉快地解决了同样的问题，但我不能使用它（即使是 Python，因为 MATLAB C API 的 Python 绑定都没有实现回调，我需要传递一个函数到 ODE 求解器）。我正在尝试 PyGSL，但它非常复杂。任何建议将不胜感激。
编辑：我在使用 PyGSL 时遇到的具体问题是选择正确的阶跃函数。其中有几个，但没有直接与 ode15s 或 ode23s 类似（bdf 公式和修改后的 Rosenbrock，如果有意义的话）。那么对于刚性系统来说，什么是好的阶跃函数呢？我必须对这个系统求解很长时间才能确保它达到稳态，而 GSL 求解器要么选择一个很小的时间步长，要么选择一个太大的时间步长。
如果你能用 Matlab 解决你的问题ode15s，你应该可以用vodescipy 的求解器来解决它。为了模拟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)。它应该可以很好地解决一个棘手的问题。
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 动力学是非线性的。不过，Aitken 加速可以与之配合使用。（如果你想要一个简短的解释，首先考虑 Y 是一个标量的简单情况。你运行系统得到 3 个 Y(T) 点。通过它们拟合指数曲线（简单代数）。然后将 Y 设置为渐近线和重复。现在只是概括为 Y 是一个向量。假设这 3 个点在一个平面上 - 如果它们不是也可以。）此外，除非你有一个强制函数（如恒定的 IV 滴注），否则 MM 消除将衰减离开，系统将接近线性。希望有帮助。
PyDSTool封装了 Radau 求解器，这是一个出色的隐式刚性积分器。这比 odeint 有更多的设置，但比 PyGSL 少得多。最大的好处是您的 RHS 函数被指定为字符串（通常，尽管您可以使用符号操作构建系统）并转换为 C，因此没有慢的 python 回调，整个事情会非常快。
我目前正在研究一点 ODE 及其求解器，所以你的问题对我来说很有趣......
根据我所听到和阅读的内容，对于棘手的问题，正确的方法是选择隐式方法作为阶跃函数（如果我错了，请纠正我，我仍在学习 ODE 求解器的奥秘）。我不能引用你在哪里读到这篇文章，因为我不记得了，但这里有一个来自 gsl-help 的线程，其中提出了类似的问题。
因此，简而言之，该方法似乎bsimp值得一试，尽管它需要一个雅可比函数。如果您无法计算雅可比行列式，我将尝试使用rk2imp、rk4imp或任何齿轮方法。