将刚性 ODE 与 Python 集成
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/2088473/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Integrate stiff ODEs with Python
提问by Chinmay Kanchi
I'm looking for a good library that will integrate stiff ODEs in Python. The issue is, scipy's odeint gives me good solutions sometimes, but the slightest change in the initial conditions causes it to fall down and give up. The same problem is solved quite happily by MATLAB's stiff solvers (ode15s and ode23s), but I can't use it (even from Python, because none of the Python bindings for the MATLAB C API implement callbacks, and I need to pass a function to the ODE solver). I'm trying PyGSL, but it's horrendously complex. Any suggestions would be greatly appreciated.
我正在寻找一个很好的库，它将在 Python 中集成僵硬的 ODE。问题是，scipy 的 odeint有时会给我很好的解决方案，但初始条件的最轻微变化都会导致它倒下并放弃。MATLAB 的刚性求解器（ode15s 和 ode23s）很高兴地解决了同样的问题，但我不能使用它（即使是从 Python 中，因为 MATLAB C API 的 Python 绑定都没有实现回调，我需要传递一个函数到 ODE 求解器）。我正在尝试 PyGSL，但它非常复杂。任何建议将不胜感激。
EDIT: The specific problem I'm having with PyGSL is choosing the right step function. There are several of them, but no direct analogues to ode15s or ode23s (bdf formula and modified Rosenbrock if that makes sense). So what is a good step function to choose for a stiff system? I have to solve this system for a really long time to ensure that it reaches steady-state, and the GSL solvers either choose a miniscule time-step or one that's too large.
编辑：我在 PyGSL 上遇到的具体问题是选择正确的步进函数。其中有几个，但没有与 ode15s 或 ode23s 的直接类似物（bdf 公式和修改后的 Rosenbrock，如果有意义的话）。那么，为刚性系统选择什么是好的阶跃函数呢？我必须解决这个系统很长时间以确保它达到稳态，而 GSL 求解器要么选择一个很小的时间步长，要么选择一个太大的时间步长。
采纳答案by Mike Dunlavey
Python can call C. The industry standard is LSODEin ODEPACK. It is public-domain. You can download the C version. These solvers are extremely tricky, so it's best to use some well-tested code.
Python可以调用C。行业标准是ODEPACK中的LSODE。它是公共领域的。您可以下载C 版本。这些求解器非常棘手，因此最好使用一些经过良好测试的代码。
Added: Be sure you really have a stiff system, i.e. if the rates (eigenvalues) differ by more than 2 or 3 orders of magnitude. Also, if the system is stiff, but you are only looking for a steady-state solution, these solvers give you the option of solving some of the equations algebraically. Otherwise, a good Runge-Kutta solver like DVERKwill be a good, and much simpler, solution.
补充：确保你真的有一个刚性系统，即如果速率（特征值）相差超过 2 或 3 个数量级。此外，如果系统是刚性的，但您只是在寻找稳态解，这些求解器可为您提供代数求解某些方程的选项。否则，像DVERK这样优秀的 Runge-Kutta 求解器将是一个很好的、更简单的解决方案。
Added here because it would not fit in a comment: This is from the DLSODE header doc:
添加在这里是因为它不适合评论：这是来自 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).
Also, yes Michaelis-Menten kinetics is nonlinear. The Aitken acceleration works with it, though. (If you want a short explanation, first consider the simple case of Y being a scalar. You run the system to get 3 Y(T) points. Fit an exponential curve through them (simple algebra). Then set Y to the asymptote and repeat. Now just generalize to Y being a vector. Assume the 3 points are in a plane - it's OK if they're not.) Besides, unless you have a forcing function (like a constant IV drip), the MM elimination will decay away and the system will approach linearity. Hope that helps.
此外，是的 Michaelis-Menten 动力学是非线性的。不过，Aitken 加速器可以与之配合使用。（如果你想要一个简短的解释，首先考虑 Y 是标量的简单情况。你运行系统以获得 3 个 Y(T) 点。通过它们拟合指数曲线（简单代数）。然后将 Y 设置为渐近线和重复。现在只是概括为 Y 是一个向量。假设 3 个点在一个平面上 - 如果它们不在也没关系。）此外，除非你有一个强制函数（比如一个恒定的 IV 点滴），否则 MM 消除会衰减离开，系统将接近线性。希望有帮助。
回答by Olivier Verdier
If you can solve your problem with Matlab's ode15s, you should be able to solve it with the vodesolver of scipy. To simulate ode15s, I use the following settings:
如果您可以使用 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)
and then you can happily solve your problem with ode15s.integrate(t_final). It should work pretty well on a stiff problem.
然后你就可以愉快地用ode15s.integrate(t_final). 它应该可以很好地解决棘手的问题。
(See also http://www.scipy.org/NumPy_for_Matlab_Users)
回答by RHC
PyDSToolwraps the Radau solver, which is an excellent implicit stiff integrator. This has more setup than odeint, but a lot less than PyGSL. The greatest benefit is that your RHS function is specified as a string (typically, although you can build a system using symbolic manipulations) and is converted into C, so there are no slow python callbacks and the whole thing will be very fast.
PyDSTool封装了 Radau 求解器，它是一个优秀的隐式刚性积分器。这比 odeint 有更多的设置，但比 PyGSL 少很多。最大的好处是您的 RHS 函数被指定为字符串（通常，尽管您可以使用符号操作构建系统）并转换为 C，因此没有缓慢的 Python 回调，整个过程会非常快。
回答by YuppieNetworking
I am currently studying a bit of ODE and its solvers, so your question is very interesting to me...
我目前正在研究一些 ODE 及其求解器，所以你的问题对我来说很有趣......
From what I have heard and read, for stiff problems the right way to go is to choose an implicit method as a step function (correct me if I am wrong, I am still learning the misteries of ODE solvers). I cannot cite you where I read this, because I don't remember, but here is a threadfrom gsl-help where a similar question was asked.
根据我所听到和阅读的内容，对于僵硬的问题，正确的方法是选择隐式方法作为阶跃函数（如果我错了，请纠正我，我仍在学习 ODE 求解器的奥秘）。我不能在我读到的地方引用你，因为我不记得了，但这里有一个来自 gsl-help的线程，其中提出了类似的问题。
So, in short, seems like the bsimpmethod is worth taking a shot, although it requires a jacobian function. If you cannot calculate the Jacobian, I will try with rk2imp, rk4imp, or any of the gear methods.
因此，简而言之，该bsimp方法似乎值得一试，尽管它需要一个雅可比函数。如果你无法计算雅可比，我会尝试rk2imp，rk4imp或任何的齿轮的方法。
