## Tuesday, November 13, 2012

### 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):

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)