solve_ivp事件条件触发时机错误：粒子出磁元件时误终止积分

我在用solve_ivp模拟粒子穿过多个圆柱形磁元件的动力学过程，自定义了一个径向终止事件，期望仅当粒子处于磁元件内部时，径向距离r达到元件内径的90%时终止积分。但实际运行中，粒子未发生碰撞、到达磁元件x_exit时，该事件也会错误触发并终止积分。

现有代码实现 #

径向终止事件函数 #

def stop_event(self, t, y, collection):
    """
    Event function to stop the integration when the particle's radial distance r
    exceeds the diameter of any magnetic element within its x-bounds.
    """
    x, y_pos, z = y[0], y[1], y[2]
    r = np.sqrt(y_pos**2 + z**2)  # radial distance r = sqrt(y^2 + z^2)
    
    # Loop through the magnetic elements to check if the particle is in any element's x-range
    for element in collection:
        mag_element = element['magnetic_element']
        x_entrance = mag_element.entrance_x  # Start x of the magnetic element
        x_exit = mag_element.exit_x      # End x of the magnetic element
        radius = mag_element.innermost_loop_diam / 2  # Critical radial boundary (90% of element radius)
            
        # Check if particle is within the x-bounds of the magnetic element
        if x_entrance <= x <= x_exit:
            return (r - 0.9 * radius) # Event triggered when this crosses 0
    
    return 1  # No event if not within the bounds of any magnetic element

积分求解简化代码 #

def solve_trajectory(self, collection, t_span, nb_points, max_step = None, h=1e-6):
    y0 = np.concatenate([self.initial_position, self.initial_velocity])  # initial state [x, y, z, vx, vy, vz]
    
    # Define the stopping event
    def radial_stop_event(t, y):
        return self.stop_event(t, y, collection)
    
    radial_stop_event.terminal = True  # The event should terminate the integration
    radial_stop_event.direction = 1    # Detect crossing in from inside to outside
    
    [...]

    # Compute dt from the number of points
    dt = (t_span[1] - t_span[0]) / (nb_points - 1)
    
    # Call solve_ivp with Runge-Kutta method (RK45)
    sol = solve_ivp(
        fun=lambda t, y: self.equations_of_motion(t, y, collection, h=h),
        t_span=t_span,
        y0=y0,
        method='RK45',
        events = [radial_stop_event, longitudinal_enter_event, longitudinal_exit_event],
        dense_output=True,
        max_step = max_step,
        rtol=1e-8, 
        atol=1e-10
    )
    
    [...]

    
    return sol, times_in_elements

问题验证数据 #

触发前参数：

t=9.203132978272565e-05, x=0.003033979317911664, r=1.4922512843731781e-15, radius=0.00135, condition = -0.0013499999999985078

触发后参数：

t=0.0002882113689884196, x=0.020175674729679623, r=7.060896498215925e-15, radius=0.00135, condition = -0.001349999999992939

问题根源 #

solve_ivp的事件检测逻辑是寻找函数值穿越0的时刻，同时结合direction参数判断穿越方向。当前代码存在两个关键问题：

• 当粒子从磁元件内部（x_entrance <= x <= x_exit）移动到外部时，事件函数的返回值会从r - 0.9*radius（负数，因为r远小于阈值）突然跳转到1，这个跳变会被solve_ivp误判为一次从负到正的穿越，正好匹配你设置的direction=1，从而触发终止事件。
• 事件函数的返回值在元件边界处不连续，导致数值积分器误识别事件。

修复方案 #

修改事件函数，确保返回值在元件边界处连续，且仅当粒子在元件内部且r超过阈值时才会产生0穿越。具体做法：

• 当粒子在元件内部时，返回r - 0.9*radius；
• 当粒子在元件外部时，返回一个始终小于0的值（比如-1），避免产生从负到正的跳变；
• 优先找到粒子当前所在的唯一磁元件（假设元件在x轴上不重叠），避免多元件判断逻辑混乱。

修改后的stop_event函数 #

def stop_event(self, t, y, collection):
    """
    Event function to stop the integration when the particle's radial distance r
    exceeds 90% of the magnetic element's inner radius, only while inside the element.
    """
    x, y_pos, z = y[0], y[1], y[2]
    r = np.sqrt(y_pos**2 + z**2)
    
    current_element = None
    # 定位粒子当前所在的磁元件（若存在）
    for element in collection:
        mag_element = element['magnetic_element']
        x_entrance = mag_element.entrance_x
        x_exit = mag_element.exit_x
        if x_entrance <= x <= x_exit:
            current_element = mag_element
            break
    
    if current_element is not None:
        radius = current_element.innermost_loop_diam / 2
        # 仅在元件内部时返回r与阈值的差值，触发条件为r超过阈值（从负到正穿越0）
        return r - 0.9 * radius
    else:
        # 粒子在所有元件外部时，返回稳定负值，避免误触发
        return -1

优化说明 #

1. 连续返回值：外部返回-1而非1，确保粒子离开元件时，返回值从r-0.9*radius（负数）平滑过渡到-1，彻底消除误判的跳变。
2. 单一元件定位：通过循环找到当前粒子所在的唯一元件，避免多元件判断导致的逻辑冲突。
3. direction参数兼容：保持direction=1不变，仅响应r从小于阈值到大于阈值的正向穿越，符合碰撞终止的需求。

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