Dynamically ending ODE integration in SciPy

Question

I have a light ray moving through space-time, i.e. a curve in $\mathbb{R}^4$, parametrized by some variable λ. The exact trajectory, i.e. the coordinate functions $x^μ(λ)$ of the curve are given by some ODE $\frac{dx^μ}{dλ} = f(λ)$.

Now, I've got two requirements for my code to satisfy:

1. The results $x^μ(λ)$ must be given on a regular grid, i.e. the step size dλ must be the same everywhere.

2. I don't know the interval λ runs over beforehand, i.e. integration must end dynamically as soon as a certain condition $g(λ, x^μ(λ))$ is met (basically I'm testing whether the coordinates $x^μ(λ)$ are outside of some coordinate interval where I don't know the space-time's properties anymore).

Requirement 2 obviously excludes scipy.odeint(). This leaves me with the scipy.ode class which allows me to pass the break condition to the integrator via a solout callback if I choose the dopri5 or the dop853 algorithm.

It's a bit more complicated, though: Since these (and all other) integrators dynamically adjust their internal step size (which is kind of the point of not doing the integration manually) I cannot simply store the values being passed to the callback if I want to fulfill requirement 1. (Update: In fact, it doesn't make sense to do so, anyway. See my remark (*) below.) Instead, I still need to loop over λ and call the scipy.ode.integrate() method for each step $λ + dλ$. This, however, means that I must somehow stop the loop:

from scipy.integrate import ode

def deriv(λ, x, params):
    # …
    return dx_dλ

end_reached = False

def solout(λ, x):
    if g(λ, x):
        # Stop integration
        # Returning -1 will NOT make solver.successful() return False, so
        # we need to inform the integration loop (see below) separately.
        nonlocal end_reached
        end_reached = True
        return -1

    # Everything's ok. Move on.
    return 0

λ0 = 0
x0 = …
params = …
dλ = 0.1

solver = ode(deriv)
solver.set_integrator('dopri5', nsteps=500)
solver.set_solout(solout)
solver.set_initial_value(x0, λ0)
solver.set_f_params(params)

λ = [λ0]
x = [x0]

while True:
    solver.integrate(solver.t + dλ)

    if not solver.successful() or end_reached:
        break

    λ.append(solver.t)
    x.append(solver.y)

All in all, this is rather ugly. Not only is the logic basically all over the place (for instance I cannot test the condition in the loop where it would belong). But due to the solout callback (I suppose) it's also at least one order of magnitude slower than without one. Update: In fact, even without the callback, the dopri5 integrator is orders of magnitude slower than vode (at least for another ODE where I wanted to use this approach).

Am I missing something or is this really the only way to meet the above requirements? I also tried raising an exception in deriv() (and subsequently catching it in the loop) in order to stop the integrator. While this works for all integrators, it will still lead to a bunch of errors being printed to stdout.

(*) As far as I understand, the values passed to deriv() are not meant to be a reliable approximation but only the final value returned by the solver is. For instance, the solver might do overshooting to choose a sensible step size in view of (rapidly) changing slopes and to thereby arrive at a best estimate.

Answer 1

The feature that you demand is called event location in Matlab ODE solvers pack, or rootfinding in SUNDIALS solvers suite terminology. Essentially this feature allows to stop integration exactly at the point where some vector function of free and dependent variables has a root. Namely, for system $$ \frac{d\boldsymbol y}{dt} = f(t, \boldsymbol y), \quad \boldsymbol y \in \mathbb R^n, $$ an ODE solver with the ability to do event location will be able to find when some function $ g(t, \boldsymbol y)$ crosses zero—and this process of zero-crossing is called an event.

Unfortunately, scipy.integrate lacks good ODE solver with the support of event location. You can use solout feature of dopri5 solver, but you will overshoot in this case, because solout is called only after step is done. In contrast, when ODE solver is able to do event location, it will work hard on finding the exact location of the root of $g(t, \boldsymbol y)$, so that the integration will be stopped much more precisely.

I also need this feature for my computations. I use the following package https://github.com/bmcage/odes. This package is quite difficult to compile, because it depends on BLAS, LAPACK and SUNDIALS. However, interface (API) is quite good, and it will work much better than using solout feature of dopri5 from scipy.integrate.

Answer 2

The solution is https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

From the documentation : ‘RK45’ or ‘RK23’ method for non-stiff problems and ‘Radau’ or ‘BDF’ for stiff problems

The documentation taken from scipy:

scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, **options)[source] Solve an initial value problem for a system of ODEs.

Parameters: . . . . .

events: callable, or list of callables, optional Events to track. If None (default), no events will be tracked. Each event occurs at the zeros of a continuous function of time and state. Each function must have the signature event(t, y) and return a float. The solver will find an accurate value of t at which event(t, y(t)) = 0 using a root-finding algorithm. By default, all zeros will be found. The solver looks for a sign change over each step, so if multiple zero crossings occur within one step, events may be missed. Additionally each event function might have the following attributes:

terminal: bool, optional Whether to terminate integration if this event occurs. Implicitly False if not assigned.

direction: float, optional Direction of a zero crossing. If direction is positive, event will only trigger when going from negative to positive, and vice versa if direction is negative. If 0, then either direction will trigger event. Implicitly 0 if not assigned.

You can assign attributes like event.terminal = True to any function in Python.

Answer 3

I ended up raising an exception in the deriv() function (and catching it in the loop) because it turned out to be the only reliable way to end integration in every single case and for every ODE.

(While I considered this option already in my original question, it seems there is no other way. If anybody presents a better option I will gladly mark her/his answer as accepted. Update: I accepted Dmitry's answer.)

Answer 4

A simpler answer, IMHO, is to solve your system using PyDSTool. It has built-in support for accurate event finding. In fact, it's easier to set up and more efficient (including faster) than Matlab's approach if you use the C-based solvers. I am the author of this package.

A tutorial using events is given here.

Answer 5

Not sure if the solve_ivp package went through an update, but I just found a workaround that gives it the capability to detect events not at the root of the function. In the "limit" function, a limiting case can be described that depends on the value of the integrated function.

However, the solver expects the return to be a list the first time this is evaluated, and a numpy array every time after that. Typically, the input to the function is passed directly back to solve_ivp, and this is not a problem. However, if a custom condition is desired, that has to be accounted for. I did it with a global variable, as shown below. If the value "0" is returned from the function "limit", the solver will stop immediately.

 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.integrate import solve_ivp

 # Define derivative function 
 def f(t, y, c):
     dydt = (0.5*t) 
     return dydt

 # stop integration when y-1 = 0 (arbitrary condition)
 def limit(t,y):
     global flag
     if flag == 1:
         test = np.array(y)-np.array(1)
         test = np.ndarray.tolist(test)
         flag = 0
     else:
         test = np.array(y) - np.array(1)
     return test

 # Define terminal condition and type-change flag
 limit.terminal = True                                                                                                                                                                                   
 global flag
 flag = 1

 # Define time spans, initial values, and constants
 tspan = np.linspace(0.5, 3, 1000)
 yinit = [-1]
 c = []

 # Solve differential equation
 sol = solve_ivp(lambda t, y: f(t, y, c),[tspan[0], tspan[-1]], yinit, t_eval=tspan,events=limit)

 t = sol.t
 y = np.transpose(sol.y)
 plt.plot(t,y)
 plt.show()