The N-body Problem

A numerical simulation in MATLAB

Image by author.

Newton’s discovery of universal attraction in the 17th Century dramatically changed our understanding of the motion of celestial bodies. His law masterly reconciles two profound physical principles: the principle of inertia, put forward by Galileo and Descartes in terrestrial mechanics, and the laws of Kepler, written down in Astronomia Nova.

Newton, in computing the orbital quantities of planets, realized that his law of gravitation wasn’t entirely accurate. The unforeseen consequence of Newton’s discovery was to question the belief that the solar system is stable: it was no longer obvious that planets kept moving immutably, without collisions or ejections. This started a two-century long competition started between astronomers, who made more and more precise observations, and geometers, who had the status and destiny of Newton’s law in their hand. At the heart of this, was the N-body problem.

In what follows, we will create a simulation of N particles interacting through a common gravitational potential. Poincare and Bruns proved that the differential equations of this system cannot be generally be integrated in closed form (that is, in terms of elementary functions instead of infinite series). So, if we want any meaningful attempt at a simulation, we must rely on numerical methods.

The system is set up as follows. We will assume N particles indexed i = 1, 2, …, N, with masses m_i, positions rᵢ = [ xᵢ, yᵢ, zᵢ ], and velocities vᵢ = [vxᵢ, vyᵢ, vzᵢ ]. Following Newton’s law of gravitation, each particle feels a force

where G=6.67×10^-11 m³/kg/s² is the Gravitational constant. To obtain the acceleration of our masses, we will compute getAcceleration as follows.

function [a] = getAcceleration(pos, mass, G, softening)
%   pos  is an N x 3 matrix of positions
%   mass is an N x 1 vector of masses
%   softening is the softening length
%   a is N x 3 matrix of accelerations

x = pos(:,1);
y = pos(:,2);
z = pos(:,3);
dx = x' - x;
dy = y' - y;
dz = z' - z;

% matrix that stores 1/r^3 for all particle pairwise particle separations
inverse = (dx.^2 + dy.^2 + dz.^2 + softening.^2).^(-3/2);

ax = G * (dx .* inverse) * mass;
ay = G * (dy .* inverse) * mass;
az = G * (dz .* inverse) * mass;

% pack together the acceleration components
a = [ax ay az];
end

The softening term is added to prevent two particles being abysmally close to each other, in which case the acceleration from the gravitational law goes to infinity.

The above code is vectorized. Although storing intermediate calculations in matrices requires a lot of memory, it’s incredibly useful with interpreted languages; this can sometimes lead to a speedup in the order of tens.

To validate our code, we rely on energy conservation — we know that this quantity should be constant through the time evolution.

The first term is the kinetic energy, defined as the momentum squared over twice the mass. The second term is the gravitational potential energy. Our code computes these quantities and keeps track of the total energy, ensuring validation of our approximations.

function [Ek, Ep] = getEnergy(pos, vel, mass, G)

% Kinetic Energy:
KE = 0.5 * sum(sum(mass.* vel.^2));

% Potential Energy:
% positions r = [x,y,z] for all particles
x = pos(:,1);
y = pos(:,2);
z = pos(:,3);

% matrix that stores all pairwise particle separations: r_j - r_i
dx = x' - x;
dy = y' - y;
dz = z' - z;

% matrix that stores r for all particle pairwise particle separations 
r = sqrt(dx.^2 + dy.^2 + dz.^2);

% sum over upper triangle, to count each interaction only once
PE = G *  sum(sum(triu(-(mass*mass')./r,1)));
end

To specify the initial conditions to our system, we’ll randomly choose values from a Gaussian distribution. You can set this up using the randn function in MATLAB.

For numerical integration, we use a leapfrog scheme which employs a kick-drift-kick form.

The evolution is performed in the code using a for-loop.

for i = 1:Nt
    vel = vel + acc * dt/2;
    pos = pos + vel * dt
    acc = getAcceleration(pos, mass, G, softening);
    vel = vel + acc * dt/2;
    t = t + dt;
end

The separation of the acceleration calculation onto the beginning and end of a step means that if time resolution is increased by a factor of two (Δt →Δt/2), then only one extra (computationally expensive) acceleration calculation is required.

There are two primary strengths to leapfrog integration when applied to mechanics problems. The first is the time-reversibility of the Leapfrog method. One can integrate forward n steps, and then reverse the direction of integration and integrate backwards n steps to arrive at the same starting position. The second strength is its symplectic nature, which sometimes allows for the conservation of a (slightly modified) energy of a dynamical system (only true for certain simple systems). This is especially useful when computing orbital dynamics, as many other integration schemes, such as the Runge–Kutta method, do not conserve energy and allow the system to drift substantially over time.

Finally, we use another for loop, for i = 1:ceil(end_t/dt), we run the leapfrog integrator and save the energies and positions for a plotting trail. Below is a result of the our simulation.

The N-body simulation for N = 10, dt = 0.01, softening = 0.1.

Find the GitHub repository here. That’s all from me. Thanks for reading and have a great day!