使用Python的ODEINT模拟行星轨道：耦合ODE绘图异常求助

问题背景 #

尝试用scipy.integrate.odeint求解四耦合一阶ODE模拟地球绕太阳的轨道，但生成的图像呈线性（不符合预期的椭圆轨道），或出现超时错误。

现有代码 #

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

m = 6*(10**24)  # 地球质量
M = 2*(10**30)  # 太阳质量
G = 6.8*(10**(-11))  # 引力常数（值错误）

def orbit(f, t):
    x, y, vx, vy = f
    r = np.sqrt((x**2)+(y**2))  # 需用numpy的sqrt
    # v变量未使用，可删除
    dx = vx
    dy = vy
    # 加速度公式错误：应为-G*M*x/r³，而非-G*m*x/r³
    dfx = -G*m*x/(r**3)
    dfy = -G*m*y/(r**3)
    dFdt = [dx, dy, dfx, dfy]
    return dFdt

t = np.linspace(0, 3.2*(10**7))  # 约1年的时间范围
x0 = 0
y0 = 1.5*(10**11)  # 初始轨道半径（日地距离）
vx0 = 3*(10**4)  # 初始切向速度
vy0 = 0

f = odeint(orbit, (x0,y0, vx0, vy0), t)

plt.plot(t, f[:, 0])
plt.title("X随时间变化")
plt.figure()
plt.plot(f[:,1], f[:,0])
plt.title("Y-X相图（轨道）")
plt.show()

当前错误现象 #

• X随时间变化图呈线性：
  
• Y-X相图呈直线：
  

问题根源与修正方案 #

1. 引力加速度公式错误 #

核心错误：地球绕太阳运动的加速度由太阳引力提供，公式应为：
$$a_x = -\frac{G M x}{r^3}, \quad a_y = -\frac{G M y}{r^3}$$
原代码中误用了地球质量m代替太阳质量M，导致加速度量级差了约5个数量级（$M/m \approx 3.3\times10^5$），地球几乎不受引力，做匀速直线运动，因此图像呈线性。

2. 引力常数值错误 #

标准引力常数为$G=6.674\times10^{{-11} \text{N·m}}2/\text{kg}^2$，原代码使用6.8e-11存在偏差。

3. 缺少必要依赖导入与细节问题 #

• 未导入numpy和matplotlib.pyplot，实际运行会报错；
• 计算r时需用np.sqrt而非原生sqrt（避免标量运算限制）；
• 相图需设置plt.axis('equal')，否则椭圆会被拉伸成直线。

修正后的代码 #

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

m = 6e24  # 地球质量
M = 2e30  # 太阳质量
G = 6.674e-11  # 标准引力常数

def orbit(f, t):
    x, y, vx, vy = f
    r = np.sqrt(x**2 + y**2)
    dx = vx
    dy = vy
    # 修正加速度公式，使用太阳质量M
    dfx = -G * M * x / (r**3)
    dfy = -G * M * y / (r**3)
    return [dx, dy, dfx, dfy]

# 时间范围：0到1年（约3.15e7秒）
t = np.linspace(0, 3.15e7, 1000)  # 增加采样点，使轨道更平滑
x0 = 0
y0 = 1.5e11  # 日地距离
vx0 = 3e4  # 地球公转速度
vy0 = 0

f = odeint(orbit, (x0, y0, vx0, vy0), t)

# 绘制X随时间变化
plt.figure(figsize=(8,4))
plt.plot(t/86400, f[:, 0])  # 转换为天为单位
plt.xlabel("时间（天）")
plt.ylabel("X位置（m）")
plt.title("地球X位置随时间变化")

# 绘制轨道相图
plt.figure(figsize=(6,6))
plt.plot(f[:,1], f[:,0], label="地球轨道")
plt.scatter(0, 0, color='orange', s=100, label="太阳")
plt.xlabel("Y位置（m）")
plt.ylabel("X位置（m）")
plt.title("地球绕太阳轨道")
plt.axis('equal')
plt.legend()
plt.show()

修正后效果 #

修正后会生成符合预期的椭圆轨道，X随时间变化呈正弦曲线。

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