The 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] < 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 set_solout
correctly?
See Question&Answers more detail:
os 与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…