scipy.integrate.solve_ivp中的第n次交叉事件检测

Question

我想知道有没有办法在检测到特定数量的事件后停止常微分方程（ODE）的计算。我知道scipy.integrate.solve_ivp可以让事件决定是否终止，但我想知道是否可以在第n个事件时终止，其中n是一个大于1的整数。

举个例子，是否可以只计算前面三个事件的交互，然后就停止计算？如果不行，有没有其他的库可以做到这一点？

我看过scipy关于ODE积分的官方文档，但没有找到相关的信息。

Answer 1

这里有一个简单的技巧。正如评论中提到的，wrapped_event 在达到所需事件数量之前会返回一个任意的值，然后才开始返回真正的事件函数。（在检测事件交叉时，由于符号问题会有一点额外的复杂性。）要在自己的代码中使用这个技巧，你只需要用 wrap_event 来包裹你的事件函数，然后把结果传递给 solve_ivp。虽然不要求遵循事件的 terminal 属性，但这很容易处理。

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

def wrap_event(event, n_events, t0, y0):
    def wrapped_event(t, y, n_events=n_events):
        new_event = event(t, y)
        if wrapped_event.count < n_events:
            new_sign = np.sign(new_event)
            if new_sign == -wrapped_event.sign:
                wrapped_event.count += 1
                wrapped_event.sign = new_sign
            return -1
        return new_event * (-1)**(n_events)
    wrapped_event.count = 1
    wrapped_event.sign = np.sign(event(t0, y0))
    wrapped_event.terminal = event.terminal
    return wrapped_event

def f(t, y):
    return [y[1], -y[0]]

def event(t, y):
    return y[0]
event.terminal = True

t0 = 0
tf = 100
y0 = [1, 0]
n_events = 3

sol = integrate.solve_ivp(f, (t0, tf), y0, dense_output=True, 
                          events=wrap_event(event, n_events, t0, y0))
t = np.linspace(sol.t[0], sol.t[-1], 300)
x = sol.sol(t)
plt.plot(t, x[0])