行星轨道模拟为何输出直线而非预期椭圆轨道？

我尝试用compute_orbit方法模拟行星轨道，但绘制计算出的位置时，结果总是直线，不是预期的椭圆轨道。我试过不同初始条件，问题依然存在。

最小复现代码 #

from astroquery.jplhorizons import Horizons
import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def get_initial_conditions(planet_id):
        obj = Horizons(id=planet_id, location='@sun', epochs=2000.0)
        eph = obj.vectors()
        position = np.array([eph['x'][0], eph['y'][0], eph['z'][0]]) 
        velocity = np.array([eph['vx'][0], eph['vy'][0], eph['vz'][0]])   
        scale_factor_position = 1  
        scale_factor_velocity = 1  
        return {
            "position": np.array(position) / scale_factor_position,
            "velocity": np.array(velocity) / scale_factor_velocity
        }

def compute_orbit(central_mass=1.989e30,rotating_mass=5.972e24, dt=10000, total_time=31536000):
        initial_conditions_data = get_initial_conditions(399)
        initial_conditions = [
            initial_conditions_data['position'][0], initial_conditions_data['velocity'][0],
            initial_conditions_data['position'][1], initial_conditions_data['velocity'][1],
            initial_conditions_data['position'][2], initial_conditions_data['velocity'][2]
                            ]

        def f(t, state):
            x, vx, y, vy, z, vz = state
            r = np.sqrt(x**2 + y**2 + z**2) + 1e-5  
            G = 6.67430e-11
            Fx = -G * central_mass * rotating_mass * x / r**3
            Fy = -G * central_mass * rotating_mass * y / r**3
            Fz = -G * central_mass * rotating_mass * z / r**3
            return [vx, Fx / rotating_mass, vy, Fy / rotating_mass, vz, Fz / rotating_mass]


        t_span = (0, total_time)
        t_eval = np.arange(0, total_time, dt)  
        sol = scipy.integrate.solve_ivp(
        f, 
        t_span, 
        initial_conditions, 
        t_eval=t_eval, 
        rtol=1e-3, 
        atol=1e-6)
        x = sol.y[0]  
        y = sol.y[2]  
        z = sol.y[4]  
        positions = np.column_stack((x, y, z))
        return positions

def plot_orbit(positions):
    x = positions[:, 0]
    y = positions[:, 1]
    z = positions[:, 2]
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot(x, y, z, label='Orbit Path', color='blue')
    ax.set_xlabel('X Position (m)')
    ax.set_ylabel('Y Position (m)')
    ax.set_zlabel('Z Position (m)')
    ax.set_title('Planet Orbit Simulation')
    ax.scatter(0, 0, 0, color='yellow', s=100, label='Central Mass (Sun)')
    ax.legend()
    ax.set_box_aspect([1, 1, 1]) 
    plt.show()

positions = compute_orbit()
plot_orbit(positions)

问题根源 #

核心问题是单位不匹配：从Horizons获取的位置单位是天文单位(AU)，速度单位是AU/天，但计算时用的是国际单位制（米、秒）的引力常数G=6.67430e-11，两者未统一，导致加速度计算完全错误，轨道自然变成直线。

修复步骤 #

1. 统一单位体系 #

• 将Horizons返回的位置（AU）转换为米：1 AU = 1.496e11 米
• 将速度（AU/天）转换为米/秒：1 天 = 86400 秒

2. 修正初始条件转换 #

修改get_initial_conditions函数，加入单位转换：

def get_initial_conditions(planet_id):
    obj = Horizons(id=planet_id, location='@sun', epochs=2000.0)
    eph = obj.vectors()
    # 单位转换：AU -> 米，AU/天 -> 米/秒
    au_to_m = 1.496e11
    day_to_sec = 86400
    position = np.array([eph['x'][0], eph['y'][0], eph['z'][0]]) * au_to_m
    velocity = np.array([eph['vx'][0], eph['vy'][0], eph['vz'][0]]) * au_to_m / day_to_sec
    return {
        "position": position,
        "velocity": velocity
    }

3. 简化加速度计算（可选） #

行星质量在加速度计算中会被约掉，可直接简化公式减少计算量：

def f(t, state):
    x, vx, y, vy, z, vz = state
    r = np.sqrt(x**2 + y**2 + z**2) + 1e-5  
    G = 6.67430e-11
    # 加速度a = -G*M/r³ * r向量，省略行星质量
    ax = -G * central_mass * x / r**3
    ay = -G * central_mass * y / r**3
    az = -G * central_mass * z / r**3
    return [vx, ax, vy, ay, vz, az]

4. 调整时间参数（可选） #

将dt设为86400（1天），能让轨道曲线更平滑，total_time保持1年（≈3.154e7秒）即可模拟完整公转周期。

完整修复代码 #

from astroquery.jplhorizons import Horizons
import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def get_initial_conditions(planet_id):
    obj = Horizons(id=planet_id, location='@sun', epochs=2000.0)
    eph = obj.vectors()
    # 单位转换：AU -> 米，AU/天 -> 米/秒
    au_to_m = 1.496e11
    day_to_sec = 86400
    position = np.array([eph['x'][0], eph['y'][0], eph['z'][0]]) * au_to_m
    velocity = np.array([eph['vx'][0], eph['vy'][0], eph['vz'][0]]) * au_to_m / day_to_sec
    return {
        "position": position,
        "velocity": velocity
    }

def compute_orbit(central_mass=1.989e30, rotating_mass=5.972e24, dt=86400, total_time=3.154e7):
    initial_conditions_data = get_initial_conditions(399)
    initial_conditions = [
        initial_conditions_data['position'][0], initial_conditions_data['velocity'][0],
        initial_conditions_data['position'][1], initial_conditions_data['velocity'][1],
        initial_conditions_data['position'][2], initial_conditions_data['velocity'][2]
    ]

    def f(t, state):
        x, vx, y, vy, z, vz = state
        r = np.sqrt(x**2 + y**2 + z**2) + 1e-5  
        G = 6.67430e-11
        # 简化加速度计算
        ax = -G * central_mass * x / r**3
        ay = -G * central_mass * y / r**3
        az = -G * central_mass * z / r**3
        return [vx, ax, vy, ay, vz, az]

    t_span = (0, total_time)
    t_eval = np.arange(0, total_time, dt)  
    sol = scipy.integrate.solve_ivp(
        f, 
        t_span, 
        initial_conditions, 
        t_eval=t_eval, 
        rtol=1e-3, 
        atol=1e-6
    )
    x = sol.y[0]  
    y = sol.y[2]  
    z = sol.y[4]  
    positions = np.column_stack((x, y, z))
    return positions

def plot_orbit(positions):
    x = positions[:, 0]
    y = positions[:, 1]
    z = positions[:, 2]
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot(x, y, z, label='Orbit Path', color='blue')
    ax.set_xlabel('X Position (m)')
    ax.set_ylabel('Y Position (m)')
    ax.set_zlabel('Z Position (m)')
    ax.set_title('Planet Orbit Simulation')
    ax.scatter(0, 0, 0, color='yellow', s=100, label='Central Mass (Sun)')
    ax.legend()
    ax.set_box_aspect([1, 1, 1]) 
    plt.show()

positions = compute_orbit()
plot_orbit(positions)

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