Basins of convergence of equilibrium points in the generalized Henon-Heiles system
BBasins of convergence of equilibrium points in the generalized H´enon-Heiles system
Euaggelos E. Zotos a, ∗ , A. Ria˜no-Doncel b , F. L. Dubeibe b a Department of Physics, School of Science, Aristotle University of Thessaloniki, GR-541 24, Thessaloniki, Greece b Facultad de Ciencias Humanas y de la Educaci´on, Universidad de los Llanos, Villavicencio, Colombia
Abstract
We numerically explore the Newton-Raphson basins of convergence, related to the libration points (which act as attractorsof the convergence process), in the generalized H´enon-Heiles system (GHH). The evolution of the position as well as ofthe linear stability of the equilibrium points is determined, as a function of the value of the perturbation parameter.The attracting regions, on the configuration ( x, y ) plane, are revealed by using the multivariate version of the classicalNewton-Raphson iterative algorithm. We perform a systematic investigation in an attempt to understand how theperturbation parameter affects the geometry as well as of the basin entropy of the attracting domains. The convergenceregions are also related with the required number of iterations.
Keywords:
Generalized H´enon-Heiles system, Equilibrium points, Basins of convergence, Fractal basins boundaries
1. Introduction
In every two-dimensional Hamiltonian system the de-termination of the fixed points is a necessary step in orderto diagnose the global dynamical behavior of the system.The usual way to do so is by solving the algebraic systemof equations ∂V∂x = ∂V∂y = 0 , (1)where V denotes the effective potential function, associ-ated to the physical system. In the majority of cases, it iscommonly necessary to solve numerically the system (1)via root finding algorithms, which use iterations and re-quire one or more initial guesses of the root as starting val-ues. In this realm, it is important to consider the possibleappearance of two issues: the non-convergence of a giveninitial condition and the number of iterations needed toreach a given fixed point. These two issues can be treatedby means of the so called basin of convergence, which givesus information about: (i) the total number of fixed pointsin the system, (ii) the non-converging initial conditionsand (iii) the number of iterations needed to reach a givenfixed point.In a general sense, the basins of convergence can be un-derstood as the set of initial guesses (conditions) that af-ter iterations, via a root-finding algorithm, tend to a fixedpoint (see e.g., [14] for the case of complex polynomials)and must not be confused with the concept of basin of at-traction, defined as the set of initial conditions that lead toa specific attractor [17], because the last one is meaning-less in the context of conservative Hamiltonian systems. ∗ Corresponding author
Email address: [email protected] (Euaggelos E. Zotos)
As in the case of basins of attraction, the basins of conver-gence could be composed by more than one basin (associ-ated to each fixed point) separated by boundaries whichcan be smooth or fractal curves. During the last yearsthe basins of convergence have been investigated in manydifferent dynamical systems, using the Newton-Raphsoniterative method, e.g. for the Hill’s problem with radi-ation and oblateness [9], the Sitnikov problem [10], therestricted four-body problem [4], the Copenhagen prob-lem with magnetized primaries [15], the restricted four-body problem with oblate primaries [16], the restrictedfour-body problem with radiation pressure [3], the planarcircular restricted three-body problem with oblateness andradiation pressure [21], the pseudo-Newtonian planar cir-cular restricted three-body problem [22], the ring problemof N + 1 bodies [5, 12], or even the restricted 2 + 2 bodyproblem [6].The H´enon-Heiles Hamiltonian is a two-dimensionaltime-independent dynamical system, originally proposedas a simplified version of the gravitational potential ex-perimented by a star orbiting around an axially symmet-ric galaxy (see e.g., [13]). An extension of this potentialup to the fourth-order was performed by Verhulst [19],while in [11] we generalized the H´enon-Heiles potential upto the fifth-order. In the present paper, we examine thebasins of convergence in the fifth-order generalization ofthe H´enon-Heiles Hamiltonian (in all that follows GHH).The GHH is highly nonlinear, so an adequate tool to an-alyze the convergence properties of iteration functions isthe Newton-Raphson algorithm. On the other hand, thefractality of the basins will be analyzed through the basinentropy, a new measure introduced recently to quantifythe uncertainty of basins (e.g. of escape, convergence orattraction) [7]. Published in International Journal of Non-Linear Mechanics March 31, 2018 a r X i v : . [ n li n . C D ] M a r he present paper has the following structure: themain properties of the dynamical system are presented inSection 2. The parametric evolution of the position aswell as the stability of the equilibrium points is investi-gated in Section 3. Section 4 contains the most relevantresults regarding the evolution of the Newton-Raphsonbasins of convergence. In Section 5 we demonstrate howthe basin entropy of the configuration ( x, y ) convergenceplanes evolves, as a function of the perturbation parame-ter. Our paper ends with Section 6, where we summarizethe main conclusions of this work.
2. Properties of the dynamical system
Let us briefly recall the main properties of the gen-eralized H´enon-Heiles system (GHH). The correspondingpotential is given by U ( x, y ) = 12 (cid:0) x + y (cid:1) + x y − y δ (cid:104) x y + x y − y − (cid:0) x + y (cid:1) (cid:105) , (2)where δ is a perturbation parameter. Note that when δ = 0the potential (2) is automatically reduced to that of theclassical H´enon-Heiles system.Potential (2) is derived as a Taylor expansion up to the5-th order of a general potential with axial and reflectionsymmetries. More information about the exact expansionand the derivation of the generalized potential is given in[11].The equations of motion read¨ x = − ∂U∂x , ¨ y = − ∂U∂y , (3)where U x ( x, y ) = ∂U∂x = x (1 + 2 y )+ 2 δx (cid:0) x ( y −
1) + y ( y − (cid:1) ,U y ( x, y ) = ∂U∂x = x + y (1 − y )+ δ (cid:0) x + x y (3 y − − y (5 y − (cid:1) . (4)In the same vein, the second order derivatives of the poten-tial U ( x, y ), which will be needed later for the computationof the multivariate Newton-Raphson iterative scheme areas follows U xx ( x, y ) = ∂ U∂x = 1 + 2 y + 2 δ (cid:0) x ( y −
1) + y ( y − (cid:1) ,U xy ( x, y ) = ∂ U∂x∂y = 2 x + δ (cid:0) x + 2 xy (3 y − (cid:1) ,U yx ( x, y ) = ∂ U∂y∂x = U xy ( x, y ) ,U yy ( x, y ) = ∂ U∂y = 1 − y + 2 δ (cid:0) x (3 y − − y (5 y − (cid:1) . (5) The Hamiltonian, which dictates the motion of the testparticle, is given by H ( x, y, ˙ x, ˙ y ) = U ( x, y ) + 12 (cid:0) ˙ x + ˙ y (cid:1) = E, (6)where ˙ x and ˙ y are the velocities, associated to the coordi-nates x and y , respectively, while E is the numerical valueof the total orbital energy of the test particle, which isconserved.
3. Parametric evolution of the equilibrium points
The necessary and sufficient conditions that must besatisfied for the existence of coplanar equilibrium points,are ˙ x = ˙ y = ¨ x = ¨ y = 0 . (7)On the other hand, the positions of the libration pointscan be determined by numerically solving the system ofequations U x ( x, y ) = U y ( x, y ) = 0 . (8)From Eqs. (2) and (8), it can be easily noticed that thetotal number of libration points in the GHH is a functionof the perturbation parameter δ . In particular • When δ ∈ [0 , . • When δ ∈ [0 . ,
4) there exist eight equilib-rium points (see panel (b) of Fig. 1). • When δ = 4 there are nine libration points (see panel(c) of Fig. 1). • When δ > δ for which six equilibrium points are present is an irrationalnumber and therefore it cannot be exactly determined.As shown in panels (a)-(d) of Fig. 1 for (a): δ = 0 . δ = 2 .
5, (c): δ = 4, (d): δ = 7, the intersection pointsof the curves corresponding to the first order derivatives(8), shall denote the position of the libration points L i , i = 1 , ...,
10. Furthermore, in Fig. 2 we present how thenumber and the exact positions of the equilibrium pointsevolve as a function of the value of the perturbation pa-rameter.In Fig. 3 we present the parametric evolution of thepositions of the equilibrium points, on the configuration( x, y ) plane, when δ ∈ [0 , L remains unperturbed at theorigin (0 , δ ≥ . δ >
4. Our computations indicate that inthe limit δ → ∞ L , L , L , L , L , and L tend to collidewith L , while on the other hand L , L , and L movefar away from the center.2 igure 1: Positions (red dots) and numbering of the equilibrium points ( L i , i = 1 , ...,
10) through the intersections of U x = 0 (green) and U y = 0 (blue), when (a-upper left): δ = 0 . δ = 2 . δ = 4(nine equilibrium points) and (d-lower right): δ = 7 (ten equilibrium points). (Color figure online). Following the procedure outlined in [22], once we knowthe positions of the libration points, we can determine theirlinear stability by means of the characteristic equation. Todo so, we defined a dense, uniform sequence of 10 values of δ in the interval [0 ,
10] and numerically solved the system (8), thus determining the coordinates ( x , y ) of the equi-librium points. The last step was to insert the coordinatesof the equilibria into the characteristic equation and deter-mine the nature of the four roots. The above-mentionednumerical analysis reveals that all the equilibrium points3 igure 2: The variation of the positions of the equilibrium points (red dots) and the contours defined by the equations U x = 0, U y = 0, as afunction of the perturbation parameter δ . (Color figure online). L i , i = 2 , ...,
10 are always linearly unstable, while L isalways Lyapunov stable, when δ >
4. The basins of convergence
One of the most well-known numerical methods forfinding successive approximations to the roots of nonlinearequations is the Newton-Raphson method. This methodis applicable to systems of multivariate functions f ( x ) = 04 igure 3: The evolution of the equilibrium points, L i , i = 1 , ..., δ ∈ [0 , δ increases is indicated by thearrows. The black dots (points A, B, and C) correspond to δ → δ = 0 . δ = 4, respectively. (Color figure online). through the iterative scheme x n +1 = x n − J − f ( x n ) , (9)where f ( x n ) denotes the system of equations, while J − is the corresponding inverse Jacobian matrix. In our caseEqs. (8) describe the system of the differential equations.The above-mentioned iterative scheme can be decom-posed for each coordinate x and y , as follows x n +1 = x n − (cid:18) U x U yy − U y U xy U yy U xx − U xy (cid:19) ( x n ,y n ) ,y n +1 = y n + (cid:18) U x U yx − U y U xx U yy U xx − U xy (cid:19) ( x n ,y n ) , (10)where x n , y n are the values of the x and y coordinates atthe n -th step of the iterative process.The philosophy behind the Newton-Raphson methodis the following: An initial condition ( x , y ), on the con-figuration plane activates the code, while the iterative pro-cedure continues until an equilibrium point (attractor) isreached, with the desired predefined accuracy. If the par-ticular initial condition leads to one of the libration pointsof the system it means that the numerical method con-verges for that particular initial condition. At this point,it should be emphasized that in general terms the methoddoes not converge equally well for all the available initialconditions. The sets of the initial conditions which leadto the same root compose the so-called Newton-Raphsonbasins of convergence or even attracting domains/regions. Nevertheless, as pointed out in the Introduction section,it should be clarified that the Newton-Raphson basins ofconvergence should not be confused, by no means, withthe basins of attractions which are present in systems withdissipation.From the iterative formulae of Eqs. (10) it becomesevident that they should reflect some of the most intrinsicproperties of the Hamiltonian system. This is true if wetake into account that they contain the derivatives of bothfirst and second order of the effective potential U ( x, y ).A double scan of the configuration ( x, y ) plane is per-formed for revealing the structures of the basins of conver-gence. In particular, a dense uniform grid of 1024 × x , y ) nodes is defined which shall be used as initial con-ditions of the iterative scheme. The number N of theiterations, required for obtaining the desired accuracy, isalso monitored during the classification of the nodes. Forour computations, the maximum allowed number of iter-ations is N max = 500, while the iterations stop only whenan attractor is reached, with an accuracy of 10 − .The Newton-Raphson basins of convergence when δ =0 (which correspond to the classical H´enon-Heiles system)are presented in panel (a) of Fig. 4. Different colors areused for each basin of convergence, while the positions ofthe four equilibrium points (attractors) are indicated byblack dots. It is seen that the 2 π/ x, y ) plane, by thegeometry of the attracting domains. The distribution ofthe corresponding number N of iterations is given in panel(b) of the same figure, using tones of blue.In the following subsections, we will determine howthe perturbation parameter δ affects the structure of theNewton-Raphson basins of convergence in the GHH, byconsidering two cases regarding the total number of the at-tractors (equilibrium points). For the classification of thenodes in the configuration ( x, y ) plane we will use color-coded diagrams (CCDs), in which each pixel is assigned adifferent color, according to the final state (attractor) ofthe corresponding initial condition. Our investigation begins with the first case where fourequilibrium points are present, that is when 0 < δ ≤ . δ > π/ x = 0 axis. As we pro-ceed to higher values of the perturbation parameter thevast majority of the CCDs remains almost unperturbed.The only significant changes appear mainly in the vicinityof the basin boundaries. More precisely, the unpredictabil-ity in these regions increases rapidly, since they display ahighly fractal geometry . When δ = 0 . With the term fractal we refer to the fractal-like geometry of theregion, without conducting any additional calculations as in [1, 2]. igure 4: (a-left): The Newton-Raphson basins of convergence on the configuration ( x, y ) plane for the classical H´enon-Heiles system, where δ = 0. The four equilibrium points are indicated by black dots. The initial conditions, leading to a certain equilibrium point, are markedusing the following color code: L (green); L (red); L (blue); L (magenta); non-converging points (white). (b-right): The distribution ofthe corresponding number N of required iterations for obtaining the Newton-Raphson basins of attraction shown in panel (a). (Color figureonline). of Fig. 5) and also when δ = 0 . N (cid:29) N of it-erations is provided, using tones of blue, in Fig. 6(a-f).We see that initial conditions inside the basins of attrac-tion converge relatively quickly ( N <
N >
20) are those in the vicinityof the basin boundaries.
In the case where 0 . ≤ δ <
4, there are eightequilibrium points on the configuration ( x, y ) plane. InFig. 7(a-i) we present the Newton-Raphson basins of con-vergence for nine values of the perturbation parameter.In panels (a) and (b) of Fig. 7, which correspond to δ = 0 . δ = 0 . δ = 0 . δ = 0 . L and L . We believe that the complete absence of converging points to attractors L and L is somehowrelated with the presence of non-converging points. For δ = 3 .
99 (see panel (h) of Fig. 7) and δ = 3 . N max = 50000 all these initial con-ditions converged, sooner or later, to one of the first twoattractors of the system ( L or L ).The distribution of the corresponding number N of it-erations needed to achieve the desired accuracy in our com-putations is presented in Fig. 8(a-i). Looking at panels (g)and (h) of Fig. 8 it is evident that the required number ofiterations, for initial conditions inside the basins locatedat the lower part of the CCDs, increases rapidly as weapproach the second critical value of the perturbation pa-rameter. The last case under consideration corresponds to theregion δ ≥
4, where there are nine or ten equilibriumpoints. In Fig. 9 we present, through the correspondingCCDs, the Newton-Raphson basins of convergence for ninevalues of the perturbation parameter. When δ = 4 . L . However, as the value of δ increases, thus moving away from the critical level, theamount of non-converging points decreases rapidly and6 igure 5: The Newton-Raphson basins of convergence on the configuration ( x, y ) plane for the first case, where four equilibrium points arepresent. (a): δ = 0 .
40; (b): δ = 0 .
45; (c): δ = 0 .
60; (d): δ = 0 .
75; (e): δ = 0 . δ = 0 . L (green); L (red); L (blue); L (magenta);non-converging points (white). (Color figure online). by δ = 0 .
77 (see panel (d) of Fig. 9) they have com-pletely disappeared. Our numerical experiments indicatethat these initial conditions must be true non-converging points, since their portion, in each CCD, remains unper-turbed by the shift of the number of allowed iterations upto
Nmax = 50000. In general terms it is seen that the con-7 igure 6: The distribution of the corresponding number N of required iterations for obtaining the Newton-Raphson basins of convergenceshown in Fig. 5(a-f). The non-converging points are shown in red. (Color figure online). figuration ( x, y ) plane is dominated by well-formed unified basins of convergence (corresponding mainly to libration8 igure 7: The Newton-Raphson basins of convergence on the configuration ( x, y ) plane for the second case, where eight equilibrium pointsare present. (a): δ = 0 . δ = 0 . δ = 0 .
77; (d): δ = 1; (e): δ = 1 .
5; (f): δ = 2; (g): δ = 3; (h): δ = 3 .
99; (i) δ = 3 . L (green); L (red); L (blue); L (magenta); L (cyan); L (yellow); L (olive); L (teal); non-converging points (white). (Color figure online). points L , L , and L ). Moreover, the most noticeablechange with increasing value of δ is the fact that basinboundaries become less noisy (fractal).The corresponding distributions of the required num- ber N of iterations are given in Fig. 10(a-i). It is inter-esting to note in panel (a) of Fig. 10 that the requirediterations number, corresponding to the equilibrium point(attractor) L , is high N >
25. We believe that this be-9 igure 8: The distribution of the corresponding number N of required iterations for obtaining the Newton-Raphson basins of convergenceshown in Fig. 7(a-f). The non-converging points are shown in red. (Color figure online). igure 9: The Newton-Raphson basins of convergence on the configuration ( x, y ) plane for the third case, where nine or ten equilibriumpoints are present. (a): δ = 4; (b): δ = 4 . δ = 4 . δ = 4 .
5; (e): δ = 5; (f): δ = 5 .
15; (g): δ = 5 .
5; (h): δ = 6; (i) δ = 10.The positions of the equilibrium points are indicated by black dots. The color code, denoting the eight attractors, is as follows: L (green); L (red); L (blue); L (magenta); L (cyan); L (yellow); L (olive); L (teal); L (purple); L (orange); non-converging points (white).(Color figure online). havior can be explained, in a way, if we take into accountthat δ = 4 is a critical value and, as we have seen so far,near the critical values of the perturbation parameter the iterations number grow significantly.11 igure 10: The distribution of the corresponding number N of required iterations for obtaining the Newton-Raphson basins of convergenceshown in Fig. 9(a-f). The non-converging points are shown in red. (Color figure online). . Parametric evolution of the basin entropy As noted in the introduction, the basin entropy con-cept was recently introduced in [7] as a new quantitativemeasure of the uncertainty of a given basin (e.g. of escape,convergence or attraction). The idea behind the methodis to subdivide the phase space into N small cells, each ofwhich containing at least one of the total number of finalstates N A . Since the probability to find a state j in the i − th cell is defined as p i,j , the entropy for the i − th cellcan be calculated by means of the Gibbs entropy as follows S i = N A (cid:88) j =1 p i,j log (cid:18) p i,j (cid:19) . (11)Then, the basin entropy is calculated as the average en-tropy for the total number of cells N , i.e., S b = 1 N N (cid:88) i =1 S i = 1 N N (cid:88) i =1 N A (cid:88) j =1 p i,j log (cid:18) p i,j (cid:19) . (12)Finally, it is worth mentioning that there is a strong de-pendence between the total number of cells N and theresult for basin entropy, such that for larger values of N amore precise value of S b is obtained. This technical prob-lem can be tackled by randomly selecting small cells inphase space through a Monte Carlo procedure (see e.g.,[7, 8]). Following this approach, we find that in our par-ticular problem, the final value for the basin entropy doesnot change for N > × cells, so in all cases, we haveused N = 2 . × cells.In Fig. 11, we present the parametric evolution of thebasin entropy for different values of the perturbation pa-rameter δ , with δ ∈ [0 , δ ∈ [0 , δ the basin entropy decreases almost mono-tonically. Particular attention deserve the critical values δ ≈ . δ = 4, where small variations in the pertur-bation parameter give place to huge changes in the basinentropy. This result can be explained by considering thatclose to δ ≈ . N A is mod-ified and hence the value of Eq. (12) is also significantlymodified. Similarly, close to δ = 4 the number of rootschange from 8 to 10 modifying the value of the basin en-tropy. On the other hand, it is important to note that thesmaller and larger values of the basin entropy correspondto δ approximately 0 and 4, respectively. It means thatthe unpredictability associated to the NR basins of con-vergence for the classical H´enon-Heiles system is smallerin comparison with the Generalized H´enon-Heiles system. Note that the differences for δ = 0 and δ ≈ . δ = 0 and δ ≈ S b , as a function of theperturbation parameter δ . The vertical dashed red lines indicate thecritical values of δ , where the total number of equilibrium pointschanges.
6. Discussion and conclusions
We numerically explored the basins of convergence, re-lated to the equilibrium points, in the generalized H´enon-Heiles system. More precisely, we demonstrated how theperturbation parameter δ influences the position as well asthe linear stability of the libration points. The multivari-ate version of the Newton-Raphson iterative scheme wasused for revealing the corresponding basins of convergenceon the configuration ( x, y ) plane. These attracting do-mains play a significant role, since they explain how eachpoint of the configuration plane is attracted by the libra-tion points of the system, which act, in a way, as attrac-tors. We managed to monitor how the Newton-Raphsonbasins of convergence evolve as a function of the pertur-bation parameter. Another important aspect of this workwas the relation between the basins of convergence and thecorresponding number of required iterations.To our knowledge, this is the first time that the Newton-Raphson basins of convergence, in the generalized H´enon-Heiles system, are numerically investigated in a systematicmanner. On this basis, the presented results are novel andthis is exactly the contribution of our work.The following list contains the most important conclu-sions of our numerical analysis.1. The stability analysis suggests that most of the equi-librium points of the system are always linearly un-stable, when δ ≥
0. Only the central libration point L is Lyapunov stable, for the same values of theperturbation parameter.13. The attracting domains, associated to the librationpoints, extend to infinity, in all studied cases. Al-ways the convergence diagrams, on the configuration( x, y ) plane are symmetrical, with respect to the ver-tical x = 0 axis. In the case of the classical H´enon-Heiles system ( δ = 0) there is an additional 2 π/ δ these initial conditions do converge to one of theattractors but only after a considerable amount of it-erations ( N (cid:29) δ the correspond-ing initial conditions must be true non-convergingpoints, since they do not display any numerical signof convergence, not even after 50000 iterations.4. As expected, the multivariate Newton-Raphson methodwas found to converge very fast (0 ≤ N <
10) forinitial conditions close to the equilibrium point andvery slow ( N ≥
25) for initial conditions of dispersedpoints lying either in the vicinity of the basin bound-aries, or between the dense regions of the librationpoints.5. The lowest value of the basin entropy was foundnear δ = 0, while on the other hand the highestvalue of S b was measured near δ = 4. This im-plies that the unpredictability, regarding the attract-ing regions, in the classical H´enon-Heiles system isconsiderably smaller with respect to the generalizedH´enon-Heiles system.A double precision numerical code, written in standard FORTRAN 77 [18], was used for the classification of the ini-tial conditions into the different basins of convergence. Inaddition, for all the graphical illustration of the paper weused the latest version 11.2 of Mathematica (cid:114) [20]. Usingan Intel (cid:114)
Quad-Core TM i7 2.4 GHz PC the required CPUtime, for the classification of each set of initial conditions,was about 5 minutes.In the future, it would be very interesting to use othertypes of iterative schemes and compare the similarities aswell as the differences on the corresponding basins of at-traction. In particular, using iterative methods of higherorder, with respect to the classical Newton-Raphson method,would be an ideal starting point. This would certainly leadto useful, and perhaps unexpected, results in the very ac-tive field of attracting domains of equilibrium points indynamical systems. Acknowledgments
FLD and ARD acknowledge financial support from Universidadde los Llanos, under Grant No. CDP 2478. FLD gratefully acknowl-edges the financial support provided by COLCIENCIAS, Colombia,under Grant No. 8840. The authors would like to express their warmest thanks to the two anonymous referees for the careful read-ing of the manuscript and for all the apt suggestions and commentswhich allowed us to improve both the quality and the clarity of thepaper.
References [1] Aguirre, J., Vallejo, J.C., Sanju´an, M.A.F., Wada basins andchaotic invariant sets in the H´enon-Heiles system. Phys. Rev. E64 (2001) 066208.[2] Aguirre, J., Viana, R.L., Sanju´an, M.A.F., Fractal Structures innonlinear dynamics. Rev. Mod. Phys. 81 (2009) 333-386.[3] Asique, Md.Ch., Prasad, U., Hassan, M.R., Suraj, Md.S., Onthe photogravitational R4BP when the third primary is a triaxialrigid body. Astrophys. Space Sci. 361 (2016) 379.[4] Baltagiannis, A.N., Papadakis, K.E., Equilibrium points andtheir stability in the restricted four-body problem. Int. J. Bifurc.Chaos 21 (2011) 2179-2193.[5] Croustalloudi, M.N., Kalvouridis, T.J., Attracting domains inring-type N-body formations. Planet. Space Science 55 (2007) 53-69.[6] Croustalloudi, M.N., Kalvouridis, T.J., The Restricted 2+2 bodyproblem: Parametric variation of the equilibrium states of theminor bodies and their attracting regions. ISRN Astronomy andAstrophysics (2013) Article ID 281849.[7] Daza, A., Wagemakers, A., Georgeot, B., Gu´ery-Odelin, D.,Sanju´an, M.A.F., Basin entropy: a new tool to analyze uncer-tainty in dynamical systems. Scientific reports, 6, (2016) 31416.[8] Daza, A., Georgeot, B., Gu´ery-Odelin, D., Wagemakers, A.,Sanju´an, M.A.F., Chaotic dynamics and fractal structures in ex-periments with cold atoms. Physical Review A, 95 (2017) 013629.[9] Douskos, C.N., Collinear equilibrium points of Hill’s problemwith radiation and oblateness and their fractal basins of attrac-tion. Astrophys. Space Sci. 326 (2010) 263-271.[10] Douskos, C.N., Kalantonis, V., Markellos, P., Perdios E., OnSitnikov-like motions generating new kinds of 3D periodic orbitsin the R3BP with prolate primaries. Astrophys. Space Sci. 337(2012) 99-106.[11] Dubeibe, F.L., Ria˜no-Doncel, A., Zotos, E.E., Dynamical analy-sis of bounded and unbounded orbits in a generalized H´enon-Heilessystem. arXiv:1712.01873.[12] Gousidou-Koutita, M., Kalvouridis, T.J., On the efficiency ofNewton and Broyden numerical methods in the investigation ofthe regular polygon problem of ( N + 1) bodies. Appl. Math. Com-put. 212 (2009) 100-112.[13] H´enon, M., Heiles, C., The applicability of the third integral ofmotion: some numerical experiments. Astron. J. 69 (1964) 73-79.[14] Kalantari, B. Polynomiography and applications in art, educa-tion, and science. Computers & Graphics 28 (2004) 417-430.[15] Kalvouridis, T.J., Gousidou-Koutita, M.C. Basins of attractionin the Copenhagen problem where the primaries are magneticdipoles. Applied Mathematics, 3 (2012) 541-548.[16] Kumari, R., Kushvah, B.S., Stability regions of equilibriumpoints in restricted four-body problem with oblateness effects. As-trophys. Space Sci. 349 (2014) 693-704.[17] Nusse, H.E., Yorke, J.A., Basins of attraction. Science 271(1996) 1376-1380.[18] Press, H.P., Teukolsky, S.A, Vetterling, W.T., Flannery, B.P.,Numerical Recipes in FORTRAN 77, 2nd Ed., Cambridge Univ.Press, Cambridge, (1992) USA.[19] Verhulst, F., Discrete symmetric dynamical systems at the mainresonances with applications to axi-symmetric galaxies. Philo-sophical Transactions of the Royal Society of London A: Math-ematical, Physical and Engineering Sciences 290 (1979) 435-465.[20] Wolfram, S., The Mathematica Book, Fifth Edition. WolframMedia, (2003) Champaign.[21] Zotos, E.E., Fractal basins of attraction in the planar circu-lar restricted three-body problem with oblateness and radiationpressure. Astrophys. Space Sci. 361 (2016) 181.
22] Zotos, E.E., Basins of convergence of equilibrium points in thepseudo-Newtonian planar circular restricted three-body problem.Astrophys. Space Sci. 362, (2017) 195.22] Zotos, E.E., Basins of convergence of equilibrium points in thepseudo-Newtonian planar circular restricted three-body problem.Astrophys. Space Sci. 362, (2017) 195.