火卫一（Phobos）绕火星（Mars）轨道运动异常：呈直线运动求助

问题描述 #

模拟中火卫一未绕火星做轨道运动，而是沿直线运动，尝试过调整初始速度但未解决，采用线性化方法从力推导位置，无法定位错误。

代码示例 #

import math
import matplotlib.pyplot as plt

G = 6.6743*(10**-11)
MMars = 0.64169 * (10**24)     #kg
MPhobos = 10.6 * (10**15)      #kg
#Vxphobos = 2138 
#Vxphobos = 7696700 
Vxphobos = 0
#Vxphobos = 3520
Vyphobos = 4100 
mars = [100,100]
phobos = [mars[1]-(6*(10**3)),0]
t = 0
phobosunit = [0,0]
force = []  # 补充初始化列表，避免报错

def distt(distx, disty):
    totaldist = math.sqrt(distx**2 + disty**2)
    return totaldist 

while t < 57552: # 火卫一公转周期（秒）
    tchange = 300
    t += tchange

    phobosr = distt(mars[0]-phobos[0], mars[1]-phobos[1])
    # 冗余计算：phobosdist与phobosr完全相同，直接复用即可
    phobosdist = phobosr
    for i in range(2):
        phobosunit[i] = -(mars[i]-phobos[i])/phobosdist
    # plt.quiver(...) 可以注释，避免绘图混乱
    axphobos = (G*MMars/(phobosr**2))*phobosunit[0]
    ayphobos = (G*MMars/(phobosr**2))*phobosunit[1]
    
    Vxphobos += axphobos*tchange
    Vyphobos += ayphobos*tchange
    phobos[0] += Vxphobos*tchange
    phobos[1] += Vyphobos*tchange
    plt.plot(phobos[0], phobos[1], 'go') 
    # print(phobosunit)  # 调试用，可注释
    force.append(G*MMars*MPhobos/(phobosr**2))

e = 10**8
plt.plot(mars[0], mars[1], 'ro')
plt.xticks([-1*e,0,1*e,2*e,3*e])
plt.yticks([-1*e,0,1*e,2*e,3*e])
plt.show()

关键错误排查与解决思路 #

• 初始位置与轨道半径不匹配：火卫一实际轨道半径约9376km（9.376×10⁶米），但代码中初始位置仅偏移6×10³米，数量级偏差极大，导致引力和速度的比例完全错误。建议将初始位置设置为phobos = [mars[0] - 9.376e6, mars[1]]，保证初始距离符合真实轨道参数。
• 初始速度计算错误：近圆轨道的速度公式为v = √(GM/r)，代入火星质量和轨道半径，计算得火卫一的轨道速度约为2138m/s（对应你注释掉的Vxphobos=2138），但当前代码中Vxphobos=0，缺少横向速度，火卫一自然会沿直线落向火星。需启用正确的初始横向速度。
• 时间步长过大：tchange=300秒的步长对于欧拉法来说太大，线性近似的误差会快速累积，导致轨道完全偏离。建议将步长缩小到1~10秒，提升模拟精度。
• 欧拉法的局限性：基础欧拉法无法很好地保持能量守恒，长期模拟会导致轨道发散。推荐改用半隐式欧拉法或Verlet积分，例如更新位置时使用平均速度：phobos[0] += (Vxphobos + Vxphobos + axphobos*tchange)/2 * tchange，减少误差。
• 冗余计算与逻辑错误：代码中phobosdist的计算完全冗余，直接复用phobosr即可；初始位置的x坐标误用了火星的y坐标（mars[1]），逻辑错误导致初始位置偏离预期。

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