[SciPy-User] Speeding up odeint for big problems

Justin Bois justinbois@gmail....
Tue Sep 28 12:15:16 CDT 2010

Thanks Ron and Sebastian for your help.  The code to calculate the RHS of
the ODEs is really large so I can't really post it.  I am solving a
reaction-diffusion-advection system in 2D.  The velocity (which contributes
to the advection) is a complicated function of the concentrations of the
species in the reaction-diffusion system.  At each time step, I have to
solve a separate, time-independent ODE for the velocity, which is then
incorporated into the RHS.  Because this ODE happens to be linear, I can
solve it using a spectral method (for my boundary conditions of interest,
regardless of my differencing scheme in the reaction-diffusion-advection
system).  The presence of this velocity, manifest in the advection term,
makes the Jacobian dense, even if I use low-order finite difference for my
differencing.  Since I necessarily have a dense Jacobian, I typically use
pseudospectral methods for computing the derivatives to get high accuracy
with fewer grid points, though I have also implemented finite difference

I have pasted below the first few lines of output of a profile I did using
the magic function %prun with ipython -pylab.

      36447356 function calls in 1212.986 CPU seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1  903.674  903.674 1212.976 1212.976
   678300  157.372    0.000  239.023    0.000 basic.py:358(_raw_fftnd)
  1356600   24.225    0.000   35.504    0.000 basic.py:109(_fix_shape)

So, aside from odeint itself, most of the computational effort is spent
computing derivatives and solving for the velocity, all of which is just
using 2D FFTs.  Curiously, odeint spends a lot of time doing something that
does not require evaluation of the RHS of the ODEs.  As I watch the
integration proceed, it takes several time steps very rapidly, and then
pauses at a particular time point, presumably to adjust the step size.  My
guess is that Sebastian is right: this is slow because the Jacobian is large
and dense.  Do you think just doing an explicit RK4 time stepping or
something might help?  Or are there tricks we can play to get odeint to roll
along more quickly?  I am passing an object with all sorts of attributes as
an argument to the RHS function, though only some are used.  I am assuming
the RHS function is all executed in Python and then the return values are
all that undergo type conversion, right?  I'm guessing it's not so much type
conversion that's the big slowdown, or do you think I'm mistaken?
odeinthas been a nice little workhorse in my work so far!
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20100928/28161f83/attachment.html 

More information about the SciPy-User mailing list