# Adaptive ODE algorithm in Python

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?

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.

• Thank you for your answer. I will try to look more at this documentation, but I must admit I have never heard of this module. Is there a reason why scipy does not have this functionality? It would seem this is a pretty standard requirement for most problems of practical interest. – Dipole Apr 4 '14 at 18:42
• Lack of developer time, mainly. SUNDIALS is a big suite of solvers to wrap. It is BSD licensed, so it is eligible for inclusion into SciPy, and I think it would be a good fit. There's been a lot of discussion on scipy-dev about it, but the volunteer work hasn't been there. Benny Malengier works on scikits-odes and does a good job. Johan Akesson (and others) do a good job with Assimulo, too, but the license terms are a bit weird; it's LGPL, but with some custom terms, so it's not an OSI-approved license. – Geoff Oxberry Apr 4 '14 at 20:30
• I think I have overcomplicated things and therefore I may have mislead you as to what I want. The second part of your answer is essentially what I want: I want to integrate my system between start state and a specified end state. I want the integrator to choose internal time steps automatically, depending on how quickly my system is changing at a given time. However I want the integrator to choose the largest internal step sizes that can be spared, such that any value of the state I query at intermediate times will be accurate enough (to a given tolerance I guess) – Dipole Apr 5 '14 at 11:59
• @Jack: That is exactly what I stated in my answer. – Geoff Oxberry Apr 5 '14 at 19:15

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.

• PyDSTool looks interesting (and the code is readable, which is a big plus). I'm not sure what you mean about SUNDIALS being "legacy". Its last bugfix release was in August, and the last major features added (e.g., interfaces to sparse direct linear solvers, ARKIMEX methods) were in March of this year. I'm also not sure what you mean here by "arbitrary precision"; at least in the event detection code I've looked at in PyDSTool, I don't see any use of arbitrary precision arithmetic in Python or C, so presumably, you're talking about a user-defined tolerance? – Geoff Oxberry Nov 4 '15 at 7:18
• Yes, I meant "arbitrary" up to the limits of floating point precision based on a user tolerance, as opposed to having to accept the last point from a fixed-step solver before the event threshold is crossed. – RHC Nov 4 '15 at 15:54
• Also, I should have been clearer regarding "legacy". I am referring to a style of interface to a legacy solver that does not add high level capabilities like named fields in the output data rather than a raw array, convenient ways to define vector fields, etc. Their python interface simply exposes the legacy solvers directly. We use legacy solvers too (i.e. classic solvers written in C and Fortran, not brand new implementations) but we wrap them substantially with high-level features. – RHC Nov 4 '15 at 15:55