将刚性 ODE 与 Python 集成答案

【问题标题】：Integrate stiff ODEs with Python将刚性 ODE 与 Python 集成
【发布时间】：2011-01-06 12:07:00
【问题描述】：

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

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

【问题讨论】：

我想帮助你处理 PyGSL。我从未使用过它，但我有使用 GSL 的经验。我刚刚查看了 pygsl (odeiv.py) 中提供的示例，它看起来与 C 中的几乎相同。您是否认为 PyGSL 是由于 python 接口或 GSL 本身而非常复杂？
好吧，可怕的复杂也许是夸大其词:)。它比 MATLAB 或 scipy 复杂一个数量级。需要澄清的是，python 接口与 C 接口几乎相同，因此库本身很复杂。另外，PyGSL 没有记录 odeiv，所以我必须使用 C 文档来弄清楚在 Python 中要做什么。不好玩。
OT：我不记得 t0:tn 在 Matlab 中做了什么，但 numpy.r_[0:1:11j] 可能会有所帮助。

标签： python scipy pygsl

【解决方案1】：

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

根据我所听到和阅读的内容，对于棘手的问题，正确的方法是选择隐式方法作为阶跃函数（如果我错了，请纠正我，我仍在学习 ODE 求解器的奥秘）。我无法引用你在哪里读到这篇文章，因为我不记得了，但这里有一个来自 gsl-help 的thread，其中提出了类似的问题。

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

【讨论】：

【解决方案2】：

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

添加：确保您确实有一个僵硬的系统，即，如果速率（特征值）差异超过 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 消除将衰减离开，系统将接近线性。希望对您有所帮助。

【讨论】：

感谢 LSODE 提示。我发现有人已经在cs.uiuc.edu/homes/mrgates2/ode 上为其编写了 Python 接口。明天会尝试一下，如果可行，请接受您的回答。
我知道我确实有一个僵硬的系统，但查看 LSODE 源代码，我意识到我的程序中可能存在错误。 scipy.integrate.odeint 求解器基于 LSODA，并且与所有 LS* 求解器一样，它需要一个时间点向量来计算解。事实证明，我计算的这个时间向量太粗糙了。我从 MATLAB 切换到 Python，并认为 numpy.arange 的工作方式与 MATLAB 中的 i=t0:tn 类似。它没有。所以我一直只解决整数次......
@Chinmay：需要一个时间点向量？根据我的经验，您从第一次运行到第二次，此时您可以根据需要对状态进行不连续的更改，然后运行到第三次，依此类推。这就是我们进行所有药物计量建模的方式。顺便说一句，我忘了提到如果您的 ODE 系统是线性的（常数雅可比行列式），您可以使用 DGPADM - 它速度很快，而且从来没有刚度问题。此外，如果系统是线性的，您可以通过代数方式求解稳态解。
...我忘了说，因为你说你在寻找稳态，我使用了一种叫做 Aitken Acceleration 的方法（系统不必是线性的）。基本上，这个想法是您假设它将以某种大致指数的方式收敛到渐近线。所以你在轨迹上取 3 个点，预测渐近线的位置，将你的系统设置在那里，然后重复。它收敛得非常快。
我担心系统是非常非线性的（Michaelis-Menten 型方程经过修改后会使它们更多非线性）。将研究 Aitken 加速。来自 odeint 文档：t : array A sequence of time points for which to solve for y. The initial value point should be the first element of this sequence.

【解决方案3】：

如果你能用 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) 愉快地解决您的问题了。它应该可以很好地解决棘手的问题。

（另见http://www.scipy.org/NumPy_for_Matlab_Users）

【讨论】：

不错！感谢那。最后，我使用vode 得出了或多或少类似的解决方案。
bdf 的最大订单数为 5。无论如何，设置订单数没有必要高于 12。因为 12 是 Adams 的最大值。

【解决方案4】：

PyDSTool 封装了 Radau 求解器，这是一个出色的隐式刚性积分器。这比 odeint 有更多的设置，但比 PyGSL 少得多。最大的好处是你的 RHS 函数被指定为一个字符串（通常，虽然你可以使用符号操作来构建一个系统）并被转换成 C，所以没有慢的 python 回调，整个事情会非常快。

【讨论】：