四阶Yoshida积分模拟卫星轨道出现螺旋外扩问题咨询

本人尝试使用4阶Yoshida积分技术对绕地运行的圆轨道卫星轨道进行建模。

模拟得到的轨道很快出现向外螺旋偏移的问题。以下为类月卫星的模拟实现代码：使用欧拉法时粒子运动表现正常，但希望采用精度更高的积分方法，因此推测问题可能出在算法本身的实现逻辑上。

已尝试直接使用引力参数代替手动计算G*M，问题并未解决；此外还尝试了减小时间步长、调整单位制、打印校验各积分步骤的数值等方法，均未定位问题原因，核心疑问为对该算法的使用方式是否正确。

G = 6.674e-20 # km^3 kg^-1 s^-2
day = 60.0 * 60.0 * 24.0 # length of day
dt = day / 10.0
M = 5.972e24 # kg

N = 1
delta = np.random.random(1) * 2.0 * np.pi / N
angles = np.linspace(0.0, 2.0 * np.pi, N) + delta
rad = np.random.uniform(low = 384e3, high = 384e3, size = (N))

x, y = rad * np.cos(angles), ringrad * np.sin(angles)
vx, vy = np.sqrt(G*M / rad) * -np.sin(angles), np.sqrt(G*M / rad) * np.cos(angles)

def update(frame):
    global x, y, vx, vy, dt, day

    positions.set_data(x, y)

    # coefficients 
    q = 2**(1/3)
    w1 = 1 / (2 - q)
    w0 = -q * w1
    d1 = w1
    d3 = w1
    d2 = w0
    c1 = w1 / 2
    c2 = (w0 + w1) / 2
    c3 = c2
    c4 = c1

    # Step 1

    x1 = x + c1*vx*dt
    y1 = y + c1*vy*dt
    dist1 = np.hypot(x1, y1)
    acc1 = -(G*M) / (dist1**2.0)
    dx1 = x1 - x
    dy1 = y1 - y
    accx1 = (acc1*dx1)/(x1)
    accy1 = (acc1*dy1)/(y1)
    vx1 = vx + d1*accx1*dt
    vy1 = vy + d1*accy1*dt

    # Step 2

    x2 = x1 + c2*vx1*dt
    y2 = y1 + c2*vy1*dt
    dist2 = np.hypot(x2, y2)
    acc2 = -(G*M) / (dist2**2.0)
    dx2 = x2 - x1
    dy2 = y2 - y1
    accx2 = (acc2*dx2)/(x2)
    accy2 = (acc2*dy2)/(y2)
    vx2 = vx1 + d2*accx2*dt
    vy2 = vy1 + d2*accy2*dt

    # Step 3

    x3 = x2 + c3*vx2*dt
    y3 = y2 + c3*vy2*dt
    dist3 = np.hypot(x3, y3)
    acc3 = -(G*M) / (dist3**2.0)
    dx3 = x3 - x2
    dy3 = y3 - y2
    accx3 = (acc3*dx3)/(x3)
    accy3 = (acc3*dy3)/(y3)
    vx3 = vx2 + d3*accx3*dt
    vy3 = vy2 + d3*accy3*dt

    # Full step

    x = x3 + c4*vx3*dt
    y = y3 + c4*vy3*dt
    vx = vx3
    vy = vy3

    return positions

问题原因与修正方案 #

你的Yoshida积分系数和整体步骤框架是对的，存在2个核心错误和1个低级笔误：

• 引力加速度计算完全错误
  引力是指向坐标原点（地心）的，加速度分量应该直接按当前位置矢量计算：

  # 正确加速度计算，r_i为当前位置到原点的距离
  a_x = - G * M * x_i / (r_i ** 3)
  a_y = - G * M * y_i / (r_i ** 3)
  

  你代码里用两次位置的差dx_i = x_i - x_{i-1}（即上一小步的漂移位移）计算加速度方向，甚至除以当前坐标值x_i做归一化，得到的加速度方向根本不指向地心，大小也完全错误，是轨道偏移的核心原因。
• 代码笔误直接导致运行错误
  初始化y坐标时你写的是ringrad * np.sin(angles)，但前面根本没有定义ringrad变量，应该替换为你定义的轨道半径rad。
• （非致命优化项）Yoshida系数是固定值，不需要每帧在update函数里重复计算，提到全局初始化一次即可，不影响精度但能减少无意义计算。

修正后的核心步进代码参考 #

# 提前算好固定系数，不用每帧重复计算
q = 2**(1/3)
w1 = 1 / (2 - q)
w0 = -q * w1
d1, d2, d3 = w1, w0, w1
c1 = w1 / 2
c2 = (w0 + w1) / 2
c3, c4 = c2, c1

def update(frame):
    global x, y, vx, vy
    positions.set_data(x, y)

    # Step 1: 漂移c1*dt
    x1 = x + c1*vx*dt
    y1 = y + c1*vy*dt
    # 计算当前位置引力加速度
    dist1_sq = x1**2 + y1**2
    dist1 = np.sqrt(dist1_sq)
    accx1 = - G*M * x1 / (dist1_sq * dist1)
    accy1 = - G*M * y1 / (dist1_sq * dist1)
    # 踢d1*dt更新速度
    vx1 = vx + d1*accx1*dt
    vy1 = vy + d1*accy1*dt

    # Step 2: 漂移c2*dt
    x2 = x1 + c2*vx1*dt
    y2 = y1 + c2*vy1*dt
    dist2_sq = x2**2 + y2**2
    dist2 = np.sqrt(dist2_sq)
    accx2 = - G*M * x2 / (dist2_sq * dist2)
    accy2 = - G*M * y2 / (dist2_sq * dist2)
    vx2 = vx1 + d2*accx2*dt
    vy2 = vy1 + d2*accy2*dt

    # Step3: 漂移c3*dt
    x3 = x2 + c3*vx2*dt
    y3 = y2 + c3*vy2*dt
    dist3_sq = x3**2 + y3**2
    dist3 = np.sqrt(dist3_sq)
    accx3 = - G*M * x3 / (dist3_sq * dist3)
    accy3 = - G*M * y3 / (dist3_sq * dist3)
    vx3 = vx2 + d3*accx3*dt
    vy3 = vy2 + d3*accy3*dt

    # 最后一步漂移c4*dt得到新位置
    x = x3 + c4*vx3*dt
    y = y3 + c4*vy3*dt
    vx, vy = vx3, vy3

    return positions

修正后轨道不会再出现螺旋偏移，4阶Yoshida积分的辛特性可以保证轨道长期稳定，能量误差不会随时间累积。

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