High-Precision Continuation of Periodic Orbits

Abstract

Obtaining periodic orbits of dynamical systems is the main source of information about how the orbits, in general, are organized. In this paper, we extend classical continuation algorithms in order to be able to obtain families of periodic orbits with high-precision. These periodic orbits can be corrected to get them with arbitrary precision. We illustrate the method with two important classical Hamiltonian problems.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Ángeles Dena 1 and Marcos Rodríguez 1 and Sergio Serrano 2 and Roberto Barrio 2

Recommended by Muhammad Aslam Noor

1, Centro Universitario de la Defensa and IUMA, 50090 Zaragoza, Spain
2, Departamento de Matemática Aplicada and IUMA, Universidad de Zaragoza, 50009 Zaragoza, Spain

Received 30 December 2011; Accepted 14 February 2012

1. Introduction

Obtaining periodic orbits in dynamical systems is the main source of information about how the orbits in general are organized, and, in Poincaré's words, they offer "the only opening through which we might try to penetrate the fortress (chaos) which has the reputation of being impregnable." The stable periodic orbits explain the dynamics of bounded regular motion [1], while the unstable periodic orbits in a chaotic set (attractor or saddle) determine its structure [2] and, thus, the chaotic behavior of the system.

Therefore, it is not surprising the high number of works that use the detection of periodic orbits in the analysis of dynamical systems in many different fields. The analysis and control of chaotic dynamical systems [3, 4], the "scar" theory in quantum mechanics [5], electrical and magnetic fields [6], hydrodynamical flows [7], electrical circuits and optical systems [8], and astrodynamics [1, 9, 10] are just a few examples.

There are several works developing different numerical algorithms to compute periodic orbits [11-19], but most of them do not work with high-precision. This is a serious difficulty if our problem requires obtaining periodic orbits with many precision digits. Some important examples of this kind of problems are the complex pole location in physics [20-23], the study of the splitting of separatrices (numerical difficulties arise because this phenomenon can exhibit exponentially small splitting) [24], and the study of the fractal properties of the chaotic attractors [25]. Currently, the only algorithms capable of computing high-precision periodic orbits are based on Lindstedt-Poincaré series [14] or Taylor series methods [19]. The first one uses very clever algorithms, but it is quite difficult to use in a general setting. Based on the second method, we extend classical continuation algorithms in order to be able to continue families of periodic orbits with high-precision, generating with arbitrary precision the values of those orbits that are considered of special interest. This algorithm is based on the pseudo-arclength continuation method [26, 27], an optimized shooting method with the Newton method [11, 15], the Taylor series method (implemented by means of the free software TIDES [28, 29]), and the singular value decomposition, SVD [30, 31].

The organization of the paper is as follows. In Section 2, we review the correction method developed in [19] and we show an example to illustrate its use with arbitrary precision. Section 3 describes the continuation method explaining some technical aspects and examples. Finally, in Section 4, we present the combined use of methods described above to take advantage of each one of them.

2. Correction of Periodic Orbits with Arbitrary Precision

Let us consider the differential system: [figure omitted; refer to PDF] with y=(_y1 ,_y2 ,...,_ys^)T and _y0 =(^y10 ,^y20 ,...,^ys0)T ∈^...s , and t∈... . And let [figure omitted; refer to PDF] be the solution of (2.1).

If we have an estimation of the initial conditions _y0 of a periodic orbit with estimated period _T0 , then y(_t0 +_T0 ;_y0 )[approximate]_y0 . If we call Δy and ΔT the corrections to the initial conditions and the period, respectively, then we have the periodicity condition: [figure omitted; refer to PDF] To obtain an approximation of ΔT and Δy , we follow the strategy of the Newton method: we take the multivalued Taylor expansion up to the first order in (2.3) to get s linear equations (with s+1 unknowns) to satisfy [figure omitted; refer to PDF] with _Is ∈^...s×s the identity matrix, Δ_T0 and Δ_y0 approximations of ΔT and Δy , respectively, and _M0 ∈^...s×s given by [figure omitted; refer to PDF]

At this point, we have to obtain numerically y(_t0 +_T0 ;_y0 ) and _M0 . Usually, to compute _M0 , the variational equations of the differential system (2.1) are formulated and integrated. However, this task is not always easy and frequently cumbersome. Free software TIDES [28, 29] allows to compute simultaneously the solution of (2.1) and its partial derivatives automatically. Besides, this software is programmed implementing multiple-precision arithmetic (via the MPFR library [32]), which allows us to compute the former integrations with arbitrary precision.

To optimize each iteration of the method, we can add a condition to have a correction (Δ_y0 ,Δ_T0 ) orthogonal to the vector field (2.1): [figure omitted; refer to PDF]

If our differential system (2.1) has an integral of the motion ( IM(t;_y0 )...1;_im , being _im ∈... a constant value) and we want to preserve it in the correction, we have to add an extra equation to the system. With the same idea described above, we take a multivalued Taylor expansion up to the first order: [figure omitted; refer to PDF] Note that, when (2.1) has k independent integrals of the motion, the last equation becomes a set of k equations with _im ∈^...k .

Putting together (2.4), (2.6), and (2.7), we arrive to the system [figure omitted; refer to PDF]

Linear system (2.8) is not necessary squared. Thus, to find a solution that minimizes the residual error, we can use, for example, the Singular Value Decomposition (SVD) method [30, 31], which leads to a robust solver of the minimization problem. After solving the system, we define _y1 =_y0 +Δ_y0 and _T1 =_T0 +Δ_T0 . Now, we can apply the above steps again for a new iteration (with _y1 and _T1 ) that provides us a better approximation and iterate the correction algorithm to correct the initial approximation up to any precision level.

Particularizing now to the case of autonomous Hamiltonian systems with two degrees of freedom ( y=(x,y,X,Y^)T ), the (i+1) th step is defined by the system [figure omitted; refer to PDF] where ...A5; is the symplectic matrix, h the desired value of Hamiltonian constant, and ^yi* =y(_t0 +_Ti ;_yi ) .

Sometimes, we are interested in a subset of the full phase space. Then we have to correct only some variables, fixing the remaining variables and removing their corresponding columns.

As an example, to show the applicability of the method described in this section, we show some results obtained for the classical Hénon-Heiles Hamiltonian [33]: [figure omitted; refer to PDF] In this system, the energy E...1;[Hamiltonian (script capital H)](x,y,X,Y) is a first integral of the problem. This problem was introduced in the study of galactic dynamics to describe the motion of stars around a galactic center, assuming the motion is restricted to the xy plane.

The study of the chaotic and regular regions is traditionally done by means of the Poincaré surfaces of section (PSS) [34] and the Lyapunov exponents. To accelerate the analysis, some fast chaos indicators have been introduced the last few years. The one used in this paper, OFLI2 [35], is defined at the final time _tf by [figure omitted; refer to PDF] where ξ and δξ are the first- and second-order sensitivities with respect to carefully chosen initial vectors and ^{x[perpendicular]} stands for the orthogonal component to the flow of the vector x . In this case, the variational equations up to second order and the initial conditions are [figure omitted; refer to PDF] Note that the last line of (2.12) is written for a single j -th component δ_ξj to simplify the notation. The OFLI2 tends to a constant value for the periodic orbits, behaves linearly (in log-log scale) for regular initial conditions, and grows exponentially for chaotic orbits [35].

Figure 1(a) shows the PSS (dark blue points) and the chaos indicator OFLI2 (different color regions, red color for chaotic behavior and blue color for regular one) for x=0 and E=1/8 . Both techniques give complementary information about the behavior of the system. Black circles sign two periodic orbits selected for the correction process. The left point corresponds to an unstable periodic orbit and the right point to a stable one. On the right side of Figure 1, we show the CPU time in seconds (on a personal computer PC Intel Core i7 CPU 860, 2.80 GHz under 2.6.32-29-generic SMP x86_64 GNU/Linux system) versus the number of digits for obtaining periodic orbits up to one thousand precision digits. We observe that the behavior is quite similar for the stable and the unstable orbits. The computational complexity behaves like ...AA;(^d4 ) being d=-_{log 10} (TOL) the requested number of digits. We can see that the CPU time for 1000 digits of precision is around 3 hours, a very acceptable time for this very high level of precision. The possibility of following orbits with high-precision enables to compute the Poincaré map and its first- and higher-order derivatives also with high-precision.

(a): OFLI2 and PSS section of the Hénon-Heiles problem for x=0 and Energy E=1/8 . (b): Computational relative error versus CPU time for the correction method applied to the two marked periodic orbits on the left picture.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

3. Continuation of Periodic Orbits with High-Precision

If F:^...n+k [arrow right]^...n (with k>0 ) is a smooth mapping and ^x* ∈^...n+k is a solution of F(x)=0 then, by applying the implicit function theorem near ^x* , the solution set of former equation is a smooth k -dimensional submanifold of ^...n+k . Therefore, when k>1 , if we want just a curve of solutions, we have to add further restrictions.

For autonomous Hamiltonian systems, once we have one elementary periodic orbit [36], it can be continued and so the periodic orbits appear in families.

The great advantage of the continuation methods is that, once we have one periodic orbit of a family of periodic orbits, the method gives us the complete family and bifurcations of it. There are several techniques to continue a family of periodic orbits. The simplest one consists of performing a prediction of the next point of the family and, then, a correction to have the new point up to the desired precision level. This technique is used in a great variety of research papers [15, 37] and programs [38-40]. In this work, we use the pseudo-arclength continuation method [26, 27, 41], which allows the continuation of the solution curve regardless of its direction. Therefore, there is no problem with folds or with "vertical" branches. This is an important feature when working in the calculation of periodic orbits of conservative systems [42, 43].

Given a point _xi on the solution curve, the prediction of the next step is a point _x...i+1 =_xi +δ·_ui , where _ui is the unit tangent vector to the curve at _xi and δ is the fixed step size. To compute the unit tangent vector _ui , we have to solve the system _Fx (_xi )_ui =0 , with the restriction ||_ui ||=1 . For _x0 , we have two valid solutions, with the same direction but opposite senses, each of them will give us a side of the curve. For i>0 , we can approximate _ui by choosing a unit secant vector, that is, the normalized vector of the difference _xi -_xi-1 . For practical implementations, this approach leads to a faster and easier to program algorithm.

Once we have _x...i+1 , we can apply the correction method described in the previous section by adding the condition ^uiT ·(^xi+1j -_x...i+1 )=0 , where ^xi+1j is the j th approximation of _xi+1 obtained by the correction algorithm. That is, in the correction process, we impose that _xi+1 is on the hyperplane HP, perpendicular to _ui passing through _x...i+1 (see Figure 2). Here, again (see Section 2), the use of the free software TIDES and SVD allows us to automate the process and to work with any arbitrary precision.

Figure 2: Sketch of one step in pseudo-arclength continuation method.

[figure omitted; refer to PDF]

Note that the phase space where we look for the curve of periodic orbits can be given by variables of the system or by functions ...A2;(y) of them (the Hamiltonian, e.g.). In the latter case, the value of Δ...A2;(y) depends on other corrected variables: [figure omitted; refer to PDF]

Following the example of Section 2, we are going to look for families of periodic orbits of the classical Hénon-Heiles Hamiltonian with _x0 =_Y0 =0 and _X0 given by _y0 and [Hamiltonian (script capital H)] . In this case, we have three variables to continue {_y0 ,[Hamiltonian (script capital H)],T} (we also have to predict next period T ) and [figure omitted; refer to PDF] In Figure 3, we show two curves corresponding to two different families of periodic orbits. The blue curve corresponds to orbits of multiplicity m=1 , while the red curve corresponds to orbits of multiplicity m=5 . The stability of the orbits [18, 44] is discriminated on the figure by using continuous or dashed lines for stable and unstable periodic orbits, respectively. Both curves have been obtained with δ=0.01 and 15 digits of precision. In Table 1, we show CPU times (in seconds) for the computation of the unstable branch of the family with multiplicity 5 (m5u) in Figure 3 for different values of δ and precision level. These CPU times are still acceptable, although obviously we cannot take values of delta of the size of the precision considered in the correction since the number of points calculated on the curve would soar.

Table 1: CPU times for the computation of the unstable branch of the family with multiplicity 5 (m5u) in Figure 3.

	CPU time (seconds)
Precision	δ=0.01	δ=0.001
Double	7.002675	62.871128
Multiple (16 digits)	131.818665	1187.724976
Multiple (32 digits)	304.062012	2980.975586
Multiple (100 digits)	2741.602051	28029.414062

Figure 3: Continuation of two families of periodic orbits (multiplicities m=1 and 5) of the classical Hénon-Heiles Hamiltonian with _x0 =_Y0 =0 depending on the coordinate y and the energy E .

[figure omitted; refer to PDF]

Summarizing the method, we start from a point in the family of periodic orbits (we can obtain it by correcting an approximation with the algorithm of Section 2). Now, if we want to compute the new initial conditions (_ynew ,_{[Hamiltonian (script capital H)]new} ,_Tnew ) at a distance δ from point (_yold ,_{[Hamiltonian (script capital H)]old} ,_Told ) , we define the first approximation, (y...,[Hamiltonian (script capital H)]...,T...)=(_yold ,_{[Hamiltonian (script capital H)]old} ,_Told )+δ·u , of the new point as the starting point. The correction of the provisional conditions (y...+Δy,[Hamiltonian (script capital H)]...+Δ[Hamiltonian (script capital H)],T...+ΔT) will be in the previously defined hyperplane HP: [figure omitted; refer to PDF] Then,

(i) we compute corrections in position Δy , momentum ΔX and period ΔT . The corrected orbit will satisfy

(a) correction conditions (2.9) with Δx=ΔY=0 ,

(b) restriction of being in the hyperplane HP (3.3);

(ii) we output Δ[Hamiltonian (script capital H)] (according to (3.2)) from Δy and ΔX ;

(iii): we iterate previous stages until the desired precision level.

As a second example, we consider the Copenhagen problem, a simplification of the three-body problem supposing that the mass of one of the three bodies is negligible and the other two masses are equal [10, 44]. The equations of motion of this problem are [figure omitted; refer to PDF] with ^r12 =(x+1/2⁾² +^y2 and ^r22 =(x-1/2⁾² +^y2 . The only known integral of motion of this problem is the Jacobian constant C [figure omitted; refer to PDF] that can be used to define the effective energy as _EJ =-C/2 . Note that now we do not use the Hamiltonian formalism and so the function ...A2; (3.1) is the Jacobi constant.

In [44], Hénon made an exhaustive analysis of classes of periodic orbits of the problem. Several families of such orbits were drawn. Here, we use the necessary precision to redraw one of them and to check that the end of the curve is a spiral (see Figure 4). This situation cannot be seen from the pictures in [44] and provides a clear application of the present approach. Note that with this approach we may obtain the curves with any desired precision level, up to the hardware and CPU time restrictions, whereas the classical continuation programs, like AUTO or MATCONT, are limited to double-precision (enough for most of the situations).

Figure 4: Family of periodic orbits for the Copenhagen problem depending on the coordinate x and the Jacobi constant C . Picture (a) is extracted from Exploration numérique du problème restreint. I. Masses égales; orbites périodiques , Hénon [44]; Picture (b) shows the same curve obtained with the method described above. Bottom: magnifications of the squared regions showing the spiral structure.

[figure omitted; refer to PDF]

4. Continuation of Periodic Orbits with Arbitrary Precision

As stated in the previous section, the CPU time for the families of periodic orbits with high-precision is acceptable as long as we do not take a too small value of δ or a too high-precision level. If you also want to make an exhaustive continuation of the different families of periodic orbits that can appear in a region of the phase space, restrictions must be more severe. However, in most cases, these restrictions are acceptable because we do not need to get the data of all the orbits of the different families with arbitrary precision. Double- or quadruple-precision is usually sufficient for most of them. Only a few of these orbits (being of special interest) have to be further corrected with arbitrary precision.

This idea has been implemented in the calculations of Figure 5 where the two families of periodic orbits have been computed in double-precision. Subsequently, the orbits of the family of multiplicity 5, corresponding to values of E=0.12 and 0.14, or y=0.5 , (values taken just as an example) have been corrected to obtain 100 digits of precision. On the right plots, the obtained periodic orbits are shown: unstable orbits above and below the stable ones. In Table 2, we show the corrected values up to 50 digits for these selected orbits.

Table 2: Corrected values up to 50 digits for the selected orbits in Figure 5.

y0	1.8817343638528059241776135906704656237302870299392 e-3
[Hamiltonian (script capital H)]	1.2 e-1
T	3.1884755881458007139894029240879646745967964214245 e1
y0	5.4846395794245312569019288801973048034848306007155 e-1
[Hamiltonian (script capital H)]	1.2 e-1
T	3.2132619977067518675898102052374230853344441909611 e1
y0	-6.9222934138872141467629060227049920239548916485414 e-2
[Hamiltonian (script capital H)]	1.4 e-1
T	3.2463680584143904506022667600424054135989589752729 e1
y0	6.4236679697082661947506127379024223882593232233963 e-1
[Hamiltonian (script capital H)]	1.4 e-1
T	3.3190847175337893134522290150974531243339047630215 e1
y0	5.0 e-1
[Hamiltonian (script capital H)]	1.1064108429643783371307847713912555690466002998356 e-1
T	3.1704681349442488934633816972756993961160952189802 e1

(a) continuation of two families of periodic orbits of the classical Hénon-Heiles Hamiltonian with _x0 =_Y0 =0 . Orbits with some selected values of the energy and coordinate y are pointed out. (b) corrected arbitrary precision periodic orbits corresponding to surrounded points following the color code; top, unstable orbits, and bottom, the stable ones.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

5. Conclusions

In this paper, we introduce a numerical scheme able to continue families of periodic orbits in high-precision and to correct selected orbits up to any arbitrary precision level. The method is based on the combined use of the pseudo-arclength continuation method, an optimized shooting method with the Newton method, the Taylor series method (implemented by means of the free software TIDES), and the singular value decomposition (SVD). We illustrate the scheme on the paradigmatic Hénon-Heiles Hamiltonian and the Copenhagen problem. On the Copenhagen problem, we have been able to show that a classical family of periodic orbits found by Hénon ends in a spiral structure, which cannot be observed without high-precision.

Acknowledgments

The authors have been supported for this research by the Spanish Research Grant MTM2009-10767. A. Dena thanks Professor G. Lord (Heriot-Watt University, Edinburgh) for his hospitality and his help during her stay at his department.

References

[1] R. Barrio, F. Blesa, S. Serrano, "Bifurcations and safe regions in open Hamiltonians," New Journal of Physics , vol. 11, 2009.

[2] R. Gilmore, M. Lefranc The Topology of Chaos , pp. xxiv+495, John Wiley & Sons, New York, NY, USA, 2002.

[3] K. Pyragas, "Control of chaos via an unstable delayed feedback controller," Physical Review Letters , vol. 86, no. 11, pp. 2265-2268, 2001.

[4] W. M. Zheng, "Predicting orbits of the Lorenz equation from symbolic dynamics," Physica D , vol. 109, no. 1-2, pp. 191-198, 1997.

[5] D. A. Wisniacki, E. Vergini, R. M. Benito, F. Borondo, "Signatures of homoclinic motion in quantum chaos," Physical Review Letters , vol. 94, no. 5, 2005.

[6] A. D. Peters, C. Jaffé, J. B. Delos, "Closed-orbit theory and the photodetachment cross section of H- in parallel electric and magnetic fields," Physical Review A , vol. 56, no. 1, pp. 331-344, 1997.

[7] E. Kazantsev, "Sensitivity of the attractor of the barotropic ocean model to external influences: approach by unstable periodic orbits," Nonlinear Processes in Geophysics , vol. 8, no. 4-5, pp. 281-300, 2001., [email protected]

[8] J. Aguirre, R. L. Viana, M. A. F. Sanjuán, "Fractal structures in nonlinear dynamics," Reviews of Modern Physics , vol. 81, no. 1, pp. 333-386, 2009.

[9] R. P. Russell, "Global search for planar and three-dimensional periodic orbits near Europa," Journal of the Astronautical Sciences , vol. 54, no. 2, pp. 199-226, 2006.

[10] R. Barrio, F. Blesa, S. Serrano, "Periodic, escape and chaotic orbits in the Copenhagen and the (n + 1)-body ring problems," Communications in Nonlinear Science and Numerical Simulation , vol. 14, no. 5, pp. 2229-2238, 2009.

[11] S. C. Farantos, "Methods for locating periodic orbits in highly unstable systems," Journal of Molecular Structure , vol. 341, no. 1-3, pp. 91-100, 1995.

[12] P. Schmelcher, F. K. Diakonos, "Detecting unstable periodic orbits of chaotic dynamical systems," Physical Review Letters , vol. 78, no. 25, pp. 4733-4736, 1997.

[13] R. L. Davidchack, Y. C. Lai, "Efficient algorithm for detecting unstable periodic orbits in chaotic systems," Physical Review E , vol. 60, no. 5 B, pp. 6172-6175, 1999.

[14] D. Viswanath, "The Lindstedt-Poincaré technique as an algorithm for computing periodic orbits," SIAM Review , vol. 43, no. 3, pp. 478-495, 2001.

[15] M. Lara, J. Peláez, "On the numerical continuation of periodic orbits: an intrinsic, 3-dimensional, differential, predictor-corrector algorithm," Astronomy and Astrophysics , vol. 389, no. 2, pp. 692-701, 2002.

[16] J. H. B. Deane, L. Marsh, "Unstable periodic orbit detection for ODEs with periodic forcing," Physics Letters, Section A , vol. 359, no. 6, pp. 555-558, 2006.

[17] Y. Saiki, "Numerical detection of unstable periodic orbits in continuous-time dynamical systems with chaotic behaviors," Nonlinear Processes in Geophysics , vol. 14, no. 5, pp. 615-620, 2007., [email protected]

[18] R. Barrio, F. Blesa, "Systematic search of symmetric periodic orbits in 2DOF Hamiltonian systems," Chaos, Solitons and Fractals , vol. 41, no. 2, pp. 560-582, 2009.

[19] A. Abad, R. Barrio, A. Dena, "Computing periodic orbits with arbitrary precision," Physical Review E , vol. 84, no. 1, Article ID 016701, 2011., [email protected]; [email protected]; [email protected]

[20] J. N. L. Connor, D. Farrelly, "Uniform semiclassical and quantum calculations of Regge pole positions and residues for complex optical nuclear heavy-ion potentials," Physical Review C , vol. 48, no. 5, pp. 2419-2432, 1993.

[21] D. Sokolovski, "Complex-angular-momentum analysis of atom-diatom angular scattering: zeros and poles of the S matrix," Physical Review A , vol. 62, no. 2, Article ID 024702, 2000.

[22] G. L. Alfimov, D. Usero, L. Vázquez, "On complex singularities of solutions of the equation H_ux -u+^up =0 ," Journal of Physics A , vol. 33, no. 38, pp. 6707-6720, 2000.

[23] D. Viswanath, S. Sahuto[gcaron]lu, "Complex singularities and the lorenz attractor," SIAM Review , vol. 52, no. 2, pp. 294-314, 2010.

[24] V. Gelfreich, C. Simó, "High-precision computations of divergent asymptotic series and homoclinic phenomena," Discrete and Continuous Dynamical Systems - Series B , vol. 10, no. 2-3, pp. 511-536, 2008.

[25] D. Viswanath, "The fractal property of the Lorenz attractor," Physica D , vol. 190, no. 1-2, pp. 115-128, 2004.

[26] H. B. Keller Numerical Solution of Bifurcation and Nonlinear Eigenvalue Problems , Academic Press, 1977.

[27] J. Broeckhove, P. Klosiewicz, W. Vanroose, "Applying numerical continuation to the parameter dependence of solutions of the Schrödinger equation," Journal of Computational and Applied Mathematics , vol. 234, no. 4, pp. 1238-1248, 2010.

[28] A. Abad, R. Barrio, F. Blesa, M. Rodríguez, "TIDES Software," http://gme.unizar.es/software/tides

[29] A. Abad, R. Barrio, F. Blesa, M. Rodríguez, "TIDES: a Taylor series Integrator for Differential Equations," ACM Transactions on Mathematical Software . In press

[30] J. W. Demmel Applied Numerical Linear Algebra , pp. xii+419, SIAM, Philadelphia, Pa, USA, 1997.

[31] L. N. Trefethen, D. Bau III Numerical Linear Algebra , pp. xii+361, SIAM, Philadelphia, Pa, USA, 1997.

[32] L. Fousse, G. Hanrot, V. Lefèvre, P. Pélissier, P. Zimmermann, "MPFR: a multiple-precision binary floating-point library with correct rounding," ACM Transactions on Mathematical Software , vol. 33, no. 2, article 13, pp. 15, 2007.

[33] M. Hénon, C. Heiles, "The applicability of the third integral of motion: some numerical experiments," The Astronomical Journal , vol. 69, pp. 73-79, 1964.

[34] H. Poincaré Les Méthodes Nouvelles de la Mécanique Céleste , Dover, New York, NY, USA, 1957.

[35] R. Barrio, "Sensitivity tools vs. Poincaré sections," Chaos, Solitons and Fractals , vol. 25, no. 3, pp. 711-726, 2005.

[36] K. R. Meyer, "Generic bifurcation of periodic points," Transactions of the American Mathematical Society , vol. 149, pp. 95-107, 1970.

[37] C. Wulff, A. Schebesch, "Numerical continuation of symmetric periodic orbits," SIAM Journal on Applied Dynamical Systems , vol. 5, no. 3, pp. 435-475, 2006.

[38] E. Doedel, "AUTO: a program for the automatic bifurcation analysis of autonomous systems," in Proceedings of the 10th Manitoba Conference on Numerical Mathematics and Computing, Vol. I (Winnipeg, Man., 1980) Congressus Numerantium , vol. 30, pp. 265-284, 1981.

[39] E. J. Doedel, R. C. Paffenroth, A. R. Champneys, T. F. Fairgrieve, B. Sandstede, X. Wang, "AUTO 2000: continuation and bifurcation software for ordinary differential equations (with HomCont),", California Institute of Technology

[40] A. Dhooge, W. Govaerts, Y. A. Kuznetsov, "MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs," ACM Transactions on Mathematical Software , vol. 29, no. 2, pp. 141-164, 2003.

[41] B. Krauskopf, H. M. Osinga, J. Galán-Vioque Numerical Continuation Methods for Dynamical Systems , Springer, Dordrecht, The Netherlands, 2007.

[42] F. J. Muñoz-Almaraz, E. Freire, J. Galán, E. Doedel, A. Vanderbauwhede, "Continuation of periodic orbits in conservative and Hamiltonian systems," Physica D , vol. 181, no. 1-2, pp. 1-38, 2003.

[43] E. J. Doedel, V. A. Romanov, R. C. Paffenroth, H. B. Keller, D. J. Dichmann, J. Galán-Vioque, A. Vanderbauwhede, "Elemental periodic orbits associated with the libration points in the circular restricted 3-body problem," International Journal of Bifurcation and Chaos in Applied Sciences and Engineering , vol. 17, no. 8, pp. 2625-2677, 2007.

[44] M. Hénon, "Exploration numérique du problème restreint. I. Masses égales ; orbites périodiques," Annales d'Astrophysique , vol. 28, pp. 499, 1965.