我正在模拟一个由三个物体在重力作用下相互作用的程序，但代码不起作用，我不知道为什么。我的计算基于速度verlet 方法。我实施对了吗？ （我很抱歉问这个问题，但我刚刚开始了解它并且有点不确定）。感谢您的任何意见或建议！
import pygame
import math
G = 6.67430e-11 # gravitational constant
M = [(5.9721910**(24)), (7.34910**(22)), (0.017.34910**(22))] # masses of the planets
R = [(450, 325), (450,(575/11)), (450,170)] # initial positions of the planets
V = [(-24.947719394204714/2,0), (1023,0), (1000,0)] # initial velocities of the planets
# Define simulation parameters
t_start = 0 # start time
t_end = 31536000 # end time
dt = 0.01 # time step
N = int((t_end - t_start) / dt) # number of time steps
scale = 1/(1409466.667) #converts meters to pixels
pygame.init()
screen = pygame.display.set_mode((900, 650))
pygame.display.set_caption("Three-Body Problem")
def acceleration(r):
r12 = [r[1][0] - r[0][0], r[1][1] - r[0][1]]
r23 = [r[2][0] - r[1][0], r[2][1] - r[1][1]]
r31 = [r[0][0] - r[2][0], r[0][1] - r[2][1]]
a = [[0, 0], [0, 0], [0, 0]]
r12_norm = math.hypot(*r12)
if r12_norm > 0:
a[0][0] = G * M[1] * r12[0] / r12_norm**3 + G * M[2] * r31[0] / math.hypot(*r31)**3
a[0][1] = G * M[1] * r12[1] / r12_norm**3 + G * M[2] * r31[1] / math.hypot(*r31)**3
r23_norm = math.hypot(*r23)
if r23_norm > 0:
a[1][0] = G * M[0] * -r12[0] / r12_norm**3 + G * M[2] * r23[0] / r23_norm**3
a[1][1] = G * M[0] * -r12[1] / r12_norm**3 + G * M[2] * r23[1] / r23_norm**3
r31_norm = math.hypot(*r31)
if r31_norm > 0:
a[2][0] = G * M[0] * r31[0] / r31_norm**3 + G * M[1] * -r23[0] / r23_norm**3
a[2][1] = G * M[0] * r31[1] / r31_norm**3 + G * M[1] * -r23[1] / r23_norm**3
return a
def velocity_verlet(r, v, a, dt):
r_new = [[0, 0], [0, 0], [0, 0]]
v_new = [[0, 0], [0, 0], [0, 0]]
for i in range(3):
r_new[i][0] = r[i][0] + v[i][0] * dt + 0.5 * a[i][0] * (dt**2)
r_new[i][1] = r[i][1] + v[i][1] * dt + 0.5 * a[i][1] * (dt**2)
a_new = acceleration(r_new)
for i in range(3):
v_new[i][0] = v[i][0] + 0.5 * (a[i][0] + a_new[i][0]) * dt
v_new[i][1] = v[i][1] + 0.5 * (a[i][1] + a_new[i][1]) * dt
return r_new, v_new, a_new
t = t_start
r = R
v = V
a = acceleration(r)
for i in range(N): # Clear the screen
screen.fill((255, 255, 255))
# Draw the planets
for j in range(3):
pygame.draw.circle(screen, (0, 0, 0), (int(r[j][0])*scale, int(r[j][1])*scale), 10)
pygame.display.update()
# Update the positions and velocities
r, v, a = velocity_verlet(r, v, a, dt)
t += dt
for event in pygame.event.get():
if event.type == pygame.QUIT:
pygame.quit()
quit()
pygame.quit()