Python模拟地球椭圆轨道运动遇问题：速度未按预期变化

问题分析 #

你的代码无法形成椭圆轨道的核心原因：

• 用固定大小的速度指向太阳，本质是直线运动，完全不符合轨道运动的物理规律
• 没有引入万有引力作为向心力，轨道运动的核心是引力提供向心力，迫使物体做曲线运动
• 忽略了开普勒第二定律（面积速度守恒），椭圆轨道的速度必须随日地距离变化：近地点速度快，远地点速度慢

实现思路 #

要准确模拟椭圆轨道，需要做到三点：

1. 用万有引力公式计算加速度：$a = \frac{GM}{r^2}$，方向始终指向太阳（G为引力常数，M为太阳质量，r为日地距离）
2. 采用数值积分方法更新速度与位置（推荐Verlet法，比基础欧拉法精度更高、稳定性更好）
3. 设置符合物理规律的初始条件：给地球一个垂直于日地连线的初始速度，结合初始距离生成椭圆轨道

代码示例（Verlet法，高精度轨道模拟） #

import math
import matplotlib.pyplot as plt

# 简化物理参数（便于计算，数值接近真实值）
G = 6.67430e-11  # 引力常数
M_SUN = 1.989e30  # 太阳质量（kg）
DT = 86400  # 时间步长：1天（秒）
SIM_DAYS = 365  # 模拟时长：1年

# 初始位置：地球位于远日点（距离太阳约1.521e11米）
sun_pos = [0.0, 0.0]
earth_pos = [1.521e11, 0.0]
# 初始速度：垂直于日地连线，远日点速度约29290 m/s（符合开普勒定律）
earth_vel = [0.0, 29290.0]

# 存储轨迹用于可视化
x_track = []
y_track = []

for _ in range(SIM_DAYS):
    # 计算当前日地距离与单位方向向量
    dx = earth_pos[0] - sun_pos[0]
    dy = earth_pos[1] - sun_pos[1]
    r = math.hypot(dx, dy)
    unit_x = dx / r
    unit_y = dy / r

    # 计算当前引力加速度（指向太阳，所以取负）
    accel_mag = G * M_SUN / (r ** 2)
    accel_x = -unit_x * accel_mag
    accel_y = -unit_y * accel_mag

    # Verlet法更新位置
    earth_pos[0] += earth_vel[0] * DT + 0.5 * accel_x * DT**2
    earth_pos[1] += earth_vel[1] * DT + 0.5 * accel_y * DT**2

    # 计算新位置下的加速度
    dx_new = earth_pos[0] - sun_pos[0]
    dy_new = earth_pos[1] - sun_pos[1]
    r_new = math.hypot(dx_new, dy_new)
    unit_x_new = dx_new / r_new
    unit_y_new = dy_new / r_new
    accel_mag_new = G * M_SUN / (r_new ** 2)
    accel_x_new = -unit_x_new * accel_mag_new
    accel_y_new = -unit_y_new * accel_mag_new

    # Verlet法更新速度
    earth_vel[0] += 0.5 * (accel_x + accel_x_new) * DT
    earth_vel[1] += 0.5 * (accel_y + accel_y_new) * DT

    # 记录轨迹点
    x_track.append(earth_pos[0])
    y_track.append(earth_pos[1])

# 可视化轨道
plt.figure(figsize=(8, 8))
plt.plot(x_track, y_track, label='Earth Orbit')
plt.scatter(sun_pos[0], sun_pos[1], color='orange', s=200, label='Sun')
plt.xlabel('X Position (m)')
plt.ylabel('Y Position (m)')
plt.title('Earth\'s Elliptical Orbit Around the Sun')
plt.legend()
plt.axis('equal')
plt.show()

关键细节说明 #

• Verlet法优势：通过前后两次加速度的平均值更新速度，能有效减少长期模拟的误差，避免轨道逐渐偏离或崩溃
• 初始速度设置：椭圆轨道的初始速度需满足公式$v = \sqrt{GM(\frac{2}{r} - \frac{1}{a})}$（a为椭圆半长轴），速度过小会导致地球撞向太阳，过大则会脱离轨道
• 开普勒定律自动满足：由于引力是有心力，系统角动量守恒，自然会出现近快远慢的速度变化，符合开普勒第二定律

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