# Approximate Analytical Methodology for the Restricted Three-Body and Four-Body Models Based on Polynomial Series

Authors: Ying-Jing Qian; Xiao-Dong Yang; Lei-Yu Yang; Wei Zhang
Journal: International Journal of Aerospace Engineering
 (2016)

---

## Abstract

The restricted three-body problem (R3BP) and restricted four-body problem (R4BP) are modeled based on the rotating frame. The conservative autonomous system for the R3BP and nonautonomous system with period parametric resonance due to the fourth body are derived. From the vibrational point of view, the methodology of polynomial series is proposed to solve for these problems analytically. By introducing the polynomial series relations among the three directions of motion, the three-degree-of-freedom coupled equations are transferred into one degree-of-freedom containing the full dynamics of the original autonomous system for the R3BP. As for the R4BP case, the methodology of polynomial series combined with the iterative approach is proposed. During the iterative approach, the nonautonomous system can be treated as pseudoautonomous equation and the final polynomial series relations and one-degree-of-freedom system can be derived iteratively.


---

## Body

## 1. Introduction

In our solar system, the planets move around the Sun, while the moons move around their host planets, which in turn also move around the Sun. They all follow the two-body theory or the perturbed two-body theory which has well-known closed form solution given in terms of elementary functions. However, to study the motion of artificial satellites in the solar system, the problem is too complicated to obtain similar types of solutions. Compared to the mass of either the Sun, the planet, or the moon, the mass of the artificial satellite is always infinite, which means that the gravitational influence of the artificial satellite on the primaries can be omitted. According to the different mission purposes, the system can be modeled correspondingly as therestricted three-body problem (R3BP) or therestricted four-body problem (R4BP). For example, if the motion of the satellite is limited in the planet-moon-moon system, such as the Saturn-Tethys-Telesto-spacecraft, Saturn-Tethys-Calypso-spacecraft, and Saturn-Dione-Helen-spacecraft systems, they are typical R4BP systems. If the motion of the satellite is limited in the planet-moon system while the Sun’s gravitational influence cannot be neglected, such as ARTEMIS mission [1], the system needs to be modeled as the R4BP as well. Meanwhile, for the mission exploring the Sun-planet system, such as GALA and PLANCK [2], the R3BP model is suitable when we consider the Sun and the barycenter of the planet and its moon as two dominant masses move around their center of mass along either circular or elliptic orbits.At first glance, the difficulties of the R3BP and R4BP are not obvious, especially when considering that the two-body problem can be analytically solved. In the past, many physicists, astronomers, and mathematicians attempted unsuccessfully to find closed form solutions to the three-body problem. Such solutions do not exist because motions of the bodies are in general unpredictable, which makes the problem one of the most challenging problems in the history of science.The analytical techniques developed for studying the dynamics around the CR3BP and the R4BP models are similar. Although the exact solutions for R3BP and R4BP cannot be derived, classical perturbation methods are available to construct the analytical approximate solutions of the periodic motions. The analytical approximate solutions propose a way to inspect the trajectory and exhibit possible parameter design. The most representative work was successfully demonstrated by Richardson [3] who derived the third-order analytical solution of halo orbits with small parametric expansions in CR3BP with Lindstedt-Poincaré (L-P) method. The third-order analytical solutions of halo orbits are extensively used as initial conditions in mission trajectories designing progress for numerical integration [4]. Similarly, Gómez and Marcote [5] solved a reduced case of the CR3BP, namely, the Clohessy-Wiltshire equation with L-P method. Jorba and Masdemont [6] expanded the Lissajous and halo orbits of CR3BP as formal series of two amplitudes and semianalytically constructed high-order solutions with the assistance of computer.Similar to CR3BP, it has been analytically proved that the R4BP has some relative equilibrium solutions like the RTBP, and there are families of periodic orbits in the vicinity of equilibrium points [7]. The quasi-halo orbits in the Sun-Earth-Moon system were computed from the bicircular model in R4BP by Andreu [8]. Papadakis [9] computed the invariant stable and unstable manifolds in R4BP. For the Sun-Jupiter-Trojan asteroid-spacecraft system, the evolution of the families of periodic orbit and their stability around the asteroid and/or the Jupiter in R4BP were calculated by Baltagiannis and Papadakis [7, 10]. The periodic orbits in the R4BP with two equal masses were explored by Burgos-García and Delgado [11]. Recently, Lei and Xu [12] investigated the motions around the equilibrium points of R4BP, in which the motions are expressed as formal series of long, short, and vertical periodic amplitudes. By means of L-P method, series solutions are constructed up to a certain order.In this study, we model the CR3BP and bicircular case of R4BP by the conservative autonomous system and nonautonomous parametric excited system, respectively, from the vibration point of view, which may shed a new light on the study of periodic orbit investigations.The classical celestial perturbation methods can give an approximate analytical solution by correcting the amplitude and the period of the linear solutions with the nonlinear terms in the three-dimensional space. However, less attentions have been paid to the relationship of the motions in each dimension of the major and minor axes and their contributions to the nonlinear dynamics. The nonlinear complex modes analysis based on the invariant manifold method by Shaw and Pierre [13–15] and Avramov [16] could provide a chance to focus on the nonlinear relations of the motions in different directions. Inspired by Shaw and Pierre and Avramov, the polynomial expansion method for autonomous system and iterative polynomial expansion method for nonautonomous system are innovatively proposed in this paper. By introducing polynomial expansion relations, we transfer the three-degree-of-freedom periodic motions for the nonlinear governing equations into one-degree-of-freedom system for the autonomous CR3BP system. By generating the pseudoautonomous equation for the R4BP system, periodic motions can be obtained by iterative polynomial expansion method.The remainder of this paper is structured as follows. Section2 introduces the dynamical model R4BP and the CR3BP model. Section 3.1 discusses how to construct high-order series expansions of periodic solution in the CR3BP with polynomial expansion method. In Section 3.2, the iterative polynomial expansion method by transferring the parameter resonance nonautonomous three-degree-of-freedom equations into a pseudoautonomous system is introduced. At last, the conclusions are drawn in Section 4.


## 2. Equations of Motion

### 2.1. Restricted Four-Body Problem

The geometry of the R4BP is shown in Figure1. The finite masses P1, P2, and P4 and the infinite mass body P3 are denoted as m1, m2, m4, and m3, respectively. The motion of P3 is influenced by the attraction of P1, P2, and P4, but the motions of the finite masses P1, P2, and P4 are assumed not to be affected by P3.Figure 1
The configuration of R4BP.In the R4BP,P1 and P2 move about their barycenter perturbed by P4. The masses P1, P2, and P4 all move about the mass center of the entire system. It is assumed that all bodies are modeled as point masses, and the primaries P1, P2, and P4 move in a plane, with the distance from P4 to the P1-P2 system being much greater than the distance from P1 to P2. However, P3 is free to move in all three dimensions.We define the inertial coordinate system,O-ξηζ, with origin at the center of mass of the four-body system, O. The coordinate system, O-XYZ, rotates at an instantaneous rate of Ω1 about the mass center of the four-body system, O. The X-axis points to the barycenter of P1 and P2, B. The synodic coordinate system is B-xyz with origin at B. The x-axis points from the barycenter, B, to P2. B-xyz moves with respect to the coordinate system, O-XYZ, at an instantaneous rate of Ω. θ is the angle between the X-axis and the x-axis, counterclockwise measured. Meanwhile, B-xyz rotates with respect to the inertial coordinate system, O-ξηζ, at an instantaneous rate of Ω2=Ω+Ω1.Nondimensionalization is applied to simplify the form of the equation. The unit length is set as the distance betweenP1 and P2. The unit time is set as T/2π, where T is the period of the P1-P2 system. We introduce μ as the mass parameter of the system; then(1)µ=m2m1+m2,1-µ=m1m1+m2,μ4=m4m1+m2.The massive P1, with mass 1-µ, is located at (-µ,0,0) in B-xyz, and P2, with mass µ, is located at (1-µ,0,0) in B-xyz.The motion ofP3 is governed by the following equations with dimensionless time and space with respect to the coordinate system, O-ξηζ.(2)ρ¨=-1-μR13R1-μR23R2-μ4R43R4,in which R1, R2, and R4 are the position vectors from P1, P2, and P4 to P3, respectively. And ρ is the position vector from the system barycenter, O, to P3, with geometrical relation as (3)ρ=ROB+r.Note that ROB and r are presented in O-XYZ and B-xyz, respectively, and can be expressed as(4)ROB=ROBX,r=xx+yy+zz.Apparently, we have time derivative ofρ with respect to the inertial frame as(5)ρ˙I=R˙OBF+Ω1×ROB+r˙R+Ω2×r,(6)ρ¨I=R¨OBF+2Ω1×R˙OBF+Ω˙1×ROB+Ω1×Ω1×ROB+r¨R+2Ω2×r˙R+Ω˙2×r+Ω2×Ω2×r,in which the superscripts “F” and “R” mean the derivatives are measured with respect to the rotating coordinates O-XYZ and B-xyz. Combining (2) and (6), the motion of P3 in synodic coordinate system B-xyz is obtained as(7)r¨R+2Ω2×r˙R+Ω˙2×r+Ω2×Ω2×r+R¨OBF+2Ω1×R˙OBF+Ω˙1×ROB+Ω1×Ω1×ROB=-1-μR13R1-μR23R2-μ4R43R4,where(8)R1=x+μx+yy+zz,R2=x-1+μx+yy+zz,R4=xx+yy+zz+ROB+ROP4X.Then, we introduce the transformation matrix betweenO-XYZ and B-xyz:(9)XYZ=cos⁡θsin⁡θ0-sin⁡θcos⁡θ0001xyz;R4 can be expressed as(10)R4=ROB+ROP4cos⁡θ+xx+ROB+ROP4sin⁡θ+yy+zz.Thus, we obtain the completed equation of motion of R4BP in B-xyz as(11)x¨-2Ω2y˙-Ω˙2y=∂U∂x-R¨OB-ROBΩ12cos⁡θ-2R˙OBΩ1+ROBΩ˙1sin⁡θ,y¨+2Ω2x˙+Ω˙2x=∂U∂y-R¨OB-ROBΩ12sin⁡θ-2R˙OBΩ1+ROBΩ˙1cos⁡θ,z¨=∂U∂z,where U is the pseudopotential function of the four-body problem and is expressed as(12)U=12Ω22x2+y2+1-μR1+μR2+μ4R4.


### 2.2. Bicircular Model in Restricted Four-Body Problem

The bicircular model is a special case of the R4BP and built on the hypothesis that the motions of three primary bodies are restricted. Specifically,P1 and P2 are revolving in circular orbits around their barycenter and the P1-P2 barycenter, B, is moving in a circular orbit around the center of mass of the system O. Besides that, the two orbital planes are assumed to coincide with each other [17].With the aforementioned hypothesis, we consider(13)Ω2=T2π=1,Ω˙2=0,R¨OB=R˙OB=0,Ω˙1=0,ROB+ROP4≈ROB.Equation (11) is simplified as(14)x¨-2y˙=1-μR13x+μ-μR23x-1+μ-μ4R43ROB+ROP4cos⁡θ+x+ROBΩ12cos⁡θ,y¨+2x˙=-1-μR13y-μR23y-μ4R43ROB+ROP4sin⁡θ-y+ROBΩ12sin⁡θ,z¨=-1-μR13z-μR23z-μ4R43z.With expansion of 1/R4,(15)1R43≈ROB+ROP4-31-3xcos⁡θROB+ROP4-3ysin⁡θROB+ROP4≈ROB-31-3xcos⁡θ+ysin⁡θROB,equation of motion (14) can be further reduced to(16)x¨-2y˙=x-1-μR13x+μ-μR23x-1+μ-μ4ROB3x-μ4ROB2cos⁡θ-3μ4xcos2⁡θ+ysin⁡θcos⁡θROB3+3μ4x2cos⁡θ+xysin⁡θROB4+ROBΩ12cos⁡θy¨+2x˙=y-1-μR13y-μR23y-μ4ROB3y-μ4ROB2sin⁡θ-3μ4xsin⁡θcos⁡θ+ysin2⁡θROB3+3μ4xycos⁡θ+y2sin⁡θROB4+ROBΩ12sin⁡θ,z¨=-1-μR13z-μR23z-μ4ROB31-3xcos⁡θ+ysin⁡θROBz.Since P4 does not affect the relative orbital positions of the other two finite masses, its motion is subjected to Kepler’s third law, shown as(17)μ4=ROB3Ω12.Therefore, (16) equals(18)x¨-2y˙=x-1-μR13x+μ-μR23x-1+μ-μ4ROB3x+3μ4xcos2⁡θ+ysin⁡θcos⁡θROB3+3μ4x2cos⁡θ+xysin⁡θROB4,y¨+2x˙=y-1-μR13y-μR23y-μ4ROB3y+3μ4xsin⁡θcos⁡θ+ysin2⁡θROB3+3μ4xycos⁡θ+y2sin⁡θROB4,z¨=-1-μR13z-μR23z-μ4ROB3z+3μ4xzcos⁡θ+zysin⁡θROB4.It is clear thatP4’s gravitational force consists of two parts: the direct constant gravity part due the average gravitational force and the parametric excitation part due to the periodic variation of the position.


### 2.3. Circular Restricted Three-Body Problem

The circular restricted three-body problem (CR3BP) model is a simplification of the R4BP without considering the fourth body’s influence,μ4=0. P1 and P2 still revolve in circular orbits around their barycenter under the influence of their mutual gravitational attractions. P3 moves in the plane defined by the two revolving primaries. The motion of the particle is governed by the following equations with dimensionless time and space units [18]:(19)x-2y˙=∂Ω∂x,y¨+2x˙=∂Ω∂y,z¨=∂Ω∂z,where Ω is the pseudopotential function of the three-body problem and is expressed as(20)Ω=12x2+y2+1-μR1+μR2+12μ1-μ.Apparently, all the governing equations of the motions around the libration points in R4BP model, bicircular model, and the CR3BP model belong to the class of gyroscopic systems [19], which always behave as the massless body reaches an extremum in one coordinate as it passes through zero in the other coordinate, such that they are π/2 out of phase.


## 2.1. Restricted Four-Body Problem

The geometry of the R4BP is shown in Figure1. The finite masses P1, P2, and P4 and the infinite mass body P3 are denoted as m1, m2, m4, and m3, respectively. The motion of P3 is influenced by the attraction of P1, P2, and P4, but the motions of the finite masses P1, P2, and P4 are assumed not to be affected by P3.Figure 1
The configuration of R4BP.In the R4BP,P1 and P2 move about their barycenter perturbed by P4. The masses P1, P2, and P4 all move about the mass center of the entire system. It is assumed that all bodies are modeled as point masses, and the primaries P1, P2, and P4 move in a plane, with the distance from P4 to the P1-P2 system being much greater than the distance from P1 to P2. However, P3 is free to move in all three dimensions.We define the inertial coordinate system,O-ξηζ, with origin at the center of mass of the four-body system, O. The coordinate system, O-XYZ, rotates at an instantaneous rate of Ω1 about the mass center of the four-body system, O. The X-axis points to the barycenter of P1 and P2, B. The synodic coordinate system is B-xyz with origin at B. The x-axis points from the barycenter, B, to P2. B-xyz moves with respect to the coordinate system, O-XYZ, at an instantaneous rate of Ω. θ is the angle between the X-axis and the x-axis, counterclockwise measured. Meanwhile, B-xyz rotates with respect to the inertial coordinate system, O-ξηζ, at an instantaneous rate of Ω2=Ω+Ω1.Nondimensionalization is applied to simplify the form of the equation. The unit length is set as the distance betweenP1 and P2. The unit time is set as T/2π, where T is the period of the P1-P2 system. We introduce μ as the mass parameter of the system; then(1)µ=m2m1+m2,1-µ=m1m1+m2,μ4=m4m1+m2.The massive P1, with mass 1-µ, is located at (-µ,0,0) in B-xyz, and P2, with mass µ, is located at (1-µ,0,0) in B-xyz.The motion ofP3 is governed by the following equations with dimensionless time and space with respect to the coordinate system, O-ξηζ.(2)ρ¨=-1-μR13R1-μR23R2-μ4R43R4,in which R1, R2, and R4 are the position vectors from P1, P2, and P4 to P3, respectively. And ρ is the position vector from the system barycenter, O, to P3, with geometrical relation as (3)ρ=ROB+r.Note that ROB and r are presented in O-XYZ and B-xyz, respectively, and can be expressed as(4)ROB=ROBX,r=xx+yy+zz.Apparently, we have time derivative ofρ with respect to the inertial frame as(5)ρ˙I=R˙OBF+Ω1×ROB+r˙R+Ω2×r,(6)ρ¨I=R¨OBF+2Ω1×R˙OBF+Ω˙1×ROB+Ω1×Ω1×ROB+r¨R+2Ω2×r˙R+Ω˙2×r+Ω2×Ω2×r,in which the superscripts “F” and “R” mean the derivatives are measured with respect to the rotating coordinates O-XYZ and B-xyz. Combining (2) and (6), the motion of P3 in synodic coordinate system B-xyz is obtained as(7)r¨R+2Ω2×r˙R+Ω˙2×r+Ω2×Ω2×r+R¨OBF+2Ω1×R˙OBF+Ω˙1×ROB+Ω1×Ω1×ROB=-1-μR13R1-μR23R2-μ4R43R4,where(8)R1=x+μx+yy+zz,R2=x-1+μx+yy+zz,R4=xx+yy+zz+ROB+ROP4X.Then, we introduce the transformation matrix betweenO-XYZ and B-xyz:(9)XYZ=cos⁡θsin⁡θ0-sin⁡θcos⁡θ0001xyz;R4 can be expressed as(10)R4=ROB+ROP4cos⁡θ+xx+ROB+ROP4sin⁡θ+yy+zz.Thus, we obtain the completed equation of motion of R4BP in B-xyz as(11)x¨-2Ω2y˙-Ω˙2y=∂U∂x-R¨OB-ROBΩ12cos⁡θ-2R˙OBΩ1+ROBΩ˙1sin⁡θ,y¨+2Ω2x˙+Ω˙2x=∂U∂y-R¨OB-ROBΩ12sin⁡θ-2R˙OBΩ1+ROBΩ˙1cos⁡θ,z¨=∂U∂z,where U is the pseudopotential function of the four-body problem and is expressed as(12)U=12Ω22x2+y2+1-μR1+μR2+μ4R4.


## 2.2. Bicircular Model in Restricted Four-Body Problem

The bicircular model is a special case of the R4BP and built on the hypothesis that the motions of three primary bodies are restricted. Specifically,P1 and P2 are revolving in circular orbits around their barycenter and the P1-P2 barycenter, B, is moving in a circular orbit around the center of mass of the system O. Besides that, the two orbital planes are assumed to coincide with each other [17].With the aforementioned hypothesis, we consider(13)Ω2=T2π=1,Ω˙2=0,R¨OB=R˙OB=0,Ω˙1=0,ROB+ROP4≈ROB.Equation (11) is simplified as(14)x¨-2y˙=1-μR13x+μ-μR23x-1+μ-μ4R43ROB+ROP4cos⁡θ+x+ROBΩ12cos⁡θ,y¨+2x˙=-1-μR13y-μR23y-μ4R43ROB+ROP4sin⁡θ-y+ROBΩ12sin⁡θ,z¨=-1-μR13z-μR23z-μ4R43z.With expansion of 1/R4,(15)1R43≈ROB+ROP4-31-3xcos⁡θROB+ROP4-3ysin⁡θROB+ROP4≈ROB-31-3xcos⁡θ+ysin⁡θROB,equation of motion (14) can be further reduced to(16)x¨-2y˙=x-1-μR13x+μ-μR23x-1+μ-μ4ROB3x-μ4ROB2cos⁡θ-3μ4xcos2⁡θ+ysin⁡θcos⁡θROB3+3μ4x2cos⁡θ+xysin⁡θROB4+ROBΩ12cos⁡θy¨+2x˙=y-1-μR13y-μR23y-μ4ROB3y-μ4ROB2sin⁡θ-3μ4xsin⁡θcos⁡θ+ysin2⁡θROB3+3μ4xycos⁡θ+y2sin⁡θROB4+ROBΩ12sin⁡θ,z¨=-1-μR13z-μR23z-μ4ROB31-3xcos⁡θ+ysin⁡θROBz.Since P4 does not affect the relative orbital positions of the other two finite masses, its motion is subjected to Kepler’s third law, shown as(17)μ4=ROB3Ω12.Therefore, (16) equals(18)x¨-2y˙=x-1-μR13x+μ-μR23x-1+μ-μ4ROB3x+3μ4xcos2⁡θ+ysin⁡θcos⁡θROB3+3μ4x2cos⁡θ+xysin⁡θROB4,y¨+2x˙=y-1-μR13y-μR23y-μ4ROB3y+3μ4xsin⁡θcos⁡θ+ysin2⁡θROB3+3μ4xycos⁡θ+y2sin⁡θROB4,z¨=-1-μR13z-μR23z-μ4ROB3z+3μ4xzcos⁡θ+zysin⁡θROB4.It is clear thatP4’s gravitational force consists of two parts: the direct constant gravity part due the average gravitational force and the parametric excitation part due to the periodic variation of the position.


## 2.3. Circular Restricted Three-Body Problem

The circular restricted three-body problem (CR3BP) model is a simplification of the R4BP without considering the fourth body’s influence,μ4=0. P1 and P2 still revolve in circular orbits around their barycenter under the influence of their mutual gravitational attractions. P3 moves in the plane defined by the two revolving primaries. The motion of the particle is governed by the following equations with dimensionless time and space units [18]:(19)x-2y˙=∂Ω∂x,y¨+2x˙=∂Ω∂y,z¨=∂Ω∂z,where Ω is the pseudopotential function of the three-body problem and is expressed as(20)Ω=12x2+y2+1-μR1+μR2+12μ1-μ.Apparently, all the governing equations of the motions around the libration points in R4BP model, bicircular model, and the CR3BP model belong to the class of gyroscopic systems [19], which always behave as the massless body reaches an extremum in one coordinate as it passes through zero in the other coordinate, such that they are π/2 out of phase.


## 3. Approximate Analytical Methodology

The governing equations of R4BP and CR3BP in (11) and (19) are widely used in the applications of orbital design, mostly solved by the traditional perturbation methods. In this section, we propose an approximate analytical methodology to derive the polynomial series relations among different directions of the three coupled equations, with which the three-degree-of-freedom problem may be transferred to simple system with one degree-of-freedom. The proposed polynomial relation method can deal with the equation of motions for both autonomous system and the nonautonomous system, which may describe the overall characteristics of the whole system with general rules.


### 3.1. Polynomial Expansion Method for Autonomous System

First, we describe the polynomial expansion methodology in the CR3BP case. The three second-order nonlinear ordinary differential equations in (19) account for the nonlinear vibrations of the three directions of the conservative system. Since the three equations are coupled by both the gyroscopic terms and nonlinear potential terms, there must exist relations among the three directions. If we can locate the relations, then the three coupled equations can be replaced by one equivalent equation.Before we implement the routine to find the relations, we need to transfer the reference frame to the principle coordinate system with the origin exact on the libration point, by which the (0, 0, 0) is the equilibrium position and the linear potential terms are decoupled. Now, we transfer the three second-order equations on the principle coordinate system into six first-order equations as(21)q˙k=pk,p˙k=fkq1,q2,q3,p1,p2,p3k=1,2,3,where q1, q2, and q3 denote the displacements on the x, y, and z directions, respectively, and p1, p2, and p3 denote the corresponding velocities. According to the invariant manifold method proposed by Shaw and Pierre [13, 20, 21], during periodic motion, the positions and velocities of all the directions are functionally related to a single position-velocity pair of the base direction, which we choose arbitrarily, for example, the displacement and velocity on the x-axis, q1 and p1. The position and velocity along the x-axis can be assumed as u=q1 and v=p1, and the positions and velocities on the y- and z-axis are related to u and v. Since we are treating the nonlinear differential equations, nonlinear relations of Qk and Pk should be considered as(22)qk=Qku,v=a1ku+a2kv+a3kuv+a4ku2+a5kv2+a6ku2v+a7kuv2+a8ku3+a9kv3+⋯,pk=Pku,v=b1ku+b2kv+b3kuv+b4ku2+b5kv2+b6ku2v+b7kuv2+b8ku3+b9kv3+⋯,where Qk and Pk are real functions of u and v. aik and bik are the coefficients to be determined. For the weak nonlinear system, the Taylor series expansion of Qk and Pk is used to obtain the relations by determining the coefficients of the power series.In order to determine relations (22), the time derivatives of Qk and Pk need to be expanded using the chain rule, as(23)Q˙k=∂Qk∂ududt+∂Qk∂vdvdt=∂Qk∂uv+∂Qk∂vf1,P˙k=∂Pk∂ududt+∂Pk∂vdvdt=∂Pk∂uv+∂Pk∂vf1k=2,3.By comparing (21) to (23), we can obtain the following equations:(24)∂Qk∂uv+∂Qk∂vf1=Pk,∂Pk∂uv+∂Pk∂vf1=fk,k=2,3,where any occurrence of qk or pk in fk is systematically replaced by Qk or Pk, respectively. To be detailed, for the libration cases, substituting the Taylor series expansions of Qk and Pk into (24), the following equations are obtained:(25)a1kv+a3kv2+2a4kuv+2a6kuv2+a7kv3+3a8ku2v+⋯+a2k+a3ku+2a5kv+a6ku2+2a7kuv+3a9kv2+⋯f1q1,q2,q3,p1,p2,p3=b1ku+b2kv+b3kuv+b4ku2+b5kv2+b6ku2v+b7kuv2+b8ku3+b9kv3+⋯,b1v+b3v2+2b4uv+2b6uv2+b7v3+3b8u2v+b2+b3u+2b5v+b6u2+2b7uv+3b9v2f1q1,q2,q3,p1,p2,p3=fkq1,q2,q3,p1,p2,p3.Gathering the like terms of the coordinates and velocities provides algebraic equations from which the Taylor series coefficients ofQk and Pk can be determined. Finally, substituting the analytical relations of (22) to the x-direction equation in (21), we transfer the three-degree-of-freedom problem into one-degree-of-freedom problem, which is a function of u. The resulting nonlinear one-degree-of-freedom function of u and its time derivation v can be treated easily by any perturbation methods or numerical methods. The remaining two directions of motion can be expressed as the nonlinear relations of u and v. In the potential applications, the resulting analytical expressions may be as both initial points of computations and constraint conditions for the numerical orbits design.


### 3.2. Iterative Polynomial Expansion Method for Nonautonomous System

By inspection of (11) or (18), we can concluded that the bicircular model in R4BP is nonautonomous system in the sense of parametric resonance. From the view point of perturbation, when the parametric period is a half of the period of the corresponding autonomous system, the parametric resonance may occur. Since the bicircular model is nonconservative and nonautonomous, the methodology used in the last subsection can not be used directly to the system with time-varying parameters. However, we can substitute the time-varying terms by the base displacement-velocity pair (q1,p1), from the mathematical point of view. The iterative idea by Avramov [16] should be considered to fulfill this methodology.Similarly, we transfer the three second-order equations on the principle coordinate system into six first-order equations.q1, q2, and q3 denote the displacements on the x, y, and z directions, respectively, and p1, p2, and p3 denote the corresponding velocities. The first step is to derive the nonlinear relations of the y- and z-axis with x-axis by determining the coefficients of the Taylor series of (22) by neglecting the time-varying parameters in the governing equations. This step can be carried out as exactly described in the last subsection. However, here the coefficients we derived for the first step are not the final values: we need to optimize these values by iterations.With the relations based on the autonomous system, which can be expressed as(26)qk=a1kq1+a2kp1+a3kq1p1+a4kq12+a5kp12+a6kq12p1+a7kq1p12+a8kq13+a9kp13+⋯,pk=b1kq1+b2kp1+b3kq1p1+b4kq12+b5kp12+b6kq12p1+b7kq1p12+b8kq13+b9kp13+⋯,we can transfer the three-dimensional system of xyz (18) into one-dimensional system on x direction,(27)q¨1+ω2q1+F1q1,q˙1,q¨1+a1q1sin⁡2θ+φ1+a2q˙1sin⁡2θ+φ2+a3q12sin⁡θ+φ3+a4q1q˙1sin⁡θ+φ4+a5q˙12sin⁡θ+φ5+⋯=0,where higher-order terms other than two have been neglected.The second step is to transfer the time-varying terms as expressions of (q1,p1). To study the case of parametric resonance, we only consider the period of the system is approximately a half of the parametric pulsating period,(28)θ˙≈2ω,orθ˙≈ω.Hence, the solution to (27) can be presented in the following form:(29)q1=A0+A1cos⁡ωt+B1sin⁡ωt+A2cos⁡2ωt+B2sin⁡2ωt.Substituting (29) back into (27) and using the harmonic balance method, we can obtain a set of five nonlinear algebraic equations with respect to the six parameters (A0, A1, B1, A2, B2, and ω). By introducing the notations(30)z1=sin⁡ωt,z2=cos⁡ωt,the following equations are obtained from (29):(31)q1=A0+A1z2+B1z1+A2z22-z12+2B2z1z2,p1=-A1ωz2+B1ωz1+4ωA2z1z2+2ωB2z22-z12.The solution of (31) with respect to z1 and z2 is presented in the form of polynomial series:(32)z1=α0+α1q1+α2p1+α3q12+α4p12+α5q1p1+⋯,z1=β0+β1q1+β2p1+β3q12+β4p12+β5q1p1+⋯,where α0,α1,…,β0,β1,…, are unknown coefficients. Substituting (31) into (32) and collecting the same terms of q1mp1n and letting the coefficients equal zero yield the following linear algebraic equations:(33)MΑ=R1,MB=R2,where R1=0,0,1,0,0,0T, R2=0,0,0,0,1,0T, A=α0,α1,α2,α3,α4,α5T, B=β0,β1,β2,β3,β4,β5T, and M is the 6×6 matrix, the elements of which are expressed in the six parameters (A0,A1,B1,A2,B2,ω).Solving the two systems of linear algebraic equation (33), polynomial series (32) can be calculated and the following two equations are obtained from these series:(34)sin⁡ωt=Sq1,p1,cos⁡ωt=Cq1,p1.Now, we express the time-varying terms by the functions of (q1,p1). Replacing all the time-varying terms by the functions S and C in the governing equations, we can transfer the parameter resonance nonautonomous three-degree-of-freedom equations into a pseudoautonomous system.The last step is to apply the method described in the last subsection again to this set of pseudoautonomous equations by determining the relations of thex- and y-axis with the x-axis and transfer the three-degree-of-freedom system into that of one degree-of-freedom. The three-step procedure can be iterated further until satisfactory results are obtained.


## 3.1. Polynomial Expansion Method for Autonomous System

First, we describe the polynomial expansion methodology in the CR3BP case. The three second-order nonlinear ordinary differential equations in (19) account for the nonlinear vibrations of the three directions of the conservative system. Since the three equations are coupled by both the gyroscopic terms and nonlinear potential terms, there must exist relations among the three directions. If we can locate the relations, then the three coupled equations can be replaced by one equivalent equation.Before we implement the routine to find the relations, we need to transfer the reference frame to the principle coordinate system with the origin exact on the libration point, by which the (0, 0, 0) is the equilibrium position and the linear potential terms are decoupled. Now, we transfer the three second-order equations on the principle coordinate system into six first-order equations as(21)q˙k=pk,p˙k=fkq1,q2,q3,p1,p2,p3k=1,2,3,where q1, q2, and q3 denote the displacements on the x, y, and z directions, respectively, and p1, p2, and p3 denote the corresponding velocities. According to the invariant manifold method proposed by Shaw and Pierre [13, 20, 21], during periodic motion, the positions and velocities of all the directions are functionally related to a single position-velocity pair of the base direction, which we choose arbitrarily, for example, the displacement and velocity on the x-axis, q1 and p1. The position and velocity along the x-axis can be assumed as u=q1 and v=p1, and the positions and velocities on the y- and z-axis are related to u and v. Since we are treating the nonlinear differential equations, nonlinear relations of Qk and Pk should be considered as(22)qk=Qku,v=a1ku+a2kv+a3kuv+a4ku2+a5kv2+a6ku2v+a7kuv2+a8ku3+a9kv3+⋯,pk=Pku,v=b1ku+b2kv+b3kuv+b4ku2+b5kv2+b6ku2v+b7kuv2+b8ku3+b9kv3+⋯,where Qk and Pk are real functions of u and v. aik and bik are the coefficients to be determined. For the weak nonlinear system, the Taylor series expansion of Qk and Pk is used to obtain the relations by determining the coefficients of the power series.In order to determine relations (22), the time derivatives of Qk and Pk need to be expanded using the chain rule, as(23)Q˙k=∂Qk∂ududt+∂Qk∂vdvdt=∂Qk∂uv+∂Qk∂vf1,P˙k=∂Pk∂ududt+∂Pk∂vdvdt=∂Pk∂uv+∂Pk∂vf1k=2,3.By comparing (21) to (23), we can obtain the following equations:(24)∂Qk∂uv+∂Qk∂vf1=Pk,∂Pk∂uv+∂Pk∂vf1=fk,k=2,3,where any occurrence of qk or pk in fk is systematically replaced by Qk or Pk, respectively. To be detailed, for the libration cases, substituting the Taylor series expansions of Qk and Pk into (24), the following equations are obtained:(25)a1kv+a3kv2+2a4kuv+2a6kuv2+a7kv3+3a8ku2v+⋯+a2k+a3ku+2a5kv+a6ku2+2a7kuv+3a9kv2+⋯f1q1,q2,q3,p1,p2,p3=b1ku+b2kv+b3kuv+b4ku2+b5kv2+b6ku2v+b7kuv2+b8ku3+b9kv3+⋯,b1v+b3v2+2b4uv+2b6uv2+b7v3+3b8u2v+b2+b3u+2b5v+b6u2+2b7uv+3b9v2f1q1,q2,q3,p1,p2,p3=fkq1,q2,q3,p1,p2,p3.Gathering the like terms of the coordinates and velocities provides algebraic equations from which the Taylor series coefficients ofQk and Pk can be determined. Finally, substituting the analytical relations of (22) to the x-direction equation in (21), we transfer the three-degree-of-freedom problem into one-degree-of-freedom problem, which is a function of u. The resulting nonlinear one-degree-of-freedom function of u and its time derivation v can be treated easily by any perturbation methods or numerical methods. The remaining two directions of motion can be expressed as the nonlinear relations of u and v. In the potential applications, the resulting analytical expressions may be as both initial points of computations and constraint conditions for the numerical orbits design.


## 3.2. Iterative Polynomial Expansion Method for Nonautonomous System

By inspection of (11) or (18), we can concluded that the bicircular model in R4BP is nonautonomous system in the sense of parametric resonance. From the view point of perturbation, when the parametric period is a half of the period of the corresponding autonomous system, the parametric resonance may occur. Since the bicircular model is nonconservative and nonautonomous, the methodology used in the last subsection can not be used directly to the system with time-varying parameters. However, we can substitute the time-varying terms by the base displacement-velocity pair (q1,p1), from the mathematical point of view. The iterative idea by Avramov [16] should be considered to fulfill this methodology.Similarly, we transfer the three second-order equations on the principle coordinate system into six first-order equations.q1, q2, and q3 denote the displacements on the x, y, and z directions, respectively, and p1, p2, and p3 denote the corresponding velocities. The first step is to derive the nonlinear relations of the y- and z-axis with x-axis by determining the coefficients of the Taylor series of (22) by neglecting the time-varying parameters in the governing equations. This step can be carried out as exactly described in the last subsection. However, here the coefficients we derived for the first step are not the final values: we need to optimize these values by iterations.With the relations based on the autonomous system, which can be expressed as(26)qk=a1kq1+a2kp1+a3kq1p1+a4kq12+a5kp12+a6kq12p1+a7kq1p12+a8kq13+a9kp13+⋯,pk=b1kq1+b2kp1+b3kq1p1+b4kq12+b5kp12+b6kq12p1+b7kq1p12+b8kq13+b9kp13+⋯,we can transfer the three-dimensional system of xyz (18) into one-dimensional system on x direction,(27)q¨1+ω2q1+F1q1,q˙1,q¨1+a1q1sin⁡2θ+φ1+a2q˙1sin⁡2θ+φ2+a3q12sin⁡θ+φ3+a4q1q˙1sin⁡θ+φ4+a5q˙12sin⁡θ+φ5+⋯=0,where higher-order terms other than two have been neglected.The second step is to transfer the time-varying terms as expressions of (q1,p1). To study the case of parametric resonance, we only consider the period of the system is approximately a half of the parametric pulsating period,(28)θ˙≈2ω,orθ˙≈ω.Hence, the solution to (27) can be presented in the following form:(29)q1=A0+A1cos⁡ωt+B1sin⁡ωt+A2cos⁡2ωt+B2sin⁡2ωt.Substituting (29) back into (27) and using the harmonic balance method, we can obtain a set of five nonlinear algebraic equations with respect to the six parameters (A0, A1, B1, A2, B2, and ω). By introducing the notations(30)z1=sin⁡ωt,z2=cos⁡ωt,the following equations are obtained from (29):(31)q1=A0+A1z2+B1z1+A2z22-z12+2B2z1z2,p1=-A1ωz2+B1ωz1+4ωA2z1z2+2ωB2z22-z12.The solution of (31) with respect to z1 and z2 is presented in the form of polynomial series:(32)z1=α0+α1q1+α2p1+α3q12+α4p12+α5q1p1+⋯,z1=β0+β1q1+β2p1+β3q12+β4p12+β5q1p1+⋯,where α0,α1,…,β0,β1,…, are unknown coefficients. Substituting (31) into (32) and collecting the same terms of q1mp1n and letting the coefficients equal zero yield the following linear algebraic equations:(33)MΑ=R1,MB=R2,where R1=0,0,1,0,0,0T, R2=0,0,0,0,1,0T, A=α0,α1,α2,α3,α4,α5T, B=β0,β1,β2,β3,β4,β5T, and M is the 6×6 matrix, the elements of which are expressed in the six parameters (A0,A1,B1,A2,B2,ω).Solving the two systems of linear algebraic equation (33), polynomial series (32) can be calculated and the following two equations are obtained from these series:(34)sin⁡ωt=Sq1,p1,cos⁡ωt=Cq1,p1.Now, we express the time-varying terms by the functions of (q1,p1). Replacing all the time-varying terms by the functions S and C in the governing equations, we can transfer the parameter resonance nonautonomous three-degree-of-freedom equations into a pseudoautonomous system.The last step is to apply the method described in the last subsection again to this set of pseudoautonomous equations by determining the relations of thex- and y-axis with the x-axis and transfer the three-degree-of-freedom system into that of one degree-of-freedom. The three-step procedure can be iterated further until satisfactory results are obtained.


## 4. Conclusions

In this study, the restricted three-body problem (R3BP) and restricted four-body problem (R4BP) are classified as conservative autonomous system and nonconservative nonautonomous parametric excited system, respectively. From the vibrational point of view, the methodology of polynomial series is proposed to solve for these problems analytically. By introducing the polynomial series relations among thex, y, and z directions, the three-degree-of-freedom coupled equations are transferred into one degree-of-freedom containing the full dynamics of the original autonomous system for the R3BP. On the other hand, the methodology of polynomial series is combined with the iterative approach to treat the nonautonomous system of bicircular R4BP. The pseudoautonomous equation and the final polynomial series relations and one-degree-of-freedom system can be derived iteratively.The methodology proposed in this study may be used in the computations of periodic orbits and further be applied in the numerical orbit design as initial conditions and constraints.



---

*Source: 9747289-2016-10-19.xml*