python 三体运动
时间: 2023-10-20 20:35:48 浏览: 72
Python中可以使用科学计算库NumPy来模拟三体运动。以下是一个简单的示例代码:
```python
import numpy as np
import matplotlib.pyplot as plt
# 定义常数
G = 6.67430e-11 # 万有引力常数
# 定义初始条件
m1 = 5.972e24 # 地球质量
m2 = 7.348e22 # 月球质量
m3 = 1.989e30 # 太阳质量
r1 = np.array([0, 0]) # 地球位置
r2 = np.array([384400000, 0]) # 月球位置
r3 = np.array([0, 0]) # 太阳位置
v1 = np.array([0, 0]) # 地球速度
v2 = np.array([0, np.sqrt(G * m3 / np.linalg.norm(r2 - r3))]) # 月球速度
v3 = np.array([0, -np.sqrt(G * m2 / np.linalg.norm(r1 - r2))]) # 太阳速度
# 定义时间步长和总时长
dt = 3600 # 时间步长(秒)
total_time = 86400 * 30 # 总时长(秒)
# 初始化轨迹数组
r1_track = [r1]
r2_track = [r2]
r3_track = [r3]
# 进行模拟
for t in range(0, total_time, dt):
# 计算加速度
a1 = G * (m2 * (r2 - r1) / np.linalg.norm(r2 - r1)**3 + m3 * (r3 - r1) / np.linalg.norm(r3 - r1)**3)
a2 = G * (m1 * (r1 - r2) / np.linalg.norm(r1 - r2)**3 + m3 * (r3 - r2) / np.linalg.norm(r3 - r2)**3)
a3 = G * (m1 * (r1 - r3) / np.linalg.norm(r1 - r3)**3 + m2 * (r2 - r3) / np.linalg.norm(r2 - r3)**3)
# 更新位置和速度
r1 += v1 * dt
v1 += a1 * dt
r2 += v2 * dt
v2 += a2 * dt
r3 += v3 * dt
v3 += a3 * dt
# 记录轨迹
r1_track.append(r1.copy())
r2_track.append(r2.copy())
r3_track.append(r3.copy())
# 绘制轨迹
r1_track = np.array(r1_track)
r2_track = np.array(r2_track)
r3_track = np.array(r3_track)
plt.plot(r1_track[:, 0], r1_track[:, 1], label='Earth')
plt.plot(r2_track[:, 0], r2_track[:, 1], label='Moon')
plt.plot(r3_track[:, 0], r3_track[:, 1], label='Sun')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Three-Body Motion')
plt.legend()
plt.show()
```
该代码使用欧拉方法进行数值模拟,通过计算每个物体的加速度并更新位置和速度。最后使用Matplotlib库将轨迹绘制出来。
请注意,这只是一个简单的示例,并且在实际的三体运动中可能存在更多复杂性和数值稳定性问题。如果需要更精确的模拟,可以考虑使用更高阶的数值积分方法,例如四阶龙格-库塔法(RK4)等。