scipy.integrate.solve_ivp能否拒绝步长以避免用无效状态评估RHS?

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)

在这个例子中，求解器在计算时超出了事件的范围，导致调用了 funct 时传入了 y < 0 的值。正如大家所预料的那样，funct 会发出警告，并返回一个无效值（NaN）。在这个特定的例子中，solve_ivp 还是能恢复过来，找到一个解决方案。然而，我真正的 funct 使用了一个第三方的包，当 funct 在无效参数下被调用时会抛出错误，这样 solve_ivp 就会失败。

现在，我知道 event 是通过一种迭代的方法来计算的，所以我想问的是：solve_ivp 是否有可能检测到 y 的值是否低于某个阈值，从而避免出现错误或警告？也就是说，能不能像“如果 y < 0 就拒绝这一步”这样处理？

Answer 1

在处理关于如何在使用scipy.integrate模块中的solve_ivp函数时，避免在事件处理过程中出现负值 � y的问题时，我们需要找到一些方法来防止这些负值出现，以免引发错误或警告。通过讨论一些策略来绕过这个问题，我们希望能确保依赖于solve_ivp的第三方包能够顺利运行，同时增强涉及微分方程的数值计算的稳定性。

一个有效的办法是使用自定义事件处理器，当 � y从负值接近零时，停止迭代。通过设置一个机制，当y低于某个特定的阈值时拒绝继续计算，我们可以避免因负值而导致的错误或警告。接下来，我们将深入探讨如何实现这一点，并通过示例和见解来展示它的好处。

> 
> 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 def custom_event(t, y):
>     return y[0] + 1e-6  # Threshold to avoid negative values if __name__ == '__main__':
>     custom_event.terminal = True
>     sol = solve_ivp(funct, t_span=[0, 3], y0=[1], events=custom_event)

在这个修改后的代码版本中，我们引入了一个custom_event函数，作为一个修改过的事件处理器。我们不是直接将y与零进行比较，而是设置一个小的正阈值（在这个例子中是 1×10−6 1×10−6 ），确保当y从负值接近零时，迭代会停止，从而有效地防止负值的出现。通过调整这个阈值，我们可以根据具体需求来微调事件检测机制的灵敏度。

使用自定义事件处理器有几个好处。首先，它让我们对事件检测过程有更大的控制权，可以针对特定情况（比如避免y出现负值）制定解决方案。这种灵活性增强了数值计算的稳健性，减少了在迭代过程中由于意外行为而产生的潜在错误或警告。

此外，自定义事件处理器还使用户能够将特定领域的知识融入事件检测机制，从而在迭代过程中做出更智能的决策。通过根据问题领域的特征设置合适的阈值或条件，用户可以确保数值解的可靠性和准确性，同时避免出现不希望的结果，比如负值。

而且，通过主动解决在迭代过程中可能出现的负值问题，自定义事件处理有助于提高数值计算的整体效率和稳定性。通过防止不必要的迭代或错误，它简化了计算过程，提高了依赖于solve_ivp函数的算法性能，从而使科学模拟和分析的执行更加顺畅。

总的来说，利用自定义事件处理器来防止在事件处理过程中出现负值y，提供了一种稳健而灵活的解决方案，增强了涉及微分方程的数值计算的可靠性、准确性和效率。通过融入领域特定的知识并微调事件检测机制，用户可以确保第三方包的顺利执行，并优化科学模拟和分析的性能。

Answer 2

我理解的情况是，这里出现的问题是因为 funct 在连续的迭代中被调用时，前后传入的 y 值分别是 [0.14553212]（事件发生前）和 [-0.01227461]（事件发生后），而这个变化的步长太大了。在这个简化的例子中，当 y 为负值时，funct 返回的是 NaN（不是一个数字），而 solve_ivp 能够处理这个情况，找到事件，结束计算，并返回事件发生之前的解。但在你的实际代码中，funct 在接收到无效输入时会抛出错误，这样 solve_ivp 就无法处理了。

你可以考虑给 funct 加个保护，检测导致错误的情况，然后返回一个合理的值（这个值不会改变结果），而不是直接调用 funct。或者，你也可以尝试调用 funct，如果出现错误，就用 except 返回一个合理的值？

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

如果这样做还是不行，你可以选择“拒绝这一步”，然后尝试用更小的步长，如果你愿意直接使用 OdeSolver 对象，就像 solve_ivp 那样。不过，solve_ivp 本身并没有实现这种逻辑的机制。