使用Python集成僵硬的ODE

Question

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.

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.

Answer 1

如果你可以用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).它应该在一个僵硬的问题上运作得很好.

(另见http://www.scipy.org/NumPy_for_Matlab_Users)

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动力学是非线性的.不过,Aitken加速可以使用它.(如果你想要一个简短的解释,首先考虑Y是标量的简单情况.你运行系统得到3个Y(T)点.通过它们拟合指数曲线(简单代数).然后将Y设置为渐近线和重复.现在只是概括为Y是一个向量.假设3个点在一个平面上 - 如果它们不是那就没关系.)此外,除非你有一个强制函数(如恒定的静脉滴注),MM消除会衰减离开,系统将接近线性.希望有所帮助.