Python：scipysolve_ivp 和事件的问题

Question

我的代码要复杂得多，但我可以通过以下示例重现该问题：

import numpy as np
from scipy.integrate import solve_ivp


def funct(t, y):
    return -np.sqrt(y)


def event(t, y):
    return y[0]-0.1


if __name__ == '__main__':
    event.terminal = True
    sol = solve_ivp(funct, t_span=[0, 3], y0=[1], events=event)

正如人们所预料的那样，如果

y < 0

，则会发出警告，但对于这个特定示例，仍然找到了解决方案。但是，我的代码停止工作，因为我使用的是第三方包，当发生类似情况时，该包会引发错误。

现在，我知道

event

是通过迭代方法计算的，所以我的问题是：是否可以避免在迭代过程中，

y

的值低于给定阈值并避免错误/警告？那么，类似“如果

y < 0

则拒绝该步骤”或类似的内容？

Answer 1

IIUC，这里的问题是

funct

在连续迭代中被

y=[0.14553212]

（事件之前）和

y=[-0.01227461]

（事件之后）调用。在此示例中，当

funct

为负数时，

y

返回 NaN，并且

solve_ivp

能够恢复、找到事件、终止并返回直到事件发生为止的解决方案。然而，在您的代码中，当使用无效输入调用时

funct

会引发错误，因此

solve_ivp

无法恢复。

如何包装可调用对象，检测会导致错误的条件，并返回一个合理的值（不会改变解决方案）而不是调用函数。

def wrapped(t, y):
    # alternatively, `try` and `except` if `funct` raises an error
    if y[0] <= 0:
        return [0]  # or maybe NaN; since that works in this example
    return funct(t, y)

sol = solve_ivp(wrapped, t_span=[0, 3], y0=[1], events=event)
sol  # no warning is raised

#  message: A termination event occurred.
#  success: True
#   status: 1
#        t: [ 0.000e+00  1.000e-01  9.784e-01  1.368e+00]
#        y: [[ 1.000e+00  9.025e-01  2.610e-01  1.000e-01]]
#      sol: None
# t_events: [array([ 1.368e+00])]
# y_events: [array([[ 1.000e-01]])]
#     nfev: 32
#     njev: 0
#      nlu: 0