太阳系八大行星模拟：使用 Python 代码模拟行星运动

import numpy as np
import matplotlib.pyplot as plt

# 常量（SI 单位）
G = 6.6743e-11    # 万有引力常数
m_s = 1.9885e30    # 太阳质量
m_m = 7.342e22    # 月球质量
m_mer = 3.301e23    # 水星质量
m_ven = 4.867e24    # 金星质量
m_ear = 5.972e24    # 地球质量
m_mar = 6.39e23    # 火星质量
m_jup = 1.898e27    # 木星质量
m_sat = 5.683e26    # 土星质量
m_ura = 8.681e25    # 天王星质量
m_nep = 1.024e26    # 海王星质量
r_s = 6.96e8    # 太阳半径
r_m = 1.737e6    # 月球半径
r_mer = 2.44e6    # 水星半径
r_ven = 6.052e6    # 金星半径
r_ear = 6.371e6    # 地球半径
r_mar = 3.39e6    # 火星半径
r_jup = 7.149e7    # 木星半径
r_sat = 6.03e7    # 土星半径
r_ura = 2.556e7    # 天王星半径
r_nep = 2.477e7    # 海王星半径

# 初始条件
# 注意：距离以米为单位，速度以米/秒为单位，时间以秒为单位。
# 数据来源：https://nssdc.gsfc.nasa.gov/planetary/planetfact.html
r0_s = np.array([0, 0])
v0_s = np.array([0, 0])
r0_mer = np.array([0.4667 * 1.496e11, 0])
v0_mer = np.array([0, 47.36e3])
r0_ven = np.array([0.7233 * 1.496e11, 0])
v0_ven = np.array([0, 35.02e3])
r0_ear = np.array([1.0 * 1.496e11, 0])
v0_ear = np.array([0, 29.78e3])
r0_mar = np.array([1.524 * 1.496e11, 0])
v0_mar = np.array([0, 24.08e3])
r0_jup = np.array([5.203 * 1.496e11, 0])
v0_jup = np.array([0, 13.07e3])
r0_sat = np.array([9.539 * 1.496e11, 0])
v0_sat = np.array([0, 9.69e3])
r0_ura = np.array([19.18 * 1.496e11, 0])
v0_ura = np.array([0, 6.81e3])
r0_nep = np.array([30.07 * 1.496e11, 0])
v0_nep = np.array([0, 5.43e3])
r0_m = r0_ear + np.array([3.844e8, 0])
v0_m = v0_ear + np.array([0, 1022])
m_solar_system=[m_mer,m_ven,m_ear,m_mar,m_jup,m_sat,m_ura,m_nep]
r0_solar_system=[r0_mer,r0_ven,r0_ear,r0_mar,r0_jup,r0_sat,r0_ura,r0_nep]
v0_solar_system=[v0_mer,v0_ven,v0_ear,v0_mar,v0_jup,v0_sat,v0_ura,v0_nep]

# 模拟参数
t0 = 0
tf = 60 * 60 * 24 * 365.25 * 20    # 模拟 20 年
dt = 60 * 60 * 24    # 时间步长为 1 天
num_steps = round((tf-t0)/dt) + 1

# 用于存储解的数组
times = np.linspace(t0, tf, num_steps)
rs = np.zeros((num_steps, 2))
rs[0] = r0_s
vs = np.zeros((num_steps, 2))
vs[0] = v0_s

# 计算重力加速度的函数
def acceleration(r, m):
    # 计算太阳和其它行星的重力加速度
    a = np.zeros(2)
    for i in range(8):
        if not np.array_equal(r, r0_solar_system[i]):
            a += G * m_s * (r0_s - r) / np.linalg.norm(r0_s - r)**3
        if not np.array_equal(r, r0_solar_system[i] + np.array([0, r_m])):
            a += G * m_m * (r0_solar_system[i] + np.array([0, r_m]) - r) / np.linalg.norm(r0_solar_system[i] + np.array([0, r_m]) - r)**3
        if not np.array_equal(r, r0_solar_system[i]+np.array([0, r0_solar_system[i][1]])):
            a += G * m_solar_system[i] * (r0_solar_system[i]+np.array([0, r0_solar_system[i][1]]) - r) / np.linalg.norm(r0_solar_system[i]+np.array([0, r0_solar_system[i][1]]) - r)**3
    return a

def acceleration_sun(r):
    # 计算太阳的重力加速度
    a = np.zeros(2)
    a += G * m_s * (r0_s - r) / np.linalg.norm(r0_s - r)**3
    return a

# 主模拟循环
for i in range(1, num_steps):
    # 使用 Verlet 方法计算新的位置和速度
    r_old = rs[i-1]
    v_old = vs[i-1]
    a_old = acceleration(r_old, m_solar_system)
    r_new = r_old + v_old * dt + 0.5 * a_old * dt**2
    a_new = acceleration(r_new, m_solar_system)
    v_new = v_old + 0.5 * (a_old + a_new) * dt
    v_new_s = v_old + 0.5 * (acceleration_sun(r_old) + acceleration_sun(r_new)) * dt
    
    # 存储新的解
    rs[i] = r_new
    vs[i] = v_new_s

# 绘制结果
fig, ax = plt.subplots(figsize=(8, 8))
ax.plot(rs[:,0], rs[:,1], 'k-', linewidth=1)
ax.plot(rs[0,0], rs[0,1], 'yo', markersize=10)
ax.plot(rs[-1,0], rs[-1,1], 'bo', markersize=10)
ax.plot(r0_mer[0], r0_mer[1], 'o', color='tab:gray', markersize=10)
ax.plot(r0_ven[0], r0_ven[1], 'o', color='tab:red', markersize=10)
ax.plot(r0_ear[0], r0_ear[1], 'o', color='tab:blue', markersize=10)
ax.plot(r0_mar[0], r0_mar[1], 'o', color='tab:orange', markersize=10)
ax.plot(r0_jup[0], r0_jup[1], 'o', color='tab:brown', markersize=10)
ax.plot(r0_sat[0], r0_sat[1], 'o', color='tab:pink', markersize=10)
ax.plot(r0_ura[0], r0_ura[1], 'o', color='tab:purple', markersize=10)
ax.plot(r0_nep[0], r0_nep[1], 'o', color='tab:green', markersize=10)
ax.set_xlim([-4e12, 4e12])
ax.set_ylim([-4e12, 4e12])
ax.set_aspect('equal', adjustable='box')
plt.axis('off')
plt.show()

代码错误原因：invalid index to scalar variable。

错误原因是：r0_m是一个二维数组，而r0_m[i]试图访问一个标量变量，而实际上r0_m是一个二维数组，不能这样访问。

如何修改：

将r0_m[i]改为r0_solar_system[i]。因为r0_solar_system是一个包含所有行星初始位置的列表，每个元素对应一个行星的二维位置，这样就可以正确访问每个行星的初始位置了。

修改后的acceleration函数如下：

def acceleration(r, m):
    # 计算太阳和其它行星的重力加速度
    a = np.zeros(2)
    for i in range(8):
        if not np.array_equal(r, r0_solar_system[i]):
            a += G * m_s * (r0_s - r) / np.linalg.norm(r0_s - r)**3
        if not np.array_equal(r, r0_solar_system[i] + np.array([0, r_m])):
            a += G * m_m * (r0_solar_system[i] + np.array([0, r_m]) - r) / np.linalg.norm(r0_solar_system[i] + np.array([0, r_m]) - r)**3
        if not np.array_equal(r, r0_solar_system[i]+np.array([0, r0_solar_system[i][1]])):
            a += G * m_solar_system[i] * (r0_solar_system[i]+np.array([0, r0_solar_system[i][1]]) - r) / np.linalg.norm(r0_solar_system[i]+np.array([0, r0_solar_system[i][1]]) - r)**3
    return a