N-body problem

Shared by Antoine Rideau

On this page you will find how to simulate using Castor the n-body problem with 3 celestial bodies : the Sun, Jupiter and Saturn.

In physics, the n-body problem is the problem of predicting the individual motions of a group of celestial objects interacting with each other gravitationally.

First of all, Newton’s law of gravity says that the gravitational force felt on a planet P by the Sun S is given by
\[\overrightarrow{F}_{S\rightarrow P} = - \overrightarrow{F}_{P\rightarrow S} = -\frac{Gm_{S}m_{P}}{d^2}\overrightarrow{u} ,\]
where :
\(G\) : gravitational constant,
\(m_{S}\), \(m_{P}\) : masses of the Sun and the planet,
\(d\) : distance between the Sun and the planet,
\(\overrightarrow{u}\) : monic vector directed from the Sun to the planet.

The system is composed of the three heavier bodies in the solar system : the Sun, Jupiter and Saturn which will be referred as S, Ju and Sa. This system is supposed isolated and only gravitational forces are applied on the three celestial bodies.
According to Newton’s second principle, with \(m\) the masses, \(\overrightarrow{a}\) the acceleration and \(\overrightarrow{F}\) the forces :
\[\begin{split}\begin{matrix} m_{Ju}\overrightarrow{a}_{Ju} = \overrightarrow{F}_{S\rightarrow Ju} + \overrightarrow{F}_{Sa\rightarrow Ju} , \\ m_{Sa}\overrightarrow{a}_{Sa} = \overrightarrow{F}_{S\rightarrow Sa} + \overrightarrow{F}_{Ju\rightarrow Sa} , \\ m_{S}\overrightarrow{a}_{S} = \overrightarrow{F}_{Ju\rightarrow S} + \overrightarrow{F}_{Sa\rightarrow S} . \end{matrix}\end{split}\]
Bodies | Masses
(relatively to the Sun)

Sun and inner planet

\(m_{0}\) = 1.00000597682

Jupiter

\(m_{1}\) = 0.000954786104043

Saturn

\(m_{2}\) = 0.000285583733151

Gravitational constant \(G = 2.95912208286 \times 10^{-4}\)

// Parameters
matrix<> m = {{1.00000597682}, {9.54786104043e-4}, {2.85583733151e-4}}; // Masses : Sun, Jupiter and Saturn
double G = 2.95912208286e-4;                                            //Gravitation's constant

The position of a object over time is given by \(q(t)\).

Bodies

Initial position (AU)

Sun

0

0

0

Jupiter

−3.5023653

−3.8169847

−1.5507963

Saturn

9.0755314

−3.0458353

−1.6483708

where AU stands for astronomical unit and \(1 AU = 1.495 978 707 \times 10^{11}\) m.

matrix<> qini = {0, 0, 0,                               // Sun's initial position
                -3.5023653, -3.8169847, -1.5507963,     // Jupiter's initial position
                9.0755314, -3.0458353, -1.6483708};     // Saturn's initial position
Moreover, using the momentum \(\overrightarrow{p} = m\overrightarrow{v} = m\overrightarrow{\dot{q}}(t)\) instead of the speed is more practical.
Indeed, the system’s total momentum is preserved over time :
\[\overrightarrow{p}_{S}(t) + \overrightarrow{p}_{Ju}(t) + \overrightarrow{p}_{Sa}(t) = Constant\]

Bodies

Initial velocity (AU/day)

Sun

0

0

0

Jupiter

0.00565429

0.00565429

−0.00190589

Saturn

0.00168318

0.00483525

0.00192462

matrix<> vini = {0, 0, 0,                                                                   // Sun's initial velocity
0.00565429, -0.00412490, -0.00190589,                                                       // Jupiter's initial velocity
0.00168318, 0.00483525, 0.00192462};                                                        // Saturn's initial velocity
matrix<> pini = reshape(mtimes(transpose(m), ones(1, numel(m))), 1, numel(vini)) * vini;    // Initial momentums

See reshape , mtimes, transpose, ones, numel.

Time is discretized into nt steps

\[t_{i} = it \times \delta t \text{ for } it = \left [ \! \left [ 0, nt-1 \right ] \! \right ]\]
// Disretization
int nt = 1501;
double dt = (tend - tini) / (nt - 1);
auto T = linspace(tini, tend, nt);

See linspace.

Scheme

A symplectic Euler scheme is used in this simulation because it preserved the energy of the system unlike either forward and backward Euler scheme.

As the system is conservative the Hamiltonian can be separated in a cinetical part \(K(p)\) and a potential part \(V(q)\) :
\[H(q,p) = K(p) + V(q) ,\]

where

\[\begin{matrix} \displaystyle K(p) = \frac{1}{2}\frac{p^2}{m} & \text{ and } & \displaystyle V(q_{i}) = \sum_{j\neq i}- \frac{Gm_{j}m_{i}}{\left | q_{i}-q_{j} \right |} . \end{matrix}\]

With such a separation, Hamilton equation are given by

\[\begin{matrix} \displaystyle \frac{\mathrm{d} q}{\mathrm{d} t} = + \frac{\mathrm{d} K}{\mathrm{d} p} & \text{ and } & \displaystyle \frac{\mathrm{d} p}{\mathrm{d} t} = - \frac{\mathrm{d} V}{\mathrm{d} q} , \end{matrix}\]

where

\[\begin{matrix} \displaystyle \frac{\mathrm{d} K(p)}{\mathrm{d} p} = \frac{p}{m} & \text{ and } & \displaystyle \frac{\mathrm{d} V(q_{i})}{\mathrm{d} q} = \sum_{j\neq i} \frac{Gm_{j}m_{i}\left ( q_{i}-q_{j} \right )}{\left | q_{i}-q_{j} \right |^3} \end{matrix}\]

which result to the symplectic Euler scheme :

\[\begin{matrix} \displaystyle q_{n+1} = q_{n} + \frac{\mathrm{d} K}{\mathrm{d} p}(p_{n}) & \text{ and } & \displaystyle p_{n+1} = p_{n} - \frac{\mathrm{d} V}{\mathrm{d} q}(q_{n+1}) \end{matrix}\]
// Scheme
auto Q = zeros(nt, numel(qini));
Q(0, col(Q)) = qini;
auto P = zeros(nt, numel(pini));
P(0, col(P)) = pini;
// Symplectic Euler
for (int it = 0; it < nt - 1; it++)
{
    matrix<> q_n = eval(Q(it, col(Q)));
    matrix<> p_n = eval(P(it, col(P)));
    Q(it + 1, col(Q)) = q_n + dt * H_p(G, m, p_n);
    P(it + 1, col(P)) = p_n - dt * H_q(G, m, eval(Q(it + 1, col(Q))));
}

See zeros, numel, col , view.

In the code, \(\displaystyle \frac{\mathrm{d} K}{\mathrm{d} p}(p)\) is represented by the function H_p

matrix<> H_p(double G, matrix<> m, matrix<> p)
{
    auto Hp = zeros(1, numel(p));
    m = reshape(mtimes(transpose(m), ones(1, numel(m))), 1, numel(p));
    Hp = p / m;
    return Hp;
}

See zeros, reshape , mtimes, transpose, ones, numel.

and \(\displaystyle \frac{\mathrm{d} V}{\mathrm{d} q}(q)\) by the function H_q

matrix<> H_q(double G, matrix<> m, matrix<> q)
{

    auto q0 = eval(q(range(0, 3)));
    auto q1 = eval(q(range(3, 6)));
    auto q2 = eval(q(range(6, 9)));
    auto Hq = zeros(1, 9);
    Hq(range(0, 3)) = (G * m(0) * m(1) * ((q0 - q1) / pow(norm(q0 - q1), 3)) + G * m(0) * m(2) * ((q0 - q2) / pow(norm(q0 - q2), 3)));
    Hq(range(3, 6)) = (G * m(1) * m(0) * ((q1 - q0) / pow(norm(q1 - q0), 3)) + G * m(1) * m(2) * ((q1 - q2) / pow(norm(q1 - q2), 3)));
    Hq(range(6, 9)) = (G * m(2) * m(0) * ((q2 - q0) / pow(norm(q2 - q0), 3)) + G * m(2) * m(1) * ((q2 - q1) / pow(norm(q2 - q1), 3)));
    return Hq;
}

See range, view, zeros, norm.

Visualisation

Simple figure with Castor

The simplest method to visualize the results is to plot them using plot or plot3, here plot3 is used to show the motion in 3 dimensions.

For each coordinates x,y and z, the Sun’s positions are subtracted in order to keep it still in the center.
Moreover, transpose is needed because of matrix Q ‘s dimensions.
// Visu
figure fig;
plot3(fig, transpose(eval(Q(row(Q), 3)) - eval(Q(row(Q), 0))), transpose(eval(Q(row(Q), 4)) - eval(Q(row(Q), 1))), transpose(eval(Q(row(Q), 5)) - eval(Q(row(Q), 2))), {"c"});
plot3(fig, transpose(eval(Q(row(Q), 6)) - eval(Q(row(Q), 0))), transpose(eval(Q(row(Q), 7)) - eval(Q(row(Q), 1))), transpose(eval(Q(row(Q), 8)) - eval(Q(row(Q), 2))), {"b"});
plot3(fig, zeros(1, nt), zeros(1, nt), zeros(1, nt), {"y"});
_images/3body.png

Orbits of Jupiter (cyan) and Saturn (blue) around the Sun (yellow) in the center.

Video output with VTK

A way to visualize the results through a video is by using C++ VTK video writer.
However, C++ VTK video writer’s behavior is very OS dependent, it works perfectly fine on MAC but can cause some issues on Linux. (Windows ?)

First of all, the source and the writer need to be initialized with name of the output file, quality, framerate and connnected together.
// Initialize source and movie
vtkNew<vtkWindowToImageFilter> source;
vtkNew<vtkOggTheoraWriter> movie;
movie->SetInputConnection(source->GetOutputPort());
movie->SetFileName("nbody.avi");
movie->SetQuality(2); // in [0,2]
movie->SetRate(25);   // frame per seconds
int Nplot = 150;      // < 200

Afterwards, the writer is initiated before the loop over time.

movie->Start();
for (int it = 0; it < nt - 1; it++){...}

Then each wanted frame is plotted and added to the final movie.

// Visu
if (it % (nt / Nplot) == 0)
{
    figure fig;
    matrix<> limits = {-10, 10};
    plot(fig, {Q(it, 3) - Q(it, 0)}, {Q(it, 4) - Q(it, 1)}, {"c"});
    plot(fig, {Q(it, 6) - Q(it, 0)}, {Q(it, 7) - Q(it, 1)}, {"b"});
    plot(fig, zeros(1), zeros(1), {"y"});
    xlim(fig, limits);
    ylim(fig, limits);
    source->SetInput(fig.GetView()->GetRenderWindow());
    source->SetInputBufferTypeToRGB();
    source->ReadFrontBufferOff();
    movie->Write();
}

See plot, label-xlim, label-ylim.

Finally, the writer is closed after the loop over time.

for (int it = 0; it < nt - 1; it++){...}
movie->End();
Orbit of Jupiter (orange) and Saturn (green) around the Sun in the center.

Video animation with Python

Another way to get an animation of the orbiting planets is to output our data, and then post processing those with Python.
To do so, first the positions in the matrix Q are stored in a .txt file using writetxt .
So as to keep the Sun still in the center, its positions are subtracted of Jupiter’s and Saturn’s positions.
// Output
writetxt("./", "dataJu.txt", cat(2, eval(Q(row(Q), 3)) - eval(Q(row(Q), 0)), eval(Q(row(Q), 4)) - eval(Q(row(Q), 1))));
writetxt("./", "dataSa.txt", cat(2, eval(Q(row(Q), 6)) - eval(Q(row(Q), 0)), eval(Q(row(Q), 7)) - eval(Q(row(Q), 1))));

See writetxt, row .

Then the following Python code shows the beautiful animation.

import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np
from collections import deque

# Data input
dataJu = np.loadtxt("./build/dataJu.txt")
dataSa = np.loadtxt("./build/dataSa.txt")

# Parameters extraction
nt = int(dataJu[0, 0])

# Data processing
dataJu = np.delete(dataJu, 0, 0)
dataSa = np.delete(dataSa, 0, 0)

# Visu initialization
fig = plt.figure(figsize=(5, 4))
ax = fig.add_subplot(autoscale_on=False, xlim=(-10, 10), ylim=(-10, 10))
ax.set_aspect('equal')

line, = ax.plot([], [], 'o', lw=2)
traceJu, = ax.plot([], [], ',-', lw=1)
traceSa, = ax.plot([], [], ',-', lw=1)
historyJu_x, historyJu_y = deque(maxlen=nt), deque(maxlen=nt)
historySa_x, historySa_y = deque(maxlen=nt), deque(maxlen=nt)


def animate(i):
    # Get planets' current positions
    thisx = [0, dataJu[i, 0], dataSa[i, 0]]
    thisy = [0, dataJu[i, 1], dataSa[i, 1]]

    # Clear the trace when the animation loops
    if i == 0:
        historyJu_x.clear()
        historyJu_y.clear()
        historySa_x.clear()
        historySa_y.clear()

    # Add the current position to the trace
    historyJu_x.appendleft(thisx[1])
    historyJu_y.appendleft(thisy[1])
    historySa_x.appendleft(thisx[2])
    historySa_y.appendleft(thisy[2])

    line.set_data(thisx, thisy)  # Update planets' positions
    # Update planets' traces
    traceJu.set_data(historyJu_x, historyJu_y)
    traceSa.set_data(historySa_x, historySa_y)
    return line, traceJu, traceSa


# Creating the Animation object
ani = animation.FuncAnimation(
    fig, animate, nt, interval=10, blit=True)
plt.show()
Orbit of Jupiter (orange) and Saturn (green) around the Sun in the center.

Code

Here is all the code at once, without the functions H_q and H_p written above :

#include "castor/matrix.hpp"
#include "castor/graphics.hpp"
#include "castor/linalg.hpp"

using namespace castor;

    int main(int argc, char const *argv[])
{
    // Parameters
    matrix<> m = {1.00000597682, 9.54786104043e-4, 2.85583733151e-4}; // Masses : Sun, Jupiter and Saturn
    double G = 2.95912208286e-4;                                      //Gravitation's constant

    matrix<> qini = {0, 0, 0,                                                                // Sun's initial position
                     -3.5023653, -3.8169847, -1.5507963,                                     // Jupiter's initial position
                     9.0755314, -3.0458353, -1.6483708};                                     // Saturn's initial position
    matrix<> vini = {0, 0, 0,                                                                // Sun's initial velocity
                     0.00565429, -0.00412490, -0.00190589,                                   // Jupiter's initial velocity
                     0.00168318, 0.00483525, 0.00192462};                                    // Saturn's initial velocity
    matrix<> pini = reshape(mtimes(transpose(m), ones(1, numel(m))), 1, numel(vini)) * vini; // Initial momentums

    double tini = 0.;
    double tend = 12500.;

    // Disretization
    int nt = 1501;
    double dt = (tend - tini) / (nt - 1);
    auto T = linspace(tini, tend, nt);

    // Initialize source and movie
    vtkNew<vtkWindowToImageFilter> source;
    vtkNew<vtkOggTheoraWriter> movie;
    movie->SetInputConnection(source->GetOutputPort());
    movie->SetFileName("nbody.avi");
    movie->SetQuality(2); // in [0,2]
    movie->SetRate(25);   // frame per seconds
    int Nplot = 150;      // < 200

    // Scheme
    auto Q = zeros(nt, numel(qini));
    Q(0, col(Q)) = qini;
    auto P = zeros(nt, numel(pini));
    P(0, col(P)) = pini;
    // Symplectic Euler
    tic();
    movie->Start();
    for (int it = 0; it < nt - 1; it++)
    {

        matrix<> q_n = eval(Q(it, col(Q)));
        matrix<> p_n = eval(P(it, col(P)));
        Q(it + 1, col(Q)) = q_n + dt * H_p(G, m, p_n);
        P(it + 1, col(P)) = p_n - dt * H_q(G, m, eval(Q(it + 1, col(Q))));

        // Visu
        if (it % (nt / Nplot) == 0)
        {
            figure fig;
            matrix<> L({-10, 10, -10, 10}); // Axis dimensions
            plot(fig, Q(it, 3) - Q(it, 0) * ones(1), Q(it, 4) - Q(it, 1) * ones(1), L, {"c"});
            plot(fig, Q(it, 6) - Q(it, 0) * ones(1), Q(it, 7) - Q(it, 1) * ones(1), L, {"b"});
            plot(fig, zeros(1), zeros(1), {"y"});
            source->SetInput(fig.GetView()->GetRenderWindow());
            source->SetInputBufferTypeToRGB();
            source->ReadFrontBufferOff();
            movie->Write();
        }
    }
    movie->End();
    toc();

    // Output for Python post-processing and visualisation
    // writetxt("./", "dataJu.txt", cat(2, eval(Q(row(Q), 3)) - eval(Q(row(Q), 0)), eval(Q(row(Q), 4)) - eval(Q(row(Q), 1))));
    // writetxt("./", "dataSa.txt", cat(2, eval(Q(row(Q), 6)) - eval(Q(row(Q), 0)), eval(Q(row(Q), 7)) - eval(Q(row(Q), 1))));

    // Visu in native Castor
    // figure fig;
    // plot3(fig, transpose(eval(Q(row(Q), 3)) - eval(Q(row(Q), 0))), transpose(eval(Q(row(Q), 4)) - eval(Q(row(Q), 1))), transpose(eval(Q(row(Q), 5)) - eval(Q(row(Q), 2))), {"c"});
    // plot3(fig, transpose(eval(Q(row(Q), 6)) - eval(Q(row(Q), 0))), transpose(eval(Q(row(Q), 7)) - eval(Q(row(Q), 1))), transpose(eval(Q(row(Q), 8)) - eval(Q(row(Q), 2))), {"b"});
    // plot3(fig, zeros(1, nt), zeros(1, nt), zeros(1, nt), {"y"});

    // drawnow(fig);

    return 0;
}

References