Coroutines for ODE solvers

Question

Are there any ODE solver packages that use coroutines and yield their results instead of functions and returning?

Briefly, a subroutine in a programming language does some computations, returns a result to whoever called it, and then it's done. A coroutine does some computations, yields a result to whoever called it, and can then be resumed where it last yielded a value. In Python, coroutines are used to implement generators.

For some scientific applications that involve solving systems of ODE, the logic about how long to integrate for can be very complex. For example, if you're writing some code to create a stream plot of a vector field, you'll continue creating an integral curve if (1) it isn't long enough yet and (2) it isn't too close to all of the other integral curves you've already calculated. The second stopping condition is much more complicated and involves some kind of spatial data structure like a KD-tree. Most ODE solving routines – for example SciPy’s solve_ivp function – present you with a very limited number of options for when to stop integrating and you have to decide what your criterion is at the outset.

By contrast, it's easy to write an adaptive integration scheme that yields the next point on the integral curve, the value of the vector field at that point, and the timestep it used. The ODE integration routine doesn't care about how long the caller wants it to go for. It could be part of a screensaver and it'll keep going until the user wiggles the mouse.

Answer 1

Not quite a package but the example at the end of the Revised6 Report on the Algorithmic Language Scheme implements a classical Runge-Kutta integrator as a function that returns 'an infinite stream of states'.

This is done in a few lines of code and is easy to port to for example Python using its coroutines, generators, itertools, &c.

I think that this is a compelling idea. Since first seeing it in one of the earlier revisions of the Scheme report towards the end of the last millenium, I have always used something based on this for integrating ODEs as well as systems of implicit ODEs (e.g. arising from finite element discretizations without 'mass lumping' and so having nondiagonal 'mass' matrices) and differential-algebraic equations which cannot be handled by standard ODE packages like SciPy referred to above. Thus far, this has all been in-house implementations, too ad hoc to publish.

I have several times over the last two decades looked for a package (particularly a Python package) implementing something like the R6RS approach to generating trajectories for dynamical systems (particularly more general than ODEs) but not yet found one. There was GarlicSim which was very nice but it had been abandoned by the time I found it. More recently, I thought that the stream interface of Scikit FiniteDiff, based on Streamz looking promising, particularly for still more complicated (e.g. coupled or hierarchical) problems.

Answer 2

You don't actually need co-routines for this -- you can achieve the same by just using regular callbacks into user code. For example, in the C programming language, the ODE solver routine would have an extra argument that denotes a pointer to a function. User code would call the integrator with a pointer-to-function for a user function that the ODE integrator then calls at the end of each time step (sending the current solution along as an argument). If the user function returns to the ODE integrator with a true value, then integration resumes. If the user function returns false, then the integrator terminates, releases memory, and returns to the caller of the ODE integrator.

In languages such as C++, you have more options: overloading base classes, or using boost-style signals in the ODE integrator object that user code can attach to. Point being, all of this is easily possible with traditional languages. It just inverts the logic of who drives whom compared to your example.

The example you show is interesting in that you want lots of integrators and some top-level logic that decides which of these integrators gets to go next. That's most easily done through the co-routines you mention, but not all that difficult either with callbacks. It would probably easiest to implement with callbacks if each integrator runs on its own thread, and the callbacks access some shared state that they modify whenever one of the threads calls a callback. That seems radically different from the co-routine approach, but it's worth remembering that co-routines typically are also implemented with one stack per co-routine -- not actually so different from separate threads. The only difference between co-routines and threads is really that only one co-routine can run at a time -- but that really just comes down to having one mutex per thread in the callback picture.

Answer 3

There are integrators which expose a more fine-grained control to the user, e.g., scipy.integrate.ode:

from scipy.integrate import ode

def f(t,y):
    return y[0]

initial = [1]
ODE = ode(f)
ODE.set_integrator("dopri5")
ODE.set_initial_value(initial)

for time in range(1,10):
    print(ODE.integrate(time))

Instead of coroutines, this provides the integrator as an instance (ODE) of a class (ode), which is equipped with a method that performs a small integration into the future (integrate) and returns the result.

I am not sure whether you gain anything from using coroutines here, but if you really wish, it is straightforward to wrap a coroutine around this:

from scipy.integrate import ode

def f(t,y):
    return y[0]

def integrator(f,initial):
    ODE = ode(f)
    ODE.set_integrator("dopri5")
    ODE.set_initial_value(initial)

    target_time = 0
    while True:
        target_time = yield ODE.integrate(target_time)

my_coroutine = integrator(f,[1.0])
next(my_coroutine)
for time in range(1,10):
    print(my_coroutine.send(time))