# [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.

More Related questions

We use cookies to improve the performance of the site. By staying on our site, you agree to the terms of use of cookies.