2

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 import ode
y0, t0 = [1.0j, 2.0], 0
def f(t, y, arg1):
    return [1j*arg1*y[0] + y[1], -arg1*y[1]**2]
def jac(t, y, arg1):
    return [[1j*arg1, 1], [0, -arg1*2*y[1]]]

# The integration:
r = ode(f, jac).set_integrator('zvode', method='bdf', with_jacobian=True)
r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0)
t1 = 10
dt = 1
while r.successful() and r.t < t1:
    r.integrate(r.t+dt)
    print("%g %g" % (r.t, r.y))

My problem is: it is using a lot of the python loops (the while-loop) which essentially slows the program down. I can try to write a C code and use ctypes to make it faster, but I won't be able to access the nice algorithms like dopri5 in scipy (unless I implement it myself).

Is there any vectorized way of coding to make this faster?

Thank you!

Yuxiang Wang
  • 8,095
  • 13
  • 64
  • 95

2 Answers2

4

'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.

hpaulj
  • 221,503
  • 14
  • 230
  • 353
3

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).

subhacom
  • 868
  • 10
  • 24
  • Never too late man! This is what I ended up using, and others looking at this post in the future may find it helpful. Thanks a lot! – Yuxiang Wang Jun 21 '14 at 05:18