地球绕太阳Python仿真力值计算不准确问题求助

你的引力矢量计算是正确的，轨道半径持续变化的核心问题出在数值积分方法的选择和时间步长设置上，具体分析和修复方案如下：

问题根源 #

1. 显式欧拉法的固有缺陷 #

你当前使用的是最简单的显式欧拉积分，步骤为：先通过当前力更新速度，再用新速度更新位置。这种方法无法保持保守系统（如引力场）的能量守恒——每一步计算都会引入微小的能量误差，随着循环次数增加，误差不断积累，最终表现为轨道半径持续膨胀或收缩。

2. 时间步长设置不合理 #

TIME_STEP=0.01是秒级步长，模拟地球公转一周需要运行约30亿次循环，不仅效率极低，而且即使是小步长，欧拉法的误差依然会逐步累积，导致轨道明显偏离圆周。

3. 太阳位置固定（次要） #

实际物理系统中，地球对太阳也有引力作用，太阳会产生微小位移。虽然因太阳质量巨大，该位移可忽略，但完全固定太阳位置会引入微小的系统误差。

修复方案 #

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

辛积分方法能近似保持系统的能量守恒，是天体轨道模拟的最优选择。以下是两种实现方案：

半隐式欧拉法（简单高效）

修改update_position方法，调整计算顺序，先更新位置再计算新位置下的力来更新速度：

TIME_STEP = 86400  # 改用1天为时间步长，平衡精度与效率
G = 6.67430e-11

class Planet:
    def __init__(self, x_position, y_position, x_velocity, y_velocity, mass):
        self.r_position = np.array([x_position, y_position], dtype='float64')
        self.velocity = np.array([x_velocity, y_velocity], dtype='float64')
        self.mass = mass

    def update_position(self, other):
        # 1. 用当前速度更新位置
        self.r_position += self.velocity * TIME_STEP
        
        # 2. 基于新位置计算引力矢量
        delta_r = self.r_position - other.r_position
        r_norm = np.linalg.norm(delta_r)
        force = -G * self.mass * other.mass * delta_r / (r_norm ** 3)
        
        # 3. 用新引力更新速度
        self.velocity += force / self.mass * TIME_STEP

Verlet积分法（更高精度）

Verlet积分通过前后位置的差值计算速度，能量守恒效果更优：

TIME_STEP = 86400
G = 6.67430e-11

class Planet:
    def __init__(self, x_position, y_position, x_velocity, y_velocity, mass):
        self.r_position = np.array([x_position, y_position], dtype='float64')
        # 初始化上一时刻位置（Verlet积分需要）
        self.r_prev = self.r_position - self.velocity * TIME_STEP
        self.mass = mass

    def update_position(self, other):
        # 计算当前位置的引力加速度
        delta_r = self.r_position - other.r_position
        r_norm = np.linalg.norm(delta_r)
        acceleration = -G * other.mass * delta_r / (r_norm ** 3)
        
        # Verlet核心公式：计算下一时刻位置
        r_next = 2 * self.r_position - self.r_prev + acceleration * TIME_STEP ** 2
        
        # 推导当前速度（用于可视化或后续计算）
        self.velocity = (r_next - self.r_prev) / (2 * TIME_STEP)
        
        # 更新位置历史
        self.r_prev = self.r_position
        self.r_position = r_next

2. 合理设置时间步长 #

将时间步长改为86400秒（1天），既能保证模拟效率，又能有效控制误差积累。若需要更高精度，可尝试14400秒（4小时）或更短，但避免使用秒级步长。

3. 同步更新太阳位置（可选） #

若追求完全真实的物理模拟，每次循环同时更新太阳和地球的位置：

while running:
    earth.update_position(sun)
    sun.update_position(earth)

由于太阳质量远大于地球，其位移会非常微小，但该调整能让系统更符合牛顿万有引力的相互作用规律。

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