[ACCEPTED]-Integrate stiff ODEs with Python-pygsl

If you can solve your problem with Matlab's 6 ode15s, you should be able to solve it with the 5 vode solver of scipy. To simulate ode15s, I use the 4 following settings:

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 3 solve your problem with ode15s.integrate(t_final). It should work 2 pretty well on a stiff problem.

(See also 1 http://www.scipy.org/NumPy_for_Matlab_Users)

Python can call C. The industry standard 27 is LSODE in ODEPACK. It is public-domain. You 26 can download the C version. These solvers are extremely 25 tricky, so it's best to use some well-tested 24 code.

Added: Be sure you really have a stiff 23 system, i.e. if the rates (eigenvalues) differ 22 by more than 2 or 3 orders of magnitude. Also, if 21 the system is stiff, but you are only looking 20 for a steady-state solution, these solvers 19 give you the option of solving some of the 18 equations algebraically. Otherwise, a good 17 Runge-Kutta solver like DVERK will be a good, and 16 much simpler, solution.

Added here because 15 it would not fit in a comment: This is from 14 the DLSODE header doc:

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 13 kinetics is nonlinear. The Aitken acceleration 12 works with it, though. (If you want a short 11 explanation, first consider the simple case 10 of Y being a scalar. You run the system 9 to get 3 Y(T) points. Fit an exponential 8 curve through them (simple algebra). Then 7 set Y to the asymptote and repeat. Now just 6 generalize to Y being a vector. Assume the 5 3 points are in a plane - it's OK if they're 4 not.) Besides, unless you have a forcing 3 function (like a constant IV drip), the 2 MM elimination will decay away and the system 1 will approach linearity. Hope that helps.

PyDSTool wraps the Radau solver, which is an excellent 8 implicit stiff integrator. This has more 7 setup than odeint, but a lot less than PyGSL. The 6 greatest benefit is that your RHS function 5 is specified as a string (typically, although 4 you can build a system using symbolic manipulations) and 3 is converted into C, so there are no slow 2 python callbacks and the whole thing will 1 be very fast.

I am currently studying a bit of ODE and 13 its solvers, so your question is very interesting 12 to me...

From what I have heard and read, for 11 stiff problems the right way to go is to 10 choose an implicit method as a step function 9 (correct me if I am wrong, I am still learning 8 the misteries of ODE solvers). I cannot 7 cite you where I read this, because I don't 6 remember, but here is a thread from gsl-help where 5 a similar question was asked.

So, in short, seems 4 like the bsimp method is worth taking a shot, although 3 it requires a jacobian function. If you 2 cannot calculate the Jacobian, I will try 1 with rk2imp, rk4imp, or any of the gear methods.