Using Velocity Verlet to simulate 3-body motion in 3D

Question

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

Answer 1

Here is an example of Verlet algorithm implementation for Lennard Jones potential in 2D. There is no problem to go from 2 dimensions to 3.

import numpy as np

def lennard_jones_potential(r, sigma: float, eps: float) -> np.ndarray:
    """
    Lennard-Jones potential
    """
    return 4 * eps * np.power(sigma, 12) / np.power(r, 12) - 4 * eps * np.power(sigma, 6) / np.power(r, 6)


def lennard_jones_force(r, sigma: float, eps: float) -> np.ndarray:
    """
    Lennard-Jones force (Derivative of LJ potential)
    """
    return 24 * eps * np.power(sigma, 6) / np.power(r, 7) - 48 * eps * np.power(sigma, 12) / np.power(r, 13)


def get_kinetic_energy(m, v):
    """
    Get kinetic energy of particle
    """
    return m * (v ** 2) / 2


def get_verlet_next_x(x, v, f, m, dt):
    """
    Get next position of particle
    """
    a = f / m
    return x + v * dt + a * (dt ** 2) / 2


def get_verlet_next_f(r, rx, sigma: float, eps: float):
    """
    Get next LJ interaction force of particle
    """
    return np.sum(lennard_jones_force(r, sigma, eps) * rx / r)


def get_verlet_next_u(rs, sigma: float, eps: float):
    """
    Get next LJ interaction potential of particle
    """
    return np.sum(lennard_jones_potential(rs, sigma, eps)) / 2


def get_verlet_next_v(v, f, f_next, m, dt):
    """
    Get next velocity of particle using Verlet algorithm
    """
    a, a_next = f / m, f_next / m
    return v + ((a + a_next) * dt) / 2


def get_verlet_next(p: np.ndarray, rx: np.ndarray, ry: np.ndarray, dt: float, sigma: float, eps: float):
    """
    Get next state values of particle using Verlet algorithm
    (force, velocity, potential, kinetic energy)
    """
    r = np.sqrt(np.power(rx, 2) + np.power(ry, 2))

    fx_next = get_verlet_next_f(r, rx, sigma, eps)
    fy_next = get_verlet_next_f(r, ry, sigma, eps)

    vx_next = get_verlet_next_v(p[P_VX], p[P_FX], fx_next, p[P_M], dt)
    vy_next = get_verlet_next_v(p[P_VY], p[P_FY], fy_next, p[P_M], dt)
    v_next = np.sqrt(np.power(vx_next, 2) + np.power(vy_next, 2))

    u_next = get_verlet_next_u(r, sigma, eps)
    ek_next = get_kinetic_energy(p[P_M], v_next)

    return fx_next, fy_next, vx_next, vy_next, v_next, u_next, ek_next

# Indexes of particle properties
P_ID, P_M, P_X, P_Y, P_VX, P_VY, P_FX, P_FY = 0, 1, 2, 3, 4, 5, 6, 7
P_X_NEXT, P_Y_NEXT, P_VX_NEXT, P_VY_NEXT, P_FX_NEXT, P_FY_NEXT = 8, 9, 10, 11, 12, 13

# Indexes of result properties
R_ID, R_M, R_X, R_Y, R_VX, R_VY, R_FX, R_FY, R_U, R_EK, R_ET, R_T = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11

def verlet_algo(particles: np.ndarray, simulation_time: float, dt: float, sigma: float, eps: float) -> np.ndarray:
    """
    Verlet algorithm for N particles in 2D space
    """
    time = np.arange(0, simulation_time + dt, dt)

    result = np.zeros(shape=(particles.shape[0], time.shape[0], 12))

    for i in range(particles.shape[0]):
        result[i, :, R_ID] = np.repeat(particles[i, P_ID], result.shape[1])
        result[i, :, R_M] = np.repeat(particles[i, P_M], result.shape[1])

    result[:, 0, R_X] = particles[:, P_X]
    result[:, 0, R_Y] = particles[:, P_Y]
    result[:, 0, R_VX] = particles[:, P_VX]
    result[:, 0, R_VY] = particles[:, P_VY]
    result[:, 0, R_FX] = particles[:, P_FX]
    result[:, 0, R_FY] = particles[:, P_FY]
    result[:, 0, R_U] = 0
    result[:, 0, R_EK] = 0
    result[:, 0, R_ET] = 0
    result[:, :, R_T] = time

    # Calculate the initial forces and energies
    for i in range(particles.shape[0]):
        p = particles[i]
        p_id = p[P_ID].astype(np.int64)
        neighbours = particles[particles[:, P_ID] != p[P_ID]]

        rx = neighbours[:, P_X] - p[P_X]
        ry = neighbours[:, P_Y] - p[P_Y]

        fx, fy, vx, vy, v, u, ek = get_verlet_next(p, rx, ry, dt, sigma, eps)

        particles[i, P_FX] = fx
        particles[i, P_FY] = fy

        result[p_id, 0, R_FX] = fx
        result[p_id, 0, R_FY] = fy
        result[p_id, 0, R_U] = u
        result[p_id, 0, R_EK] = ek
        result[p_id, 0, R_ET] = ek + u

    for curr in range(len(time) - 1):
        next = curr + 1

        # Get following coordinates of all particles:
        particles[:, P_X_NEXT] = get_verlet_next_x(particles[:, P_X], particles[:, P_VX], particles[:, P_FX],
                                                   particles[:, P_M], dt)
        particles[:, P_Y_NEXT] = get_verlet_next_x(particles[:, P_Y], particles[:, P_VY], particles[:, P_FY],
                                                   particles[:, P_M], dt)

        for i in range(particles.shape[0]):
            p = particles[i]
            p_id = p[P_ID].astype(np.int64)
            neighbours = particles[particles[:, P_ID] != p[P_ID]]

            # Get distances between current particle and all other particles:
            rx = neighbours[:, P_X_NEXT] - p[P_X_NEXT]
            ry = neighbours[:, P_Y_NEXT] - p[P_Y_NEXT]

            # Get next force, velocity and energies for current particle using verlet algorithm:
            fx_next, fy_next, vx_next, vy_next, v_next, u_next, ek_next = get_verlet_next(p, rx, ry, dt, sigma, eps)

            # Save next values of current particle:
            p[P_VX_NEXT] = vx_next
            p[P_VY_NEXT] = vy_next
            p[P_FX_NEXT] = fx_next
            p[P_FY_NEXT] = fy_next

            # Save next values of current particle to result array:
            result[p_id, next, R_X] = p[P_X_NEXT]
            result[p_id, next, R_Y] = p[P_Y_NEXT]
            result[p_id, next, R_VX] = vx_next
            result[p_id, next, R_VY] = vy_next
            result[p_id, next, R_FX] = fx_next
            result[p_id, next, R_FY] = fy_next
            result[p_id, next, R_U] = u_next
            result[p_id, next, R_EK] = ek_next
            result[p_id, next, R_ET] = ek_next + u_next

        # Apply next values to all particles:
        particles[:, P_X] = particles[:, P_X_NEXT]
        particles[:, P_Y] = particles[:, P_Y_NEXT]
        particles[:, P_VX] = particles[:, P_VX_NEXT]
        particles[:, P_VY] = particles[:, P_VY_NEXT]
        particles[:, P_FX] = particles[:, P_FX_NEXT]
        particles[:, P_FY] = particles[:, P_FY_NEXT]

    return result
    
particles = np.array([
    #id,m, x,  y,  vx,vy,fx,fy,ẋ, ẏ, vẋ,vẏ,fẋ,fẏ (ẋ == x_next, ...)
    [0, 2, 11, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [1, 1, 10, 13, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [2, 3, 12, 13, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
], dtype=np.float64)

result = verlet_algo(particles, simulation_time=50, dt=0.002, sigma=1, eps=0.7, animate=True)

Energies plot:

More detailes you can see on my own GitHub repository.

Answer 2

The basic example orbit, or one variant of it, is a circle of radius 1 with the central mass in the center. If the time scale is chosen so that the orbital speed is 1, then the constant mu = GM in the equation q'' = -mu*q/|q|^3 is mu=1. Let's distribute that as M=1, so that also G=1.

Now change to the scales of the linked post, radius about 10, total mass assumed as central mass of M=60, G=9.8, time scale T. In the new positions and masses the equation reads as

T^2*(q/10)'' = -(G/9.8)*(M/60)*(q/10)/|q/10|^3
q'' = -G*M*q/|q|^3 * 50/(3*9.8)/T^2

Fixing T=1.3 almost removes the last factor, so the orbital speed of the new configuration is about R/T=10/1.3=7.7, which is in the scale of the velocities on the linked page. The trajectories there are not circles, the bodies "fall into each other".

You want to scale the system to the dimensions of the solar system and use the natural gravity constant. So the radius becomes R=100e9 m = 100e6 km or 2/3 AU. The total mass is M=60e30 kg or 30 sun masses, the time unit changes to some unspecified T. The same dance as above changes the gravity formula to

T^2*(q/1e11)'' = -(G/6.67e-11)*(M/6e31)*(q/1e11)/|q/1e11|^3

q'' = -G*M*q/|q|^3 * 1e13/(6*6.67)/T^2

To set the last factor to one requires about T=5e5, which gives an orbital speed of R/T=10e10/5e5=2e5 or about 200 km/s. So to get similar plots like on the web page you need velocities in this range, for instance by scaling the velocities there by a factor of 2e5/8=25000, so that a velocity component of 3 becomes 75000. Your initial velocities are much higher, so that the bodies will very rapidly move away from each other without influencing each other significantly.

The above velocities give reasonable but not very interesting results. Reducing the velocity components to size 3e4 gives more interactions in the picture