Funktioniert scipy.integrate.ode.set_solout?

Question

Die Schnittstelle scipy.integrate.ode zu Integrationsroutinen bietet eine Methode zum Stoppen der Integration, wenn in einem Schritt, set_solout, eine Einschränkung verletzt wird. Allerdings kann ich diese Methode nicht einmal in den einfachsten Beispielen zum Laufen bringen. Hier ist ein Versuch:

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) 

Die Integration von solout eindämmen sollte, wenn t = log(2) = 0.693..., sondern weiterhin glücklich bis t = 5, wenn y = 0.007.

Ist das ein Fehler in scipy, oder verwende ich set_solout nicht richtig?

Answer 1

Es stellt sich heraus, Sie müssen set_soloutvor aufrufen set_initial_value anrufen. (Ich fand das heraus, indem ich die set_solout tests in der scipy Testsuite studierte.) Also, die Reihenfolge der zwei Anrufe in meinem Fragecode umzukehren erzeugt das korrekte Ergebnis. Wenn dieses Verhalten korrekt ist, sollte es in der Dokumentation für set_solout erwähnt werden. Ich habe an issue with SciPy on GitHub gepostet.

UPDATE: Dieses Problem ist in SciPy 0.17.0 behoben; set_solout funktioniert auch wenn nach set_initial_value aufgerufen wird, und der Fragecode wird das richtige Ergebnis erzeugen.