基于Pygame的绕日地球模拟轨迹偏移问题排查与解决求助

问题描述 #

使用Pygame开发的绕日地球运动模拟程序运行后，地球轨道随时间逐渐偏移，核心代码如下：

核心类定义 #

import pygame  
from pygame import *  
import math  

class orb:  
    def __init__(self, mass=1, s0=[0,0], v0=[0,0], a=[0,0], r=10, color=(255, 255, 255), fix = False, F=0, Fx=0, Fy=0):  
        self.a= a  
        self.mass=mass  
        self.r=r  
        self.t = 0  
        self.v0=v0  
        self.v=[self.v0[0], self.v0[1]]  
        self.s0=s0  
        self.s=[self.s0[0], self.s0[1]]  
        self.color=color  
        self.fix = fix  
        self.F=F  
        self.Fx = Fx  
        self.Fy = Fy  

位置计算原代码 #

def comPos(self, dt = 0.01):
    if not self.fix:  
        self.v[0] = self.v0[0] + self.a[0] * dt
        self.v[1] = self.v0[1] + self.a[1] * dt
        ds_x = self.v0[0] * dt + (1/2) * self.a[0] * dt**2
        ds_y = self.v0[1] * dt + (1/2) * self.a[1] * dt**2
        self.s[0] = self.s0[0] + ds_x
        self.s[1] = self.s0[1] + ds_y
        self.v0[0] = self.v[0]
        self.v0[1] = self.v[1]
        self.s0[0] = self.s[0]
        self.s0[1] = self.s[1]
        self.t += dt

引力计算原代码 #

def comForce(self, other, G = 10000.0, F = 0.0, dist_cr = 0.0, Fx = 0.0, Fy = 0.0):
    self.F = G * self.mass * other.mass / ( (self.s[0]- other.s[0])**2 + (self.s[1] - other.s[1])**2 )
    self.Fx = (self.F * math.fabs(math.cos(math.atan((other.s[1] - self.s[1]) / (other.s[0] - self.s[0])))))
    self.Fy = (self.F * math.fabs(math.sin(math.atan((other.s[1] - self.s[1]) / (other.s[0] - self.s[0])))))
    if ((other.s[0] - self.s[0])**2 + (other.s[1] - self.s[1])**2)**(1/2) >= dist_cr:
        if other.s[0] - self.s[0] == 0:
            self.a[0] = 0
        else:
            self.a[0] = (( (other.s[0] - self.s[0]) / math.fabs(other.s[0] - self.s[0])) * (self.Fx/self.mass))
        if other.s[1] - self.s[1] == 0:
            self.a[1] = 0
        else:
            self.a[1] = (( (other.s[1] - self.s[1]) / math.fabs(other.s[1] - self.s[1])) * self.Fy / self.mass)

主循环原代码 #

screen = pygame.display.set_mode((800, 600))
clock = pygame.time.Clock()
earth = orb(s0=[150, 300], v0=[0,60], r=10, mass= 3)
sun = orb(s0=[400, 300], r=70, color = [255, 0, 0], fix=True, mass=100)
orb = [sun, earth]

running = True
while running:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False
    screen.fill((0, 0, 0))
    pygame.draw.circle(screen, earth.color, earth.s, earth.r, 0)
    pygame.draw.circle(screen, sun.color, sun.s, sun.r, 0)
    earth.comForce(sun)
    earth.comPos()
    pygame.display.update()
    clock.tick(720)  

轨道偏移原因分析 #

• 欧拉法积分误差：当前使用的显式欧拉法是一阶数值积分，每一步用初始时刻的加速度更新速度和位置，未考虑积分过程中加速度的变化，长时间模拟会积累能量误差，导致轨道膨胀或收缩。
• 三角函数计算精度损失：通过math.atan转余弦、正弦的方式易引入精度误差，还需额外处理x位移差为0的特殊场景，增加出错概率。
• 初始参数不匹配：初始速度不符合圆周运动的向心力要求，初始状态就存在能量偏差，加速轨道偏移。

解决方案 #

1. 改用半隐式欧拉法（辛积分） #

半隐式欧拉法是辛积分方法，能保持保守系统的能量守恒，适合天体运动模拟，长时间运行轨道不会漂移。修改comPos方法：

def comPos(self, dt = 0.01):
    if not self.fix:
        # 先通过当前加速度更新速度
        self.v[0] += self.a[0] * dt
        self.v[1] += self.a[1] * dt
        # 用更新后的速度更新位置
        self.s[0] += self.v[0] * dt
        self.s[1] += self.v[1] * dt
        # 同步初始状态变量
        self.v0[0], self.v0[1] = self.v[0], self.v[1]
        self.s0[0], self.s0[1] = self.s[0], self.s[1]
        self.t += dt

2. 优化引力加速度计算 #

直接通过向量分量计算引力方向，消除三角函数误差和特殊分支：

def comForce(self, other, G = 10000.0, dist_cr = 0.0):
    dx = other.s[0] - self.s[0]
    dy = other.s[1] - self.s[1]
    dist_sq = dx**2 + dy**2
    dist = math.sqrt(dist_sq)
    
    if dist >= dist_cr:
        # 计算万有引力大小
        F = G * self.mass * other.mass / dist_sq
        # 直接通过向量分量计算加速度（方向指向天体）
        self.a[0] = F * dx / (self.mass * dist)
        self.a[1] = F * dy / (self.mass * dist)
    else:
        self.a[0] = 0
        self.a[1] = 0

3. 匹配圆周运动初始速度 #

根据万有引力提供向心力的公式v = sqrt(G*M/r)，计算合适的初始速度：

# 计算地球到太阳的初始距离
initial_dist = math.sqrt((400 - 150)**2 + (300 - 300)**2)
G_val = 10000.0
# 计算圆周运动所需初始速度
v_initial = math.sqrt(G_val * sun.mass / initial_dist)
# 初始化地球（y方向速度对应左侧位置的圆周运动方向）
earth = orb(s0=[150, 300], v0=[0, v_initial], r=10, mass=3)

修正后主循环 #

screen = pygame.display.set_mode((800, 600))
clock = pygame.time.Clock()
sun = orb(s0=[400, 300], r=70, color=[255, 0, 0], fix=True, mass=100)

# 计算匹配的初始速度
initial_dist = math.sqrt((sun.s[0] - 150)**2 + (sun.s[1] - 300)**2)
G_val = 10000.0
v_initial = math.sqrt(G_val * sun.mass / initial_dist)
earth = orb(s0=[150, 300], v0=[0, v_initial], r=10, mass=3)

running = True
while running:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False
    screen.fill((0, 0, 0))
    # 转换为整数避免Pygame绘图报错
    pygame.draw.circle(screen, earth.color, (int(earth.s[0]), int(earth.s[1])), earth.r, 0)
    pygame.draw.circle(screen, sun.color, (int(sun.s[0]), int(sun.s[1])), sun.r, 0)
    earth.comForce(sun, G=G_val)
    earth.comPos(dt=0.005)  # 减小时间步长进一步降低误差
    pygame.display.update()
    clock.tick(720)  

效果说明 #

修改后，半隐式欧拉法保持系统能量守恒，优化后的引力计算消除了三角函数误差，匹配的初始速度让地球进入稳定圆周轨道，长时间模拟轨道不会偏移。

内容的提问来源于stack exchange，提问作者연수루