Create Your Own N-body Simulation (With Python)

For today’s recreational coding exercise, we will look at the gravitational N-body problem. We will create a simulation of a dynamical system of particles interacting with each other gravitationally. Such a system may describe the orbits of planets in the Solar System or stars in our Galaxy.

You may find the accompanying Python code on github. (And if you prefer to use Matlab, please see my Matlab version of the article)

But first, below is a gif of what running our simulation looks like:

Force Calculation

We will assume a system of N point particles, indexed by i=1,…,N. Each particle has a:

• mass mᵢ,
• position rᵢ = [ xᵢ, yᵢ, zᵢ ] ,
• velocity vᵢ = [ vxᵢ, vyᵢ, vzᵢ ]

Each particle feels the gravitational attraction of all the other particles according to Newton’s law of universal gravitation (the famous ‘inverse square-law’). That is, each particle feels an acceleration:

where G=6.67×10^-11 m³/kg/s² is the Gravitational constant.

We can write a Python function to perform the calculation on an input N×3 matrix of positions: