如何在N体模拟中合并天体并修正单向引力问题？

一、单向引力问题的原因与修复 #

问题根源 #

当前代码中gravity方法使用set_ax/set_ay直接赋值加速度，而非累加。当两个天体互相调用gravity时，第二次调用会覆盖第一次计算的加速度值，导致最终只有最后一次的引力生效，表现为单向引力。同时，每次计算引力前未重置加速度，会导致加速度持续累积，进一步加剧错误。

修复步骤 #

1. 修改加速度更新方式：将gravity方法中的set_ax/set_ay替换为change_ax/change_ay，实现加速度的累加（每个天体的加速度是所有其他天体引力的总和）。
2. 重置初始加速度：在每次计算引力前，将所有天体的加速度重置为0，避免累积错误。
3. 调整主循环逻辑：先统一处理所有引力计算，再统一更新天体位置，避免重复更新导致的逻辑混乱。

关键代码修改 #

修改gravity方法中的加速度更新逻辑

# 原代码（错误）
if self.y_pos > otherbody.y_pos:
    self.set_ay(-m.sin(theta)*a)
else:
    self.set_ay(m.sin(theta)*a)

if self.x_pos > otherbody.x_pos:
    self.set_ax(-m.cos(theta)*a)
else:
    self.set_ax(m.cos(theta)*a)

# 修改后（正确）
if self.y_pos > otherbody.y_pos:
    self.change_ay(-m.sin(theta)*a)
else:
    self.change_ay(m.sin(theta)*a)

if self.x_pos > otherbody.x_pos:
    self.change_ax(-m.cos(theta)*a)
else:
    self.change_ax(m.cos(theta)*a)

修改主循环逻辑

while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
            exit()
    WINDOW.blit(BACKGROUND, BACKGROUND_RECT)
    
    # 重置所有天体的加速度为0
    for body in body_group:
        body.set_ax(0)
        body.set_ay(0)
    
    # 遍历body_pairs副本，避免合并时列表变化导致迭代错误
    for body, otherbody in list(body_pairs):
        if body not in body_group or otherbody not in body_group:
            continue
        body.gravity(otherbody)
        otherbody.gravity(body)
    
    # 统一更新所有天体位置
    for body in body_group:
        body.animate()
    
    body_group.draw(WINDOW)
    pygame.display.update()
    CLOCK.tick(60)

二、天体合并功能实现 #

1. 修复碰撞判断逻辑 #

原代码中dx < self.radius*2 and dy < self.radius*2的碰撞判断完全错误，正确逻辑应为：两个天体的距离小于等于两者半径之和，即：

r = m.sqrt(dx**2 + dy**2)
if r <= self.radius + otherbody.radius:
    # 执行合并逻辑
else:
    # 计算引力

2. 合并规则（遵循物理守恒） #

• 质量：新天体质量 = 两者质量之和
• 半径：按体积守恒（均匀球体体积与半径立方成正比），新半径 = (self.radius³ + otherbody.radius³) ** (1/3)
• 位置：质心坐标 = ( (m1*x1 + m2*x2)/(m1+m2), (m1*y1 + m2*y2)/(m1+m2) )
• 速度：动量守恒 = (m1*vx1 + m2*vx2)/(m1+m2)（vy同理）
• 颜色：可直接沿用其中一个天体的颜色，或取RGB分量平均值

3. 合并功能代码实现 #

修改gravity方法，加入合并逻辑：

def gravity(self, otherbody):
    if self == otherbody:
        return
    dx = self.x_pos - otherbody.x_pos
    dy = self.y_pos - otherbody.y_pos
    r = m.sqrt(dx**2 + dy**2)

    # 碰撞合并判断
    if r <= self.radius + otherbody.radius:
        total_mass = self.mass + otherbody.mass
        # 计算质心位置
        new_x = (self.mass * self.x_pos + otherbody.mass * otherbody.x_pos) / total_mass
        new_y = (self.mass * self.y_pos + otherbody.mass * otherbody.y_pos) / total_mass
        # 计算新速度（动量守恒）
        new_vx = (self.mass * self.vx + otherbody.mass * otherbody.vx) / total_mass
        new_vy = (self.mass * self.vy + otherbody.mass * otherbody.vy) / total_mass
        # 计算新半径（体积守恒）
        new_radius = (self.radius**3 + otherbody.radius**3) ** (1/3)
        # 获取原天体颜色
        new_color = self.image.get_at((self.radius, self.radius))

        # 创建新天体并添加到组中
        new_body = Body(total_mass, new_radius, (new_x, new_y), new_color, new_vx, new_vy)
        body_group.add(new_body)
        # 移除原天体
        self.kill()
        otherbody.kill()
        # 更新全局天体列表与配对
        global body_list, body_pairs
        body_list = list(body_group)
        body_pairs = list(itertools.combinations(body_list, 2))
        return  # 合并后不再计算引力

    # 正常计算引力逻辑
    try:
        a = G * otherbody.mass / r**2
        theta = m.asin(dy / r)

        if self.y_pos > otherbody.y_pos:
            self.change_ay(-m.sin(theta)*a)
        else:
            self.change_ay(m.sin(theta)*a)
        
        if self.x_pos > otherbody.x_pos:
            self.change_ax(-m.cos(theta)*a)
        else:
            self.change_ax(m.cos(theta)*a)
    except ZeroDivisionError:
        pass

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