m3 = 1.989e30 # 太阳质量

# 初始位置 (m)

r1_0 = np.array([0, 0])

r2_0 = np.array([3.844e8, 0]) # 月球距地球的平均距离

r3_0 = np.array([1.496e11, 0]) # 地球距太阳的平均距离

小主，

# 初始速度 (m/s)

v1_0 = np.array([0, 0])

v2_0 = np.array([0, 1022]) # 月球轨道速度

v3_0 = np.array([0, ]) # 地球轨道速度

initial_conditions = np.concatenate([r1_0, v1_0, r2_0, v2_0, r3_0, v3_0])

# 定义运动方程

def equations(t, y):

r1 = y[0:2]

v1 = y[2:4]

r2 = y[4:6]

v2 = y[6:8]

r3 = y[8:10]

v3 = y[10:12]

r12 = np.linalg.norm(r2 - r1)

r13 = np.linalg.norm(r3 - r1)

r23 = np.linalg.norm(r3 - r2)

dv1_dt = G * m2 * (r2 - r1) / r12**3 + G * m3 * (r3 - r1) / r13**3

dv2_dt = G * m1 * (r1 - r2) / r12**3 + G * m3 * (r3 - r2) / r23**3

dv3_dt = G * m1 * (r1 - r3) / r13**3 + G * m2 * (r2 - r3) / r23**3

dr1_dt = v1

dr2_dt = v2

dr3_dt = v3

return np.concatenate([dr1_dt, dv1_dt, dr2_dt, dv2_dt, dr3_dt, dv3_dt])

# 数值积分

t_span = (0, 3.154e7) # 积分时间为一年

t_eval = np.linspace(0, 3.154e7, 1000) # 1000个时间点

solution = solve_ivp(equations, t_span, initial_conditions, t_eval=t_eval)

# 将位置结果提取出来