Comparing the basins of attraction for several methods in the circular Sitnikov problem with spheroid primaries
Abstract
The circular Sitnikov problem, where the two primary bodies are prolate or oblate spheroids, is numerically investigated. In particular, the basins of convergence on the complex plane are revealed by using a large collection of numerical methods of several order. We consider four cases, regarding the value of the oblateness coefficient which determines the nature of the roots (attractors) of the system. For all cases we use the iterative schemes for performing a thorough and systematic classification of the nodes on the complex plane. The distribution of the iterations as well as the probability and their correlations with the corresponding basins of convergence are also discussed. Our numerical computations indicate that most of the iterative schemes provide relatively similar convergence structures on the complex plane. However, there are some numerical methods for which the corresponding basins of attraction are extremely complicated with highly fractal basin boundaries. Moreover, it is proved that the efficiency strongly varies between the numerical methods.
Keywords:
Sitnikov problem Equilibrium points Oblateness Fractal basin boundaries∎
1 Introduction
A special version of the classical restricted threebody problem is the socalled Sitnikov problem (Sitnikov, 1960). In this case, two equally massed primary bodies move in circular or elliptic orbits, while the test particle oscillates along the vertical axis, perpendicular to the orbital plane of the primaries. The simplest case, where the two primary bodies move in circular orbits is also known as the MacMillan problem (McMillan, 1911). Over the years, a large number of studies have been devoted on the Sitnikov problem, while introducing various types of perturbation such as the radiation pressure (e.g., Perdios & Kalantonis, 2006), the prolateness of the primaries (e.g., Douskos et al., 2012), as well as the oblateness of the primaries (e.g., Rahman et al., 2015).
In dynamical astronomy and celestial mechanics the equilibrium points of a system play a role of great importance since at these locations the test particle is able to maintain its relative position, with respect to the primary bodies. This is true because at the libration points of the system the combined gravitational attraction of the primaries provides precisely the required centripetal force. Unfortunately, in many systems, such as those of the body problem (with ), there are no explicit formulae for the positions of the libration points. Therefore, the locations of the equilibrium points can be obtained only by means of numerical methods. In other words, we need a multivariate iterative scheme for solving the system of the first order derivatives. It is well known that the results of any numerical method strongly depend on the initial conditions (staring points of the iterative procedure). Indeed, for some initial conditions the iterative formulae converge quickly, while for other starting points a considerable amount of iterations is required for reaching to a root (equilibrium point). Fast converging points usually belong to basins of attraction, while on the other hand slow converging points are located in fractal regions. On this basis, the knowledge of the basins of attraction of a dynamical system is very important because these basins reveal the optimal (regarding fast convergence) starting points for which the iterative formulae require the lowest amount of iterations, for leading to an equilibrium point. In addition, being aware of the fractal regions we know exactly which points should be avoided as initial conditions of the iterative formulae.
Furthermore, the basins of attraction, associated with the equilibrium points, contain useful information regarding the intrinsic dynamical properties of a system. This should be true if we take into account that the iterative scheme usually contains both the first and the second order derivatives of the effective potential. It is well know that the first order derivatives are directly linked with the equations of motion of the test particle, while the second order derivatives are used for computing the variational equations. In addition, the variational equations are used for the calculation of the monodromy matrix of the periodic orbits and therefore they are also linked with the stability properties of the test particle. All the abovementioned issues justify the reasons of why one should be aware of the basins of attraction in a dynamical system.
The literature is replete of papers on the basins of attraction in several types of dynamical systems. In the vast majority of them the NewtonRaphson iterative scheme (which is the simplest one) is used for revealing the convergence properties in several dynamical systems such as the Sitnikov problem (e.g., Douskos et al., 2012), the Hill problem with oblateness and radiation pressure (e.g., Douskos, 2010; Zotos, 2017b), the circular restricted threebody problem with oblateness and radiation pressure (e.g., Zotos, 2016), the Copenhagen problem with radiation pressure (e.g., Kalvouridis, 2008), the pseudoNewtonian planar circular restricted threebody problem (e.g., Zotos, 2017c), the circular restricted fourbody problem (e.g., Baltagiannis & Papadakis, 2011; Kumari & Kushvah, 2014; Zotos, 2017a, d), the circular restricted fourbody problem with radiation pressure (e.g., Asique et al., 2016), the circular restricted fourbody problem with various perturbations (e.g., Suraj et al., 2017a, b), the circular restricted fivebody problem (e.g., Zotos & Suraj, 2018), the ring problem of bodies (e.g., Croustalloudi & Kalvouridis, 2007; GousidouKoutita & Kalvouridis, 2009), or even the restricted 2+2 body problem (e.g., Croustalloudi & Kalvouridis, 2013).
In Douskos et al. (2012) the NewtonRaphson basins of attraction of the Sitnikov problem with prolate primaries have been briefly investigated. In the present paper will use a large variety of numerical methods in an attempt to reveal the corresponding basins of convergence. In addition, we will try to evaluate each iterative scheme and therefore obtain a general overview regarding the convergence speed as well as the efficiency of all the numerical methods.
The structure of the present article is as follows: the basic properties of the dynamical model are presented in Section 2. In Section 3 we discuss the parametric evolution of the roots of the system. The following Section contains all the numerical outcomes, regarding the basins of attraction of all the numerical methods. The paper ends with Section 5, where we provide the main conclusions of our analysis.
2 The mathematical model
The dynamical system is mainly composed of two primary bodies, and , which move, around their common center of gravity, in circular orbits. The dimensionless masses of the primaries are and , respectively, where is the mass parameter (Szebehely, 1967). A third body (which behaves as a test particle) is moving in the combined gravitational field of the primaries however, its motion does not perturb the orbits of the two main bodies. This is because the mass of the third body is considerable smaller with respect to the masses of the primaries. In a dimensionless rotating system of coordinates the centers of the two main bodies are located at and , where and .
It is assumed that the primary bodies do not have a spherically symmetric shape but they resemble a spheroid. For this reason we introduce the oblateness coefficient , . When the primary is a prolate spheroid, while when the primary is an oblate spheroid.
According to AbdulRaheem & Singh (2006); Douskos & Markellos (2006); Oberti & Vienne (2003); Sharma & Subba Rao (1975) the timeindependent effective potential of the circular restrictedthree body problem, where the primaries are spheroids is given by
(1) 
where
(2) 
are the distances of the test particle from the respective primary bodies, while is the mean motion of the primaries, which is defined as
(3) 
The equations which govern the motion of a test particle read
(4) 
This dynamical system admits only one integral of motion (also known as the Jacobi integral). The corresponding Hamiltonian is
(5) 
where of course , , and are the velocities, while is the numerical value of the Jacobi constant which is conserved.
In order to obtain the potential function of the circular Sitnikov problem all we have to do is to set , , and in Eq. (1). Then we obtain
(6) 
where . Looking at Eq. (6) we realize that it describes the motion of a massless test particle, oscillating along a straight line which is perpendicular to the orbital plane of the two equally massed primaries. The spatial geometry of circular Sitnikov problem is presented in Fig. 1.
The motion of the test particle, along the vertical axis, is described by the equation
(7) 
while the corresponding energy (Jacobi) integral, regarding the vertical motion, has the form
(8) 
3 The equilibrium points (roots) of the system
Following the approach successfully used in Douskos et al. (2012) (see Section 3), from now on the coordinate is considered as a complex variable and it is denoted by . The transition to complex numbers is imperative because all the impressive fractal basin structures appear only on the complex plane, as it was discussed in Douskos (2010).
In order to locate the positions of the equilibrium points (roots) we have to set the right hand side of Eq. (7) equal to zero as
(9) 
which is reduced to
(10) 
Looking at Eq. (10) we observe that the root is always present, regardless the value of the oblateness coefficient. This root corresponds to the inner collinear equilibrium point of the circular restricted threebody problem. However since the left hand side of Eq. (10) is a fifth order polynomial it means that there are four additional roots, given by
(11) 
The nature of these four roots strongly depends on the numerical value of the oblateness coefficient. Our analysis reveals that, along with the root

When there are two real and two imaginary roots.

When there are two imaginary roots.

When there are four imaginary roots.

When only the root exists.

When there are four complex roots.

When there are two real roots.

When there are four real roots.
It is seen, that the values are in fact critical values of the oblateness coefficient, since they determine the change on the nature of the four roots.
It would be very interesting to determine how the positions of the four roots, on the complex plane, evolve as a function of the oblateness coefficient. Fig. 2 shows the parametric evolution of the four roots , , on the complex plane, when , with and . When the two real roots tend to , while the two imaginary roots tend to infinity. As we proceed to higher values of all four roots tend to the central region. When the two real roots collide at the origin which increases the multiplicity of the root from 1 to 3. At the same time, the two imaginary roots are located at on the vertical axis. As soon as a new pair of imaginary roots emerge from the origin . As the value of increases approaching 0, all four imaginary roots move on collision courses. The collision occurs when , while the roots are exactly at . For positive values of the oblateness coefficient (or in other words for oblate primaries) four complex roots emerge, one at each of the quadrants of the complex plane. As long as lies in the interval the combined traces of the four complex roots create an oval shape. When the four complex roots collide, in two pairs, on the horizontal axis, thus resulting to two real roots of multiplicity 2. For two pairs of real roots emerge, while the roots of each pair move away from each other. Specifically, as the outer roots and tend to infinity, while the roots and tend to .
4 The basins of attraction
A plethora of methods for numerically solving an equation with one variable parameter have been developed over the years. In this article we will consider and compare sixteen methods, whose order of convergence is varying from 2 to 16. In particular the methods under consideration are the following

The NewtonRaphson’s optimal method of second order (Conte & de Boor, 1973).

The Halley’s method of third order (Halley, 1964).

The Chebyshev’s method of third order (Truab, 1964)

The super Halley’s method of fourth order (Gutiérrez & Hernádez, 2001).

The modified super Halley’s optimal method of fourth order (Chun & Ham, 2008).

The King’s method of fourth order (King, 1973).

The Jarratt’s method of fourth order (Jarratt, 1966).

The KungTraub’s optimal method of fourth order (Kung, 1974).

The Maheshwari’s optimal method of fourth order (Maheshwari, 2009).

The Murakami’s method of fifth order (Murakami, 1978).

The Neta’s method of sixth order (Neta, 1979).

The ChunNeta’s method of sixth order (Chun & Neta, 2012).

The NetaJohnson’s method of eighth order (Neta & Johnson, 2008).

The NetaPetkovic’s optimal method of eighth order (Neta & Petković, 2010).

The Neta’s method of fourteenth order (Neta, 1981).

The Neta’s method of sixteenth order (Neta, 1981).
The analytical expressions of all the abovementioned iterative schemes are given in the Appendix of Zotos (2017b).
All the computational methodology that we are going to use in order to classify the initial conditions on the complex plane are described in detail in Section 3 of Zotos (2017b).
4.1 Case I: Two real and two imaginary roots, along with (0,0)
We begin with the first case, that is when , where there are two real and two imaginary roots, along with the root. The basins of attraction on the complex plane, for , using the sixteen iterative schemes are presented in Fig. 3. Our calculations suggest that the convergence structure on the complex plane for the majority of the numerical methods is, in general terms, very similar. In particular, for all the numerical methods except for the super Halley and the King methods we observe the following aspects:

The extent of all the basins of attractions, associated with the five roots, is finite. Furthermore, the shape of all the attracting regions resemble the shape of a lobe.

Around the attracting domains the complex plane is covered by a unified sea of initial conditions (yellow regions) which do not converge to any of the five roots of the system. In fact for these initial conditions all the iterative schemes (except for the super Halley and the King methods) lead, sooner or later, to extremely large numbers, which is a numerical indication that these particular initial conditions lead to infinity.

It is interesting to note that the convergence structure, corresponding to Neta iterative scheme, is a bit different, with respect to the convergence structures of the other thirteen methods.
On the other hand, the convergence structures of both the super Halley and the King methods have significant differences comparing to the rest of the numerical methods. In panel (4) of Fig. 3 it is seen that all the basins of attraction, corresponding to the five roots, extend to infinity. Apart from the usual basins of convergence there exist additional basins (orange regions) in which the initial conditions display a strange behavior. Our analysis indicates that the initial conditions which form these basins converge to , which are not roots of the system. This means that the super Halley method for a considerable amount of initial conditions exhibits a false convergence.
According to panel (6) of Fig. 3 all the basins of attraction also extend to infinity in the case of the King method. Between the convergence regions we can identify several areas (white) in which the initial conditions do not converge to any of the five roots of the system. Initially we suspected that maybe these initial conditions are just extremely slow converging nodes. To check this we increased the maximum allowed number of iterations from 500 to 50000 and we reclassified these initial conditions. We found that for all these initial conditions, during the iterative procedure, the final state smoothly oscillates between two complex numbers of the form . Therefore, we argue that for these initial conditions we have strong numerical evidence that they do not converge to any of the roots of the system.
Another interesting aspect, shown in the convergence diagram of the King method, concerns the geometry of the CCD. More precisely, one can observe that the overall pattern, especially in the vicinity of the basin boundaries, is very noisy or in other words highly fractal. This directly implies that for the initial conditions inside these chaotic domains it is extremely difficult (or even impossible) to know beforehand their final state (root). At this point, it should be noted that when we state that an area is fractal we simply mean that it displays a fractallike geometry, without using any quantitative arguments, such as the fractal dimension, as in Aguirre et al. (2001, 2009).
In Fig. 4 we provide, using tones of blue, the corresponding distributions of the number of the required iterations for obtaining the basins of attraction shown in Fig. 3. In almost all cases we observe the expected behavior according to which the fastest converging nodes are those with initial conditions inside the basins of attraction, while the slowest initial conditions are those located in the vicinity of the basin boundaries. However in panel (4) it becomes evident that a considerable amount of initial conditions require more than 50 iterations for converging to one of the two imaginary roots and , while in all other cases, the vast majority of the initial conditions converge within the first 20 iterations. We suspect that the phenomenon of the extremely slow converging points is directly related with the existence of false converging points. Moreover, it should be emphasized that the regions on the complex plane, in which the extremely slow converging nodes are located, are highly fractal.
The corresponding probability distributions of iterations are given in Fig. 5. The definition of the probability is the following: assume that after iterations initial conditions on the complex plane converge to one of the roots of the system. Then , where is the total number of initial conditions in every CCD. Our results suggest that for almost all the numerical methods more than 98% of the initial conditions converge within the first 30 iterations, while only for the super Halley method the tail of the corresponding histogram extends to 60 iterations. The most probable number of iterations (see the vertical, dashed, red lines in the histograms) seems, in general terms, to decrease as we proceed to numerical methods of higher order.
4.2 Case II: Four imaginary roots, along with (0,0)
In this case, where the system admits four imaginary roots, along with the classical root. In Fig. 6 we provide the CCDs for the sixteen numerical methods, when . Once more, the convergence properties of most of the iterative schemes are very similar, while the only two cases which display complete different patterns are those corresponding to super Halley and King methods. For both these methods the structure of the corresponding CCDs is highly complex.
For the super Halley method it is seen in panel (4) of Fig. 6 that apart from the basins of attraction, associated with the five roots, there exist two additional types of basins. The first type corresponds to initial conditions which display a false convergence to (orange regions), while the second type corresponds to initial conditions that exhibit a false convergence to complex numbers of the form (pink regions). In the case of the King method (see panel (6) of Fig. 6)we have again the appearance of nonconverging initial conditions which infinitely oscillate between two complex numbers of the form .
The corresponding distributions of the number of iterations and the probability are given in Figs. 7 and 8, respectively. It should be mentioned that the most smooth distributions of the probability histograms correspond to numerical methods where either the initial conditions converge to one of the five roots or tend to infinity. The most noisy histograms on the other hand, in which multiple peaks are present (see e.g., panel (4) of Fig. 8) correspond to problematic numerical methods in which the phenomenon of false convergence occurs.
4.3 Case III: Four complex roots, along with (0,0)
The next case under consideration corresponds to , when there are four complex roots, along with the central root . The basins of attraction on the complex plane for the sixteen methods, when are presented in Fig. 9. It is seen that only the CCDs corresponding to super Halley and King methods are different, while almost all the other CCDs have similar convergence patterns.
For the super Halley method, shown in panel (4) of Fig. 9, we observe that the vast majority of the complex plane (orange regions) is covered by initial conditions for which the super Halley iterative method displays a false convergence to . In panel (6), regarding the King method, we encounter again a substantial amount of nonconverging initial conditions for which the iterative scheme oscillates between two imaginary roots. Looking carefully at panel (11) of Fig. 9 one can also identify a small portion of nonconverging initial conditions in the case of the Neta method. Our calculations indicate that for these nodes the iterative scheme oscillates between two complex numbers, which do not coincide with the complex roots of the system.
In Fig. 10 we provide the corresponding distributions of the required number of iterations, while the corresponding probability distributions are shown in Fig. 11. Combining the results of both these figures we may conclude that the highest numbers of required iterations are observed in numerical methods with problematic behavior, where false and nonconverging points are present.
4.4 Case IV: Four real roots, along with (0,0)
We close our numerical investigation with the last case, that is for , when the system admits four real roots, along with the universal root . Fig. 12 depicts the basins of attraction of the sixteen numerical methods, when . As in the previous subsections, the convergence properties of most of the examined numerical methods are, in general terms, very similar.
Once more, for the super Halley method (see panel (4) of Fig. 12) we encountered the phenomenon of initial conditions with false convergence to . In the same vein, the phenomenon of nonconverging initial conditions is also observed for the King method (see panel (6)). At this point we should note the complicated basin structures (with the highly fractal basin boundaries) on the complex plane which are produced by the King numerical method. Nonconverging initial conditions are also present in the case of the Neta method (see panel (11) of Fig. 12). However in this case the corresponding iterative procedure oscillates between two complex numbers, while in the King method there is an infinite oscillation between two imaginary roots.
In Fig. 13 one can observe how the corresponding numbers of the required iteration are distributed on the complex plane, for the numerical methods presented in Fig. 12. In panel (6) of Fig. 13 it is clearly seen how the extremely slow converging points (when ) are located in the vicinity of the fractal basin boundaries. Indeed, in panels (4) and (6) of Fig. 14 we see how the tails of the corresponding histograms are much more extended with respect to histograms of all the other numerical methods.
Before ending this section we would like to state that mainly for saving space we did not present any results regarding the three critical values of the oblateness coefficient . In fact, for these particular values of the system has either one or three roots and therefore the overall structure of the basins of attraction is less interesting.
5 Concluding remarks
In this work we used a large variety of numerical methods in order to reveal the basins of attraction on the complex plane in the circular Sitnikov problem, when the two primary bodies are either prolate or oblate spheroids. All the magnificent basin structures on the complex plane were identified by classifying dense grids of initial conditions, using the corresponding iterative schemes. In particular, we managed to determine how the geometry of the convergence structures changes as a function of the order of the applied numerical methods. Furthermore, all the correlations between the attracting domains and the corresponding distributions of the probability as well as the required number of iterations have been successfully established.
It should be emphasized that this is the first time that the basins of attraction in the circular Sitnikov problem with spheroid primaries are numerically investigates in such a thorough and systematic manner, using a plethora of numerical methods. On this basis, we argue that all the presented numerical results are novel, while they add considerably to our existing knowledge on the field of basins of attraction.
The following list contains the most important conclusions of our numerical analysis:

For almost all the numerical methods (except for the super Halley and the King methods) and for all the examined cases (regarding the nature of the five roots of the system) all basins of attractions are finite. On the other hand, for the super Halley and the King methods all the basins of convergence extend to infinity.

For all the methods where the basins of attractions are finite, it was found that the rest of the complex plane is covered by initial conditions for which the corresponding iterative schemes lead to extremely large numbers (in other words, we have a numerical indication that these nodes lead to infinity).

We observed that not all numerical methods display the same degree of efficiency. In particular, for the super Halley method we identified a nonzero amount of problematic initial conditions which display false convergence to final states which are different from the roots of the system.

Our experiments indicated that for the King and the Neta methods there exist basins of initial conditions which do not converge not even after 50000 iterations. Additional numerical computations revealed that for these initial conditions the corresponding iterative schemes infinitely oscillate between two (imaginary or complex) numbers, which directly implies that these nodes are true nonconverging points.

The most complicated convergence structures on the complex plane, full of highly fractal basin boundaries, correspond to iterative methods with problematic behavior, that is when false or nonconverging points are present.

For the majority of the numerical methods more than 98% of the classified initial conditions converge, to one of the five roots, within the first 25 iterations. Only for the problematic methods (super Halley and King) the required number of iterations extends to more than 30 iterations.

If we exclude the two most problematic methods (super Halley and King) then we may conclude that there is a clear relation between the convergence speed of the iterative schemes and the order of the methods. More specifically, the most probable number of iterations seems to decrease, with increasing order of the method.
A double precision numerical routine, written in standard FORTRAN 77
(Press et al., 1992), was used for the classification of the initial conditions on the complex plane. Using a QuadCore i7 2.4 GHz PC we needed about 4 minutes of CPU time, for performing the classification in each grid of initial conditions. Moreover, all the graphical illustration of the paper has been created using the latest version 11.3 of Mathematica (Wolfram, 2003).
We hope that the presented numerical outcomes, regarding the convergence properties of the circular Sitnikov problem with spheroid primaries, to be useful in the active field of numerical methods and the associated basins of attraction.
Acknowledgments
I would like to express my warmest thanks to the two anonymous referees for the careful reading of the manuscript and for all the apt suggestions and comments which allowed us to improve both the quality and the clarity of the paper.
References
 AbdulRaheem & Singh (2006) AbdulRaheem, A., Singh, J.: Combined effects of perturbations, radiation, and oblateness on the stability of equilibrium points in the restricted threebody problem. Astron. J. 131, 18801885 (2006)
 Aguirre et al. (2001) Aguirre, J., Vallejo, J.C., Sanjuán, M.A.F.: Wada basins and chaotic invariant sets in the HénonHeiles system. Phys. Rev. E 64, 066208 (2001)
 Aguirre et al. (2009) Aguirre, J., Viana, R.L., Sanjuán, M.A.F.: Fractal Structures in nonlinear dynamics. Rev. Mod. Phys. 81, 333386 (2009)
 Asique et al. (2016) Asique, Md.Ch., Prasad, U., Hassan M.R., Suraj Md.S.: On the photogravitational R4BP when the third primary is a triaxial rigid body. Astrophys. Space Sci. 361, 379 (2016)
 Baltagiannis & Papadakis (2011) Baltagiannis, A.N., Papadakis, K.E.: Equilibrium points and their stability in the restricted fourbody problem. Int. J. Bifurc. Chaos 21, 21792193 (2011)
 Chun & Ham (2008) Chun, C., Ham, Y.: Some secondderivativefree variants of superHalley method with fourthorder convergence. Appl. Math. Comput. 195, 537541 (2008)
 Chun & Neta (2012) Chun, N., Neta, B.: A new sixthorder scheme for nonlinear equations. Appl. Math. Lett. 25, 185189 (2012)
 Conte & de Boor (1973) Conte, S.D., de Boor, C.: Elementary Numerical Analysis: An Algorithmic Approach, McGraw Hill Co., New York, 1973
 Croustalloudi & Kalvouridis (2007) Croustalloudi, M.N., Kalvouridis, T.J.: Attracting domains in ringtype Nbody formations. Planet. Space Sci. 55, 5369 (2007)
 Croustalloudi & Kalvouridis (2013) Croustalloudi, M.N., Kalvouridis, T.J.: The Restricted 2+2 body problem: Parametric variation of the equilibrium states of the minor bodies and their attracting regions. ISRN Astronomy and Astrophysics Article ID 281849 (2013)
 Douskos (2010) Douskos, C.N.: Collinear equilibrium points of Hill’s problem with radiation and oblateness and their fractal basins of attraction. Astrophys. Space Sci. 326, 263271 (2010)
 Douskos & Markellos (2006) Douskos, C.N., Markellos, V.V.: Outofplane equilibrium points in the restricted threebody problem with oblateness. Astron. Astrophys. 446, 357360 (2006)
 Douskos et al. (2012) Douskos, C., Kalantonis, V., Markellos, P., Perdios, E.: On Sitnikovlike motions generating new kinds of 3D periodic orbits in the R3BP with prolate primaries. Astrophys. Space Sci. 337, 99106 (2012)
 GousidouKoutita & Kalvouridis (2009) GousidouKoutita, M., Kalvouridis, T.J.: On the efficiency of Newton and Broyden numerical methods in the investigation of the regular polygon problem of bodies. Appl. Math. Comput. 212, 100112 (2009)
 Gutiérrez & Hernádez (2001) Gutiérrez, J.M., Hernández, M.A.: An acceleration of Newton’s method: superHalley method. Appl. Math. Comput. 117, 223239 (2001)
 Halley (1964) Halley, E.: A new, exact and easy method of finding the roots of equations generally and that without any previous reduction. Phil. Trans. Roy. Soc. London 18, 136148 (1964)
 Jarratt (1966) Jarratt, P.: Multipoint iterative methods for solving certain equations. Comput. J. 8, 398400 (1966)
 Kalvouridis (2008) Kalvouridis, T.J.: On some new aspects of the photogravitational Copenhagen problem. Astrophys. Space Sci. 317, 107117 (2008)
 King (1973) King, R.F.: A family of fourthorder methods for nonlinear equations. SIAM Numer. Anal. 10, 876879 (1973)
 Kumari & Kushvah (2014) Kumari, R., Kushvah, B.S.: Stability regions of equilibrium points in restricted fourbody problem with oblateness effects. Astrophys. Space Sci. 349, 693704 (2014)
 Kung (1974) Kung, H.T., Traub, J.F.: Optimal order of onepoint and multipoint iterations. J. Assoc. Comput. Mach. 21, 643651 (1974)
 Maheshwari (2009) Maheshwari, A.K.: A fourth order iterative method for solving nonlinear equations. Appl. Math. Comput. 211, 383391 (2009)
 McMillan (1911) McMillan, W.D.: An integrable case in the restricted problem of three bodies. Astron. J. 27, 1113 (1911)
 Murakami (1978) Murakami, T.: Some fifth order multipoint iterative formulae for solving equations. J. Inform. Process. 1, 138139 (1978)
 Neta (1979) Neta, B.: A sixth order family of methods for nonlinear equations. Int. J. Comput. Math. 7, 157161 (1979)
 Neta (1981) Neta, B.: On a family of multipoint methods for nonlinear equations. Int. J. Comput. Math. 9, 353361 (1981)
 Neta & Johnson (2008) Neta, B., Johnson, A.N.: High order nonlinear solver. J. Comput. Methods Sci. Eng. 8, 245250 (2008)
 Neta & Petković (2010) Neta, B., Petković, M.S.: Construction of optimal order nonlinear solvers using inverse interpolation. Appl. Math. Comput. 217, 24482455 (2010)
 Oberti & Vienne (2003) Oberti, P., Vienne, A.: An upgraded theory for Helene, Telesto, and Calypso. Astron. Astrophys. 397, 353359 (2003)
 Perdios & Kalantonis (2006) Perdios, E.A., Kalantonis, V. S.: Sitnikov motions in the photogravitational restricted threebody problem. In Recent Advances in Astronomy and Astrophysics, 848, 743747 (2006)
 Press et al. (1992) Press, H.P., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P.: Numerical Recipes in FORTRAN 77, 2nd Ed., Cambridge Univ. Press, Cambridge, USA, 1992
 Rahman et al. (2015) Rahman, M.A., Garain, D.N., Hassan, M.R.: Stability and periodicity in the Sitnikov threebody problem when primaries are oblate spheroids. Astrophys. Space Sci. 357, 64 (2015)
 Sharma & Subba Rao (1975) Sharma, R.K., Subba Rao, P.V.: Collinear equilibria and their characteristic exponents in the restricted threebody problem when the primaries are oblate spheroids. Celest. Mech. 12, 189201 (1975)
 Sitnikov (1960) Sitnikov, K.: Existence of oscillating motions for the threebody problem. Dokl. Akad. Nauk. USSR, 133, 303306 (1960)
 Suraj et al. (2017a) Suraj, M.S., Aggarwal, R., Arora, M.: On the restricted fourbody problem with the effect of small perturbations in the Coriolis and centrifugal forces. Astrophys Space Sci. 362, 159 (2017a)
 Suraj et al. (2017b) Suraj, M.S., Asique, M.C., Prasad, U. Hassan, M.R., Shalini, K.: Fractal basins of attraction in the restricted fourbody problem when the primaries are triaxial rigid bodies. Astrophys Space Sci. 362, 211 (2017b)
 Szebehely (1967) Szebehely, V.: Theory of Orbits, Academic Press, New York, 1967
 Truab (1964) Traub, J.F.: Iterative Methods for Solution of Equations. PrenticeHall, Englewood Cliffs, NJ, 1964
 Wolfram (2003) Wolfram, S.: The Mathematica Book, Fifth Edition. Wolfram Media, Champaign, 2003
 Zotos (2016) Zotos, E.E.: Fractal basins of attraction in the planar circular restricted threebody problem with oblateness and radiation pressure. Astrophys. Space Sci. 361, 181 (2016)
 Zotos (2017a) Zotos, E.E.: Revealing the basins of convergence in the planar equilateral restricted fourbody problem. Astrophys. Space Sci. 362, 2 (2017a)
 Zotos (2017b) Zotos, E.E.: Comparing the fractal basins of attraction in the Hill problem with oblateness and radiation. Astrophys. Space Sci. 362, 190 (2017b)
 Zotos (2017c) Zotos, E.E.: Basins of convergence of equilibrium points in the pseudoNewtonian planar circular restricted threebody problem Astrophys. Space Sci. 362, 195 (2017c)
 Zotos (2017d) Zotos, E.E.: Equilibrium points and basins of convergence in the linear restricted fourbody problem with angular velocity. Chaos, Solitons & Fractals 101, 819 (2017d)
 Zotos & Suraj (2018) Zotos, E.E., Suraj, Md.S.: Basins of attraction of equilibrium points in the planar circular restricted fivebody problem. Astrophys. Scpace Sci. 363, 20 (2018)