穹頂天魂提示您：看後求收藏（奇妙書庫www.qmshu.tw），接著再看更方便。
# 質量 (kg)
m1 = 5.972e24 # 地球質量
m2 = 7.348e22 # 月球質量
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)
# 將位置結果提取出來
r1_solution = solution.y[0:2, :]
r2_solution = solution.y[4:6, :]
r3_solution = solution.y[8:10, :]
# 軌跡繪圖
plt.plot(r1_solution, r1_solution, label='Earth')
plt.plot(r2_solution, r2_solution, label='moon')
plt.plot(r3_solution, r3_solution, label='Sun')
plt.legend plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.title('three-body problem Simulation')
plt.savefig('three_body_simulation.png')
這段程式碼執行後，將會生成一個以地球
