scipy.integrate.ode.set_solout有效吗？

Question

scipy.integrate.ode这个接口提供了一种方法，可以在积分过程中如果有约束条件被违反时停止计算，这个方法叫做set_solout。不过，我在使用这个方法时遇到了问题，甚至在最简单的例子中也无法正常工作。以下是我尝试的一个例子：

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)

按照我的预期，当t = log(2) = 0.693...时，积分应该被solout停止，但实际上它却一直计算到t = 5，这时y = 0.007。

这是不是scipy的一个bug，还是我没有正确使用set_solout？

Answer 1

原来你需要先调用 set_solout，再调用 set_initial_value。我是通过研究 set_solout 的测试 在 scipy 的测试套件中发现的。所以，如果把这两个调用的顺序调换，代码就能得到正确的结果。

即使这个行为是正确的，也应该在 set_solout 的文档中提到这一点。我在 GitHub 上向 SciPy 提出了一个问题。

更新：这个问题在 SciPy 0.17.0 中已经修复；现在即使在 set_initial_value 之后调用 set_solout，也能正常工作，代码会得到正确的结果。