Adaptive ODE algorithm in Python

Question

I want to integrate a particle path in 2D using the integrate.ode module. Things that are a bit different in my case are that, I only want to integrate up to a certain position, determined by the maximum allowed x coordinate of the particle: x_max.

The main issue I have is that the particle may first move very slowly and then gather more speed later on. Hence I don't want to waste effort with small time steps in this region. The algorithm should be able to adjust such that smaller time steps are used when the particle velocity becomes high.

In the end result if I plot the particle trajectory in the "phase-space" I should have a smooth line.

I have some rough pseudo-code below for this purpose:

backend = "dopri5"
x_max = 1
solver = ode(f)
solver.set_integrator(backend)  
solver.set_initial_value(y0, t0)
t, y = [t0], [y0]
k = 1.2
while solver.successful() and solver.y[0] < x_max
   solver.integrate(solver.t+dt)
   t.append(solver.t)
   y.append(solver.y)
   v_current = numpy.linalg.norm(y[-1])
   v_previous = numpy.linalg.norm(y[-2])
   if numpy.abs( v_current-v_previous ) > k * v_previous:
       dt = 0.8*dt
       del y[-1]
   else:
       dt = dt*1.2

Trouble is this algorithm may not be that robust, as choosing the values k, 1.2, 0.8 is somewhat arbitrary and may cause some stability issues with the algorithm.

EDIT: I also want to be able to plot points that are equally spaced in time on the trajectory, to give an indication of how the speed of the particle changes.

Can anyone suggest a better way of doing this?

Answer 1

Yes, the event detection capability in SUNDIALS should work for you. I do not know if it is implemented in the unofficial Python interfaces to SUNDIALS; you should look at Assimulo, scikits-odes, CasADi, and PySUNDIALS, roughly in that order.

As for equally spaced times, you should be able to specify times to which you want to integrate the system in your integrator; the integrator should store the internal state of its time steps, and interpolate its current solution array to give you solutions at the times you query.

Answer 2

SUNDIALS is OK, IMHO, but as an external library of legacy components it does not have the front-end interface capabilities that a python-centric environment does. I built the PyDSTool system to do this kind of work. It has a high-level set of constructions for dynamical systems and modeling beyond just the solvers. This includes a sophisticated and robust way to specify stopping points that are state or time dependent (called "events"). The events are solved for internally, automatically, and to arbitrary precision. You can then sample the solution very naturally according to any mesh you desire for plotting purposes. See this example for instance, and the other tutorials and supplied demos to see how to sample and plot with it.

P.S. PyDSTool can auto-generate C code to generate much faster-running solvers than any integrator that uses callbacks to python-defined functions, which is pretty much any other solution out there.