pycharm实现三体运动
时间: 2023-05-16 12:07:36 浏览: 98
可以使用Python的matplotlib库来实现三体运动的可视化。首先,需要定义三个质点的初始位置、速度和质量等参数,然后使用牛顿万有引力定律计算它们之间的相互作用力,再根据牛顿第二定律计算它们的加速度,最后使用欧拉法或者其他数值积分方法来更新它们的位置和速度。以下是一个简单的示例代码:
```python
import numpy as np
import matplotlib.pyplot as plt
# 定义三个质点的初始参数
m1, m2, m3 = 1, 1, 1
x1, y1, vx1, vy1 = -1, 0, 0, 0.5
x2, y2, vx2, vy2 = 1, 0, 0, -0.5
x3, y3, vx3, vy3 = 0, np.sqrt(3), -0.5, 0
# 定义时间步长和总时间
dt = 0.01
t = np.arange(0, 10, dt)
# 定义牛顿万有引力定律
def force(x1, y1, x2, y2, m1, m2):
r = np.sqrt((x1-x2)**2 + (y1-y2)**2)
f = m1 * m2 / r**2
fx = f * (x2-x1) / r
fy = f * (y2-y1) / r
return fx, fy
# 定义欧拉法更新位置和速度的函数
def update(x1, y1, vx1, vy1, x2, y2, vx2, vy2, x3, y3, vx3, vy3, dt):
fx12, fy12 = force(x1, y1, x2, y2, m1, m2)
fx13, fy13 = force(x1, y1, x3, y3, m1, m3)
fx23, fy23 = force(x2, y2, x3, y3, m2, m3)
ax1, ay1 = (fx12+fx13) / m1, (fy12+fy13) / m1
ax2, ay2 = (-fx12+fx23) / m2, (-fy12+fy23) / m2
ax3, ay3 = (-fx13-fx23) / m3, (-fy13-fy23) / m3
x1, y1, vx1, vy1 = x1+vx1*dt, y1+vy1*dt, vx1+ax1*dt, vy1+ay1*dt
x2, y2, vx2, vy2 = x2+vx2*dt, y2+vy2*dt, vx2+ax2*dt, vy2+ay2*dt
x3, y3, vx3, vy3 = x3+vx3*dt, y3+vy3*dt, vx3+ax3*dt, vy3+ay3*dt
return x1, y1, vx1, vy1, x2, y2, vx2, vy2, x3, y3, vx3, vy3
# 循环更新位置和速度,并绘制轨迹
for i in range(len(t)):
x1, y1, vx1, vy1, x2, y2, vx2, vy2, x3, y3, vx3, vy3 = update(x1, y1, vx1, vy1, x2, y2, vx2, vy2, x3, y3, vx3, vy3, dt)
plt.plot(x1, y1, 'ro', markersize=3)
plt.plot(x2, y2, 'bo', markersize=3)
plt.plot(x3, y3, 'go', markersize=3)
plt.show()
```
这段代码使用了欧拉法来更新位置和速度,因此可能存在数值误差和能量不守恒等问题。如果需要更精确的模拟,可以使用其他数值积分方法,比如四阶龙格-库塔法。
阅读全文