模拟中月球为何向地球螺旋靠近？求问题排查与建议

问题概述 #

我正在编写一个模拟地月间引力相互作用的半真实程序，但运行时月球持续向地球螺旋靠近，无法维持稳定轨道。以下是我的代码：

from math import sin,cos,sqrt,atan2,pi
import pygame
pygame.init()

class Planet:
    dt = 1/100
    G = 6.67428e-11 #G constant
    scale = 1/(1409466.667) #1 m = 1/1409466.667 pixels
    def __init__(self,x=0,y=0,radius=0,color=(0,0,0),mass=0,vx=0,vy=0):
        self.x = x #x-coordinate pygame-window
        self.y = y #y-coordinate pygame-window
        self.radius = radius
        self.color = color
        self.mass = mass
        self.vx = vx #velocity in the x axis
        self.vy = vy #velocity in the y axis
        
    def draw(self,screen):
        pygame.draw.circle(screen, self.color, (self.x, self.y), self.radius)
    
    def orbit(self,trace):
        pygame.draw.rect(trace, self.color, (self.x, self.y, 2, 2))
        
    def update_vel(self,Fnx,Fny):
        ax = Fnx/self.mass #Calculates acceleration in x- and y-axis for body 1.
        ay = Fny/self.mass
        self.vx -= ((ax * Planet.dt)/Planet.scale)
        self.vy -= ((ay * Planet.dt)/Planet.scale)
        self.update_pos()
     
    def update_pos(self):
        self.x += ((self.vx * Planet.dt)) #changes position considering each body's velocity.
        self.y += ((self.vy * Planet.dt))
        
    def move(self,body):
        dx = (self.x - body.x) #Calculates difference in x- and y-axis between the bodies
        dy = (self.y - body.y)
        r = (sqrt((dy**2)+(dx**2))) #Calculates the distance between the bodies
        angle = atan2(dy, dx) #Calculates the angle between the bodies with atan2!
        if r < self.radius: #Checks if the distance between the bodies is less than the radius of the bodies. Uses then Gauss gravitational law to calculate force.
            F = 4/3 * pi * r
            Fx = cos(angle) * F
            Fy = sin(angle) * F
        else:  
            F = (Planet.G*self.mass*body.mass)/((r/Planet.scale)**2) #Newtons gravitational formula.
            Fx = cos(angle) * F
            Fy = sin(angle) * F
        return Fx,Fy

def motion():
    for i in range(0,len(bodies)):
        Fnx = 0 #net force
        Fny = 0
        for j in range(0,len(bodies)):
            if bodies[i] != bodies[j]:
                Fnx += (bodies[i].move(bodies[j]))[0]
                Fny += (bodies[i].move(bodies[j]))[1]
            elif bodies[i] == bodies[j]:
                continue
        bodies[i].update_vel(Fnx,Fny)
        bodies[i].draw(screen)
        bodies[i].orbit(trace)
        Fnx,Fny=0,0 

screen = pygame.display.set_mode([900,650]) #width - height
trace = pygame.Surface((900, 650))
pygame.display.set_caption("Moon simulation")
FPS = 150 #how quickly/frames per second our game should update.

earth = Planet(450,325,30,(0,0,255),5.97219*10**(24),-24.947719394204714/2) #450= xpos,325=ypos,30=radius
luna = Planet(450,(575/11),10,(128,128,128),7.349*10**(22),1023)
bodies = [earth,luna]

running = True
clock = pygame.time.Clock()

while running: #if user clicks close window
    clock.tick(FPS)    
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False
            
    screen.fill((0,0,0))
    pygame.Surface.blit(screen, trace, (0, 0))
    motion()

    pygame.display.flip() 

pygame.quit()

我计划在修复地月系统后扩展为三体系统，欢迎各类建议与指导。

问题根源分析 #

1. 引力方向完全错误 #

move方法中，dx = self.x - body.x计算的是当前行星到目标行星的反向距离，导致引力方向与实际相反——本该吸引的力变成了排斥或错误的加速力，直接破坏轨道平衡。

2. 速度缩放逻辑颠倒 #

update_vel方法中，将真实加速度转换为像素速度时用了除法而非乘法：

self.vx -= ((ax * Planet.dt)/Planet.scale)

根据scale = 1/(1409466.667)的定义（1米对应scale像素），真实速度变化量（m/s）需乘以scale才能转换为像素/秒，除法会导致速度变化被错误放大，轨道迅速衰减。

3. 积分顺序引入累积误差 #

当前代码计算单个行星受力后立即更新其位置，后续行星的受力计算会基于已改变的位置，每帧都引入误差，长期累积后轨道必然崩溃。

4. 近距离引力公式错误 #

当行星距离小于半径时，使用的F = 4/3 * pi * r完全不符合高斯引力定律，正确的内部引力公式应为F = G*self.mass*body.mass*r/(R³)（R为行星半径），当前公式会导致近距离力计算完全失效。

5. 初始速度不符合物理规律 #

月球的初始速度1023未基于真实地月距离和引力公式计算，地球的初始速度也未匹配系统质心运动，导致整体动量不守恒，轨道天然不稳定。

修正方案 #

1. 修正引力方向 #

修改move方法的距离计算，让方向指向目标行星：

def move(self,body):
    dx = body.x - self.x  # 改为目标行星减当前行星，得到指向目标的方向
    dy = body.y - self.y
    r = sqrt(dy**2 + dx**2)
    angle = atan2(dy, dx)
    # 修正碰撞判断为两行星半径之和
    if r < self.radius + body.radius:
        # 替换为正确的内部引力公式
        r_m = r / Planet.scale
        R_m = self.radius / Planet.scale
        F = (Planet.G * self.mass * body.mass * r_m) / (R_m**3)
    else:
        r_m = r / Planet.scale  # 转换为真实距离（米）
        F = (Planet.G * self.mass * body.mass) / (r_m**2)
    Fx = cos(angle) * F
    Fy = sin(angle) * F
    return Fx,Fy

2. 修复速度缩放逻辑 #

调整update_vel的速度更新公式，修正符号和缩放方式：

def update_vel(self,Fnx,Fny):
    ax = Fnx / self.mass
    ay = Fny / self.mass
    # 真实加速度转速度变化（m/s），再转换为像素/秒，方向与引力一致用加号
    self.vx += ax * Planet.dt * Planet.scale
    self.vy += ay * Planet.dt * Planet.scale
    # 移除此处的update_pos调用，统一后续更新

3. 改进积分顺序 #

修改motion函数，先批量计算所有行星的合力，再统一更新速度和位置，避免顺序误差：

def motion():
    # 第一步：批量计算所有行星的合力
    forces = []
    for i in range(len(bodies)):
        Fnx = 0
        Fny = 0
        for j in range(len(bodies)):
            if i != j:
                fx, fy = bodies[i].move(bodies[j])
                Fnx += fx
                Fny += fy
        forces.append((Fnx, Fny))
    
    # 第二步：统一更新所有行星的速度
    for i in range(len(bodies)):
        bodies[i].update_vel(*forces[i])
    
    # 第三步：统一更新位置并绘制轨迹
    for body in bodies:
        body.update_pos()
        body.draw(screen)
        body.orbit(trace)

4. 设置符合物理规律的初始速度 #

基于真实地月距离计算月球的环绕速度，同时将地球初始速度设为0（或匹配质心运动）：

# 地月真实距离（3.84e8米）转换为像素
moon_distance_m = 3.84e8
moon_distance_px = moon_distance_m * Planet.scale
# 月球环绕速度（m/s）：sqrt(G*M_earth / r)
moon_v_m = sqrt(Planet.G * earth.mass / moon_distance_m)
# 转换为像素/秒
moon_v_px = moon_v_m * Planet.scale

earth = Planet(450, 325, 30, (0,0,255), 5.97219e24, 0, 0)
# 月球初始位置在地球正上方，速度沿x轴横向环绕
luna = Planet(450, 325 - moon_distance_px, 10, (128,128,128), 7.349e22, moon_v_px, 0)
bodies = [earth, luna]

扩展三体系统的建议 #

• 替换稳定的积分算法：欧拉法误差累积快，建议改用Verlet积分或Runge-Kutta（RK4）法，大幅提升长期模拟稳定性。
• 优化受力计算：N体系统中，i对j的力与j对i的力大小相等方向相反，可只计算一次并复用，减少一半计算量。
• 切换到质心参考系：将所有行星的位置和速度转换到系统质心参考系，避免整体漂移，让模拟更直观稳定。
• 完善碰撞处理：若需处理碰撞，需实现动量守恒和能量损失计算，而非简单的内部引力公式。
• 引入近似算法：当行星数量增多时，可使用Barnes-Hut算法降低计算复杂度，提升模拟效率。

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