Approximation Solution of the Nonlinear Circular Sitnikov Restricted Four–Body Problem

Abstract

In this paper, the approximated periodic solutions of the circular Sitnikov restricted four–body problem (RFBP) were constructed using the Lindstedt–Poincaré method, by removing the secular terms, and compared with numerical solution. It can be observed that, in the numerical as well as approximated solutions patterns, the initial conditions are important. In the sense of a numerical solution, the motion is periodic in a certain interval, but beyond this interval, the motion is not periodic. But, the Lindstedt–Poincaré method constantly gives regular and periodic motion all time. Finally, we observed that the solution obtained by the Lindstedt–Poincaré method gives the true motion of the circular Sitnikov RFBP and the fourth approximate solution has more accuracy than the first, second, and third approximate solutions.

1. Introduction

The dynamical system of the restricted three-body problem (RTBP) has a significant role in celestial mechanics. It has many applications in the astrodynamics and stellar dynamics fields [1,2,3,4]. The RTBP is a special case of three-body problem whereas the restricted four-body (RFBP) problem is a generalization of the RTBP. In the RTBP, both primaries circumambulate around their common center of mass while the infinitesimal body does not have gravitational influence on the primaries bodies.

Several studies have been carried out on the RTBP to analyze infinitesimal body motion [5,6,7,8]. Further considerable work in the frame of perturbed RTBP was undertaken in [9,10] to explore the equilibrium points, linear stability and feature of motion about these points.

One of the major features of the RTBP is that it can be reduced to some simpler models, which also has massive significance in celestial mechanics—including but not limited to Hill’s system [11,12,13], Robe’s model [14,15] and the Sitnikov problem [6,16]. However, the Sitnikov problem is considered the simplest model as it is a reduced model of the three-body problem which can be obtained as a sub-case of the circular or elliptic RTBP where the infinitesimal body oscillates along the perpendicular to the plane of the primaries—which is the Z axis. Moreover, it is considered as the simplest sub-case of the $N-$body problem, and in many cases it can be used as a first approximation for astronomical problems in real situations.

In the Sitnikov problem, the existence of oscillation motion was first proven by [17]. Recently, the effect of the gravitational radius in the sense of chaotic scattering phenomenon and the escape regions were studied under the relativistic effect. Furthermore, inflection points in a quantitative behavior were analyzed by using the basin entropy [18]. The Sitnikov problem is also generalized in the sense that the primary bodies’ motions are homographic and form over a long time following a permissible central configuration [19]. The author has studied the continuation of the global periodic solutions for all eccentricity values in the interval $[0,1)$.

The major motivation of the current study was that the Sitnikov problem is a simple model. Although it is broadly studied in celestial mechanics, it is still an effective model which can be used to explore periodic, symmetric and chaotic motions [20,21]. The perturbation techniques used to find periodic orbits in the Sitnikov problem can be applied to some similar real stellar systems. The aim of this paper was to find an approximated analytical periodic solution for a Sitnikov RFBP using the Lindstedt—Poincaré method by removing the secular terms and comparing it with a numerical solution to verify the importance of this perturbation method.

In this article, we studied the Sitnikov problem extended to four–body problems and found the approximate nonlinear solutions. Moreover, it was a particular case of the RTBP where both primaries had equal masses and were moving around their center of mass in the elliptical or circular orbit. In the elliptical Sitnikov problem, the position of infinitesimal mass in a new analytic way is represented by [16]. Bifurcation analysis and periodic orbits analysis in the problem of the Sitnikov four-body model were carried out by [22]. The effect of radiation pressure on the Sitnikov RFBP was discussed by [23]. Several authors have carried out significant analyses of the Sitnikov three-body, four-body and N-body problems; for example, considerable work has been established in [19,20,21].

This manuscript is organized into the following sections. In Section 1, we describe a brief introduction of the periodic solution of Sitnikov restricted three and four-body problems. Furthermore, the equations of motion and dynamical characteristics of the circular Sitnikov four-body problem are described in Section 2. In Section 3, we obtained the first-, second-, third- and fourth-order approximations with the help of the Lindstedt–Poincaré method. The results of the numerical simulation and a comparison among obtained solutions are investigated in Section 4. Finally, in Section 5, we include the discussion and conclusion of this paper.

2. Equations of Motion of the Proposed Model

It is obvious that an equilateral triangular configuration is a particular solution of the restricted problem of a three- or four-body system. We considered the three primary bodies $m_1$, $m_2$, and $m_3$ with equal mass, i.e., $m_1=m_2=m_3=m=1/3$, which take positions at the vertices of an equilateral triangle of the unit side, where these masses are moving in circular orbits around the center of mass of a system, i.e., the center of the triangle. The equations of motion of the fourth body $m_4$ (infinitesimal body) in the dimensionless rotating coordinate system within the frame of the restricted four-body problem are written as [24]

$$\begin{array}{cc}\hfill \ddot{x}-2\dot{y}=& {\Omega }_x,\hfill \\ \hfill \ddot{y}+2\dot{x}=& {\Omega }_y,\hfill \\ \hfill \ddot{z}=& {\Omega }_z,\hfill \end{array}$$

(1)

where:

$$\Omega (x,y,z)={\displaystyle \frac{(x^2+y^2)}{2}}+m_i\left({\displaystyle \frac{1}{r_1}}+{\displaystyle \frac{1}{r_2}}+{\displaystyle \frac{1}{r_3}}\right),$$

(2)

and $r_i$$(i=1,2,3)$ is given by

$$r_i=\sqrt{{(x-x_i)}^2+{(y-y_i)}^2+z^2},$$

(3)

we also remark that $r_i$ represents the distances from the infinitesimal body to the three primaries $m_i$ which are located at the following points:

$$\begin{array}{cc}\hfill & (x_1,y_1)=\left({\displaystyle \frac{1}{\sqrt{3}}},0\right),\hfill \\ \hfill & (x_2,y_2)=\left({\displaystyle \frac{-1}{2\sqrt{3}}},{\displaystyle \frac{1}{2}}\right),\hfill \\ \hfill & (x_3,y_3)=\left({\displaystyle \frac{-1}{2\sqrt{3}}},{\displaystyle \frac{-1}{2}}\right).\hfill \end{array}$$

(4)

The Sitnikov RFBP is a sub-case of the RFBP which characterizes a dynamical system as follows. Three equal bodies (which are called primary bodies) revolve around their common center of mass where the infinitesimal body moves along a line perpendicular to the orbital plane of the primaries motion [17]. The Sitnikov motions are constructed from the equations of RFBP as a special case, i.e., when $m_1=1/3=m_2=m_3$; $x\left(t\right)=y\left(t\right)=0$, and $r_1=\sqrt{1/3+z^2}=r_2=r_3$. In consideration of Equation (1), we can find equations of motion for the circular Sitnikov RFBP where the infinitesimal body moves along the vertical z axis. Therefore:

$$\ddot{z}={\Omega }_z,$$

(5)

with:

$${\Omega }_z={\displaystyle \frac{-z}{{(1/3+z^2)}^{3/2}}},$$

(6)

where $\Omega$ is the potential function:

$$\Omega \left(z\right)={\displaystyle \frac{1}{\sqrt{1/3+z^2}}},$$

(7)

From Figure 1, we can see that the potential function is symmetric at approximately $z=0$; hence, the potential function represents a harmonic oscillation motion.

In the Sitnikov problem, the motion of the infinitesimal body is described by a nonlinear differential equation given in Equation (5). Moreover, in terms of elliptic functions, only the exact analytical solution can be derived and it is not possible in structures of mathematics. Hence, the alternative possible solutions are the approximated solutions which give the characteristics of the exact solution and are developed with the help of perturbation techniques. The perturbation techniques aim to obtain a solution of nonlinear differential equations with its own limitations. This depends on a parameter which is very small and does not appear in the equation of motion such that Equation (1). The infinitesimal body moves with the Jacobian integral along the z axis as its total energy which is given by

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{2}}{\dot{z}}^2-{\displaystyle \frac{1}{\sqrt{1/3+z^2}}}=E.\end{array}$$

(8)

From (1), we observed that when the potential is symmetric and close to the region $z=0$, then $E=-\sqrt{3}$ that means the motion of the infinitesimal body is bounded in the interval $-\sqrt{3}\le E\le -1$ and the unbounded motions’ regions can be found when $E>-1$. Moreover, the zero velocity of the infinitesimal body are in three cases. When it reaches infinity ($z=\pm \infty$) with zero energy ($E=0$) in the first two cases, then there are two critical points at positive infinity as well as negative infinity. The motion is considered unbounded and the fourth body will stay stationary at infinity; whereas in the third case, the total energy takes the negative value $-\sqrt{3}$, i.e., $E=-\sqrt{3}$ and the fourth body moves to stay at the center of mass, i.e., $z=0$, hence the motion is bounded and the fourth body also remains at the center of mass. It is observed that the from-zero velocity surfaces as well as the phase portrait are shown in Figure 2.

The negative total energy means that the body has no ability to escape the gravitational pull of the central body. The negative total energy also means there is a certain position in which the kinetic energy will approximate towards zero. On the other hand, the negative energy does not mean that it is less than zero. It just implies that the orbiting body needs an amount of energy to be added, whether it comes to stable equilibrium or its own energy is zero [25].

Since the total energy is close to $-\sqrt{3}$, hence it is negative and therefore, the fourth body will be close to the center of mass, thereby $\left|z\right|\ll 1$. Therefore, with the help of Equation (1), it is up to $O\left(z^3\right)$ obtained from:

$$\ddot{z}+{\omega }_r^{2}z-{\omega }_d^2{z}^3=0,$$

(9)

where ${\omega }_r^2=3\sqrt{3}$ and ${\omega }_d^2=27\sqrt{3}/2$.

With the help of the perturbation technique, the solution of the periodic motion of Equation (9) is obtained. The term $-{\omega }_d^2{z}^3$ appears in the equation as the perturbation term.

3. Solution by Lindstedt–Poincaré Method

When regular perturbation approaches fail to obtain periodic solutions in ordinary differential equations, a method is used which was developed by Anders Lindstedt (27 June 1854–16 May 1939) and Jules Henri Poincaré (29 April 1854–17 July 1912). Since $\mid z\mid \ll 1$, we can then introduce a very small artificial parameter to the system which can be replaced by one in the final solution [6]. In this case, we interchange the variable z by $\sigma z$ where $\sigma$ is very small ($\sigma \ll 1$); then, we can write $ϵ={\sigma }^2{\omega }_d^2$. Hence, $ϵ$ is also very small ($ϵ\ll 1$) where $ϵ$ is called the strength of the perturbation. Then, the equation of motion in second-order approximation: Equation (9) can be rewritten as

$$\begin{array}{c}\hfill \ddot{z}+{\omega }_r^{2}z-ϵz^3=0,\end{array}$$

(10)

A new independent variable $s=\omega t$ was introduced where $\omega$ is an unspecified function of $ϵ$, initially as a place-keeping parameter. The coefficient of the new governing equation of the second derivative contains $\omega$. This allows the amplitude and the frequency to interact, which is a property that is observed in nonlinear systems. The function $\omega$ can be chosen in such a way as to eliminate the secular terms [26].

Using $s=\omega t$, Equation (10) becomes:

$$\begin{array}{c}\hfill {\omega }^2{z}^{″}+{\omega }_r^{2}z-ϵz^3=0,\end{array}$$

(11)

Assuming the expansion for $\omega$ as

$$\begin{array}{c}\hfill \omega ={\omega }_r+ϵ{\omega }_1+{ϵ}^2{\omega }_2+{ϵ}^3{\omega }_3+...,\end{array}$$

(12)

where ${\omega }_1,{\omega }_2,...$ are unknown constants. It is similar to straightforward expansion. The variable z can be represented by an expansion having the form:

$$\begin{array}{cc}\hfill z\left(t\right)=& z_0\left(t\right)+ϵz_1\left(t\right)+{ϵ}^2{z}_2\left(t\right)+{ϵ}^3{z}_3\left(t\right)+\cdots \hfill \\ \hfill =& \sum _{i=0}^{\infty }{ϵ}^i{z}_i\left(t\right).\hfill \end{array}$$

(13)

Substituting the values of Equations (12) and (13) into Equation (11) and equating the coefficients of ${ϵ}^0,{ϵ}^1,{ϵ}^2,{ϵ}^3$ and ${ϵ}^4$ to zero, we obtained the following differential equations:

$$\begin{array}{cc}\hfill {\omega }_r^2{\displaystyle \frac{d^2{z}_0}{d{s}^2}}+{\omega }_r^2{z}_0=& 0,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\omega }_r^2{\displaystyle \frac{d^2{z}_1}{d{s}^2}}+{\omega }_r^2{z}_1=& z_0^3+2{\omega }_1{\omega }_r{z}_0,\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill {\omega }_r^2{\displaystyle \frac{d^2{z}_2}{d{s}^2}}+{\omega }_r^2{z}_2=& 3z_0^2{z}_1+2{\omega }_1{\omega }_r\left[z_1-{\displaystyle \frac{z_0^3}{{\omega }_r^2}}\right]+\left(2{\omega }_2{\omega }_r-3{\omega }_1^2\right)z_0,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\omega }_r^2{\displaystyle \frac{d^2{z}_3}{d{s}^2}}+{\omega }_r^2{z}_3=& 3z_0{z}_1^2+3z_0^2{z}_2-3z_1{\omega }_1^2-2z_0{\omega }_1{\omega }_2\hfill \\ \hfill & +{\displaystyle \frac{3z_0^3{\omega }_1^2}{{\omega }_r^2}}-{\displaystyle \frac{2z_0(3z_0{z}_1{\omega }_1+{\omega }_1^3+z_0^2{\omega }_2)}{{\omega }_r}}\hfill \\ \hfill & +2\left[z_2{\omega }_1+z_1{\omega }_2+z_0(3{\omega }_1^3+{\omega }_3)\right]{\omega }_r-4z_0{\omega }_1{\omega }_2{\omega }_r^2,\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\omega }_r^2{\displaystyle \frac{d^2{z}_4}{d{s}^2}}+{\omega }_r^2{z}_4=& z_1^3+6z_0{z}_1{z}_2+3z_0^2{z}_3-3z_2{\omega }_1^2-9z_0{\omega }_1^4-6z_1{\omega }_1{\omega }_2+z_0{\omega }_2^2\hfill \\ \hfill & -6z_0{\omega }_1{\omega }_3-{\displaystyle \frac{4z_0^3{\omega }_1^3}{{\omega }_r^3}}+{\displaystyle \frac{9z_0^2{z}_1{\omega }_1^2}{{\omega }_r^2}}+{\displaystyle \frac{4z_0{\omega }_1^4}{{\omega }_r^2}}+{\displaystyle \frac{6z_0^3{\omega }_1{\omega }_2}{{\omega }_r^2}}\hfill \\ \hfill & -{\displaystyle \frac{6z_0{z}_1^2{\omega }_1}{{\omega }_r}}-{\displaystyle \frac{6z_0^2{z}_2{\omega }_1}{{\omega }_r}}+{\displaystyle \frac{4z_1{\omega }_1^3}{{\omega }_r}}-{\displaystyle \frac{6z_0^2{z}_1{\omega }_2}{{\omega }_r}}-{\displaystyle \frac{2z_0^3{\omega }_3}{{\omega }_r}}\hfill \\ \hfill & +2z_3{\omega }_1{\omega }_r+2z_2{\omega }_2{\omega }_r+12z_0{\omega }_1^2{\omega }_2{\omega }_r+2z_1{\omega }_3{\omega }_r\hfill \\ \hfill & +2z_0{\omega }_4{\omega }_r-4z_0{\omega }_2^2{\omega }_r^2.\hfill \end{array}$$

(18)

For simple calculation, let us assume the initial conditions:

$$\begin{array}{c}\hfill z\left(0\right)=z,\dot{z}\left(0\right)=0\end{array}$$

(19)

Utilizing conditions in Equation (19) with Equation (14), the solution of the first linear second-order initial value problem is given by

$$\begin{array}{c}\hfill z_0\left(s\right)=zcos\left(s\right).\end{array}$$

(20)

After utilizing Equations (15) and (20), then we obtain:

$$\begin{array}{c}\hfill {\displaystyle \frac{d^2{z}_1}{d{s}^2}}+z_1=\left[{\displaystyle \frac{3z^3}{4{\omega }_r^2}}+{\displaystyle \frac{2{\omega }_{1}z}{{\omega }_r}}\right]cos\left(s\right)+{\displaystyle \frac{z^3}{4{\omega }_r^2}}cos\left(3s\right).\end{array}$$

(21)

Now the first term on the right-hand side of Equation (21) will lead to a secular term in the obtained solution; if the coefficient of $cos\left(s\right)$ does not equal zero, then this coefficient must equal zero to avoid an unwanted term (secular term). Thus, we obtain:

$$\begin{array}{c}\hfill {\displaystyle \frac{3z^3}{4{\omega }_r^2}}+{\displaystyle \frac{2{\omega }_{1}z}{{\omega }_r}}=0,\end{array}$$

hence, we obtain:

$$\begin{array}{c}\hfill {\omega }_1={\displaystyle \frac{-3z^2}{8{\omega }_r}}.\end{array}$$

(22)

Therefore, Equation (21) takes the form:

$$\begin{array}{c}\hfill {\displaystyle \frac{d^2{z}_1}{d{s}^2}}+z_1={\displaystyle \frac{z^3}{4{\omega }_r^2}}cos\left(3s\right),\end{array}$$

(23)

and its general solution is:

$$\begin{array}{c}\hfill z_1\left(s\right)=c_{1}cos\left(s\right)+c_{2}sin\left(s\right)-{\displaystyle \frac{z^3}{32{\omega }_r^2}}cos\left(3s\right).\end{array}$$

(24)

Applying the initial conditions $z_1\left(0\right)=0$, ${\dot{z}}_1\left(0\right)=0$ ($z_1\left(0\right)=0$ and $z_0\left(0\right)+z_1\left(0\right)=z$, satisfying our initial conditions) into Equation (24), we obtain:

$$\begin{array}{c}\hfill c_1={\displaystyle \frac{z^3}{32{\omega }_r^2}},c_2=0.\end{array}$$

(25)

Hence, the final solution of Equation (23) is:

$$\begin{array}{c}\hfill z_1\left(s\right)={\displaystyle \frac{z^3}{32{\omega }_r^2}}\left[cos\left(s\right)-cos\left(3s\right)\right].\end{array}$$

(26)

Substituting the values of (20), (22) and (26) into (16), we obtain:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d^2{z}_2}{d{s}^2}}+z_2=& \left[{\displaystyle \frac{21z^5}{128{\omega }_r^4}}+{\displaystyle \frac{2{\omega }_{2}z}{{\omega }_r}}\right]cos\left(s\right)\hfill \\ \hfill & +{\displaystyle \frac{3z^5}{16{\omega }_r^4}}cos\left(3s\right)-{\displaystyle \frac{3z^5}{128{\omega }_r^4}}cos\left(5s\right).\hfill \end{array}$$

(27)

Again, the coefficient of $cos\left(s\right)$ must be equal to zero to eliminate the secular term in $z_2$. This leads to:

$$\begin{array}{c}\hfill {\omega }_2={\displaystyle \frac{-21z^4}{256{\omega }_r^3}}.\end{array}$$

(28)

Therefore, Equation (27) reduces to:

$$\begin{array}{c}\hfill {\displaystyle \frac{d^2{z}_2}{d{s}^2}}+z_2={\displaystyle \frac{3z^5}{16{\omega }_r^4}}cos\left(3s\right)-{\displaystyle \frac{3z^5}{128{\omega }_r^4}}cos\left(5s\right),\end{array}$$

(29)

and its general solution is:

$$\begin{array}{cc}\hfill z_2\left(s\right)=& c_{1}cos\left(s\right)+c_{2}sin\left(s\right)-{\displaystyle \frac{3z^5}{128{\omega }_r^4}}cos\left(3s\right)+{\displaystyle \frac{z^5}{1024{\omega }_r^4}}cos\left(5s\right).\hfill \end{array}$$

(30)

With the help of initial conditions $z_2\left(0\right)=0$, ${\dot{z}}_2\left(0\right)=0$ ($z_2\left(0\right)=0$ and $z_0\left(0\right)+z_1\left(0\right)+z_2\left(0\right)=z$, satisfying our initial conditions) in Equation (30), we obtain:

$$\begin{array}{c}\hfill c_1={\displaystyle \frac{3z^5}{128{\omega }_r^4}}-{\displaystyle \frac{z^5}{1024{\omega }_r^4}},c_2=0.\end{array}$$

(31)

Hence, the final solution of Equation (29) is:

$$\begin{array}{c}\hfill z_2\left(s\right)={\displaystyle \frac{z^5}{1024{\omega }_r^4}}\left[23cos\left(s\right)-24cos\left(3s\right)+cos\left(5s\right)\right].\end{array}$$

(32)

Substituting the values of (20), (22), (24), (28) and (32) into (17), we obtain:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d^2{z}_3}{d{s}^2}}+z_3=& \left[{\displaystyle \frac{531z^7}{1024{\omega }_r^6}}-{\displaystyle \frac{225z^7}{512{\omega }_r^4}}+{\displaystyle \frac{2z{\omega }_3}{{\omega }_r}}\right]cos\left(s\right)\hfill \\ \hfill & +{\displaystyle \frac{297z^7}{2048{\omega }_r^6}}cos\left(3s\right)-{\displaystyle \frac{9z^7}{256{\omega }_r^6}}cos\left(5s\right)+{\displaystyle \frac{3z^7}{2048{\omega }_r^6}}cos\left(7s\right).\hfill \end{array}$$

(33)

Again, we have to put the coefficient of $cos\left(s\right)$ equal to zero to eliminate the secular term in $z_3$. This means that:

$$\begin{array}{c}\hfill {\omega }_3={\displaystyle \frac{-531z^6}{2048{\omega }_r^5}}+{\displaystyle \frac{225z^6}{1024{\omega }_r^3}}.\end{array}$$

(34)

Therefore, Equation (33) reduces to:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d^2{z}_3}{d{s}^2}}+z_3=& {\displaystyle \frac{297z^7}{2048{\omega }_r^6}}cos\left(3s\right)-{\displaystyle \frac{9z^7}{256{\omega }_r^6}}cos\left(5s\right)+{\displaystyle \frac{3z^7}{2048{\omega }_r^6}}cos\left(7s\right),\hfill \end{array}$$

(35)

and its general solution takes the form:

$$\begin{array}{cc}\hfill z_3\left(s\right)=& c_{1}cos\left(s\right)+c_{2}sin\left(s\right)-{\displaystyle \frac{297z^7}{16384{\omega }_r^6}}cos\left(3s\right)\hfill \\ \hfill & +{\displaystyle \frac{3z^7}{2048{\omega }_r^6}}cos\left(5s\right)-{\displaystyle \frac{z^7}{32768{\omega }_r^6}}cos\left(7s\right).\hfill \end{array}$$

(36)

With the initial conditions $z_3\left(0\right)=0$, ${\dot{z}}_3\left(0\right)=0$ ($z_3\left(0\right)=0$ and $z_0\left(0\right)+z_1\left(0\right)+z_2\left(0\right)+z_3\left(0\right)=z$, satisfying our initial conditions) in Equation (36), we obtain:

$$\begin{array}{c}\hfill c_1={\displaystyle \frac{547z^7}{32768{\omega }_r^6}},c_2=0.\end{array}$$

(37)

Hence, the final solution of Equation (33) is:

$$\begin{array}{cc}\hfill z_3\left(s\right)=& {\displaystyle \frac{547z^7}{32768{\omega }_r^6}}cos\left(s\right)-{\displaystyle \frac{297z^7}{16384{\omega }_r^6}}cos\left(3s\right)\hfill \\ \hfill & +{\displaystyle \frac{3z^7}{2048{\omega }_r^6}}cos\left(5s\right)-{\displaystyle \frac{z^7}{32768{\omega }_r^6}}cos\left(7s\right).\hfill \end{array}$$

(38)

Utilizing Equations (20), (22), (24), (28), (32), (34) and (38) with Equation (18), we obtain:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d^2{z}_4}{d{s}^2}}+z_4=& \left[\begin{array}{cc}\hfill & 262144z{\omega }_4{\omega }_r\hfill \\ \hfill & +3z^9(9383-7200{\omega }_r^2)\hfill \end{array}\right]{\displaystyle \frac{cos\left(s\right)}{131072{\omega }_r^8}}+{\displaystyle \frac{z^9(31321-16200{\omega }_r^2)}{131072{\omega }_r^8}}cos\left(3s\right)\hfill \\ \hfill & -{\displaystyle \frac{2649z^9}{65536{\omega }_r^8}}cos\left(5s\right)+{\displaystyle \frac{27z^9}{8192{\omega }_r^8}}cos\left(7s\right)-{\displaystyle \frac{5z^9}{65536{\omega }_r^8}}cos\left(9s\right).\hfill \end{array}$$

(39)

We have to put the coefficient of $cos\left(s\right)$ equal to zero to eliminate the secular term in $z_4$. This means that:

$$\begin{array}{c}\hfill {\omega }_4={\displaystyle \frac{-28149z^8}{262144{\omega }_r^7}}+{\displaystyle \frac{675z^8}{8192{\omega }_r^5}}.\end{array}$$

(40)

Therefore, Equation (39) is reduced to the form:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d^2{z}_4}{d{s}^2}}+z_4=& {\displaystyle \frac{z^9(31321-16200{\omega }_r^2)}{131072{\omega }_r^8}}cos\left(3s\right)-{\displaystyle \frac{2649z^9}{65536{\omega }_r^8}}cos\left(5s\right)\hfill \\ \hfill & +{\displaystyle \frac{27z^9}{8192{\omega }_r^8}}cos\left(7s\right)-{\displaystyle \frac{5z^9}{65536{\omega }_r^8}}cos\left(9s\right).\hfill \end{array}$$

(41)

and its general solution is:

$$\begin{array}{cc}\hfill z_4\left(s\right)=& c_{1}cos\left(s\right)+c_{2}sin\left(s\right)-{\displaystyle \frac{z^9(31321-16200{\omega }_r^2)}{1048576{\omega }_r^8}}cos\left(3s\right)\hfill \\ \hfill & +{\displaystyle \frac{883z^9}{524288{\omega }_r^8}}cos\left(5s\right)-{\displaystyle \frac{9z^9}{131072{\omega }_r^8}}cos\left(7s\right)+{\displaystyle \frac{z^9}{1048576{\omega }_r^8}}cos\left(9s\right).\hfill \end{array}$$

(42)

Using initial conditions $z_4\left(0\right)=0$, ${\dot{z}}_4\left(0\right)=0$ ($z_4\left(0\right)=0$ and $z_0\left(0\right)+z_1\left(0\right)+z_2\left(0\right)+z_3\left(0\right)+z_4\left(0\right)=z$, satisfying our initial conditions) in Equation (42), we obtain:

$$\begin{array}{c}\hfill c_1={\displaystyle \frac{14813z^9}{524288{\omega }_r^8}}-{\displaystyle \frac{2025z^9}{131072{\omega }_r^6}},c_2=0.\end{array}$$

(43)

Hence, the final solution of Equation (39) is:

$$\begin{array}{cc}\hfill z_4\left(s\right)=& \left({\displaystyle \frac{14813z^9}{524288{\omega }_r^8}}-{\displaystyle \frac{2025z^9}{131072{\omega }_r^6}}\right)cos\left(s\right)-{\displaystyle \frac{z^9(31321-16200{\omega }_r^2)}{1048576{\omega }_r^8}}cos\left(3s\right)\hfill \\ \hfill & +{\displaystyle \frac{883z^9}{524288{\omega }_r^8}}cos\left(5s\right)-{\displaystyle \frac{9z^9}{131072{\omega }_r^8}}cos\left(7s\right)+{\displaystyle \frac{z^9}{1048576{\omega }_r^8}}cos\left(9s\right).\hfill \end{array}$$

(44)

Finally, substituting the values of $z_0\left(t\right),z_1\left(t\right),z_2\left(t\right)$, $z_3\left(t\right)$ and $z_4\left(t\right)$ in Equation (13) as well as ${\omega }_1,{\omega }_2,{\omega }_3$ and ${\omega }_4$ in Equation (12), we obtain:

$$\begin{array}{cc}\hfill z=& zcos\left(\omega t\right)+ϵ{\displaystyle \frac{z^3}{32{\omega }_r^2}}\left[cos\left(\omega t\right)-cos\left(3\omega t\right)\right]\hfill \\ \hfill & +{ϵ}^2{\displaystyle \frac{z^5}{1024{\omega }_r^4}}\left[23cos\left(\omega t\right)-24cos\left(3\omega t\right)+cos\left(5\omega t\right)\right]\hfill \\ \hfill & +{ϵ}^3\left[{\displaystyle \frac{547z^7}{32768{\omega }_r^6}}cos\left(\omega t\right)-{\displaystyle \frac{297z^7}{16384{\omega }_r^6}}cos\left(3\omega t\right)\right.\hfill \\ \hfill & \left.+{\displaystyle \frac{3z^7}{2048{\omega }_r^6}}cos\left(5\omega t\right)-{\displaystyle \frac{z^7}{32768{\omega }_r^6}}cos\left(7\omega t\right)\right]\hfill \\ \hfill & +{ϵ}^4\left[\left({\displaystyle \frac{14813z^9}{524288{\omega }_r^8}}-{\displaystyle \frac{2025z^9}{131072{\omega }_r^6}}\right)cos\left(\omega t\right)\right.-{\displaystyle \frac{z^9(31321-16200{\omega }_r^2)}{1048576{\omega }_r^8}}cos\left(3\omega t\right)\hfill \\ \hfill & +{\displaystyle \frac{883z^9}{524288{\omega }_r^8}}cos\left(5\omega t\right)-{\displaystyle \frac{9z^9}{131072{\omega }_r^8}}cos\left(7\omega t\right)\left.+{\displaystyle \frac{z^9}{1048576{\omega }_r^8}}cos\left(9\omega t\right)\right],\hfill \end{array}$$

(45)

with:

$$\begin{array}{cc}\hfill \omega =& {\omega }_r-ϵ{\displaystyle \frac{3z^2}{8{\omega }_r}}-{ϵ}^2{\displaystyle \frac{21z^4}{256{\omega }_r^3}}+{ϵ}^3\left({\displaystyle \frac{-531z^6}{2048{\omega }_r^5}}+{\displaystyle \frac{225z^6}{1024{\omega }_r^3}}\right)\hfill \\ \hfill & +{ϵ}^4\left({\displaystyle \frac{-28149z^8}{262144{\omega }_r^7}}+{\displaystyle \frac{675z^8}{8192{\omega }_r^5}}\right).\hfill \end{array}$$

(46)

Equations (45) and (46) are the desired solutions up to fourth-order approximation of the system, while all terms with order $O\left({ϵ}^5\right)$ and higher are ignored. At the end, the parameter $ϵ$ can be replaced by one for obtaining the final form solution according to the place-keeping parameters method.

4. Numerical Results

A comparison was carried out among the numerical: the first-, second-, third- and the fourth-order approximated solutions in the Sitnikov RFBP. The investigation includes the numerical solution of Equation (5) and the first, second, third and fourth-order approximated solutions of Equation (10) obtained using the Lindstedt–Poincaré method which are given in Equations (45) and (46), respectively.

The comparison of the solution obtained from the first-, second-, third- and fourth-order approximation with a numerical solution obtained from (1) is shown in Figure 3, Figure 4 and Figure 5, respectively. We take three different initial conditions to make the comparison. The infinitesimal body starts its motion with zero velocity in general, i.e., $\dot{z}\left(0\right)=0$ and at different positions ($z\left(0\right)=0.1$, $0.2$, $0.3$).

The investigation of motion of the infinitesimal body was divided into two groups. In a first group, three different solutions were obtained for three different initial conditions, which are shown in Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10. In these figures, the purple, green and red curves refer to the initial condition $z\left(0\right)=0.1$, $z\left(0\right)=0.2$ and $z\left(0\right)=0.3$, respectively. However, in a second group, three different solutions were obtained for the above given initial conditions. This group includes Figure 3, Figure 4 and Figure 5, in which the green, blue and red curves indicate the numerical solution (NA), third-order approximated (TA) and fourth-order approximations (FA) of the Lindstedt–Poincaré method, respectively, in these figures.

In Figure 10, we see that the motion of the infinitesimal body is periodic, and its amplitude decreases when the infinitesimal body begins moving closer to the center of mass. Moreover, in numerical simulation, the behavior of the solution is changed by the different initial conditions. Furthermore, the numerical solution is not bounded (we mean that the motion is not periodic) to the infinitesimal body which begins its motion from a large distance to the center of mass of the primaries. Perhaps, the solution has periodic behavior with regular periodicity for the proper selection of initial conditions during a certain time period.

From Figure 3, Figure 4 and Figure 5, we impose that the infinitesimal (fourth-) body starts moving from three different position with respect the center of mass of primaries. Motion is periodic in the interval $0\le t\le 102$ and its amplitude decreases when its motion starts nearer to the position because of the common center of masses of the primaries, which agree with [6]. On the other hand, when $0\le t\le 102$, the numerically obtained motion is then periodic, however, when $t>102$, it escapes from the path of the orbit.

5. Conclusions

In this paper, the periodic solutions of the Sitnikov restricted four-body problem (RFBP) were constructed by using the Lindstedt–Poincaré method by removing the unbounded terms, which are called the secular. These solutions are compared with numerical solution to verify the importance of the used perturbation method (Lindstedt–Poincaré method). One of the major significant aims of removing secular terms is determining the conditions which enforce the motion to be periodic or establish the periodicity conditions of motion. Under these conditions, we can find approximated periodic solutions in closed form. It was observed that in the numerical and the obtained periodic solutions’ patterns, the initial conditions are very important. The motion was numerically examined and compared with the approximated solution obtained by the Lindstedt–Poincaré method.

We remark that the motion is periodic in the time interval $0\le t\le 102$; however, beyond this interval, the motion is not periodic in the numerical sense. However, the Lindstedt–Poincaré method gives regular and periodic motion for a time of $t>0$. In addition, we note that the motion obtained by the first-, second-, third- and fourth-order approximate solutions is regular and periodic when the test particle starts its motion closer to the center of the mass. On the other hand, the test particle starts its motion far from the center of mass and the solutions obtained by numerical simulation may not be periodic for a long time—whereas all the first- to fourth-order approximated solutions may become regular and periodic motions. We demonstrated that the obtained solution by the Lindstedt–Poincaré method gives the true motion of the circular Sitnikov RFBP and the fourth-order approximate solution has more accuracy than the first-, second- and third-order approximate solutions.

In the future study, we aimed to present more analysis of the dynamical motion of the $N-$body problem, particularly the circular and elliptical Sitnikov motion. In this context, we will use some power tools to find periodic solutions such as the method of parameter variation and averaging or the multiple scales technique, which is considered a generalization of the Lindstedt–Poincaré method.