scipy.integrate.ode interface to integration routines provides a method for stopping the integration if a constraint is violated at any step,
set_solout. However, I cannot get this method to work, even in the simplest examples. Here’s one attempt:
import numpy as np from scipy.integrate import ode def f(t, y): """Exponential decay.""" return -y def solout(t, y): if y < 0.5: return -1 else: return 0 y_initial = 1 t_initial = 0 r = ode(f).set_integrator('dopri5') # Integrator that supports solout r.set_initial_value(y_initial, t_initial) r.set_solout(solout) # Integrate until t = 5, but stop when solout constraint violated r.integrate(5) # The time when solout should have terminated integration: intersection_time = np.log(2)
The integration should have been stopped by solout when
t = log(2) = 0.693..., but instead happily continues until
t = 5, when
y = 0.007.
Is this a bug in
scipy, or am I not using
It turns out you need to call
set_solout before calling
set_initial_value. (I figured this out by studying the
set_solout tests in the
scipy test suite.) So, reversing the order of the two calls in my question code produces the correct result.
Even if this behavior is correct, it ought to be mentioned in the documentation for
set_solout. I’ve posted an issue with SciPy on GitHub.
UPDATE: This issue is fixed in SciPy 0.17.0;
set_solout will work even if called after
set_initial_value, and the question code will produce the correct result.