Integration of Few Body Celestial Systems Implementing Explicit Numerical Methods

An $N$-body configuration, composed of $N$ point masses, is determined by the time-dependent positions, $q^{i}_{\alpha}(t)$, and the momenta, $p_{i\hskip 0.85358pt\alpha}(t)$, for each member of the system. Here, $t$ is the time, $i=1,2,3$ denotes the components and $\alpha=1,...,N$ labels each particle. Such a configuration whose members interact solely via gravitational forces, is defined as a Hamiltonian system of the from

where $G$ is the gravitational constant, which we take to be unity. This system is conservative and theoretically speaking the total energy, E, and total angular momentum, L, are conserved and constant in time. They are constants of motion. Monitoring the time evolution of the values of E and L given a particular system throughout a simulation, is of importance. Because, if E and L are not constant in time, this indicates that the simulation failed. Small deviations are expected, as long as they are in agreement with the error of the respective numerical method. Finally, the $N$-body configuration is a “stiff” problem, thus not all numerical methods give reliable solutions and we would need a special type of stability. (A differential system is said to be stiff, when the Jacobian has at least a large negative eigenvalue.)

The $N$-body problem is quite complicated for $N\geq 3$ it cannot always be solved analytically and a complete analytic solution exists only for $N=2$. For this reason, we resort to numerical solutions. In this project, we are going to use the 4-stage Runge-Kutta explicit iterative method (RK4) to solve the few body problem ($N\leq 9$) in three dimensions butcherBook ; butcherRunge ; gautchi . As it will be noted, this is not the best fitted numerical solution for this problem, mainly if one is interested in long running simulations. Due to the fact that it is explicit and not an (semi-)unconditional method, it is not suitable for stiff problems. We would need an unacceptably small step size level in order to better approximate the theoretical solution, which would make the process quite slow. So, for this reason we will mainly use a relatively small step size. But it is intriguing nonetheless to see how a method, with relatively low computational expense and no noteworthy stability, fairs with such a problem, as an alternative for methods with heavy computational needs and greater stability (e.g. A-stability).

As we can see from Appendix A, by choosing the node and weight parameters $\tau_{i},w_{i},a_{ij}$, we can produce Runge-Kutta methods butcherCoeff that are consistent and have sufficient stability in order to achieve the maximum possible order. Therefore, choosing the parameters in the form of Matrices as below,

we get the 4-stage explicit Runge-Kutta formulae. Best written as,

It has been proven by M. W. Kutta in the 1900s, that such explicit methods of stage $1\leq p\leq 4$, attain maximum order of exactly $p$. That means that, our method has global order of 4 ($\mathcal{O}(h^{4})$), since we used 4-stages. This can be proven by application of Taylor’s expansion for (vector-valued) functions, though for high orders like our method it is technically much more demanding than the lower ones. Therefore, one can use graph-theoretical tools, as established by J.C. Butcher butcherRungeTree .

In this example, we are going to simulate the celestial system of Pluto and its moon Charon, neglecting the gravitational effects from other moons, other celestial objects and planets as they are irrelevant. It is a unique pair, since Charon is such a large moon compared to its parent body that the barycenter of the system lies outside of Pluto which is the primary mass. This means that Pluto orbits about a point which lies at a distance larger than its radius. Charon orbits in an opposite fashion to Pluto about the center of mass, always respecting the conservation laws. From data known, the mass of Charon is $12.2\%$ the mass of Pluto (null , https://solarsystem.nasa.gov/moons/pluto-moons/charon/by-the-numbers/). As for the initial conditions we take those seen in Table 1 and assume that the center of mass is at the point $(0,0,0)$. Pluto’s initial position is taken to be at the arbitrary $(-1,0,0)$ and by use of momentum conservation, we determine the initial position and velocity of Charon. We present part of the trajectory (about three quarters to a complete period) in Fig. 1.

Since the total angular momentum of a bounded system of $N$-bodies is an integral of motion, we get that the third component of the total angular momentum equals

with $L_{1}=L_{2}=0$. We notice that it is constant and independent of time. As we can also see from Fig. 2, this conservation is respected in our simulation, while small changes in values are of order $10^{-11}$, which is in agreement with the error predicted by the RK4 method. The conservation of the angular momentum vector is also consistent with the plane movement seen in Fig. 3, as it is orthogonal to the subspace created by the vectors $\boldsymbol{q}$ and $\boldsymbol{p}$. Further, the energy of the system has also been checked to be constant within the errors of the method used.

Often, in $N$-body systems, one seeks for a simplistic solution of the problem for a qualitative description of the orbits. As an example, for a more precise calculation of the Earth’s orbit, it would suffice to take into consideration only the Sun and the large planets while we may ignore the gravitational effect of smaller celestial objects, like moons, dwarf planets, asteroids and comets. Motivated by such an approximation, in this sub-section we consider a system of three masses with a hierarchy in their values. This is often referred to as the restricted three-body problem pojoman . In particular, we consider a central heavy object representing a star, along with a light planet orbiting around that star and a third even lighter body, which could represent a comet. So, we have the case $m_{\text{star}}\gg m_{\text{planet}}\gg m_{\text{comet}}$. To find an approximate numerical solution to the differential system in Eq. (1) we utilize this mass hierarchy to our benefit. For instance, to calculate the force acted on the planet or the comet we only take into consideration the contribution from the star while the effect of the planet on the comet is a small perturbation which we can be neglected. Such a planet-comet interaction would only affect the orbit of the comet in case the latter encounters the planet closely, in which case the forces would be larger. The corresponding differential system is given by,

The condition for which a planet-comet interaction is not negligible is when the force of the planet on the comet is of the same order as the force of the star on the same comet; i.e. $m_{\text{star}}\left\lVert\boldsymbol{q}_{\text{\tiny comet}}-\boldsymbol{q}_{\text{\tiny planet}}\right\rVert^{2}\sim m_{\text{\tiny planet}}\left\lVert\boldsymbol{q}_{\text{\tiny comet}}-\boldsymbol{q}_{\text{\tiny star}}\right\rVert^{2}$. However, since $m_{\text{\tiny comet}}\ll m_{\text{\tiny planet}}$, the deviation of the planet’s trajectory due to a strong planet-comet interaction will not be affected significantly, and thus we neglect the force of the comet on the planet (the second term in the right-hand-side of Eq. (14b)). Nevertheless, in the time interval for which we have integrated Eq. (14) there are no strong encounters between the planet and the comet and therefore the second term in the right-hand-side of Eq. (14c) is irrelevant in our case. Initial conditions for Eq. (14) are shown in Table 2.

In this last example, we attempt to calculate the orbits of the planets in our Solar System. The Sun comprises nearly all of the matter in the Solar System, reaching up to $99.86\%$ of total mass and is expected to have a nearly zero acceleration in the Solar System frame, https://solarsystem.nasa.gov/solar-system/sun/overview/. According to Kepler’s laws of planetary motion, each object will travel along an ellipse with the Sun being approximately stationary at one of the focal points. Consequently, the closer a planet is to the Sun, the faster it will revolve around it and the smaller its orbital period. Therefore, by the time an outer planet (e.g. Jupiter) completes a full revolution, an inner planet (like Venus) will have already completed several orbits. In particular, in our simulations, we have integrated the system over such a time interval such that Uranus and Neptune orbits have not closed (see the black and purple lines in Fig. 5).

We treat the Solar System as a $N=9$ body problem, where the initial conditions and masses of the planets have all been normalized in terms of the Earth’s. We consider a system of units in which the mass of the Earth is taken to be $m_{\text{\tiny Earth}}=1$. The mass of the Sun and of the other planets are then computed according to the relative mass ratio normalized to $m_{\text{\tiny Earth}}$, https://solarsystem.nasa.gov/planets/overview/. The initial conditions are given in Table 3, where the conversion factor of $0.03$ for the mass has been used so that the ratio between distances and masses is realistic in our numerical model. In order to conserve momentum at the center of mass frame (which almost coincides with the position of the Sun), we take the initial positions and momenta of the planets to cancel out. Finally, we present the trajectories obtained from the numerical calculation in Fig. 5.

We also give the energy over time in Fig. 6, where we notice that the fluctuations displayed are within the error margin of the 4-stage Runge-Kutta method. Nevertheless, it is approximately constant as expected from a theoretical standpoint. The angular momentum was also monitored and found to be conserved within errors and fluctuates like the energy in a similar manner, but it is not depicted here.

Moreover, in Fig 7 we compare the orbits of the inner planets (Mercury, Venus, Earth and Mars) obtained from our numerical model with orbits found online (https://theskylive.com/3dsolarsystem). Looking at Fig. 7, we notice that there are the small differences in the trajectories of the planets. More generally, planets nearest to the Sun do not stray as much as planets that are further away (like Uranus, which is not depicted in Fig. 7) for which we find that the difference in their orbits is larger. Nevertheless, our trajectories were calculated using a computational inexpensive scheme with relatively small step size. Since, as we already said the initial conditions affect our solution curves, they also play a significant role in the deviation of the two sets of orbits. For that reason it was relatively a successful simulation. For better results, a smaller step size would probably not suffice in the long run and a more sophisticated analysis would be needed (see Sec. IV).