I have a few hundred lines of code that calculate a system of ODEs. The resulting system of several hundred to several thousand equations take a long time to integrate. (I'm using SciPy's integrate interface; testing on a small case suggested it's several times faster than GSL's ode_solver for my problem.)
Of course, most of the time is spent in evaluating my equations. I'm already using fast_callable to speed up the calculations. It made a wonderful difference. But it's still taking hours or even days for the larger systems. I want to put this integration inside an optimizer, so any speed gain is great.
Stealing the example from the reference manual (http://www.sagemath.org/doc/numerical_sage/scipy.html), I'm currently doing something like the following.
import scipy
from scipy import integrate
var('x, y')
mu = 10.0
dy = (y, -x + mu*y*(1-x**2))
dy_fc = tuple(fast_callable(expr, domain=float, vars=(x,y)) for expr in dy)
def f_1(y,t):
return [f(*y) for f in dy_fc]
xx=scipy.arange(0,100,.1)
yy=integrate.odeint(f_1,[1,0],xx)
I don't think I can speed up the integrate routine. Can I do anything to speed up f_1
?
Thanks!