A&A 389, 692-701 (2002)
DOI: 10.1051/0004-6361:20020598
M. Lara^{1} - J. Peláez^{2}
1 - Real Instituto y Observatorio de la Armada,
11110 San Fernando, Spain
2 -
Escuela Técnica Superior de Ingenieros Aeronáuticos.
Universidad Politécnica,
28040 Madrid, Spain
Received 17 January 2000 / Accepted 25 March 2002
Abstract
This paper deals with the analytic continuation of periodic orbits
of conservative dynamical systems with three degrees of freedom. For
variations of any parameter (or integral), it relies on
numerical analysis in order to implement a predictor-corrector
algorithm to compute the initial conditions of the periodic
orbits pertaining to the family. The method proposed here is not
restricted to symmetric problems and, since the procedure involves
the computation of the variational equations, a side effect is the
trivial computation of the linear stability of the periodic orbits.
As an illustration of the robustness of the method, several families of
periodic orbits of the Restricted Three-Body Problem are computed.
Key words: methods: numerical - celestial mechanics
In this paper we propose an algorithm for the numerical computation of families of periodic orbits of conservative dynamical systems with three degrees of freedom. Initial conditions for specific periodic orbits pertaining to a family, which is defined by the variation of a parameter, are obtained by means of an intrinsic, three dimensional, differential, predictor-corrector algorithm.
The use of differential correction algorithms for the numerical computation of either two or three dimensional periodic orbits is not a new result. The contributions (Deprit & Henrard 1967; Henrard & Lemaitre 1986; Caranicolas 1994) in reference to systems with two degrees of freedom, or (Howell 1984; Karimov & Sokolsky 1989; Belbruno et al. 1994; Contopoulos & Barbanis 1994; Scheeres 1999) with respect to systems with three degrees of freedom, can be mentioned among many others.
Often the dynamical system admits some kind of symmetry, and consequently the traditional approaches used for computing periodic orbits were based on those symmetries. When dealing with force fields without symmetries, a different approach must be used. The normal procedure is then the application of the Poincaré map, and differential corrections are obtained through the computation of the state transition matrix along the periodic orbit.
For conservative systems the monodromy matrix has one unit eigenvalue with multiplicity two - related to the time invariance of the system - thus preventing the computation of the corrections. Two approaches are normally used to compute the nontrivial eigenvalues, both of them based on the integration of the variations in Cartesian coordinates. The first computes the complete state transition matrix and uses basic techniques of matrix algebra for obtaining the eigenvalues of a singular matrix. The second eliminates the two unit eigenvalues from the system by creating a lower dimensional map.
The alternative presented here is based on the computation of the variational equations when projected onto the Frenet frame. While Cartesian variations merge periodic and secular terms, the intrinsic formulation of the variations results in the separation of the tangent displacements (Deprit 1981), which retain the secular part of the variations and cause the varied periodic orbit to have a new period. Therefore, the formulation of the variational equations in the Frenet frame directly reduces the dimension of the state transition matrix to compute, and eliminates the trivial eigenvalues. Of course, our algorithm is not restricted to symmetric problems, and is valid for the computation of families of periodic orbits for variations of any parameter or integral for a conservative dynamical system with three degrees of freedom.
A vectorial predictor-corrector algorithm for the computation of the analytic continuation of a periodic orbit is given in Sect. 2 and the formulation of the variational equations in the Frenet frame (Deprit 1981) is recalled in Sect. 3, where a re-formulation of the predictor-corrector procedure is introduced - the major contribution of this paper. Since the algorithm proposed here involves the computation of the variations, the calculation of the linear stability of the periodic orbits requires no additional effort: the nontrivial characteristic exponents are computed from the normal and binormal variations. The linear stability of the periodic orbits (Sect. 4) can then be studied by using the two stability indices first proposed by Broucke (1969).
The algorithm has been checked recalculating some of the halo orbits of the Restricted Three-Body Problem computed in Gómez et al. (1985). Although that problem is a symmetric one, the procedure described here does not make use of this fact. All the tests done showed very good agreement with the results presented in that work for either the period or any of the stability indices. In order to illustrate the robustness of the method, an application to the Restricted Three-Body Problem is provided in Sect. 5 where we compute families of three dimensional periodic orbits that bifurcate from Strömgren's (1935) planar family m of (retrograde) periodic orbits.
Let
As a consequence of the uniqueness theorem for differential equations, Eq. (3) will hold for all t as it does for a value t_{0}, say
t_{0}=0. Therefore, the new solution shall verify the periodicity condition
Rewriting Eq. (4) as
Equation (6) provides the basic scheme for implementing the
Poincaré's continuation method. The right hand side of Eq. (6)
vanishes for a periodic solution - Eq. (3) - of Eq.
(1), and the linear system
The partial derivatives of x with respect to the initial conditions
are computed from the homogeneous variational system
Finally, we note that the previous procedure is directly applicable to
nonautonomous systems
In what follows we restrict the study to dynamic systems of three degrees of freedom and change the notation for convenience. Hereafter, vectors are three dimensional.
We limit our analysis to dynamical systems determined by the Lagrangian
One passes from Cartesian to intrinsic coordinates by means of a
rotation
of the reference frame. Therefore, an intrinsic variation
(28) |
(29) |
(30) |
Q_{1}= |
Q_{2}=-R_{2} |
For details on the derivation of the Eqs. (34)-(35) the
reader is referred to Deprit (1981) where his
has been
replaced by the partial derivative
and the terms Q_{5} and R_{5} have been added,
respectively, to the right hand sides of the Eqs. (20_{1}) and (20_{2}) of his paper; these extra terms come from the right hand
side of the nonhomogeneous Eqs. (25). Note that the radius
of curvature
and the radius of torsion
that appear in the
formulas given by Deprit (1981) can be obtained from
(36) |
Lets now re-formulate the algorithm given in Sect. 2 but now
using intrinsic coordinates. First we note that a variational displacement of
Eq. (2), now written as
We formulate the problem in the Frenet frame by multiplying Eq. (39)
by the rotation matrix Eq. (31) that must be evaluated in t=T_{0}. We
get
For the predictor stage - being exactly periodic the initial orbit given by Eq.
(37) - the right hand side of Eq. (41) vanishes and
.
Therefore
Once
is obtained the corrections the initial conditions are easily
obtained from
Finally, in order to maintain the right direction of the corrections without
leaving the energy manifold, the new initial conditions are computed as in
Deprit & Henrard (1967)
The solution Eq. (37) is only approximately periodic. Thus, the
right hand side of (41) no longer vanishes. Assuming that the
difference
is of the same order as the
desired corrections we neglect
.
Then Eq. (41) is written
(54) |
Finally, Eq. (51) provides the desired corrections.
We choose initial conditions
Linear stability of periodic orbits depends on the eigenvalues of the resolvent
of the variational equations associated with the fundamental period of the
periodic orbit. In system (21) admitting Hamiltonian formulation the
eigenvalues appear in reciprocal pairs
,
(i=1,2,3)and one eigenvalue takes the value 1 with multiplicity 2. The trivial
eigenvalues correspond to the tangent displacements. Thus, the nontrivial
eigenvalues are obtained from Eq. (47) by solving the equation
Taking in account that the eigenvalues appear in reciprocal pairs, Eq. (60) is
(61) |
(62) |
(63) |
Finally, in reference to planar solutions, one should note that the
two stability indices k_{1} and k_{2} correspond to the in plane
stability and to the "vertical'' or out of plane stability
(Hénon 1973). The stability in the normal (or in plane)
direction is directly evaluated by
As illustration, the procedure described above is applied to one of the most widely studied problems in celestial mechanics: the (circular) Restricted Three-Body Problem. At this point we want to make clear that the aim of this section is not to review the problem but just to show the robustness of the intrinsic predictor-corrector algorithm by computing several families of periodic orbits.
Given a particle of infinitesimal mass that is attracted by two
primaries each moving around the other, in a synodic frame with the
x-axis is in the direction of the primaries, the z-axis in the
normal direction to the orbital plane of the primaries and the
y-axis completing a direct frame, the equations of motion are
There is an integral of these equations, the Jacobian constant, that
we write
(70) |
(71) |
For both primaries with equal masses () let us vary the Jacobi constant ( and both Q_{5} and R_{5} vanish).
We begin with a Keplerian approximation of a circular orbit in the
x-y plane far away from the primaries. For a semimajor axis
a=4, the Keplerian period in the synodic frame is
for a retrograde orbit. The corresponding initial conditions of this
very rough approximation of a periodic orbit are
,
y_{0}=4 and
.
The
maximum difference in absolute value between each of these initial
conditions and the same coordinate or velocity at the approximate
period T is
,
but three
iterations are enough for the corrector to obtain the sequence
,
and finally the improved periodic orbit with
.
The corrected initial conditions and
period for
are
Figure 1: Normal and binormal stability indices versus the Jacobi constant h for the family m of periodic orbits. The horizontal dashed lines correspond to the critical values 2. | |
Open with DEXTER |
In that way we compute the planar family m - in Strömgren's (1935) notation - of periodic orbits for variations of the Jacobi constant. The value cannot be fixed constant for all the computations and for orbits closer to the primaries, smaller values are required to obtain predictions amenable to improvement by the corrector. Figure 1 shows the behavior of the normal and binormal stability indices of this planar family; due to the large values of both indices, a scale proportional to is used. While for - orbit m_{1} in Hénon (1965) - the binormal index takes the critical values -2 for the bifurcation value and +2 for - respectively orbits m_{1v} and m_{2v} of Hénon (1973) -. The initial conditions below ( ) correspond to the critical orbits.
In the case
a bifurcated branch family of
three-dimensional periodic orbits starts from the planar,
retrograde orbit m_{1v} with period doubling and ends at the
equilateral equilibrium L_{4}. We call this family the
L_{4}-family and of course the symmetrical L_{5}-family also exists.
We remark here that we do not claim to establish a way for
computing bifurcated orbits - see, for instance
(Hénon 1973; Markellos 1981; Belbruno et
al. 1994) -, but simply to show the robustness of the
predictor-corrector algorithm. With this in mind, we compute the
bifurcated family as follows. Close to the bifurcation value, we take
a periodic orbit of the planar family
m for
h | 0.595060530821579 | -1.234038469178420 |
y | 1.174109021245709 | 0.850091106717253 |
1.972301154871868 | 0.111468952994908 | |
0.495891121768287 | 0.519684176787866 | |
T | 7.063185569369677 | 6.363331100559023 |
Figure 2: Stability indices versus the Jacobi constant for the L_{4}-family of periodic orbits. For h<-1.234039, k_{1} and k_{2} are complex numbers. | |
Open with DEXTER |
Figure 3: Periodic orbits of the L_{4}-family when k_{1}=k_{2}. Left , right . Primaries are represented by dots. | |
Open with DEXTER |
Figure 4: Stability indices k_{1} and k_{2} versus the Jacobi constant h for the L_{1}-family of periodic orbits. Both indices are complex numbers for -0.47480835<h<0.46334917. | |
Open with DEXTER |
Figure 5: Stability behavior of the L_{1}-family on the a_{1}-a_{2} plane. A scale proportional to is used. | |
Open with DEXTER |
In the case we found a transversal bifurcation. We call the bifurcated family the L_{1}-family because it ends at the L_{1} collinear point. As before, we modify the initial velocity of a planar orbit just passed the bifurcation value and use these initial conditions as starters for the corrector procedure. For h_{0}=0.4952105308215789 the initial conditions of a periodic orbit of the planar family are
Again we find two possibilities associated with the sign () of . The stability k_{1} and k_{2} indices of the L_{1}-family are presented in Fig. 4 versus the Jacobi constant h. For -0.47480835<h<0.46334917 this family has complex instability. The behavior on the a_{1}-a_{2} plane is presented in Fig. 5, where we appreciate that the family passes from Broucke's (1969) region 1 (stability) to region 2 (complex instability) through a Krein collision (Howard & MacKay 1987), then to region 4 (even-even instability) and finally to region 6 (even semi-instability). For convenience the values of a_{1} and a_{2} have been scaled proportional to their respective . Figure 6 shows the two periodic orbits of this family when k_{1}=k_{2} and the respective initial conditions ( ) are listed below.
h | -0.474808349178421 | 0.463349173821579 |
y | 0.630644393432460 | 0.772133837453028 |
1.116409595398506 | 1.913885053382574 | |
0.828732085352047 | 0.184702299430020 | |
T | 2.722825258179725 | 2.474956115380520 |
Figure 6: Periodic orbits of the L_{1}-family when k_{1}=k_{2}. Left , right . Primaries are represented by dots. | |
Open with DEXTER |
Figure 7: Stability behavior of the family for variations of . A scale proportional to is used. | |
Open with DEXTER |
We vary now the real parameter .
The partial derivative
(73) |
x= | |
y=0.8500912035204521, | |
=0.1114681281729226, | |
=0.5196824227332966, |
In this case the value can be fixed for all the family and, for instance, for a value the typical sequence of iterations of the predictor-corrector procedure is , and . Greater values can be chosen and even for the algorithm shows good convergence.
For a wide range of dynamical systems, tangent displacements are separated from normal and binormal ones when using intrinsic variations.
Differential corrections algorithms for the continuation of periodic orbits rely on the integration of the variational equations. When dealing with time-invariant systems one confronts the problem of eliminating the trivial eigenvalues of the monodromy matrix. This inconvenience does not exist when using the Frenet frame, that has the important advantage of strictly adhering to the geometry of the problem.
The analytical continuation of periodic orbits for variations of any parameter or integral of a conservative dynamical system with three degrees of freedom, has been implemented by means of a predictor-corrector algorithm based on the integration of the variational equations when projected onto the Frenet frame. The algorithm proposed here is not restricted to symmetric problems and, since the computation of the intrinsic variations is part of the procedure, the orbital stability in linear approximation can be obtained without additional calculations. Application to the computation of some families of 3-dimensional periodic orbits of the Restricted Three-Body problem that bifurcate from the Strömgren m-family show the reliability and robustness of the algorithm.
Acknowledgements
Thanks are due to Dr. Llibre for providing initial conditions of halo orbits in order to test an early stage of the software. Illustrative discussions with Drs. Deprit and Elipe took place at a 1996 seminar at the University of Zaragoza. Discussions with Dr. Scheeres during a visit of the first author to the University of Michigan at Ann Arbor in May 2001 were very fruitful. We thank him also for reading the paper and for his many suggestions for improving the English. An anonymous referee pointed to the practical application of the method (Sect. 3.3.3), avoiding unnecessary computational effort. M. Lara acknowledges partial support from the Centre Nationale d'Etudes Spatiales at Toulouse (France).