如何在scipy.solve_ivp求解过程中随事件修改self._alpha变量

我正在使用scipy.solve_ivp求解方程phi''(t) = a(t) - b * (phi'(t))^2，希望在phi'(t)符号变化时修改b的值。

该方程代表翼面的扭矩平衡（此背景对问题无关键影响），求解逻辑封装在Solver类中。核心函数为sequential_solver、phi_dot_zero_crossing_event和phi_derivatives，可运行代码如下：

# some constants to run the code ...
NUM_SAMPLES = 100  # every time window has this many samples
MASS = 1
WING_LENGTH = 0.07 
AERODYNAMIC_CENTER = 0.7 * WING_LENGTH
GYRATION_RADIUS = 0.6 * WING_LENGTH
MoI = MASS * GYRATION_RADIUS ** 2
AIR_DENSITY = 1.2 
WING_AREA = 0.5 * WING_LENGTH * (0.5 * WING_LENGTH) * np.pi  # 1/2 ellipse with minor radios ~ 1/2 major = length/2
C_D_MAX = 3.4
C_D_0 = 0.4
C_L_MAX = 1.8 
class Solver:
    def __init__(self) -> None:
        self._alpha = 45

    def get_c_drag(self) -> float:
        """
        calculates the drag coefficient based on the angle of attack
        """
        return (C_D_MAX + C_D_0) / 2 - (C_D_MAX - C_D_0) / 2 * np.cos(2 * self._alpha)

    def get_c_lift(self) -> float:
        """
        calculates the lift coefficient based on the angle of attack
        """
        return C_L_MAX * np.sin(2 * self._alpha)

    def tau_z(self, t: float) -> np.float:
        """
        the moment we apply using the motors
        """
        return np.sin(2 * np.pi * t)


    def tau_drag(self, phi_dot):
        """
        the drag moment
        """
        return 0.5 * AIR_DENSITY * WING_AREA * self.get_c_drag() * (GYRATION_RADIUS ** 2) * (phi_dot ** 2)

    def phi_derivatives(self, t, y):
        """
        A function that defines the ODE that is to be solved: I * phi_ddot = tau_z - tau_drag.
        We think of y as a vector y = [phi,phi_dot]. the ode solves dy/dt = f(y,t)
        :return:
        """
        phi, phi_dot = y[0], y[1]
        dy_dt = [phi_dot, (self.tau_z(t) - self.tau_drag(phi_dot)) / MoI]  # d[phi,phi_dot]]/dt = [phi_dot,phi_ddot]
        return dy_dt

    def phi_dot_zero_crossing_event(self, t, y):
        return y[1]

    def sequential_solver(self, use_windows):
        """
        solves the ODE
        :return:
        """
        phi_0 = 2e-3
        phi_dot_0 = 0
        prev_t = 0
        start_t, end_t, delta_t = 1e-1, 5, 1e-1
        inner_delta_t = delta_t / NUM_SAMPLES

        
        time = np.arange(start_t, end_t, inner_delta_t)
        sol = solve_ivp(self.phi_derivatives, t_span=(start_t, end_t), y0=[phi_0, phi_dot_0], t_eval=time,
                            events=self.phi_dot_zero_crossing_event)
        phi, phi_dot = sol.y
        _, phi_ddot = self.phi_derivatives(time, [phi, phi_dot])

目前仅能通过phi_dot_zero_crossing_event检测到phi_dot的过零事件，但不知道如何利用该事件修改self._alpha的值：当phi_dot>0时设为+45，否则设为-45。

要在过零事件触发时修改self._alpha，需要分段求解ODE——每次检测到过零事件后，以事件点为分界，更新参数后继续求解后续区间。具体实现步骤如下：

1. 配置过零事件的终止属性 #

修改phi_dot_zero_crossing_event，设置terminal=True让求解器在检测到过零时自动停止，并指定检测所有方向的过零：

def phi_dot_zero_crossing_event(self, t, y):
    return y[1]
# 标记事件为终止事件，触发后求解器停止
phi_dot_zero_crossing_event.terminal = True
# direction=0表示检测所有方向的过零（正→负、负→正都触发）
phi_dot_zero_crossing_event.direction = 0

2. 重构求解逻辑为循环分段求解 #

在sequential_solver中用循环不断求解到过零事件点，更新self._alpha后，以事件点的状态作为初始值继续求解下一段：

def sequential_solver(self, use_windows):
    phi_0 = 2e-3
    phi_dot_0 = 0
    current_t = 1e-1
    end_t = 5
    inner_delta_t = 1e-1 / NUM_SAMPLES

    # 存储所有分段的求解结果
    all_times = []
    all_phi = []
    all_phi_dot = []

    while current_t < end_t:
        # 定义当前求解区间
        t_span = (current_t, end_t)
        y0 = [phi_0, phi_dot_0]
        # 生成当前区间的采样时间点
        time_segment = np.arange(current_t, end_t, inner_delta_t)
        # 求解当前段ODE
        sol = solve_ivp(
            self.phi_derivatives,
            t_span=t_span,
            y0=y0,
            t_eval=time_segment,
            events=self.phi_dot_zero_crossing_event
        )

        # 保存当前段结果
        all_times.extend(sol.t)
        all_phi.extend(sol.y[0])
        all_phi_dot.extend(sol.y[1])

        # 检查是否触发过零事件
        if sol.status == 1:  # status=1表示因事件终止求解
            event_t = sol.t_events[0][0]
            event_y = sol.y_events[0][0]
            # 根据过零前后的phi_dot变化方向更新_alpha
            if len(sol.y[1]) >= 2:
                prev_phi_dot = sol.y[1][-2]
                current_phi_dot = sol.y[1][-1]
                # 从负到正，后续phi_dot>0，设为45
                if current_phi_dot > prev_phi_dot:
                    self._alpha = 45
                # 从正到负，后续phi_dot<0，设为-45
                else:
                    self._alpha = -45
            # 更新初始状态和当前时间，准备下一段求解
            phi_0, phi_dot_0 = event_y[0], event_y[1]
            current_t = event_t
        else:
            # 未触发事件，已求解到终点，退出循环
            break

    # 合并所有分段结果
    time = np.array(all_times)
    phi = np.array(all_phi)
    phi_dot = np.array(all_phi_dot)
    _, phi_ddot = self.phi_derivatives(time, [phi, phi_dot])
    return time, phi, phi_dot, phi_ddot

3. 关键说明 #

• scipy.solve_ivp的事件函数仅能触发求解终止，无法直接在求解过程中修改参数，因此必须通过分段求解的方式手动重启求解器，应用更新后的self._alpha。
• 过零方向的判断通过对比事件点前后的phi_dot值实现，确保_alpha的设置符合后续phi_dot的符号趋势。

内容的提问来源于stack exchange，提问作者Hadar