使用Python求解刚性常微分方程
我在找一个好的库，能在Python中处理刚性常微分方程（ODEs）。我的问题是，虽然scipy的odeint有时候能给我不错的结果，但只要初始条件稍微一变，它就会出问题，完全不管用。而MATLAB的刚性求解器（ode15s和ode23s）能很好地解决这个问题，但我不能使用MATLAB（即使是通过Python，因为没有任何Python的MATLAB C API绑定实现回调功能，而我需要把一个函数传给ODE求解器）。我在尝试使用PyGSL，但它复杂得令人头疼。任何建议都非常感谢。
补充说明：我在使用PyGSL时遇到的具体问题是选择合适的步长函数。虽然有好几个步长函数，但没有直接对应于ode15s或ode23s的（比如bdf公式和修改后的Rosenbrock，如果这样说有道理的话）。那么，对于一个刚性系统，选择哪个步长函数比较好呢？我需要长时间求解这个系统，以确保它达到稳态，而GSL的求解器要么选择的时间步长太小，要么太大。
4 个回答
PyDSTool 是一个工具，它使用了 Radau 解算器，这是一种非常优秀的隐式刚性积分器。虽然它的设置比 odeint 要多一些，但比 PyGSL 要简单很多。最大的好处是，你的右侧函数可以用字符串来指定（通常是这样，当然你也可以通过符号运算来构建系统），然后这个字符串会被转换成 C 语言代码，这样就没有慢慢的 Python 回调了，整个过程会非常快。
如果你能用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) 来解决你的问题了。这在处理比较复杂的问题时效果应该不错。
（也可以参考这个 链接）
Python可以调用C语言。行业标准是LSODE，这是在ODEPACK中的一个工具。它是公共领域的，你可以下载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).
另外，是的，米哈伊利斯-门腾动力学是非线性的。不过，艾特肯加速法可以用在上面。（如果你想要一个简单的解释，先考虑Y是一个标量的简单情况。你运行系统得到3个Y(T)的点。然后用简单的代数方法拟合一条指数曲线。接着把Y设置为渐近线，再重复这个过程。现在只需将Y推广为一个向量。假设这3个点在一个平面上——即使它们不在同一平面也没关系。）此外，除非你有一个强制函数（比如恒定的静脉滴注），否则米哈伊利斯-门腾的消除效应会逐渐消失，系统会趋向线性。希望这对你有帮助。