为何Numba计算到特定索引后向量出现NaN值？

可能的触发原因 #

• 数值溢出/下溢：流浪行星靠近地球时，引力加速度会陡增，导致速度、位置的数值超出当前浮点数类型的范围，最终变成NaN。另外你改object dtype完全没必要——这会废掉Numba的优化，还可能引入更多类型问题。
• 积分步长不合理：如果用固定步长的积分器（比如欧拉法），引力突变时步长太大，直接导致计算发散出NaN。
• 引力公式奇点：两个天体距离趋近于0时，G*M/r²会趋近于无穷大，直接计算会产出无穷值，后续运算就变成NaN。
• 精度不足：如果初始化数组时用了float32这类低精度类型，小范围的浮点数更容易出现溢出问题。

具体解决办法 #

1. 改用高精度浮点数 #

把数组的dtype明确设为np.float64（默认虽然是这个，但最好手动确认），别用object：

positions = np.zeros((n_steps, 3), dtype=np.float64)
velocities = np.zeros((n_steps, 3), dtype=np.float64)

float64的动态范围远大于float32，能大幅降低溢出概率。

2. 给引力计算加距离阈值 #

在计算引力时，给最小距离设一个阈值，避免分母过小或为0：

@njit
def compute_gravity(pos1, pos2, mass1, mass2, G, min_dist=1e6):
    r = pos2 - pos1
    dist_sq = np.dot(r, r)
    # 强制最小距离，避免数值爆炸
    dist_sq = max(dist_sq, min_dist**2)
    dist = np.sqrt(dist_sq)
    force_mag = G * mass1 * mass2 / dist_sq
    force = force_mag * r / dist
    return force

这里的min_dist可以设为地球半径的几倍，既不影响模拟精度，又能规避奇点问题。

3. 换成自适应步长积分器 #

固定步长在引力突变场景下很容易崩，换成自适应步长的积分方法（比如RK45），Numba能完美兼容这类逻辑：

@njit
def adaptive_rk_step(pos, vel, mass, dt, G, min_dist, tolerance=1e-8):
    # 计算当前加速度
    acc = compute_gravity(pos, rogue_pos, mass, rogue_mass, G, min_dist)
    # 简化示例：通过对比半步长和全步长的误差调整步长
    # 实际需要完整的RK45误差计算逻辑
    dt_half = dt / 2
    pos_half = pos + vel * dt_half + 0.5 * acc * dt_half**2
    vel_half = vel + acc * dt_half
    acc_half = compute_gravity(pos_half, rogue_pos, mass, rogue_mass, G, min_dist)
    pos_full = pos + vel_half * dt + 0.5 * acc_half * dt**2
    
    # 计算误差，调整步长
    error = np.linalg.norm(pos_full - (pos + vel*dt + 0.5*acc*dt**2))
    dt_new = dt * np.sqrt(tolerance / error) if error > 0 else dt
    
    vel_new = vel + 0.5*(acc + acc_half)*dt_new
    pos_new = pos + vel*dt_new + 0.5*acc*dt_new**2
    return pos_new, vel_new, dt_new

引力变化剧烈时会自动缩小步长，避免数值发散。

4. 检查初始参数合理性 #

确认流浪行星的初始质量、速度、位置是否合理——如果设成恒星级质量（1e30kg）或者初始距离过近，直接就会触发引力溢出。可以先把流浪行星设为木星级（1e27kg），初始距离放在柯伊伯带以外（比如1e13m）测试。

5. 加NaN检测调试 #

在Numba函数里加检测逻辑，出现NaN时打印状态并中断，方便定位问题：

@njit
def simulate():
    for i in range(n_steps):
        # ... 你的计算逻辑 ...
        if np.isnan(positions[i]).any():
            print(f"NaN出现在第{i}步")
            print(f"当前地球位置: {positions[i]}")
            print(f"流浪行星位置: {rogue_pos[i]}")
            break

能快速定位是哪种场景触发的NaN，针对性调整。

内容的提问来源于stack exchange，提问作者M.M. CAN