[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
rk4imp, or any of the gear methods.
More Related questions