数值积分中积分边界的实现方法探讨

我使用Metropolis-Hastings马尔可夫链蒙特卡洛算法，通过数值积分模拟非均匀磁场中粒子的轨迹，以拟合实际粒子数据。目前存在的问题是拟合过程会覆盖其他粒子的轨迹（如附图所示的正反粒子，拟合终止于另一粒子轨迹起点附近）。我希望单独积分单个粒子，例如从z≈337处开始积分，在成对产生原点z≈550处停止。尝试在积分中加入break语句，通过判断zmodel[index_max] == zmax来终止，但由于浮点精度问题该条件从未触发。现询问是否有其他可单独设置积分边界的数值积分方法。

相关代码 #

def evaluation(theta,phi,E,xi,yi,zi):  ### For creating model/experimental data

    initial_vel = BROH(E)[0]
    gamma_2 = BROH(E)[2]
    relative_mass = BROH(E)[3]

    first_x = np.zeros(len(actual_x))
    first_y = np.zeros(len(actual_y))
    first_z = np.zeros(len(actual_z))

    xmodel = np.zeros(len(actual_x))   ### Store model data here
    ymodel = np.zeros(len(actual_y))
    zmodel = np.zeros(len(actual_z))

    velocity_x = np.zeros(len(actual_x))  ### Store velocity values to calculate subsequent x,y,z model data
    velocity_y = np.zeros(len(actual_y))
    velocity_z = np.zeros(len(actual_z))

    Bx = np.zeros(len(actual_x))
    By = np.zeros(len(actual_y))
    Bz = np.zeros(len(actual_z))

    first_x[0] = xi         ### Initial guesses for x,y,z
    first_y[0] = yi
    first_z[0] = zi

    velocity_x[0] = initial_vel*np.sin(theta)*np.cos(phi)  ### Initial values for velocities
    velocity_y[0] = initial_vel*np.sin(theta)*np.sin(phi)
    velocity_z[0] = initial_vel*np.cos(theta)

    index = 0
    for i in range(len(actual_x) - 1):  ### Loop over experimental/model trajectory
        
        zbzero = zradius[2][0] #for evemt 93  # for event 71 550
        zb = abs(first_z[i] - zbzero)
        if zb > 1000:
            zb = 1000
        
        global Qcharge
        Qcharge = -1.  #positive or negative charge +1 or -1 
        Bz = 1678.5 + 0.080008*zb - 0.019289*zb**2 + 1.3946e-5*zb**3 + 3.0161e-8*zb**4
        Bz = Qcharge*Bz  #for opposite/ normal charge/positive 
        
        Rr = first_x[i]**2 + first_y[i]**2
        if Rr > 1000:
            Rr = 1000
        
        Fact = np.sqrt(Rr) / 40
        Br = Fact*(6.2674e-3 + 0.67562*zb + 1.2677e-4*zb**2 - 6.8352e-6*zb**3 + 6.6604e-9*zb**4)
        Phir = np.arctan2(first_y[i],first_x[i])
        Br = Qcharge*Br #for opposite/ normal charge/positive 
        
        Bx = -2/3*Br*np.cos(Phir)
        By = -2/3*Br*np.sin(Phir)
        
        B_field = np.array([Bx,By,Bz])
        velocity = np.array([velocity_x[i],velocity_y[i],velocity_z[i]])
        cross_product = np.cross(velocity,B_field)
        
        ### Calculate subsequent velocities for model/experimental
        velocity_x[i+1] = velocity_x[i] + const*cross_product[0]*dt / relative_mass
        velocity_y[i+1] = velocity_y[i] + const*cross_product[1]*dt / relative_mass
        velocity_z[i+1] = velocity_z[i] + const*cross_product[2]*dt / relative_mass  

        first_x[i+1] = first_x[i] + velocity_x[i]*dt + 0.5*const*cross_product[0]*dt**2 / relative_mass   
        first_y[i+1] = first_y[i] + velocity_y[i]*dt + 0.5*const*cross_product[1]*dt**2 / relative_mass  
        first_z[i+1] = first_z[i] + velocity_z[i]*dt + 0.5*const*cross_product[2]*dt**2 / relative_mass
        
        if first_x[i+1] > -150 and first_x[i+1] < 150:
            if first_y[i+1] > -150 and first_y[i+1] < 150:
                if first_z[i+1] > 0 and first_z[i+1] < 1000:
                    
                    global index_max
                    index = index + 1
                    xmodel[index] = first_x[i+1] + 0.5*const*cross_product[0]*dt**2 / relative_mass 
                    ymodel[index] = first_y[i+1] + 0.5*const*cross_product[1]*dt**2 / relative_mass  
                    zmodel[index] = first_z[i+1] + 0.5*const*cross_product[2]*dt**2 / relative_mass
                    index_max = index
                    
        if zmodel[index_max] == zmax:
            break
                
    return xmodel[1:index_max], ymodel[1:index_max], zmodel[1:index_max], index_max

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