Python实现Leapfrog蛙跳法模拟地球轨道 轨迹显示不全问题

你的代码跑不出完整轨道和运动方程、matplotlib逻辑无关，是3个非常低级的书写错误：

• 初速度笔误：地球近日点公转速度约30.287km/s，对应值是3.0287*10**4 m/s，你写成了3.0287*20**4，算出来速度接近4.85e8 m/s，比光速还快，会直接让地球沿直线飞出太阳系。
• 模拟时长单位不匹配：积分步长h=3600的单位是秒（对应1小时），但你写的t_stop=24*365*5是5年对应的总小时数，没有转换成秒（少乘了3600），实际模拟时长只有12小时左右，根本走不完一段可观测的弧段。
• 绘图没开等比例坐标轴：就算前面的参数对，默认matplotlib的x/y轴比例不一致，椭圆轨道会被拉成奇怪的形状。
• 额外冗余：你导入了scipy.integrate但全程没调用，属于无效导入。

import matplotlib.pyplot as plt
import numpy as np

# 物理参数
G = 6.6738e-11
M = 1.9891e30
h = 3600  # 积分步长：1小时，单位秒
r_perihelion = 1.4710e11  # 近日点距离，单位米
v_perihelion = 3.0287e4   # 近日点切向速度，修正原20**4的笔误

def LeapFrog(f, t_start, t_stop, z0, h):
    t_vec = np.arange(t_start, t_stop, h) 
    n = len(t_vec)
    state_dim = len(z0)  
    z_vec = np.zeros((n, state_dim))
    z_vec[0, :] = z0
    # 半步初始速度
    z_half = z_vec[0, :] + 0.5 * h * f(z0, t_vec[0]) 
    
    for i in range(n - 1):
        z_vec[i+1, :] = z_vec[i, :] + h * f(z_half, t_vec[i] + 0.5*h)
        z_half += h * f(z_vec[i+1, :], t_vec[i] + h)   
    return t_vec, z_vec

def orbit_ode(z, t):   
    x, y, vx, vy = z
    r = np.sqrt(x**2 + y**2)
    dz = np.zeros(4)
    dz[0] = vx 
    dz[1] = vy
    dz[2] = -G * M * x / r**3
    dz[3] = -G * M * y / r**3
    return dz

t_start = 0
# 模拟总时长5年，统一转换成秒单位
t_stop = 5 * 365 * 24 * 3600
# 初始状态：近日点位于(0, r_perihelion)，初速度沿x轴正方向
z0 = np.array([0, r_perihelion, v_perihelion, 0])
t_vec, z_vec = LeapFrog(orbit_ode, t_start, t_stop, z0, h)

# 绘图
plt.figure(figsize=(8,8))
plt.plot(z_vec[:,0], z_vec[:,1], 'g', linewidth=1, label='地球轨道')
plt.plot(0, 0, 'yo', markersize=10, label='太阳位置')
plt.axis('equal')  # 强制坐标轴等比例，避免轨道变形
plt.legend()
plt.show()

• 蛙跳法是辛积分，对哈密顿系统的能量守恒性很好，如果你觉得步长1小时跑5年太慢，可以把步长改成1天（h=86400），轨道依然会保持稳定闭合，不会出现普通欧拉法的轨道漂移问题。
• 如果想更清楚看到近日点远日点的差异，可以在初始条件里把初速度调小一点，轨道偏心率会更高。

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