Vectorized Scipy Ode Solver

My question is with respect to the current scipy ode solver. From the scipy doc page, their usage is: # A problem to integrate and the corresponding jacobian: from scipy.integrate

Solution 1:

'vectorization' means doing a bunch of calculations in parallel, all at once. Yes, the detailed implementation will involve iteration, but it's in C and the order does not matter to you, the Python programmer.

But an ode solution like this is essentially a serial operation. You have to solve the problem at time t before you solve it at time t+dt. You can't vectorize the solution through time. The best you can do is choose an ode solver that makes intelligent choices for the time steps (dt), big steps where possible, small ones when needed to capture rapid changes.

A good ode solver lets you vectorize the spatial dimension - i.e. solving 10 odes in parallel. You might also be able to vectorize calculating the Jacobian, returning both y+dy and y-dy at once. Basically you want to make the calculation of f and jac as fast as possible.

Solution 2:

This is a bit late, but I think what you need is scipy.integrate.odeint. It takes an array of time points and computes the solution at each of them. This makes the loop to be computed outside Python (in FORTRAN).