Multiparticle gravity system

This system was inspired by these articles

Introduction

Given two bodies $p_1$ and $p_2$ of masses $m_1$ and $m_2$; separated by a distance $r$, Newton’s law of universal gravitation establishes that the gravitational force that attracts one body to the other is:

$$F = G \frac{m_1 m_2}{r^2} \tag{1}$$

where $F$ is the gravitational force of two bodies; $r$ is the distance between them and $G$ is the gravitational constant: $G = 6.67430 \times 10^{-11} \frac{\mathrm{m}^3}{\mathrm{kg}\cdot\mathrm{s}^2}$. The gravitational force is the same in both bodies and point from one body to the other.

Given the positions of $p_1$, $\mathbf{s}_{1}$, and $p_2$, $\mathbf{s}_{2}$, the vector pointing from $p_1$ to $p_2$ is

$$\mathbf{r}_{12} = \mathbf{s}_{2} - \mathbf{s}_{1} \tag{2}$$

And the unitary vector in the direction from $p_1$ to $p_2$ is

$$\mathbf{u}_{12} = \frac {\mathbf{r}_{12}}{\sqrt {\mathbf{r}_{12} \cdot \mathbf{r}_{12}}} \tag{3}$$

The vector form of eq. 1 is:

$$\mathbf{F}_{12} = G \frac{m_1 m_2}{\mathbf{r}_{12} \cdot \mathbf{r}_{12}} \mathbf{u}_{12} \tag{4}$$

As can be seen in the figure above: $\mathbf{F}_{12} = -\mathbf{F}_{21}$.

Other equations

Newton’s second law of motion:

$$\mathbf{F} = m \mathbf{a} \Rightarrow \mathbf{a} = \frac{\mathbf{F}}{m} \tag{5}$$

Constant linear acceleration equations

$$\begin{align} \mathbf{s} & = \mathbf{s}_{0} + \mathbf{v} t + \frac12 \mathbf{a}t^2 \tag{6} \\ \mathbf{v} & = \mathbf{v}_{0} + \mathbf{a} t \tag{7} \\ \end{align}$$

Considering multi-body interaction illustrated in above figure, eq. 3 and eq. 4 become:

$$\begin{align} \mathbf{u}_{ij} & = \begin{cases} \frac {\mathbf{r}_{ij}} {r_{ij}}, & i \ne j\\ \mathbf{0}, & i = j \end{cases} \tag{8} \\ \\ \mathbf{F}_{ij} & = \begin{cases} G \frac {m_i m_j} {r_{ij}^2} \mathbf{u}_{ij}, & i \ne j \\ \mathbf{0}, & i = j \end{cases} \tag{9} \\ \\ \mathbf{F}_{ij} & = -\mathbf{F}_{ji} \tag{10} \end{align}$$

Approximations

It should be noted that, from eq. 2, $\mathbf{r} = \mathbf{r}(t)$ and consequently $\mathbf{F} = \mathbf{F}(t)$ (eq. 4). The simulation will progress in small time increments $\Delta t$ and consider that $\mathbf{F}$ is constant during these increments. The consequence of this procedure is shown in the figure below.

The dashed circle is where the body would really be after $\Delta t = t_1 - t_0$ time elapsed, but since I consider that force $\mathbf{F}$ is constant during $\Delta t$ interval, the calculated position is the one shown by the filled green disk. The deviation $\varepsilon$ has the following relations:

$$\begin{align} \varepsilon & \propto \Delta t \tag{11} \\ \varepsilon & \propto | \mathbf{v} | \tag{12} \\ \varepsilon & \propto \textrm{orbit curvature} \tag{13} \\ \end{align}$$

Of all factors affecting $\varepsilon$, $\Delta t$ is the only one that can be manipulated during the simulation.

Reducing the number of operations

In order to reduce the number of operations performed, I use some variations of equations of previous sections:

$$\begin{align} \mathbf{\Delta v}_i & = \mathbf{a}_i \Delta t \tag{14} \\ \mathbf{s}_i & = \mathbf{s}_{i_0} + (\mathbf{v}_i + \frac12 \mathbf{\Delta v}_i) \Delta t \tag{15} \\ \mathbf{v}_i & = \mathbf{v}_{i_0} + \mathbf{\Delta v}_i \tag{16} \\ \mathbf{a}_i & = \frac {G\sum_{j=1}^{n} \frac {m_i m_j} {r_{ij}^2} \mathbf{u}_{ij}}{m_i} \\ & = G\sum_{j=1}^{n} \frac {\mathbf{u}_{ij}} {r_{ij}^2} m_j \tag{17} \\ \end{align}$$

The velocity change $\mathbf{\Delta v}_i$ (eq. 14) has to be calculated to determine next velocity $\mathbf{v}_i$ (eq. 16) and can be reused in next position calculation $\mathbf{s}_i$ (eq. 15). For acceleration calculation $\mathbf{a}_i$ (eq. 17) I use the following matrix product:

$$\underbrace{ \begin{bmatrix} \mathbf{a}_1 \\ \mathbf{a}_2 \\ \vdots \\ \mathbf{a}_n \\ \end{bmatrix} }_{\mathbf{A}} = \underbrace{ \begin{bmatrix} \mathbf{0} & \frac{\mathbf{u}_{12}}{r_{12}^2} & \cdots & \frac{\mathbf{u}_{1n}}{r_{1n}^2} \\ \frac{\mathbf{-u}_{12}}{r_{12}^2} & \mathbf{0} & \cdots & \frac{\mathbf{u}_{2n}}{r_{2n}^2} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\mathbf{-u}_{1n}}{r_{1n}^2} & \frac{\mathbf{-u}_{2n}}{r_{2n}^2} & \cdots & \mathbf{0} \\ \end{bmatrix} }_{\mathbf{Ur}} \underbrace{ \left( G \begin{bmatrix} m_1 \\ m_2 \\ \vdots \\ m_n \\ \end{bmatrix} \right) }_{\mathbf{Gm}} \tag{18}$$

Using Hadamard product, eq. 18 becomes $\underbrace{ \begin{bmatrix} \mathbf{a}_1 \\ \mathbf{a}_2 \\ \vdots \\ \mathbf{a}_n \\ \end{bmatrix} }_{\mathbf{A}} = \left( \underbrace{ \begin{bmatrix} \mathbf{0} & \mathbf{u}_{12} & \cdots & \mathbf{u}_{1n} \\ \mathbf{-u}_{12} & \mathbf{0} & \cdots & \mathbf{u}_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{-u}_{1n} & \mathbf{-u}_{2n} & \cdots & \mathbf{0} \\ \end{bmatrix} }_{\mathbf{U}} \circ \underbrace{ \begin{bmatrix} 0 & r_{12}^{-2} & \cdots & r_{1n}^{-2} \\ r_{12}^{-2} & 0 & \cdots & r_{2n}^{-2} \\ \vdots & \vdots & \ddots & \vdots \\ r_{1n}^{-2} & r_{2n}^{-2} & \cdots & 0 \\ \end{bmatrix} }_{\mathbf{R}} \right) \underbrace{ \begin{bmatrix} G\,m_1 \\ G\,m_2 \\ \vdots \\ G\,m_n \\ \end{bmatrix} }_{\mathbf{Gm}} \\ \\ \\ \mathbf{A}=(\mathbf{U} \circ \mathbf{R}) \, \mathbf{Gm} \tag{19}$

Note that $\mathbf{Gm}$ is constant, $\mathbf{U}=-\mathbf{U}^\mathsf{T}$, and $\mathbf{R}=\mathbf{R}^\mathsf{T}$.

One of the properties of the Hadamard product is

$$\mathbf{x^*}(\mathbf{K} \circ \mathbf{L})\mathbf{y}=\mathrm{tr}(\mathbf{D_x^*KD_yL^\mathsf{T}}) \tag{20}$$

If we have $\mathbf{x^*}=\mathbf{1}$, then $\mathbf{D_x^*}=\mathbf{I}$, eq. 20 becomes

$$\mathbf{A}=\mathrm{tr}(\mathbf{UD_{Gm}R}) \tag{21}$$

Implementation

Since $\mathbf{Gm}$ is constant (eq. 19), it can be calculated at the beginning of the simulation and then reused in every iteration to calculate $\mathbf{A}$.

Scale

For solar system simulations:

$$10^{10}\,\mathrm{m} < r < 10^{13}\,\mathrm{m} \\ 10^{3}\,\mathrm{m/s} < v < 10^{5}\,\mathrm{m/s} \\ 10^{22}\,\mathrm{kg} < m < 10^{31}\,\mathrm{kg} \\ 10^{11}\,\mathrm{m^3/s^2} < G\,m < 10^{20}\,\mathrm{m^3/s^2} \\ $$

In order to avoid numbers that are too big or too small, we should strive to combine big numbers with small numbers (like I did with $G\,m$). Also, eq. 17 could be written as

$$\mathbf{a}_i = G\sum_{j=1}^{n} \frac {\mathbf{r}_{ij}} {r_{ij}^3} m_j$$

I could have made $\mathbf{U}_{ij}=\mathbf{r}_{ij}$ and $\mathbf{R}_{ij}=r_{ij}^{-3}$, but in that case $10^{-39}<R_{ij}<10^{-30}$, which is a very small number. So I preferred to formulate the problem as I did, because $\mathbf{U}_{ij}=\mathbf{u}_{ij}$, which is the unitary vector, and $\mathbf{R}_{ij}=r_{ij}^{-2}$ which makes $10^{-23}<R_{ij}<10^{-20}$, that are not so small as the other.

Simulations

In order to assess the effect of using a bigger $\Delta t$, I ran two simulations for the inner solar system, for 300 years with two different values for $\Delta t$: $2\mathrm{s}$ and $20\mathrm{s}$. The commands ran were:

bin/exe files/inner-solar-system 300 2 >/tmp/solar-system.2.out
bin/exe files/inner-solar-system 300 10 >/tmp/solar-system.20.out

In order to generate the graphs for these simulations I ran the following commands:

gnuplot -c files/inner-solar-system.gnuplot /tmp/solar-system.2.out 2
gnuplot -c files/inner-solar-system.gnuplot /tmp/solar-system.20.out 20

The results are shown below.

Note that, for the bigger $\Delta t$, Mercury orbit is expanding a lot. For the smaller $\Delta t$ the expansion is more limited.