Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
316 views
in Technique[技术] by (71.8m points)

python - Does scipy.integrate.ode.set_solout work?

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

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)

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.


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...