N-body simulation

In tree methods, such as a Barnes–Hut simulation, an octree is usually used to divide the volume into cubic cells and only interactions between particles from nearby cells need to be treated individually; particles in distant cells can be treated collectively as a single large particle centered at the distant cell's center of mass (or as a low-order multipole expansion). This can dramatically reduce the number of particle pair interactions that must be computed. To prevent the simulation from becoming swamped by computing particle-particle interactions, the cells must be refined to smaller cells in denser parts of the simulation which contain many particles per cell. For simulations where particles are not evenly distributed, the well-separated pair decomposition methods of Callahan and Kosaraju yield optimal O(n log n) time per iteration with fixed dimension.

Another possibility is the particle mesh method in which space is discretised on a mesh and, for the purposes of computing the gravitational potential, particles are assumed to be divided between the surrounding 2x2 vertices of the mesh. Finding the potential energy Φ is easy, because the Poisson equation

$${\displaystyle \nabla ^{2}\Phi =4\pi G{\rho },\,}$$

where G is Newton's constant and ${\displaystyle \rho }$ is the density (number of particles at the mesh points), is trivial to solve by using the fast Fourier transform to go to the frequency domain where the Poisson equation has the simple form

$${\displaystyle {\hat {\Phi }}=-4\pi G{\frac {\hat {\rho }}{k^{2}}},}$$

where ${\displaystyle {\vec {k}}}$ is the comoving wavenumber and the hats denote Fourier transforms. Since ${\displaystyle {\vec {g}}=-{\vec {\nabla }}\Phi }$, the gravitational field can now be found by multiplying by ${\displaystyle -i{\vec {k}}}$ and computing the inverse Fourier transform (or computing the inverse transform and then using some other method). Since this method is limited by the mesh size, in practice a smaller mesh or some other technique (such as combining with a tree or simple particle-particle algorithm) is used to compute the small-scale forces. Sometimes an adaptive mesh is used, in which the mesh cells are much smaller in the denser regions of the simulation.

Several different gravitational perturbation algorithms are used to get fairly accurate estimates of the path of objects in the Solar System.

People often decide to put a satellite in a frozen orbit. The path of a satellite closely orbiting the Earth can be accurately modeled starting from the 2-body elliptical orbit around the center of the Earth, and adding small corrections due to the oblateness of the Earth, gravitational attraction of the Sun and Moon, atmospheric drag, etc. It is possible to find a frozen orbit without calculating the actual path of the satellite.

The path of a small planet, comet, or long-range spacecraft can often be accurately modeled starting from the 2-body elliptical orbit around the Sun, and adding small corrections from the gravitational attraction of the larger planets in their known orbits.

Some characteristics of the long-term paths of a system of particles can be calculated directly. The actual path of any particular particle does not need to be calculated as an intermediate step. Such characteristics include Lyapunov stability, Lyapunov time, various measurements from ergodic theory, etc.

Softening is a numerical trick used in N-body techniques to prevent numerical divergences when a particle comes too close to another (and the force goes to infinity). This is obtained by modifying the regularized gravitational potential of each particle as

$${\displaystyle \Phi =-{\frac {1}{\sqrt {r^{2}+\epsilon ^{2}}}},}$$

(rather than 1/r) where ${\displaystyle \epsilon }$ is the softening parameter. The value of the softening parameter should be set small enough to keep simulations realistic.

The simplest implementation of N-body simulations where ${\textstyle n\geq 3}$ is a naive propagation of orbiting bodies; naive implying that the only forces acting on the orbiting bodies is the gravitational force which they exert on each other. In object-oriented programming languages, such as C++, some boilerplate code is useful for establishing the fundamental mathematical structures as well as data containers required for propagation; namely state vectors, and thus vectors, and some fundamental object containing this data, as well as the mass of an orbiting body. This method is applicable to other types of N-body simulations as well; a simulation of point masses with charges would use a similar method, however the force would be due to attraction or repulsion by interaction of electric fields. Regardless, acceleration of particle is a result of summed force vectors, divided by the mass of the particle:

$${\displaystyle {\vec {a}}={\frac {1}{m}}\sum {\vec {F}}}$$

An example of a programmatically stable and scalable method for containing kinematic data for a particle is the use of fixed length arrays, which in optimised code allows for easy memory allocation and prediction of consumed resources; as seen in the following C++ code:

struct Vector3
{
    double e[3] = { 0 };

    Vector3() {}
    ~Vector3() {}

    inline Vector3(double e0, double e1, double e2)
    {
        this->e[0] = e0;
        this->e[1] = e1;
        this->e[2] = e2;
    }
};

struct OrbitalEntity
{
    double e[7] = { 0 };

    OrbitalEntity() {}
    ~OrbitalEntity() {}

    inline OrbitalEntity(double e0, double e1, double e2, double e3, double e4, double e5, double e6)
    {
        this->e[0] = e0;
        this->e[1] = e1;
        this->e[2] = e2;
        this->e[3] = e3;
        this->e[4] = e4;
        this->e[5] = e5;
        this->e[6] = e6;
    }
};

Note that OrbitalEntity contains enough room for a state vector, where:

• ${\textstyle e_{0}=x}$, the projection of the objects position vector in Cartesian space along ${\displaystyle \left[1\;0\;0\right]}$
• ${\textstyle e_{1}=y}$, the projection of the objects position vector in Cartesian space along ${\displaystyle \left[0\;1\;0\right]}$
• ${\textstyle e_{2}=z}$, the projection of the objects position vector in Cartesian space along ${\displaystyle \left[0\;0\;1\right]}$
• ${\textstyle e_{3}={\dot {x}}}$, the projection of the objects velocity vector in Cartesian space along ${\displaystyle \left[1\;0\;0\right]}$
• ${\textstyle e_{4}={\dot {y}}}$, the projection of the objects velocity vector in Cartesian space along ${\displaystyle \left[0\;1\;0\right]}$
• ${\textstyle e_{5}={\dot {z}}}$, the projection of the objects velocity vector in Cartesian space along ${\displaystyle \left[0\;0\;1\right]}$

Additionally, OrbitalEntity contains enough room for a mass value.

Commonly, N-body simulations will be systems based of some type of equations of motion; of these, most will be dependent on some initial configuration to "seed" the simulation. In systems such as those dependent on some gravitational or electric potential, the force on a simulation entity is independent on its velocity. Hence, to seed the forces of the simulation, merely initial positions are needed, but this will not allow propagation- initial velocities are required. Consider a planet orbiting a star- it has no motion, but is subject to gravitational attraction to its host star. As a time progresses, and time steps are added, it will gather velocity according to its acceleration. For a given instant in time, ${\displaystyle t_{n}}$, the resultant acceleration of a body due to its neighbouring masses is independent of its velocity, however, for the time step ${\displaystyle t_{n+1}}$, the resulting change in position is significantly different due the propagation's inherent dependency on velocity. In basic propagation mechanisms, such as the symplectic euler method to be used below, the position of an object at ${\displaystyle t_{n+1}}$ is only dependent on its velocity at ${\displaystyle t_{n}}$, as the shift in position is calculated via

$${\displaystyle {\vec {r}}_{t_{n+1}}={\vec {r}}_{t_{n}}+{\vec {v}}_{t_{n}}}$$

Without acceleration, ${\textstyle {\vec {v}}_{t_{n}}}$is static, however, from the perspective of an observer seeing only position, it will take two time steps to see a change in velocity.

A solar-system-like simulation can be accomplished by taking average distances of planet equivalent point masses from a central star. To keep code simple, a non-rigorous approach based on semi-major axes and mean velocities will is used. Memory space for these bodies must be reserved before the bodies are configured; to allow for scalability, a malloc command may be used:

OrbitalEntity* orbital_entities = malloc(sizeof(OrbitalEntity) * (9 + N_ASTEROIDS));

orbital_entities[0] = { 0.0,0.0,0.0,        0.0,0.0,0.0,      1.989e30 };   // a star similar to the sun
orbital_entities[1] = { 57.909e9,0.0,0.0,   0.0,47.36e3,0.0,  0.33011e24 }; // a planet similar to mercury
orbital_entities[2] = { 108.209e9,0.0,0.0,  0.0,35.02e3,0.0,  4.8675e24 };  // a planet similar to venus
orbital_entities[3] = { 149.596e9,0.0,0.0,  0.0,29.78e3,0.0,  5.9724e24 };  // a planet similar to earth
orbital_entities[4] = { 227.923e9,0.0,0.0,  0.0,24.07e3,0.0,  0.64171e24 }; // a planet similar to mars
orbital_entities[5] = { 778.570e9,0.0,0.0,  0.0,13e3,0.0,     1898.19e24 }; // a planet similar to jupiter
orbital_entities[6] = { 1433.529e9,0.0,0.0, 0.0,9.68e3,0.0,   568.34e24 };  // a planet similar to saturn
orbital_entities[7] = { 2872.463e9,0.0,0.0, 0.0,6.80e3,0.0,   86.813e24 };  // a planet similar to uranus
orbital_entities[8] = { 4495.060e9,0.0,0.0, 0.0,5.43e3,0.0,   102.413e24 }; // a planet similar to neptune

where N_ASTEROIDS is a variable which will remain at 0 temporarily, but allows for future inclusion of significant numbers of asteroids, at the users discretion. A critical step for the configuration of simulations is to establish the time ranges of the simulation, ${\displaystyle t_{0}}$ to ${\displaystyle t_{\text{end}}}$, as well as the incremental time step ${\displaystyle dt}$ which will progress the simulation forward:

double t_0 = 0;
double t = t_0;
double dt = 86400;
double t_end = 86400 * 365 * 10; // approximately a decade in seconds
double BIG_G = 6.67e-11; // gravitational constant

The positions and velocities established above are interpreted to be correct for ${\displaystyle t=t_{0}}$.

The extent of a simulation would logically be for the period where ${\displaystyle t_{0}\leq t<t_{\text{end}}}$.

An entire simulation can consist of hundreds, thousands, millions, billions, or sometimes trillions of time steps. At the elementary level, each time step (for simulations with particles moving due to forces exerted on them) involves

• calculating the forces on each body
• calculating the accelerations of each body (${\displaystyle {\vec {a}}}$)
• calculating the velocities of each body (${\displaystyle {\vec {v}}_{n}={\vec {v}}_{n-1}+{\vec {a}}_{n}}$)
• calculating the new position of each body (${\displaystyle {\vec {r}}_{n+1}={\vec {r}}_{n}+{\vec {v}}_{n}}$)

The above can be implemented quite simply with a while loop which continues while ${\displaystyle t}$ exists in the aforementioned range:

while (t < t_end)
{
    for (size_t m1_idx = 0; m1_idx < 9 + N_ASTEROIDS; m1_idx++)
    {       
        Vector3 a_g = { 0,0,0 };

        for (size_t m2_idx = 0; m2_idx < 9 + N_ASTEROIDS; m2_idx++)
        {
            if (m2_idx != m1_idx)
            {
                Vector3 r_vector;

                r_vector.e[0] = orbital_entities[m1_idx].e[0] - orbital_entities[m2_idx].e[0];
                r_vector.e[1] = orbital_entities[m1_idx].e[1] - orbital_entities[m2_idx].e[1];
                r_vector.e[2] = orbital_entities[m1_idx].e[2] - orbital_entities[m2_idx].e[2];

                double r_mag = sqrt(
                        r_vector.e[0] * r_vector.e[0]
                      + r_vector.e[1] * r_vector.e[1]
                      + r_vector.e[2] * r_vector.e[2]);

                double acceleration = -1.0 * BIG_G * (orbital_entities[m2_idx].e[6]) / pow(r_mag, 2.0);

                Vector3 r_unit_vector = { r_vector.e[0] / r_mag, r_vector.e[1] / r_mag, r_vector.e[2] / r_mag };

                a_g.e[0] += acceleration * r_unit_vector.e[0];
                a_g.e[1] += acceleration * r_unit_vector.e[1];
                a_g.e[2] += acceleration * r_unit_vector.e[2];
            }
        }

        orbital_entities[m1_idx].e[3] += a_g.e[0] * dt;
        orbital_entities[m1_idx].e[4] += a_g.e[1] * dt;
        orbital_entities[m1_idx].e[5] += a_g.e[2] * dt;
    }

    for (size_t entity_idx = 0; entity_idx < 9 + N_ASTEROIDS; entity_idx++)
    {
        orbital_entities[entity_idx].e[0] += orbital_entities[entity_idx].e[3] * dt;
        orbital_entities[entity_idx].e[1] += orbital_entities[entity_idx].e[4] * dt;
        orbital_entities[entity_idx].e[2] += orbital_entities[entity_idx].e[5] * dt;
    }
    
    t += dt;
}

Focusing on the inner four rocky planets in the simulation, the trajectories resulting from the above propagation is shown below:

• von Hoerner, Sebastian (1960). "Die numerische Integration des n-Körper-Problemes für Sternhaufen. I". Zeitschrift für Astrophysik (in German). 50: 184. Bibcode:1960ZA.....50..184V.
• von Hoerner, Sebastian (1963). "Die numerische Integration des n-Körper-Problemes für Sternhaufen. II". Zeitschrift für Astrophysik (in German). 57: 47. Bibcode:1963ZA.....57...47V.
• Aarseth, Sverre J. (2003). Gravitational N-body Simulations: Tools and Algorithms. Cambridge University Press. ISBN 978-0-521-12153-8.
• Bertschinger, Edmund (1998). "Simulations of structure formation in the universe". Annual Review of Astronomy and Astrophysics. 36 (1): 599–654. Bibcode:1998ARA&A..36..599B. doi:10.1146/annurev.astro.36.1.599.
• Binney, James; Tremaine, Scott (1987). Galactic Dynamics. Princeton University Press. ISBN 978-0-691-08445-9.
• Callahan, Paul B.; Kosaraju, Sambasiva Rao (1992). "A decomposition of multidimensional point sets with applications to k-nearest-neighbors and n-body potential fields (preliminary version)". STOC '92: Proc. ACM Symp. Theory of Computing. ACM..