On the convergence dynamics of the Sitnikov problem with non-spherical primaries

TL;DR

This study explores how the shape of celestial bodies affects the convergence behavior of the Newton-Raphson method in the Sitnikov problem, revealing the influence of oblateness on convergence speed, efficiency, and basin fractality.

Contribution

It provides a detailed numerical analysis of the impact of oblateness on convergence dynamics and basin fractality in the Sitnikov problem, which was not previously examined.

Findings

Oblateness parameter A significantly affects convergence speed.

Convergence basins exhibit fractal structures influenced by oblateness.

Higher oblateness increases basin boundary complexity.

Abstract

We investigate, using numerical methods, the convergence dynamics of the Sitnikov problem with non-spherical primaries, by applying the Newton-Raphson (NR) iterative scheme. In particular, we examine how the oblateness parameter $A$ influences several aspects of the method, such as its speed and efficiency. Color-coded diagrams are used for revealing the convergence basins on the plane of complex numbers. Moreover, we compute the degree of fractality of the convergence basins on the complex space, as a relation of the oblateness, by using different computational tools, such the fractal dimension as well as the (boundary) basin entropy.

Figures

Equations

$$\Omega (x,y,z)=\sum _{i=1}^2\frac{m_i}{r_i}\left(1+\frac{A_i}{2r_i^2}-\frac{3A_i{z}^2}{2r_i^4}\right)+\frac{n^2}{2}\left(x^2+y^2\right),$$

$$r_1$$

$$r_2$$

$$n=\sqrt{1+3\left(A_1+A_2\right)/2}.$$

$$\ddot{x}$$

$$\ddot{y}$$

$$\ddot{z}$$

$$J=2\Omega (x,y,z)-\left({\dot{x}}^2+{\dot{y}}^2+{\dot{z}}^2\right)=C.$$

$$\Omega (z)=\frac{1}{r}+\frac{A}{2r^3}-\frac{3A{z}^2}{2r^5},$$

$$\ddot{z}=-\frac{z}{r^3}-\frac{9Az}{2r^5}+\frac{15A{z}^3}{2r^7},$$

$$J(z,\dot{z})=2\Omega (z)-{\dot{z}}^2=C_z.$$

$$\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}_{n+1}=\frac{12\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{3}\left(A\left(50-80\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{2}\right)+\left(1+4\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{2}\right)^{2}\right)}{128\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{6}+48\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{4}-6A\left(128\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{4}-96\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{2}+3\right)-1}.$$

$$P(N|l,d)=\frac{1}{2d}\{\begin{array}{ll}\exp \left(-\frac{l-N}{d}\right),& \text{if }N<l\\ \exp \left(-\frac{N-l}{d}\right),& \text{if }N\ge l\end{array},$$

Full text

††journal: International Journal of Applied and Computational Mathematics

1 Introduction

The simplest scenario, according to which the two primary bodies perform planar circular orbits (with zero eccentricity, $e=0$), is known as the MacMillan problem [9], while this notion was originally introduced in [15]. Moreover, the first qualitative results on the Sitnikov problem have been conducted in [20] and [10].

Knowing the exact coordinates of the points of equilibrium of a system is an important issue. However, this is not possible for many complicated dynamical systems for which there are no analytical equations for the positions of the equilibrium points. This automatically means that only by using numerical methods we can obtain the locations of the libration points. As we know, in all numerical methods the initial conditions are very important. Indeed, for some starting points the numerical methods may converge relatively fast to a root, while for other initial conditions they may require a considerable amount of iterations. Usually, points with fast convergence belong to the basins of convergence (Boc), while slow converging points are situated in the vicinity of the fractal basin boundaries. Therefore, it is of great importance to know the location of the Boc for a dynamical system, because then we automatically are aware of the optimal initial conditions for the numerical methods. Here, we would like to point out that the Boc of a dynamical system strongly depend on the chosen numerical method. In other words, different numerical methods yield to completely different Boc, for the same dynamical system (see e.g., [23, 24]).

The Boc define complicated geometrical structures on the complex or the configuration $(x,y)$ plane. Another important aspect is knowing the degree of fractality of the convergence regions. A quantitatively estimation of the degree of fractality can be easily achieved by computing several numerical indicators, such as the uncertainty or fractal dimension (see e.g., [2, 3]) or the basin entropy (see e.g., [4]). Both these quantities can provide safe results, regarding the degree of fractality of a dynamical system.

Our article has the following layout: in Section 2 we describe the mathematical formulation of the dynamical model. The next Section 3 contains all the numerical outcomes of our work, about the properties of the Sitnikov problem with non-spherical primaries. In the last Section 4, we emphasize the conclusions of our computational analysis.

2 Mathematical formulation of the system

The system consists of two primaries whose dimensionless masses are $m_{1}=\mu$ and $m_{2}=1-\mu$, where $\mu=m_{2}/(m_{1}+m_{2})\leq 1/2$ is the well known mass parameter [21]. The centres of the two primaries lie on the horizontal $Ox$ axis and in particular at $(x_{1},0,0)$ and $(x_{2},0,0)$, where $x_{1}=-\mu$ and $x_{2}=1-\mu$. For each primary it is assumed that its shape resembles a spheroid, according to the value of the corresponding oblateness $A_{i}$, $i=1,2$.

According to [1, 7, 19], the function of the potential of the restricted circular problem with two oblate primaries is given by

[TABLE]

where

[TABLE]

are the respective distances between the test particle (third body), and the centers of the two primaries. Moreover, the mean motion is

[TABLE]

The third body of the system moves according to the following equations

[TABLE]

For this system there is only one known motion integral, which reads

[TABLE]

By setting $x=y=0$, $\mu=1/2$, and $A_{1}=A_{2}=A$ in Eq. (1) we obtain the potential function of the Sitnikov problem

[TABLE]

where $r=\sqrt{z^{2}+1/4}$.

Thus, the vertical motion, along the $z$ axis, of the test particle is governed by the following equation

[TABLE]

while we can deduce that the corresponding Jacobi integral becomes

[TABLE]

In Paper I we demonstrated that the points of equilibrium (roots) of the Sitnikov problem can be obtained through the equation of mition (7). In addition, we seen that the value of the oblateness $A$ greatly influences the nature of the roots. Specifically:

• •

For $A<-1/18$ there exist two pairs of real and imaginary roots.

• •

For $A=-1/18$ there exists a pair of imaginary roots.

• •

For $A\in(-1/18,0)$ there exist two pairs of imaginary roots.

• •

For $A=0$ the root \ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}=0 is the only root.

• •

For $A\in(0,5/6)$ there exist two pairs of complex roots.

• •

For $A=5/6$ there exists a pair of real roots.

• •

For $A>5/6$ there exist two pairs of real roots.

We also concluded that the levels $A=\{-1/18,0,5/6\}$ are in fact critical levels of the oblateness.

3 Numerical results of the basins of convergence

The Boc on the plane of complex numbers can be determined by means of the numerical method of Newton-Raphson (NR). In Paper I we shown that the corresponding iterative scheme reads

[TABLE]

At this point, we should emphasize that from now on the coordinate $z$ is treated as a complex variable

$z$

  , while the same approach was also successfully followed in [8, 23, 24, 25]. In [6] it was demonstrated that the use of complex variables is necessary, because all the beautiful and impressive Boc, with the basin boundaries with fractal-like geometry, appear only on the plane of complex numbers.

In Paper I we presented in detail the structure of the Boc in the four intervals of $A$. In Fig. 2 we remind to the reader the structure as well as the geometry and the shape of the convergence regions, for four characteristic values of the oblateness $A$. It is seen, that in all cases the Boc have finite area. Moreover, the majority of the plane of complex numbers is occupied by initial conditions (yellow regions) for which the NR iterative method quickly diverges to very large complex numbers, thus numerically indicating divergence to infinity. The distributions of the required iterations $N$ are presented in Fig. 3. It is evident that for initial conditions $(\mathcal{R},\mathcal{I})$ near the roots the required iterations are low $(N\approx 5)$, while near the boundaries of the basins the NR iterative scheme needs more than 15 iterations for obtaining a root, with the predefined accuracy.

Fig. 4(a-d) shows the corresponding histograms with the distributions of probability. The definition of the probability is $P=N_{0}/N_{t}$, where $N_{0}$ is the number of the initial conditions $(\mathcal{R},\mathcal{I})$ which display true convergence, while $N_{t}$ is the number of the total nodes on the plane of complex numbers.

The histograms shown in panels (a-d) of Fig. 4 can be used for extracting more information, regarding the convergence properties of the NR method. For instance, we can use the Laplace distribution for obtaining the best fits of the right-hand sides (tails) of the histograms (see blue solid lines). We choose to use the Laplace distribution because this is the most natural choice particularly in systems displaying transient chaos (e.g., [11, 17, 18]).

The Laplace probability density function (PDF) is given by

[TABLE]

where the parameters $l$ and $d>0$ are the location parameter and the parameter of diversity, respectively. From the PDF we need only the $N\geq l$ part because the Laplace distributions refer only to the tails of the probability histograms.

In Paper I we investigated using numerical methods the convergence dynamics of the Sitnikov problem with non-spherical primaries however for some individual values of the oblateness $A$. In the present work, we will perform a more systematic numerical analysis in an attempt to determine how $A$ affects the convergence properties of the system. For this task, we classified 1000 grids of $1024\times 1024$ starting points $(\mathcal{R},\mathcal{I})$ inside the square region $R=[-2,2]\times[-2,2]$ on the plane of complex numbers, for the range $A\in[-0.5,1]$, thus following the pioneer works of Nagler [12, 13]. In our calculations the desired accuracy, regarding the coordinates of the attractor, was set to $10^{-15}$, while the maximum allowed number of iterations was $N_{\rm max}=500$.

The evolution of the most probable number of iterations $N^{*}$, per grid, is illustrated in panel (a) of Fig. 5. It is seen, that around the critical levels $A=-1/18$ and $A=5/6$ we have two peaks, while on the other hand around the critical level $A=0$ we observe the lowest value of $N^{*}$. In Paper I we had discussed (see bottom row of Fig. 4) that when $A=-1/18$ for a set of starting points the NR scheme requires a significant number of iterations for converging to the attractor \ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}=0. This is exactly why in part (a) of Fig. 5 we see that the highest value of $N^{*}$ is measured near $A=-1/18$.

In parts (b) and (c) of Fig. 5 we present the parametric evolution of the location parameter $l$ and the diversity $d$, respectively, as a function of the oblateness $A$. In part (b) we also included, for comparison reasons, using blue color the evolution of the $N^{*}$ of iterations. One can observe that in general terms the location parameter almost coincides with the average number of iterations (almost always $|l-<N>|\leq 2$). This implies that the Laplace probability density function (PDF) can satisfactorily fit the tails of the probability histograms. According to the plot shown in part (c) of Fig. 5 the diversity is, in most of the cases, low $(d<3)$, thus indicating the dispersion of the values of $N$ is very close to $N^{*}$. On the other hand, in the vicinity of the critical levels of the oblateness $A=-1/18$ and $A=0$ the diversity exhibits a local maximum and a local minimum, respectively. Part (d) of Fig. 5 illustrates the evolution of the differential entropy $h=1+\ln(2d)$, where $d$ is the diversity. One can see, that the evolution of both $d$ and $h$ displays similar overall patterns.

The numerical analysis presented in Paper I revealed the fractal regions on the plane of complex numbers. One of the most convenient ways of measuring the degree of fractality of a system is by computing the uncertainty or fractal dimension $D_{0}$ (see e.g., [14]), thus following the computational methodology used in [2, 3]. Fig. 6 shows the parametric evolution of the uncertainty dimension, as a function of the oblateness $A$. In the first interval of $\delta$ the fractal dimension is very close to 1, which implies zero fractality. In the first two intervals $D_{0}$ decreases, while in the third interval its value is mostly reduced. It is seen, that $D_{0}$ displays the maximum value near the value $A=0.7$, while the lowest value is measured at $A=0$.

Another efficient way for quantitatively measuring the degree of fractality of a system is by computing the so-called basin entropy [4, 5]. This method determines the fractality of a basin diagram by the process of examination of its topological properties. The evolution of the entropy of the basins $S_{b}$, as a function of $A$, is illustrated in panel (a) of Fig. 7. Once more, we note that in the vicinity of all the critical values of $A$ there are three local minima of $S_{b}$, mainly because for these values of $A$ the total number of the numerical attractors decreases from 5 to three (when $A=-1/18$ and $A=5/6$) and one (when $A=0$). The maximum value of $S_{b}$ was measured near $A=0.7$. Therefore, we may argue that two different methods (i.e., the uncertainty dimension and the basin entropy) indicate that the degree of the fractality of the Boc on the plane of complex numbers is maximum near the same value of $A$.

Apart from the basin entropy there is also the boundary basin entropy $S_{bb}$ [4], from which we can extract additional information about the degree of fractality of the Boc. The parametric evolution of $S_{bb}$ is given in part (b) of Fig. 7. From this type of plot we can also deduce information regarding the fractality of the Boc on the plane of complex numbers. More specifically, we can use the so-called “log 2 criterion”, according to which if $S_{bb}>\log 2$ then the basin boundaries are certainly fractal (here note that the converse statement is not valid). As it is seen in part (b) of Fig. 7 the basin boundaries are certainly fractal only when $0.65<A<0.78$. Once more, the lowest values of $S_{bb}$ are reported in the vicinity of the critical values of the oblateness.

At this point, we would like to briefly discuss the efficiency of the NR method. The classification of the 1000 grids of initial conditions suggested that at least for this dynamical system, the numerical method does display an ill behaviour, for some specific sets of initial conditions. In particular, for a large amount of initial conditions the iterator quickly diverges to extremely large complex numbers. In part (a) of Fig. 8 we give the parametric evolution of the area $R$ on the plane of complex numbers, covered by starting points which diverge to infinity. It is seen, that the highest value of $R$ is observed at $A=0$, while for higher values of $A$ the value of the area decreases.

Before closing this section, it would be very illuminating to discuss how the oblateness influences the geometry of the Boc. As we have already seen, the converge regions on the complex plane have finite area. Therefore, we define as $\mathcal{R}_{max}$ and $\mathcal{I}_{max}$ the maximum values of the Boc, along both axes, respectively, while the ratio $r$ is defined as $\mathcal{I}_{max}/\mathcal{I}_{max}$. Part (b) of Fig. 8 shows the evolution of $\mathcal{R}_{max}$, $\mathcal{I}_{max}$ and $r$, as a function of the oblateness. One can see, that for negative values of $A$ the overall structure of the Boc is elongated along the vertical axis. On the contrary, for $A>0$ the value of $\mathcal{R}_{max}$ increases and at about $A=0.15$ we have that $\mathcal{R}_{max}=\mathcal{I}_{max}$. Moreover, for $A>0.15$ the value of $\mathcal{R}_{max}$ increases rapidly which implies that the overall shape of the convergence regions becomes elongated along the horizontal axis.

4 Concluding remarks

The present paper can be considered as a continuation of Paper I. The scope of this work was to study the convergence dynamics of the Sitnikov problem with non-spherical primaries. Using the NR iterative method we revealed the Boc on the plane of complex numbers, by means of color-coded diagrams. Moreover, we demonstrated how the oblateness $A$ influences the speed and the accuracy of the method. At the same time, it was determined how the same parameter affects the degree of fractality of the convergence regions, by computing modern quantitative indices, such as the uncertainty dimension and the (boundary) basin entropy.

In this work we demonstrate for the first time how the oblateness of the Sitnikov problem with spheroidal primaries influences the overall properties of the system. Additionally, we relate also for the first time different techniques for measuring the fractal degree of a dynamical system. More specifically, we computed and compared the results of both the fractal dimension and the (boundary) basin entropy. On this basis, we claim that the resented outcomes of the article are interesting and novel and add new information on the field of the convergence properties of Hamiltonian systems.

The most important findings of our numerical exploration are listed here:

For the majority of the studied cases, regarding the value of $A$, the NR iterative scheme requires an average number of 8 iterations for leading to one of the numerical attractors (roots), while only near the critical values of $A$ the average number of iterations increases. 2. 2.

Our numerical analysis reported the existence of only starting points for which the numerical iterator quickly diverges to infinity. On the other hand, there is no indication of false or non-converging starting points on the plane of complex numbers. 3. 3.

The highest fractal degree (measured by using the uncertainty dimension as well as the entropy of the basins) of the convergence diagrams on the plane of complex numbers correspond to about $A=0.7$, while the lowest values of the degree of fractality were measured near the critical levels of $A$, where the number of the numerical attractors is reduced. 4. 4.

Exploiting the information of the computation of the boundary basin entropy and the “log 2 criterion” we proved that in this system the basin boundaries of the convergence regions on the plane of complex numbers are certainly fractal for $0.65<A<0.78$. 5. 5.

For negative values of the oblateness the overall shape of the convergence regions on the plane of complex numbers is elongated along the vertical axis, while for $A>0.15$ the geometry changes and the shape becomes elongated along the horizontal axis.

The numerical routine of the NR iterative method was written in the standard version of FORTRAN 77 [16]. For classifying the nodes on the plane of complex numbers, roughly about 2.5 minutes of CPU time, per grid, was required using an Quad-Core Intel i7 vPro 4.0 GHz processor. The version 11.3 of Mathematica*®* [22] has been deployed for constructing all the graphics of the paper.

Compliance with Ethical Standards

• •

Funding: The authors state that they have not received any research grants.

• •

Conflict of interest: The authors declare that they have no conflict of interest.