Problems with gravitational N-body simulation

Question

Background I am writing a simple N-body simulation code in Python, with core physics solvers such as acceleration and position integration implemented in Cython. The code uses the leapfrog method to integrate the positions.

Conditions

1. Masses: 10, 100 for 2 bodies; 1 for all bodies for number of bodies > 2.
2. Position: Randomly generated
3. Velocity: Randomly generated
4. Timestep: 1e-3

The bug

1. For number of bodies > 2: The bodies repel each other instead of attracting each other.
2. For number of bodies = 2: They do not orbit each other but instead travel in straight lines.
3. General bug (regardless of number of bodies): Repulsive forces

Expected behavior

1. Two bodies must orbit each other
2. The forces must be attractive

Efforts to resolve the bug

1. Negative sign for acceleration
2. Multiply acceleration by -1
3. Add a new temporary expression

Code Acceleration function (Cython):

def acceleration(np.ndarray[np.float64_t, ndim=2] pos, np.ndarray[np.float64_t, ndim=1] mass):
cdef int N = pos.shape[0] # Number of bodies   
# Memoryview for acceleration
cdef np.float64_t [:,:] acc = np.zeros((N,3),dtype="float64")
cdef double soft = 1e-4 # Softening length
cdef G = 6.673e-11 # Gravitational constant
# Pairwise separations
cdef double dx
cdef double dy
cdef double dz
# Total separation vectors
cdef double r
cdef double tmp
# Mass
cdef mj
# Acceleration calculation loop
for i in range(N):
    for j in range(N):
        # Remove gravitational self forces
        if i==j:
            continue
        
        # Calculate pairwise separation vectors
        dx = pos[j,0] - pos[i,0]
        dy = pos[j,1] - pos[i,1]
        dz = pos[j,2] - pos[i,2]

        # Vector magnitude of separation vector
        r = dx**2 + dy**2 + dz**2
        r = np.sqrt(r)

        # Mass
        mj = mass[j]

        tmp = G * mj * r**3

        # Calculate accelerations
        acc[i,0] += tmp * dx
        acc[i,1] += tmp * dy
        acc[i,2] += tmp * dz

return np.asarray(acc)

Position integration:

def leapfrog(np.ndarray[np.float64_t, ndim=2] pos, np.ndarray[np.float64_t, ndim=2] vel, np.ndarray[np.float64_t, ndim=2] acc, np.ndarray[np.float64_t, ndim=1] mass):
cdef double dt = 1e-3 # Timestep

# The Leapfrog integration method
# v(t+1/2) = v(t) + a(t) x dt/2
# x(t+1) = x(t) + v(t+1/2) x dt
# a(t+1) = (G * m/((dx^2 + dy^2 + dz^2)^(3/2))) * dx * x
# v(t+1) = v(t+1/2) + a(t+1) x dt/2
vel += acc * dt/2
pos += vel * dt
acc = acceleration(pos, mass)
vel += acc * dt/2

return pos, acc

Main simulation loop (Python): for _ in range(Nt): # Calculate positions and get new acceleration values pos, acc = leapfrog(pos, vel, acc, m)

Plot (Python):

plt.scatter(pos_arr[:,0], pos_arr[:,1])

Please help me solve this issue.

Refer to the images for more info regarding the bug:

Two bodies repelling each other + moving in a straight line instead of orbiting

3D View of the bug.

Edit 1: Velocity = 0 Here is the result with velocity = 0 (vel = np.zeros((N,3)), where N=2. The red point is the first element of the position array and the green point is the last point.

Edit 2: tmp * dx/2 : This is the result obtained by dividing the separation vectors by r. Updated code:

acc[i,0] += tmp * dx/r
acc[i,1] += tmp * dy/r
acc[i,2] += tmp * dz/r

Edit 3: velocity = 1e-6 This is what the result looks like for setting both velocities = 1e-6. They remain stationary.

Answer 1

The problem is solved: using a better acceleration algorithm solved the problem. The code is now available at GitHub. The only problem now is to animate the plot. Please let me know how to animate it in the comments. Thanks for the help all along!

Here are the results of a two-body simulation: