Adaptive Time Step Integrators in SciPy

For one of my projects at school I have to numericaly integrate a bunch of ODEs, something I would normally do using ode45 in Matlab or NDSolve in Mathematica. However, since the rest of my code is written in Python, and I really don't want to install the bloatware that is Matlab on my clean linux machine, I've decided I was going to try out the ODE libraries that come with numpy and scipy. As it turns out, they're not too bad, but their documentation is. So, it took me a while to figure out how to do this, and I'm guessing someone else will probably benefit from this post (I know it would have saved me about a day's work, at least).

So, it turns out that scipy doesn't really implement its own ode solver, but instead wraps a bunch of FORTRAN code, and just silently drops like 80% of the functionality of that code without documenting much of what it kept and what it ignored. Anyway, after sifting through a bunch of forums and a whole bunch of source code, I came to the conclusion that the only function that can actually do adaptive time stepping and supply dense output is the wrapper around the vode integrator. The other thing that I needed was to detect when my trajectory crosses a certain boundary defined by $\phi(y) = 0$ and trigger an event at that point. In my case I simply needed to restart integration after adding the crossing point to the solution array. The following code does what I needed:

def f(t, y):
    # define your \dot{y}=f(y,t) here

def phi(y):
    # this function defines the boundary that is crossed


from scipy.integrate import ode 

# set up the solver, with a maximum time step of 1e-1
# an initial condition of y0 = init and t0 = 0
solver = ode(f)
solver.set_integrator('vode', max_step=1e-1)
solver.set_initial_value(init, 0)
tfinal = 30  # run this for 30 seconds
t = [0] 
y = [init]

# call the solver iteratively with the 'step' option
# what this does is it returns the result of the integration
# after the next (variable) time step, and the length of that step
while solver.successful() and solver.t < tfinal :
    solver.integrate(tfinal, step=True)
    y.append(solver.y);
    t.append(solver.t);

    # check for zero crossing
    gg1 = phi(y[-1])
    gg0 = phi(y[-2]) 
    if gg0*gg1 < 0 : 
        # find zero crossing
        # simple linear interpolation will do for now
        tcross = t[-2] - gg0*(t[-1]-t[-2])/(gg1-gg0)
        ycross = y[-2] - gg0*(y[-1]-y[-2])/(gg1-gg0)
        
        # replace wrong y and t
        y[-1] = ycross
        t[-1] = tcross

        # reset integration and start again
        solver.set_initial_value(ycross, tcross)
        print("crossed the boundary at time ", tcross)