我使用scipy.solve_ivp来求解一个循环中的常微分方程，这样在每次迭代中我都会改变常微分方程中的某个项（这就像施加的力一样）。
就我的问题而言--我似乎无法理解为什么当我跟踪过零事件时，解有时会停滞，并且基本上会清零（如下图所示）。代码有点长，但实际上只有一小部分与我的问题相关
（1）初始化某些常量

MASS = 2.6e-4
WING_LENGTH = 0.07  # meters
GYRATION_RADIUS = 0.6 * WING_LENGTH
MoI = MASS * GYRATION_RADIUS**2 # Moment of Inertia
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

# drag coefficients from whitney & wood (JFM 2010):

C_D_MAX = 3.4
C_D_0 = 0.4
C_L_MAX = 1.8
ZERO_CROSSING = 1
RADIAN45 = np.pi / 4
RADIAN135 = 3 * RADIAN45

（2）定义一个RobotSimulation类的 Package 器--这是一个简单的接口，允许轻松地更改ODE的值。重要的函数是solve_dynamics（我相信其他函数不是问题的一部分）

def phi_dot_zero_crossing_event(t, y):
    """
    this event is given to solve_ivp to track if phi_dot == 0
    :param t: unused time variable
    :param y: a vector of [phi,phi_dot]
    :return:
    """
    return y[1]

class RobotSimulation:
    def __init__(self, motor_torque=lambda x: 0, alpha=RADIAN45, phi0=0.0, phi_dot0=0.01, start_t=0, end_t=0.05, delta_t=0.001):
        self.solution = 0
        self.motor_torque = motor_torque
        self.alpha = alpha  # angle of attack
        self.phi0 = phi0
        self.phi_dot0 = phi_dot0
        self.start_t = start_t
        self.end_t = end_t
        self.delta_t = delta_t

    def set_time(self, new_start, new_end) -> None:
        """
        we update the start time to be the end of the prev time and the end time
        as the last end time + the window size we wish
        """
        self.start_t = new_start
        self.end_t = new_end

    def set_init_cond(self, new_phi0, new_phi_dot0) -> None:
        """
        we update the initial conditions to be the last value of previous iterations
        """
        self.phi0 = new_phi0
        self.phi_dot0 = new_phi_dot0

    def set_motor_torque(self, new_motor_torque) -> None:
        """
        sets the moment to the new value
        """
        self.motor_torque = new_motor_torque

    def flip_alpha(self) -> None:
        """
        (Irrelevant for my question) changes the angle of attack's sign
        :return:
        """
        self.alpha = RADIAN135 if self.alpha == RADIAN45 else RADIAN45

    def 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 drag_torque(self, phi_dot):
        """
        the drag moment
        """
        return 0.5 * AIR_DENSITY * WING_AREA * self.c_drag() * GYRATION_RADIUS * phi_dot * np.abs(phi_dot)

    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.motor_torque(t) - self.drag_torque(phi_dot)) / MoI]
        return dy_dt

    def solve_dynamics(self, phiarr, phidotarr, phiddotarr, angarr, torquearr):
        """
        solves the ODE
        :return:
        """
        phi_0, phi_dot_0 = self.phi0, self.phi_dot0
        start_t, end_t, delta_t = self.start_t, self.end_t, self.delta_t
        phi_dot_zero_crossing_event.terminal = True

        ang = []
        times_between_zero_cross = []
        sol_between_zero_cross = []
        while start_t < end_t:
            torquearr.append(self.motor_torque(start_t))
            sol = solve_ivp(self.phi_derivatives, t_span=(start_t, end_t), y0=[phi_0,phi_dot_0],events=phi_dot_zero_crossing_event)
            self.solution = sol.y
            ang.append(self.alpha * np.ones(len(sol.t)))  # set alpha for every t based on solution's size
            times_between_zero_cross.append(sol.t)
            sol_between_zero_cross.append(sol.y)
            if sol.status == ZERO_CROSSING:
                start_t = sol.t[-1] + delta_t
                phi_0, phi_dot_0 = sol.y[0][-1], sol.y[1][-1]  # last step is now initial value
                self.flip_alpha()
            else:  # no zero crossing
                break
        time = np.concatenate(times_between_zero_cross)
        phi, phi_dot = np.concatenate(sol_between_zero_cross, axis=1)
        ang = np.concatenate(ang)

        _, phi_ddot = self.phi_derivatives(time, [phi, phi_dot])
        phiarr.append(phi)
        phidotarr.append(phi_dot)
        phiddotarr.append(phi_ddot)
        angarr.append(ang)

（3）最后，主函数依次调用solve_dynamics以获得不同的动作值，并更新时间窗口和初始条件：

phi_arr, phi_dot_arr, phi_ddot_arr, angarr, torquearr = [], [], [], [], []
phi0, phidot0 = 0, 0.01
start_t, end_t, delta_t = 0, 0.05, 0.001
sim = RobotSimulation()
for action in [1, -1, 1, -1, 1, -1]:
    sim.set_motor_torque(lambda x: action)
    sim.solve_dynamics(phi_arr, phi_dot_arr, phi_ddot_arr, angarr, torquearr)
    phi0, phidot0 = sim.solution[0][-1], sim.solution[1][-1]
    start_t = end_t
    end_t += 0.05
    sim.set_init_cond(phi0, phidot0)
    sim.set_time(start_t, end_t)
phi_arr = np.concatenate(phi_arr)
phi_dot_arr = np.concatenate(phi_dot_arr)
phi_ddot_arr = np.concatenate(phi_ddot_arr)
angle_arr = np.concatenate(angarr)

直观地说，由于我所应用的操作只是重复+1和-1，因此我无法理解为什么图中的蓝色点（见下图）的操作如此不同