A&A 389, 692-701 (2002)
DOI: 10.1051/0004-6361:20020598

On the numerical continuation of periodic orbits

An intrinsic, 3-dimensional, differential, predictor-corrector algorithm

M. Lara¹ - J. Peláez²

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

1 Introduction

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.

   
2 On Poincaré's continuation method: An algorithm

Let

$$\dot{\vec{x}}=\vec{F}(\vec{x};\sigma)$$ (1)

be an autonomous differential system of n equations depending on a parameter $\sigma$. The solution

$$\vec{x}=\vec{x}(t;\vec{\xi};\sigma)$$ (2)

of Eq. (1) is a function of the parameter $\sigma$ and the initial conditions $\vec{\xi}=\vec{x}(0;\vec{\xi};\sigma)$. Let us suppose that a periodic solution of system (1)

$$\vec{x}(t;\vec{\xi}_0;\sigma_0) =\vec{x}(t+T_0;\vec{\xi}_0;\sigma_0)$$ (3)

with period T₀>0 is known for a certain value $\sigma=\sigma_0$ of the parameter and for the initial conditions $\vec{\xi}=\vec{\xi}_0$. The Poincaré continuation method is concerned with the problem of computing the analytic continuation of Eq. (3) for values of the parameter close to the initial value $\sigma_0$. That is, for $\sigma=\sigma_0+\Delta\sigma$ new initial conditions $\vec\xi=\vec\xi_0+\Delta\vec\xi$ and period $T=T_0+\Delta{T}$ are required to produce a new periodic solution of Eq. (1).

As a consequence of the uniqueness theorem for differential equations, Eq. (3) will hold for all t as it does for a value t₀, say t₀=0. Therefore, the new solution shall verify the periodicity condition

$$\vec{x}(T;\vec{\xi};\sigma)-\vec{\xi}={\bf0},$$ (4)

where the existence of the implicit functions $\vec{\xi}$ is directly related to the nonvanishing of the Jacobian determinant of the left hand side of Eq. (4). In reference to this Jacobian there can be different possibilities depending on the existence of integrals of Eq. (1). The interested reader should consult Siegel & Moser (1971).

2.1 The predictor-corrector algorithm

Rewriting Eq. (4) as

$$\vec{x}(T_0+\Delta{T};\vec{\xi}_0+\Delta\vec\xi;\sigma_0+\Delta\sigma) -(\vec{\xi}_0+\Delta\vec\xi)={\bf0}$$ (5)

and expanding it around $(T_0;\vec{\xi}_0;\sigma_0)$, at first order one obtains the linear system

$$(\nabla_{\xi}\vec{x}-\vec{I})\cdot\Delta\vec\xi +\vec{F}(\v... ...l\vec{x}}{\partial\sigma}~\Delta\sigma= -(\vec{x}-\vec\xi_0),$$ (6)

that must be evaluated at $(T_0;\vec{\xi}_0;\sigma_0)$, where $\vec{I}$ is the $n\times{n}$ identity matrix.

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

$$(\nabla_{\xi}\vec{x}-\vec{I})\cdot \frac{\Delta\vec\xi}{\D... ...ta{T}}{\Delta\sigma} =-\frac{\partial\vec{x}}{\partial\sigma}$$ (7)

evaluated at $(T_0;\vec{\xi}_0;\sigma_0)$, can be used to compute the initial conditions $\vec{\xi}_1=\vec{\xi}_0+\Delta\vec{\xi}$ and period $T_1=T_0+\Delta~T$ of a new periodic solution of Eq. (1) for the value $\sigma_1=\sigma_0+\Delta\sigma$ of the parameter. Since the solution of Eq. (7) is a tangent prediction it could happen that $\vec{x}(T_1;\vec{\xi}_1;\sigma_1)-\vec\xi_1\ne{\bf0}$; in that case new corrections $\Delta\vec{\xi}_1$ and $\Delta~T_1$ must be computed such that

$$\vec{x}(T_1+\Delta{T}_1;\vec{\xi}_1+\Delta\vec\xi_1;\sigma_1) -(\vec{\xi}_1+\Delta\vec\xi_1)={\bf0}.$$ (8)

Equation (8) is formally equivalent to Eq. (5) - now with $\Delta\sigma\equiv0$ - so we can again use Eq. (6) to obtain the corrections. Then, iterative corrections could be in order from

$$(\nabla_{\xi}\vec{x}-\vec{I}) \cdot\Delta\vec\xi_i+\vec{F}(\vec{x};\sigma_1)~\Delta{T}_i =-(\vec{x}-\vec\xi_i)$$ (9)

that must be evaluated at $(T_i;\vec{\xi}_i;\sigma_1)$. Note that now the right hand side no longer vanishes.

The partial derivatives of x with respect to the initial conditions are computed from the homogeneous variational system

$$\nabla_{\xi}\vec{\dot{x}} =(\nabla_{x}\vec{F})\cdot (\nabla_{\xi}\vec{x}),$$ (10)

starting from the initial conditions

$$\Big(\nabla_{\xi}\vec{x}\Big)_{t=0}=\vec{I},$$ (11)

as follows from $\vec{x}(0;\vec\xi;\sigma)=\vec\xi$. Analogously, the partial derivative $\partial\vec{x}/\partial\sigma$ is a particular solution of the variational equations

$$\frac{\partial\vec{\dot{x}}}{\partial\sigma} =(\nabla_{x}\ve... ...\partial\sigma} +\frac{\partial\vec{F}}{\partial\sigma}\cdot$$ (12)

Equations (7) and (9) provide a predictor-corrector scheme for computing the analytic continuation of Eq. (3), but one should note that both equations are linear systems of n equations in n+1unknowns. One possibility is to look for new periodic orbits of the family with the same period ( $\Delta{T}=0$) but these new periodic solutions do not exist in general. The normal procedure is to fix one of the initial conditions and allow the variation of the period (Szebehely 1967). Because a tangent displacement is a translation along the initial solution, it is avoided by constraining the variation of the initial conditions to a (n-1)-dimensional plane that is not tangent to the solution at $\vec{\xi}_0$.

Finally, we note that the previous procedure is directly applicable to nonautonomous systems

$$\dot{\vec{x}}=\vec{F}(t;\vec{x};\sigma),$$ (13)

but only when the dependence of F on the independent variable t is through periodic functions. Then, the period of the solution is the same of those periodic functions. Applications for that case have been done, for instance, in Scheeres (1998) and Peláez & Lara (2001).

2.2 Integrals

When system (1) admits one integral

$$\psi(\vec{x};\sigma)=\gamma,$$ (14)

a new constraint must be added to the periodicity condition Eq. (4), namely

$$\psi(\vec{\xi};\sigma)=\gamma.$$ (15)

Thus, the fundamental system (6) must be augmented with the equation

$$\psi(\vec{\xi}_0+\Delta\vec{\xi};\sigma_0+\Delta\sigma)=\gamma_0+\Delta\gamma,$$ (16)

that at first order is written as

$$\nabla_{\xi}\psi\cdot\Delta\vec{\xi}+ \frac{\partial\psi}{\partial\sigma}\Delta\sigma=\Delta\gamma.$$ (17)

We append one more equation to system (6) and also introduce a new unknown $\Delta\gamma$. We have to add an additional constraint to get as many unknowns as equations. Any relation between $\Delta{T}$, $\Delta\sigma$ and $\Delta\gamma$ would be suitable for this task, but there are three natural choices on which we focus this paper:

$\Delta{T}=0$: a family of periodic orbits with constant period is generated when varying the parameter $\sigma$, where the integral changes its value from orbit to orbit. The linear system

$$(\nabla_{\xi}\vec{x}-\vec{I})\cdot\frac{\Delta\vec\xi}{\Delta\sigma} =-\frac{\partial\vec{x}}{\partial\sigma}$$ (18)

will give the corrections $\Delta\vec\xi$ to the initial conditions and Eq. (17) will provide the new value of the energy integral.

$\Delta\sigma=0$: a natural family of periodic orbits is constructed where each orbit has a different period. The system to solve is now

$$\begin{array}{rcl} \vec{F}(\vec{x};\sigma)\displaystyle\fra... ...playstyle\frac{\Delta\vec\xi}{\Delta\gamma}&=&1.\\ \end{array}$$ (19)

$\Delta\gamma=0$: is the case of the isoenergetic family of periodic orbits with different period for each value of the parameter $\sigma$. For this case we need to solve the system

$$\begin{array}{rcl} \vec{F}(\vec{x};\sigma)\displaystyle\frac... ...laystyle\frac{\partial\psi}{\partial\sigma}\cdot\\ \end{array}$$ (20)

While system (18) is generally nondegenerate, systems (19) and (20) are not and must be adequately reduced as discussed in Siegel & Moser (1971).

  
3 Intrinsic formulation

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

$${\cal L}=\frac{1}{2}\dot{\vec{x}}\cdot\dot{\vec{x}}+ \dot{\vec{x}}\cdot\vec{A}(\vec{x})+W(\vec{x}),$$ (21)

where either the scalar function W or the vector one $\vec{A}$ can depend on one or more parameters $\gamma_i$ (i>0). Such systems are conservative and admit the orbital integral

$$\gamma=\frac{1}{2}\dot{\vec{x}}\cdot\dot{\vec{x}}-W(\vec{x}).$$ (22)

The equations of motion are

$$\ddot{\vec{x}}+\vec{B}\times\dot{\vec{x}} -{\vec{\nabla}_{\vec{x}}U}={\bf {0}},$$ (23)

where the vector function $\vec{B}=\vec{B}(\vec{x})$ is the curl $\vec{B}={\vec{\nabla}_{\vec{x}}\times\vec{A}}$ and the notation

$$U\equiv W(\vec{x})+\gamma_0,$$ (24)

has been introduced. The corresponding variational equations produced by a variation of any parameter (or integral) $\sigma$ of the dynamical system are

$$\begin{array}{l} \delta\ddot{\vec{x}}+ ({\vec{\nabla}_{\vec... ...ma}-\vec{B}_\sigma\times\dot{\vec{x}})\delta\sigma, \end{array}$$ (25)

where $\vec{H}={\vec{\nabla}_{\vec{x}}}({\vec{\nabla}_{\vec{x}}U})$is the Hessian in $\vec{x}$ of the effective potential function U. Equation (25) admits the variational integral

$$\dot{\vec{x}}\cdot\delta\dot{\vec{x}} -{\nabla_{\mbox{\footn... ...ta\vec{x} -\frac{\partial U}{\partial\sigma}\delta\sigma =0.$$ (26)

The variational integral Eq. (26) reduces by one the number of dimensions of the variational system and consequently the number of equations to be integrated and the number of integration constants to be determined. When using intrinsic variations not only the order of the variational system is reduced, but the variational equations turn out to be separable. The separability of the variations was a well-known fact for dynamical systems with two degrees of freedom, where the variational system can be reduced to a differential equation of order 2 for the normal displacement and a quadrature for the tangent displacement. Deprit (1981) has shown how the separability of the variational equations carries to dynamical systems with three degrees of freedom determined by the Lagrangian (21). For these systems, the separability of the tangent displacements is a consequence of the existence of the variational integral, and Palmore (1982) extended this separability to a wide range of conservative dynamical systems with ndegrees of freedom. Apparently, that separability was found independently in Karimov & Sokolsky (1989).

  
3.1 Rotation to the Frenet frame

One passes from Cartesian to intrinsic coordinates by means of a rotation $\vec{R}$ of the reference frame. Therefore, an intrinsic variation

$$\vec{s}=\left(\begin{array}{c}p\\ q\\ r\end{array}\right),$$ (27)

where p is the component in the tangent direction, q in the normal one and r in the binormal, is obtained from the corresponding $\delta\vec{x}$ by means of

$$\vec{s}=\vec{R}\cdot\delta\vec{x}.$$ (28)

In the simplest case of a dynamical system with two degrees of freedom the rotation angle

$$\phi=\arctan\displaystyle\frac{\dot{y}}{\dot{x}}$$ (29)

corresponds to the inclination of the direction of the velocity vector with respect to the x-axis and

$$\vec{R}=\left(\begin{array}{ccc} \phantom{-}\cos\phi&\sin\phi&0\\ -\sin\phi&\cos\phi&0\\ 0&0&1 \end{array}\right)\cdot$$ (30)

In the general case of three degrees of freedom the rotation matrix $\vec{R}$is given by the components of the Frenet frame projected onto the Cartesian frame,

$$\vec{R}=(\vec{t},\vec{n},\vec{b})^\tau,$$ (31)

where $\tau$ means transposition and

$$\vec{t}=\frac{1}{V}\dot{\vec{x}},\qquad \vec{n}=\frac{1}{N}\dot{\vec{t}},\qquad \vec{b}=\vec{t}\times\vec{n},$$ (32)

where

$$V^2=\dot{\vec{x}}\cdot\dot{\vec{x}},\qquad N^2=\dot{\vec{t}}\cdot\dot{\vec{t}}.$$ (33)

The necessary components of the Frenet frame and their derivatives are listed in Table 1. Note that the derivative  ${\rm d}^3\vec{x}/{\rm d}t^3$ of the equations of motion needs to be computed, which in principle is not a problem because the analytical expression for $\ddot{\vec{x}}$ is known; of course this derivative can also be obtained from numerical approximations.

   

Table 1: From left to right and from top to bottom, sequence to compute the elements of the rotation matrix Eq. (31) and their derivatives.

$V=\displaystyle\sqrt{\dot{\vec{x}}\cdot\dot{\vec{x}}}$	$\vec{t}=\frac{1}{V}\dot{\vec{x}}$
$\dot V=\frac{1}{V}\dot{\vec{x}}\cdot\ddot{\vec{x}}$	$\dot{\vec{t}}=\frac{1}{V}(\ddot{\vec{x}}-\dot V\vec{t})$
$N=\displaystyle\sqrt{\dot{\vec{t}}\cdot\dot{\vec{t}}}$	$\vec{n}=\frac{1}{N}\dot{\vec{t}}$
$\vec{b}=\vec{t}\times\vec{n}$	$\ddot V=\frac{1}{V}(\dot{\vec{x}}\cdot\frac{{\rm d}^3}{{\rm d}t^3}\vec{x}+ \ddot{\vec{x}}\cdot\ddot{\vec{x}}-\dot V^2)$
$\ddot{\vec{t}}=\frac{1}{V}(\frac{{\rm d}^3}{{\rm d}t^3}\vec{x}- 2\dot V\dot{\vec{t}}-\ddot V\vec{t})$	$\dot N={\vec{n}}\cdot\ddot{\vec{t}}$
$\dot{\vec{n}}=\frac{1}{N}(\ddot{\vec{t}}-\dot N\vec{n})$	$\dot{\vec{b}}=\vec{t}\times\dot{\vec{n}}$

3.2 Intrinsic variations

Intrinsic variations are separated into one system of differential equations involving, exclusively, the normal and binormal displacements

 

$$\ddot q Q_1~q+Q_2~\dot r+Q_3~r+(Q_4~U_\sigma+Q_5)\delta\sigma,$$

$$\ddot r R_1~r+R_2~\dot q+R_3~q+(R_4~U_\sigma+R_5)\delta\sigma,$$ (34)

and the quadrature

$$\displaystyle\frac{{\rm d}}{{\rm d}t}\left(\frac{p}{V}\right)... ...(\vec{B}\cdot\vec{n})r+\displaystyle\frac{U_\sigma}{V}\right],$$ (35)

producing the tangent displacement. The functions Q_i, R_i, ( $i=1,\dots,5$) are given in Table 2. Thus, the general solution of the variational system (34)-(35) depends only on four arbitrary integration constants. The tangent displacement will be obtained from (35) after solving Eq. (34).

   

Table 2: Auxiliary functions for Eq. (34). The matrix $\vec{W}$ is the product of the rotation matrix $\vec{R}$ and the Hessian $\vec{H}$ of U.

$M_{\phantom{1}}=\dot{\vec{b}}\cdot\vec{n}$

$R_1=M(M-\vec{B}\cdot\vec{t})-(\vec{B}\cdot\vec{n})^2 +(\vec{W}\cdot\vec{b})\cdot\vec{b} +V[{\vec{\nabla}_{\vec{x}}\vec{B}}\cdot\vec{b}]\cdot\vec{n}$

$R_2=(M-\vec{B}\cdot\vec{t})+M$

$R_3=(N+\vec{B}\cdot\vec{b})(\vec{B}\cdot\vec{n})+\dot M +(\vec{W}\cdot\vec{n})\cdot\vec{b} +V[{\vec{\nabla}_{\vec{x}}\vec{B}}\cdot\vec{n}]\cdot\vec{n}$

$R_4=(\vec{B}\cdot\vec{n})/V$

$R_5=\vec{b}\cdot{\vec{\nabla}_{\vec{x}}U_\sigma} +V(\vec{n}\cdot\vec{B}_\sigma)$

$Q_1=M(M-\vec{B}\cdot\vec{t})-(\vec{B}\cdot\vec{b})^2 +(\vec{W}\cdot\vec{n})\cdot\vec{n} -V[{\vec{\nabla}_{\vec{x}}\vec{B}}\cdot\vec{n}]\cdot\vec{b}$

Q₁= $-3N(N+\vec{B}\cdot\vec{b})$

Q₂=-R₂

$Q_3=(2N+\vec{B}\cdot\vec{b})(\vec{B}\cdot\vec{n})-\dot M +(\vec{W}\cdot\vec{n})\cdot\vec{b} -V[{\vec{\nabla}_{\vec{x}}\vec{B}}\cdot\vec{b}]\cdot\vec{b}$

$Q_4=-(2N+\vec{B}\cdot\vec{b})/V$

$Q_5=\vec{n}\cdot{\vec{\nabla}_{\vec{x}}U_\sigma}- V(\vec{b}\cdot\vec{B}_\sigma)$

For details on the derivation of the Eqs. (34)-(35) the reader is referred to Deprit (1981) where his $\gamma$ has been replaced by the partial derivative $U_\sigma$ and the terms Q₅ and R₅ have been added, respectively, to the right hand sides of the Eqs. (20₁) and (20₂) of his paper; these extra terms come from the right hand side of the nonhomogeneous Eqs. (25). Note that the radius of curvature $\rho$ and the radius of torsion ${\cal{T}}$ that appear in the formulas given by Deprit (1981) can be obtained from

$$\displaystyle\frac{V}{\rho}=\displaystyle\sqrt{\dot{\vec{t}}\... ...uad\displaystyle\frac{V}{{\cal{T}}}=\dot{\vec{b}}\cdot\vec{n}.$$ (36)

Note, too, that the integration of the variational system Eqs. (34) needs the computation of the derivative of M, present in Q₃ and R₃.

  
3.3 The intrinsic predictor-corrector

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

$$\vec{x}=\vec{x}(t;\vec{\xi},\dot{\vec{\xi}};\sigma),$$ (37)

is given by

$$\delta\vec{x}=\nabla_{\xi}\vec{x}\cdot\delta\vec\xi +\nabla_... ...vec\xi} +\frac{\partial\vec{x}}{\partial\sigma}~\delta\sigma.$$ (38)

Therefore, in the differential approximation ( $\delta\equiv\Delta$) the linear system (6) is equivalent to

$$\Delta{T}~\dot{\vec{x}}(T_0)+\Delta\vec{x}(T_0)-\Delta\vec\xi= \vec\xi_0-\vec{x}(T_0),\\ [1.5ex]$$ (39)

$$\Delta{T}~\vec{F}(\vec{x}(T_0),\dot{\vec{x}}(T_0)) +\Delta\do... ...(T_0)-\Delta\dot{\vec\xi}= \dot{\vec\xi}_0-\dot{\vec{x}}(T_0),$$ (40)

where we omitted the dependence of x and $\dot{\vec{x}}$ on the initial conditions $\vec{\xi}_0$, $\dot{\vec{\xi}}_0$and on the parameter $\sigma$ for the sake of brevity.

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₀. We get

$$\begin{array}{l} \Delta{T}~V(T_0)~\vec{t}(T_0)+\vec{s}(T_0)-... ...*{4cm}=\vec{R}(T_0)\cdot[\vec{\xi}_0-\vec{x}(T_0)]. \end{array}$$ (41)

3.3.1 Tangent prediction

For the predictor stage - being exactly periodic the initial orbit given by Eq. (37) - the right hand side of Eq. (41) vanishes and $\vec{R}(T_0)=\vec{R}(0)$. Therefore

$$\Delta{T}~V(T_0)~\vec{t}(T_0)+\vec{s}(T_0)-\vec{s}(0)={\bf0}.$$ (42)

That is

   

$$\Delta{T}~V(T_0)+p(T_0)-p(0)$$ (43)

q(T₀)-q(0)=0, (44)

r(T₀)-r(0)=0. (45)

The particular solution $\vec{s}=\vec{s}^*$ of the variational equations that fulfills the periodicity condition (42) depends on four integration constants $\vec{\alpha}=(\alpha_1,\alpha_2,\alpha_3,\alpha_4)^\tau$. The values of these constants are computed from the linear system formed by Eqs. (44)-(45) and their respective derivatives:

$$[\vec{M}(T_0)-\vec{M}(0)]\cdot\vec{\alpha}=-\vec{\beta}(T_0),$$ (46)

where

$$\vec{M}(t)= \left(\begin{array}{cccc} q_1(t)&q_2(t)&q_3(t)&... ..._1(t)&\dot r_2(t)&\dot r_3(t)&\dot r_4(t) \end{array}\right),$$ (47)

and

$$\vec{\beta}=\Big(q_5(t),r_5(t),\dot q_5(t),\dot r_5(t)\Big)^\tau.$$ (48)

System (46) is nondegenerate except for some critical points that, as we will see below, correspond to bifurcations of the dynamical system.

Once $\vec{s}^*$ is obtained the corrections the initial conditions are easily obtained from

$$\begin{array}{ll} \Delta\vec\xi&=\vec{R}^\tau(0)\cdot\vec{s}^... ...c{s}^*(0) +\vec{R}^\tau(0)\cdot\dot{\vec{s}}^*(0). \end{array}$$ (49)

The correction to the period is obtained then from Eq. (43)

$$\Delta{T}=\frac{p(0)-p(T_0)}{V(T_0)}\cdot$$ (50)

Note that the value of $\delta\sigma$ is irrelevant in the integration of Eqs. (34)-(35). Then one can compute the variations (p,q,r) scaled by $\delta\sigma$. If that is the case, the corrections Eqs. (49)-(50) must be multiplied by  $\delta\sigma$.

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)

$$\begin{array}{ll} \vec\xi&=\vec\xi_0+\Delta\vec\xi,\\ [1ex] \... ...\vert\dot{\vec\xi}_0+\Delta\dot{\vec\xi}\vert}\cdot \end{array}$$ (51)

3.3.2 Isoenergetic corrections

The solution Eq. (37) is only approximately periodic. Thus, the right hand side of (41) no longer vanishes. Assuming that the difference $\Delta\vec{R}=\vec{R}(T_0)-\vec{R}(0)$ is of the same order as the desired corrections we neglect $\Delta\vec{R}\cdot\Delta\vec\xi$. Then Eq. (41) is written

$$\Delta{T}~V(T_0)~\vec{t}(T_0)+\vec{s}(T_0)-\vec{s}(0) =\vec{R}(T_0)\cdot[\vec{\xi}_0-\vec{x}(T_0)],$$ (52)

and the procedure for computing the displacements $\vec{s}^*$ and $\dot{\vec{s}}^*$ follows the path described for the tangent predictions. In this case we deal only with the homogeneous part of Eq. (34). The four necessary integration constants are again obtained from Eq. (46), where now

$$\vec{\beta}=\left(\begin{array}{c} \vec{n}^\tau(t)\\ \vec{b}^... ...) \end{array}\right)\cdot[\dot{\vec{x}}(t)-\dot{\vec{\xi}}_0].$$ (53)

Equations (49) provide the corrections to the initial conditions and the correction to the period is obtained from

$$\Delta{T}= \frac{[\vec{\xi}_0-\vec{x}(T_0)]\cdot\vec{t}(T_0)+p(0)-p(T_0)}{V(T_0)},$$ (54)

where the tangent displacement has been computed integrating the homogeneous part of Eq. (35).

Finally, Eq. (51) provides the desired corrections.

  
3.3.3 Practical procedure

Despite the fact that the intrinsic variational equations would be much more complicated than the Cartesian ones, one should note that there is no need to solve them. At any point, it is possible to transform the solutions of the variational equations between rectangular and intrinsic coordinates by simple rotations. In fact, for the computation of the four necessary integration constants we only need to know the intrinsic displacements $\vec{s}^*$ and $\dot{\vec{s}}^*$ at the times t=0 and t=T. Therefore, in practice we integrate the much simpler variational equations in rectangular coordinates, proceeding as follows.

We choose initial conditions

$$p_i(0)=0,\qquad i=1,\dots,4,$$ (55)

and

$$\vec{M}(0)=\vec{I}_4,$$ (56)

where $\vec{I}_4$ is the $4\times4$ identity matrix; each value $\dot{p}_i(0)$ is computed from the derivative

$$\dot{p}_i=\frac{\dot{V}}{V}~p_i+ (2N+\vec{B}\cdot\vec{b})~q_i-(\vec{B}\cdot\vec{n})~r_i,$$ (57)

of the quadrature Eq. (35), evaluated at t=0. The rotation matrices $\vec{R}(0)$, $\dot{\vec{R}}(0)$ are then computed, and using the inverse rotations

 

$$\delta\vec{x}_i(0) \vec{R}^\tau(0)\cdot\vec{s}_i(0),$$ =

$$\delta\dot{\vec{x}}_i(0) \dot{\vec{R}^\tau}(0)\cdot\vec{s}_i(0) +\vec{R}^\tau(0)\cdot\dot{\vec{s}}_i(0),$$ = (58)

we obtain the four sets of the initial Cartesian variations. Once we have integrated those Cartesian variations through a period, we apply the direct rotations

 

$$\vec{s}_i(T) \vec{R}(T)\cdot\delta\vec{x}_i(T),$$ =

$$\dot{\vec{s}}_i(T) \dot{\vec{R}}(T)\cdot\delta\vec{x}_i(T) +\vec{R}(T)\cdot\delta\dot{\vec{x}}_i(T),$$ = (59)

at t=T, and obtain $\vec{s}_i(T)$, $\dot{\vec{s}}_i(T)$. Then we are ready for the computation of the four integration constants from Eq. (46), and finally the corrections from Eq. (49).

  
4 Brief comment on stability

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 $(\lambda_i,1/\lambda_i)$, (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

$${\rm {det}}\vert\vec{M}(T_0)-\lambda\vec{I}_4\vert=0.$$ (60)

Remark here that the critical case of any of the nontrivial eigenvalues having modulus 1 is a singularity of the algorithm provided in this paper. In such case the linear system Eq. (46) associated with the periodicity condition is degenerate.

Taking in account that the eigenvalues appear in reciprocal pairs, Eq. (60) is

$$\lambda^4+a_1\lambda^3+a_2\lambda^2+a_1\lambda+1=0,$$ (61)

or

$$(\lambda^2-k_1\lambda+1)~(\lambda^2-k_2\lambda+1)=0,$$ (62)

where a₁=-(k₁+k₂) and a₂=2+k₁k₂. Thus, the two stability indices $k_1=\lambda_1+1/\lambda_1$ and $k_2=\lambda_2+1/\lambda_2$, given by

$$k_i=-{\frac{1}{2}}a_1-(-1)^i~ \displaystyle\sqrt{2+({\frac{1}{2}}a_1)^2-a_2},\qquad(i=1,2)$$ (63)

are normally used. The condition $k_i\in{\rm {I}R}$ and |k_i|<2gives linear stability while any other possibility means instability.

Finally, in reference to planar solutions, one should note that the two stability indices k₁ and k₂ 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

$$k_{\rm n}=q_1(T_0)+\dot{q}_3(T_0)$$ (64)

while the out of plane or binormal stability index simply is

$$k_{\rm b}=r_2(T_0)+\dot{r}_4(T_0).$$ (65)

The value $\vert k_{\rm n}\vert=2$ means that possibly new families of periodic orbits bifurcate in the plane and $\vert k_{\rm b}\vert=2$ means that possibly three-dimensional families bifurcate out of the plane.

  
5 An example

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.

5.1 Equations of motion and variations

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

 

$$% \ddot{x}-2\dot{y}$$ = W_x,

$$\ddot{y}+2\dot{x}$$ = W_y, (66)

$$\ddot{z}$$ = W_z,

with

$$W=\frac{1}{2}\left(x^2+y^2\right)+\displaystyle\frac{1-\mu}{r_1} +\displaystyle\frac{\mu}{r_2},$$ (67)

and

 

$$(x+\mu)^2+y^2+z^2,$$

$$(x-1+\mu)^2+y^2+z^2.$$ (68)

As usual, the units are chosen in such a way that the distance between the primaries, the mean motion of the Keplerian orbit of the primaries and the Gauss constant are equal to one.

There is an integral of these equations, the Jacobian constant, that we write

$$\frac{1}{2}\left(\dot{x}^2+\dot{y}^2+\dot{z}^2\right)-W(x,y,z;\mu)=h.$$ (69)

Taking the value h as a parameter and introducing the notation

$$U\equiv h+W(x,y,z;\mu),$$ (70)

integral Eq. (69) is written in the form

$$\frac{1}{2}\left(\dot{x}^2+\dot{y}^2+\dot{z}^2\right)-U(x,y,z;\mu,h)=0.$$ (71)

The variations are Eqs. (34) where Q₅ and R₅ depend on the parameter that varies. Other functions are computed from Tables 2 and 1 where $\vec{B}$ has modulus 2 and the z-axis direction, its gradient is zero and the elements of the Hessian are

$$\begin{array}{ll} U_{xx}=& 1+\displaystyle\frac{1-\mu}{r_1^5}... ...aystyle\frac{\mu}{r_2^5}\left(3z^2-r_2^2\right).\\ \end{array}$$ (72)

5.2 Families for variations of the Jacobi constant

For both primaries with equal masses ($\mu=1/2$) let us vary the Jacobi constant ( $U_h\equiv{1}$ and both Q₅ and R₅ 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 $T\approx5.585$for a retrograde orbit. The corresponding initial conditions of this very rough approximation of a periodic orbit are $x_0=z_0=\dot{y}_0=\dot{z}_0=0$, y₀=4 and $\dot{x}_0=4.5$. The maximum difference in absolute value between each of these initial conditions and the same coordinate or velocity at the approximate period T is $\varepsilon_0=1.1\times10^{-2}$, but three iterations are enough for the corrector to obtain the sequence $\varepsilon_1=4.1\times10^{-4}$, $\varepsilon_2=1.7\times10^{-7}$and finally the improved periodic orbit with $\varepsilon_3<10^{-13}$. The corrected initial conditions and period for $x=z=\dot{y}=\dot{z}=0$ are

y=3.96199469992294,

$$\dot{x}=4.46677589984367,$$

T=5.57243120610132.

For a variation $\Delta{h}=-0.05$ the predictor computes the initial conditions of an orbit that is approximately periodic, again with $\varepsilon_0=1.1\times10^{-2}$. The sequence of corrections is now $\varepsilon_1=1.2\times10^{-4}$, $\varepsilon_2=2.6\times10^{-11}$and finally $\varepsilon_3<10^{-14}$. Smaller values $\Delta{h}$allow of course much better predictions.

	Figure 1: Normal $k_{\rm n}$ and binormal $k_{\rm b}$ stability indices versus the Jacobi constant h for the family m of periodic orbits. The horizontal dashed lines correspond to the critical values $\pm$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 $\Delta{h}$ cannot be fixed constant for all the computations and for orbits closer to the primaries, smaller $\Delta{h}$ values are required to obtain predictions amenable to improvement by the corrector. Figure 1 shows the behavior of the normal $k_{\rm n}$ and binormal $k_{\rm b}$stability indices of this planar family; due to the large values of both indices, a scale proportional to $\sin h^{-1}k$ is used. While $k_{\rm n}=2$for $h\approx0.30818$ - orbit m₁ in Hénon (1965) - the binormal index $k_{\rm b}$ takes the critical values -2 for the bifurcation value $h\approx0.8$ and +2 for $h\approx0.5$ - respectively orbits m_1v and m_2v of Hénon (1973) -. The initial conditions below ( $x=z=\dot{y}=\dot{z}=0$) correspond to the critical orbits.

$$\begin{tabular}{lccc@{}}\hline% & \multicolumn{1}{c}{$m_1$} &... ...k_{\rm n}$ & $+2.00$ & $-1.45$ & $-1.76$\\ \hline\end{tabular}$$

In the case $k_{\rm b}=-2$ 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₄. We call this family the L₄-family and of course the symmetrical L₅-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 $x=z=\dot{y}=\dot{z}=0$

y=1.209894432634087,

$$\dot{x}=2.139024877316912,$$

for which h₀=0.791930530821579. Then, we slightly modify the initial velocity in the z-axis direction by setting $\dot{x}=2.138$and computing $\dot{z}$ from the integral Eq. (69) with h=h₀; two symmetric solutions exist for $\pm\dot{z}$. For an approximate period T=7.041614672725651 - that doubles the period of the planar periodic orbit - these modified initial conditions allow the corrector to converge to a three-dimensional periodic orbit of the bifurcated family with the sequence $\varepsilon_0=9.2\times10^{-2}$, $\varepsilon_1=1.2\times10^{-3}$, $\varepsilon_2=1.5\times10^{-4}$, $\varepsilon_3=6.4\times10^{-7}$, $\varepsilon_4=6.2\times10^{-11}$, and $\varepsilon_5<10^{-13}$. The three-dimensional corrected periodic solution is ( $x=z=\dot{y}=0$, h=h₀)

y=1.21460588387117,

$$\dot{x}=2.13950776590580,$$

T=7.06289508950945.

Once one periodic orbit of the three-dimensional family is obtained, the convergence of the algorithm is indeed much better. Figure 2 presents the stability indices of the bifurcated L₄-family that are now k₁ and k₂ because this family is not planar. The L₄-family remains with linear stability until $h\approx0.48$ where k₁>2, whereas at $h\approx -1.234039$ (k₁=k₂) the family enters a zone of complex instability remaining there until its termination at the equilateral L₄ equilibrium. Figure 3 shows the two periodic orbits of this family when k₁=k₂; the respective initial conditions are (

$x=z=\dot{y}=0$)

h	0.595060530821579	-1.234038469178420
y	1.174109021245709	0.850091106717253
$\dot{x}$	1.972301154871868	0.111468952994908
$\dot{z}$	0.495891121768287	0.519684176787866
T	7.063185569369677	6.363331100559023

	Figure 2: Stability indices versus the Jacobi constant for the L4-family of periodic orbits. For h<-1.234039, k1 and k2 are complex numbers.
Open with DEXTER

	Figure 3: Periodic orbits of the L4-family when k1=k2. Left $h\approx 0.59506$, right $h\approx -1.234039$. Primaries are represented by dots.
Open with DEXTER

	Figure 4: Stability indices k1 and k2 versus the Jacobi constant h for the L1-family of periodic orbits. Both indices are complex numbers for -0.47480835<h<0.46334917.
Open with DEXTER

	Figure 5: Stability behavior of the L1-family on the a1-a2 plane. A scale proportional to $\sinh^{-1}a_{1,2}$ is used.
Open with DEXTER

In the case $k_{\rm b}=+2$ we found a transversal bifurcation. We call the bifurcated family the L₁-family because it ends at the L₁ 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.4952105308215789 the initial conditions of a periodic orbit of the planar family are

$$x=z=\dot{y}=\dot{z}=0,$$

y=0.775707482134065,

$$\dot{x}=1.938879462050800.$$

For the period T₀=2.467597374742529 of this orbit, making $\dot{x}=1.93887$ and computing $\dot{z}$ from the integral Eq. (69) with h=h₀, the corrector algorithm converges to the three-dimensional periodic solution ( $x=z=\dot{y}=0$)

y=0.775715315248107,

$$\dot{x}=1.938876894065857,$$

T=2.4676231846873.

Again we find two possibilities associated with the sign ($\pm$) of $\dot{z}$. The stability k₁ and k₂ indices of the L₁-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₁-a₂ 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₁ and a₂ have been scaled proportional to their respective $\sin h^{-1}$. Figure 6 shows the two periodic orbits of this family when k₁=k₂ and the respective initial conditions ( $x=z=\dot{y}=0$) are listed below.

h	-0.474808349178421	0.463349173821579
y	0.630644393432460	0.772133837453028
$\dot{x}$	1.116409595398506	1.913885053382574
$\dot{z}$	0.828732085352047	0.184702299430020
T	2.722825258179725	2.474956115380520

	Figure 6: Periodic orbits of the L1-family when k1=k2. Left $h\approx 0.46334917$, right $h\approx -0.47480835$. Primaries are represented by dots.
Open with DEXTER

	Figure 7: Stability behavior of the family for variations of $\mu$. A scale proportional to $\sinh^{-1}a_{1,2}$ is used.
Open with DEXTER

5.3 Variations of $\mu$

We vary now the real parameter $\mu$. The partial derivative

$$U_\mu=\displaystyle\frac{1}{r_2} -\displaystyle\frac{1}{r_1} ... ...tyle\frac{x+\mu}{r_1^3} -\mu\displaystyle\frac{x-1+\mu}{r_2^3}$$ (73)

depends on the position vector and the functions Q₅ and R₅ in Eq. (34) are no longer zero, so the gradient of $U_\mu$ must be computed. A periodic solution of the case $\mu=1/2$ serves as starter of the algorithm and we use an orbit of the L₄-family for that purpose: the orbit that enters the area of complex instability ( h=-1.23403946917842). The initial conditions of that orbit are

$$z=\dot{y}=0,$$ y=0.8500912035204521,

$$\dot{x}$$

$$\dot{z}$$

and the period is T=6.363330478421584The new family is then computed for decreasing values of $\mu$ from $\mu=0.5$ until $\mu=0$. Figure 7 presents the stability behavior of this family, where we appreciate that starting from complex instability the family enters Broucke's region 1 of linear stability. That happens for a value $\mu_{\rm K}=0.057246492698$and the family remains stable for smaller values  $\mu<\mu_{\rm K}$ of the parameter.

In this case the value $\Delta\mu$ can be fixed for all the family and, for instance, for a value $\Delta\mu=-0.003$ the typical sequence of iterations of the predictor-corrector procedure is $\varepsilon_0<10^{-3}$, $\varepsilon_1<10^{-7}$and $\varepsilon_2<10^{-13}$. Greater values $\vert\Delta\mu\vert$ can be chosen and even for $\Delta\mu=-0.1$ the algorithm shows good convergence.

Conclusions

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).

References

 
• Belbruno, E., Llibre, J., & Ollé, M. 1994, Celest. Mech., 60, 99 In the text  
• Broucke, R. 1969, AIAA Journal, 7, 1003 In the text  
• Caranicolas, N. D. 1994, A&A, 287, 752 In the text NASA ADS  
• Contopoulos, G., & Barbanis, B. 1994, Celest. Mech., 59, 279 In the text  
• Deprit, A. 1981, Celest. Mech., 24, 185 In the text  
• Deprit, A., & Henrard, J. 1967, AJ, 72, 158 In the text NASA ADS  
• Gómez, G., Llibre, J., Martínez, R., & Simó, C. 1985, Station keeping of libration points, ESOC Con. 5648/83/D/JS(SC) In the text  
• Hénon, M. 1965, Ann. Astrophys., 28, 992 In the text  
• Hénon, M. 1973, A&A, 28, 415 In the text NASA ADS  
• Henrard, J., & Lemaitre, A. 1986, Celest. Mech., 39, 213 In the text  
• Howell, K. C. 1984, Celest. Mech., 32, 53 In the text  
• Howard, J. E., & MacKay, R. S. 1987, J. Math. Phys., 28, 1036 In the text  
• Karimov, S. R., & Sokolsky, A. G. 1989, Celest. Mech., 46, 335 In the text  
• Markellos, V. V. 1981, Celest. Mech., 25, 3 In the text  
• Palmore, J. 1982, Celest. Mech., 27, 151 In the text  
• Peláez, J., & Lara, M. 2001, Advances in the Astronautical Sciences, 108, 1189 In the text  
• Scheeres, D. J. 1998, Celest. Mech., 70, 75 In the text  
• Scheeres, D. J. 1999, Satellite Dynamics About Asteroids: Computing the Poincaré Map for the General Case, Hamiltonian Systems with Three or More Degrees of Freedom, ed. C. Simó, NATO ASI Ser. 533, 554 In the text  
• Siegel, C. L., & Moser, J. K. 1971, Lectures on Celestial Mechanics (Springer-Verlag) In the text  
• Strömgren, E. 1935, Bull. Astron., Paris 9, 87 In the text  
• Szebehely, V. 1967, Theory of Orbits (Academic Press, N. Y.) In the text