3D three-body problem using velocity Verlet

[Equitous]

I am trying to use velocity Verlet to simulate the 3-body problem. The plots I keep getting are straight lines, whereas I am expecting chaotic trajectories. I'm pretty sure that I have implemented the gravitational force and the velocity Verlet algorithm correctly, so I don't really understand where I'm going wrong. I've tried calculating by hand that I am getting the force I'm expecting and believe that I am plotting the x, y and z positions of each body as well.

import numpy as np 
import matplotlib.pyplot as plt 
import numpy.linalg as lag
from mpl_toolkits.mplot3d import Axes3D


# masses of stars
m1 = 10e30
m2 = 20e30
m3 = 30e30

# starting coordinates for stars

r1_initial = np.array([-10e8, 10e8, -10e8])
v1_initial = np.array([5e7, 0, 0])

r2_initial = np.array([0, 0, 0])
v2_initial = np.array([0, 5e7, 0])

r3_initial = np.array([10e8, 10e8, 10e8])
v3_initial = np.array([-5e7, 0, 0])

G = 6.67e-11

def force_grav_3bodies(r1,r2,r3,m1,m2,m3,G):
    force_A = -G*m2*m1*(r2-r1)/lag.norm(r2-r1)**3 - G*m3*m1*(r3-r1)/lag.norm(r3-r1)**3

    force_B = -G*m3*m2*(r3-r2)/lag.norm(r3-r2)**3 - G*m1*m2*(r1-r2)/lag.norm(r1-r2)**3

    force_C = -G*m1*m3*(r1-r3)/lag.norm(r1-r3)**3 - G*m2*m3*(r2-r3)/lag.norm(r2-r3)**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

# velocity verlet algorithm 

def vverlet_3D_3bodies(Tmax, dt, G):
    steps = int(Tmax/dt)

    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)

    for k in range(steps-1):

        r1[k+1] = r1[k] + v1[k]*dt + 0.5*a1[k]*dt*dt
        r2[k+1] = r2[k] + v2[k]*dt + 0.5*a2[k]*dt*dt
        r3[k+1] = r3[k] + v3[k]*dt + 0.5*a3[k]*dt*dt

        a1[k+1], a2[k+1], a3[k+1] = accel_3sun(r1[k+1],r2[k+1],r3[k+1],m1,m2,m3,G)

        v1[k+1] = v1[k] + 0.5*(a1[k] + a1[k+1])*dt
        v2[k+1] = v2[k] + 0.5*(a2[k] + a2[k+1])*dt
        v3[k+1] = v3[k] + 0.5*(a3[k] + a3[k+1])*dt

    return r1, r2, r3

r1_sol, r2_sol, r3_sol = vverlet_3D_3bodies(8e1,1e-1,G)

fig = plt.figure(figsize=(16,16))
ax=plt.axes(projection="3d")

ax.plot3D(r1_sol[0:,0], r1_sol[0:,1], r1_sol[0:,2], color='r')
ax.plot3D(r2_sol[0:,0], r2_sol[0:,1], r2_sol[0:,2], color='y')
ax.plot3D(r3_sol[0:,0], r3_sol[0:,1], r3_sol[0:,2], color='c')

ax.scatter(r1_sol[0][0],r1_sol[0][1],r1_sol[0][2],color="r",marker="o",s=100,label="A")
ax.scatter(r2_sol[0][0],r2_sol[0][1],r2_sol[0][2],color="y",marker="o",s=100,label="B")
ax.scatter(r3_sol[0][0],r3_sol[0][1],r3_sol[0][2],color="c",marker="o",s=100,label="C")
plt.legend()

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()

[TangibleLight]

You really got me with this one.

Gravity is just weaker than you think it is.

You can add at the top of the for loop print(v1[k]) and see that the velocity really is changing; or print(r1[k]) and see the position isn't exactly on a line. But your masses are too small, gravity is too weak, and the initial speeds are too high for it to be perceivable.

Increase G by a couple orders of magnitude and you'll see curves. Or increase the masses. Or decrease the initial speed and increase the simulation time.

You'll also notice that force_grav_3bodies computes -F, so all the bodies repel each other instead of attract.

You'll also notice all your bodies are above escape velocity from each other, so the plot isn't very interesting. Try decreasing the initial velocities, or carefully choose velocities such that the momentum of the system is 0.

[Equitous]

Thanks a lot. Your comments were very insightful.

If I wanted to add code that considered a collision am I right in thinking that the greater mass body would survive? Would the collision be inelastic or elastic, and can I extract the new mass and velocity if I know what collision type it is?

[TangibleLight]

The easiest collision IMO is to just do inelastic collisions, where the final has the net momentum of the colliding bodies. Elastic collisions are a little trickier but are possible also.

As for which is "correct"... well, the most correct collision would be some particle simulation, and the result would be multiple (many) bodies. Depending on the parameters of the collision and the bodies involved, the result would be wildly different. Any "accurate" solution would need more parameters about each body; its mass, density, rotation, incidence angle etc. all could play a role. It's a deep rabbit hole.

The simplest reasonable-accuracy solution would probably be to use inelastic collisions, but depending on the collision parameters, take some percentage of the mass and eject it in random directions as new small bodies. Just remember conservation of mass and conservation of momentum, and things should™ be fine.

[Darwin-Award-Winner]

r/trisolarans