1D重力模拟无法实现摆锤式振荡的原因及解决方案咨询

你遇到的这个问题其实是数值重力模拟里的经典坑——核心原因是普通数值积分方法的误差累积，导致系统能量不守恒，天体的运动逐渐偏离预期的振荡轨道。那些能稳定来回振荡的模拟器，主要靠这几招解决：

1. 换用更稳定的辛积分器（比如Velocity Verlet） #

欧拉法这类简单积分方法的误差会随时间快速累积，尤其是在引力突变的近距离相遇时，能量会不断“泄漏”或者凭空增加，导致天体越来越跑远。而Velocity Verlet积分是专门针对保守系统（比如引力系统）设计的辛积分器，能长期保持能量守恒，是物理模拟里的首选。

给你个简化的伪代码示例：

# 初始化位置、速度、加速度
x1, x2 = initial_positions
v1, v2 = initial_velocities
a1, a2 = calculate_accelerations(x1, x2, m1, m2)

# 积分循环
dt = time_step
for _ in range(num_steps):
    # 更新位置
    x1 += v1*dt + 0.5*a1*dt**2
    x2 += v2*dt + 0.5*a2*dt**2
    # 计算新的加速度
    new_a1, new_a2 = calculate_accelerations(x1, x2, m1, m2)
    # 更新速度（用新旧加速度的平均值）
    v1 += 0.5*(a1 + new_a1)*dt
    v2 += 0.5*(a2 + new_a2)*dt
    # 更新加速度变量
    a1, a2 = new_a1, new_a2

这种方法能极大降低能量误差，让天体的振荡长期保持稳定。

2. 添加引力软核修正 #

当两个天体距离趋近于0时，引力会趋向无穷大，这会导致加速度突变，积分器根本算不准。给引力公式加个软核项，就能避免这种极端情况：

# 原引力公式（容易出问题）
F = G*m1*m2 / r**2
# 加软核后的公式
epsilon = 0.01  # 很小的常数，根据你的模拟尺度调整
F = G*m1*m2 / (r**2 + epsilon**2)

软核项会让近距离的引力趋于一个有限值，积分过程更平滑，不会出现加速度突然爆炸的情况，自然能减少能量误差。

3. 定期校正系统能量 #

如果你的模拟对精度要求极高，还可以定期检查系统的总能量（动能+引力势能），当误差超过阈值时，手动校正速度来拉回能量：

def calculate_total_energy(x1, x2, v1, v2, m1, m2):
    kinetic = 0.5*m1*v1**2 + 0.5*m2*v2**2
    potential = -G*m1*m2 / abs(x1 - x2)
    return kinetic + potential

# 记录初始能量
initial_energy = calculate_total_energy(x1, x2, v1, v2, m1, m2)

# 每N步校正一次
if step % N == 0:
    current_energy = calculate_total_energy(x1, x2, v1, v2, m1, m2)
    # 计算能量校正因子
    correction_factor = sqrt(initial_energy / current_energy)
    # 按比例调整速度
    v1 *= correction_factor
    v2 *= correction_factor

这招相当于给系统“踩刹车”或“加油门”，强制把能量拉回初始值，保证天体不会越跑越远。

4. 针对振荡场景的特殊约束（可选） #

如果你的需求就是让两个天体严格在初始距离附近振荡，也可以直接给位置加约束（比如用弹簧力模拟等效的引力效果），不过这属于场景特定的hack，不如前几种方法通用。

总的来说，最常用的方案是Velocity Verlet积分+软核修正，这组合几乎能解决绝大多数1D引力模拟的稳定性问题，和你提到的那些模拟器用的思路基本一致。

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