As the title states, I want to use the velocity Verlet algorithm to simulate the trajectories of 3 bodies in 3D. The code is based from this. My problem is the trajectories outputted seem to be small oscillations propagating as straight lines, whereas I should be getting something chaotic like random squiggles. I'm not sure if the problem lies in my implementation of the velocity Verlet algorithm or in the mathematics defining the gravitational force between 3 bodies or whatever else.
import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg as lag
from mpl_toolkits.mplot3d import Axes3D
plt.style.use('dark_background')
# masses of stars
m1 = 10e30
m2 = 20e30
m3 = 30e30
# starting coordinates for suns
r1_initial = np.array([-10e10, 10e10, -11e10])
v1_initial = np.array([3e10, 0, 0])
r2_initial = np.array([0, 0, 0])
v2_initial = np.array([0, 3e10, 0])
r3_initial = np.array([10e10, 10e10, 12e10])
v3_initial = np.array([-3e10, 0, 0])
def force_grav_3bodies(r1,r2,r3,m1,m2,m3,G):
force_A = -G*m1*m2*(r1-r2)/lag.norm(r1-r2)**3 - G*m1*m3*(r1-r3)/lag.norm(r1-r3)**3
force_B = -G*m2*m1*(r2-r3)/lag.norm(r2-r3)**3 - G*m2*m3*(r2-r1)/lag.norm(r2-r1)**3
force_C = -G*m3*m1*(r3-r1)/lag.norm(r3-r1)**3 - G*m3*m2*(r3-r2)/lag.norm(r3-r2)**3
return force_A, force_B, force_C
def accel_3sun(r1,r2,r3,m1,m2,m3,G):
fg3bA, fg3bB, fg3bC = force_grav_3bodies(r1,r2,r3,m1,m2,m3,G)
accel_sunA, accel_sunB, accel_sunC = fg3bA/m1, fg3bB/m2, fg3bC/m3
return accel_sunA, accel_sunB, accel_sunC
dt = 0.1
Tmax = 2000
steps = int(Tmax/dt)
G = 6.67e-11
r1 = np.zeros([steps,3])
v1 = np.zeros([steps,3])
a1 = np.zeros([steps,3])
r2 = np.zeros([steps,3])
v2 = np.zeros([steps,3])
a2 = np.zeros([steps,3])
r3 = np.zeros([steps,3])
v3 = np.zeros([steps,3])
a3 = np.zeros([steps,3])
# initial conditions
r1[0], r2[0], r3[0] = r1_initial, r2_initial, r3_initial
v1[0], v2[0], v3[0] = v1_initial, v2_initial, v3_initial
a1[0], a2[0], a3[0] = accel_3sun(r1[0],r2[0],r3[0],m1,m2,m3,G)
# velocity verlet algorithm
for k in range(1,steps-1):
r1[k] = r1[k-1] + v1[k-1]*dt + 0.5*a1[k-1]*dt*dt
r2[k] = r2[k-1] + v2[k-1]*dt + 0.5*a2[k-1]*dt*dt
r3[k] = r3[k-1] + v3[k-1]*dt + 0.5*a3[k-1]*dt*dt
a1[k], a2[k], a3[k] = accel_3sun(r1[k],r2[k],r3[k],m1,m2,m3,G)
v1[k] = v1[k-1] + 0.5*(a1[k-1] + a1[k])*dt
v2[k] = v2[k-1] + 0.5*(a2[k-1] + a2[k])*dt
v3[k] = v3[k-1] + 0.5*(a3[k-1] + a3[k])*dt
fig = plt.figure(figsize=(8, 8))
ax = fig.gca(projection='3d')
plt.gca().patch.set_facecolor('black')
plt.plot([i[0] for i in r1], [j[1] for j in r1], [k[2] for k in r1] , '^', color='r', lw = 0.5, markersize = 0.1, alpha=0.5)
plt.plot([i[0] for i in r2], [j[1] for j in r2], [k[2] for k in r2] , '^', color='y', lw = 0.5, markersize = 0.1, alpha=0.5)
plt.plot([i[0] for i in r3], [j[1] for j in r3], [k[2] for k in r3] , '^', color='c', lw = 0.5, markersize = 0.1, alpha=0.5)
ax.scatter(r1[steps-2][0],r1[steps-2][1],r1[steps-2][2],color="r",marker="^",s=100)
ax.scatter(r2[steps-2][0],r2[steps-2][1],r2[steps-2][2],color="y",marker="^",s=100)
ax.scatter(r3[steps-2][0],r3[steps-2][1],r3[steps-2][2],color="c",marker="^",s=100)
ax.scatter(r1[0][0],r1[0][1],r1[0][2],color="r",marker="o",s=100,label="A")
ax.scatter(r2[0][0],r2[0][1],r2[0][2],color="y",marker="o",s=100,label="B")
ax.scatter(r3[0][0],r3[0][1],r3[0][2],color="c",marker="o",s=100,label="C")
plt.legend()
plt.axis('on')
ax.w_xaxis.set_pane_color((0.0, 0.0, 0.0, 1.0)), ax.w_yaxis.set_pane_color((0.0, 0.0, 0.0, 1.0)), ax.w_zaxis.set_pane_color((0.0, 0.0, 0.0, 1.0))
plt.show()
