N-body simulation python

Question

I am trying to code an N-body simulation code in python and have successfully managed to produce a system involving the Sun, Earth and Jupiter as below using a leapfrog approximation method.

However, when I try and extend the same code for N bodies of equal mass all with zero velocity, I don't get the expected result of a system forming. Instead, the following is produced where the bodies spread out after initially being attracted to each other.

This same pattern is replicated regardless of the number of initial particles used.

This second image is just an enlarged version of the first showing they are initially attracted to each other.

Leading me to believe the error must lie in my initial conditions:

N = 3
mass = 1e30
R = 1e10
V = np.zeros([N,3])
M = np.full([N],mass)
P = np.random.uniform(-R, R, (N,3))
epsilon = 0.1 * R

acceleration calculation:

def calc_acceleration(position, mass, softening):
    
    G = 6.67 * 10**(-11)
    
    N = position.shape[0] # N = number of rows in particle_positions array
    acceleration = np.zeros([N,3])
    
    #print(N)
    for i in range(N):
        #print(i)
        for j in range(N):
            if i != j:
                #print("j", j)
                dx = position[i,0] - position[j,0]
                dy = position[i,1] - position[j,1]
                dz = position[i,2] - position[j,2]
                
                #print(dx,dy,dz)
                
                inv_r3 = ((dx**2 + dy**2 + dz**2 + softening**2)**(-1.5))
                
                acceleration[i,0] += - G * mass[j] * dx * inv_r3
                acceleration[i,1] += - G * mass[j] * dy * inv_r3
                acceleration[i,2] += - G * mass[j] * dz * inv_r3

    return(acceleration)

leap frog functions:

def calc_next_v_half(position, mass, velocity, softening, dt):
    half_velocity = np.zeros_like(velocity)
    half_velocity = velocity + calc_acceleration(position, mass, softening) * dt/2
    return(half_velocity)
       

def calc_next_position(position, mass, velocity, dt):
    next_position = np.zeros_like(position)
    
    next_position = position + velocity * dt
    
    return(next_position)

actual program function:

def programe(position, mass, velocity, softening, time, dt):
    
    no_of_time_steps = (round(time/dt))

    all_positions = np.full((no_of_time_steps, len(mass), 3), 0.0)
    all_velocities = []
    
    kinetic_energy = []
    potential_energy = []
    total_energy = []
        
        
    for i in range(no_of_time_steps):
        all_positions[i] = position
        all_velocities.append(velocity)

        'leap frog'
        velocity = calc_next_v_half(position, mass, velocity, softening, dt)    
        position = calc_next_position(position, mass, velocity, dt)    
        velocity = calc_next_v_half(position, mass, velocity, softening, dt)
        

    return(all_positions, all_velocities, kinetic_energy, potential_energy, total_energy)

This seems like a numerical artifact, they accelerate towards eachother then the numbers get really huge from the r^2 dependence, leading it to be rocketted out of the system. I think if you set some threshold on the proximity that the bodies can have it will help and potentially stop it from blowing up. — Philip Ciunkiewicz, Jan 16 '21 at 19:19

Lutz Lehmann · Answer 1 · 2021-10-21T12:41:10.990

The problem is that the symplectic methods only have their special properties as long as the systems stays well away from any singularities. For a gravitational system this is the case if it is hierarchical like in a solar system with sun, planets and moons where all orbits have low eccentricities.

However, if you consider a "star cluster" with objects of about equal mass, you do not get Kepler ellipses and the likelihood for very close encounters becomes rather high. The more so as your initial condition of zero velocity results in an initial free fall of all stars towards the common center of gravity, as can also be seen in your detail picture.

Due to the potential energy falling down into a singularity, the kinetic energy increases as the distance decreases, so close encounters equal high speed. With a constant step size like in the leapfrog-Verlet method, the sampling rate becomes too small to represent the curve, capture the swing-by fully. Energy conservation is grossly violated and the high speed is kept beyond the close encounter, leading to the unphysical explosion of the system.

N-body simulation python

1 Answers1