# Numerical-Symbolic Methods for Searching Relative Equilibria in the Restricted Problem of Four Bodies.

1 Introduction

The many-body problem is a famous model of Celestial Mechanics proposed first by Newton, many papers were devoted to its investigation and many important results were obtained (see, for example, [17,23]). Recall that the problem is to compute and predict the motion of two or more massive particles attracting each other according to the Newton law of universal gravitation if their initial positions and velocities are given. In case of two interacting particles a general solution to the problem is known and it enables to specify main types of motion in the two-body problem (see [17]). However, adding one more particle to the system complicates the equations of motion substantially and the problem becomes non-integrable. Note that in 1912 K. Sundman (see [21]) succeeded in finding a general solution of the three-body problem in the form of power series in terms of some parameter. But this solution is so complicated and the series obtained converge so slowly that it is useless for real applications. In case of four and more interacting particles the problem is even more complicated and this stimulates development of different approximate analytical and numerical methods for its investigation (see [4,5]).

Another approach to the many-body problem is based on the search for exact particular solutions of the equations of motion and studying their properties (see [17,23]. The first such solutions in the three-body problem found by Euler (1767) and Lagrange (1772) describe planar motion of the three particles on circular orbits around their common center of mass. These solutions are known as the homographic ones and correspond to the points of equilibria (or relative equilibria) of the particles in the rotating coordinate system [23]. In one case the three particles are located at the vertices of an equilateral triangle (triangular configuration); in the other case they are located at the same straight line and form a collinear configuration. It has been proven that there exist only three collinear and two triangular equilibrium configurations in the three-body problem. However, in case of four and more interacting particles of given masses the number of different equilibrium configurations is not known and the problem of searching for such configurations continues to be a challenge to mathematicians and remains open for further research (see [1,18,20]).

In many particular cases of the many-body problem which are important in the astronomical respect a mass of one body is negligibly small in comparison to the masses of the others. A typical example of great interest is the Earth, the Moon and a spacecraft constituting the three-body system. Due to the smallness of mass of the spacecraft it is quite reasonable to assume that it does not affect the motion of the Earth and the Moon. Consequently, one can formulate a problem of studying the motion of particle [P.sub.3] of negligible mass in the gravitational field created by particles [P.sub.0], [P.sub.1] of masses [m.sub.0], [m.sub.1], respectively, moving on circular orbits about their common center of mass according to the corresponding solution of the two-body problem. Such a model was first proposed by L. Euler and is known as the circular restricted three-body problem (see [22]). As in case of a general three-body problem, there exist five exact particular solutions of the differential equations of motion determining the equilibrium positions of particle [P.sub.3] in the rotating coordinate system; these equilibria are called the points of libration [L.sub.j] (j = 1, 2, ..., 5). The libration points are of great interest for applications and so their stability was a subject of many papers during the past two hundred years. As a result it was proven that three points [L.sub.1], [L.sub.2], [L.sub.3] situated at the line [P.sub.0][P.sub.1] (collinear equilibrium positions) are unstable while the libration points [L.sub.4], [L.sub.5] (triangular equilibrium positions) may be stable if the mass ratio [[mu].sub.1] = [m.sub.1]/[m.sub.0] is sufficiently small (see [12,22]).

The systems of particular interest in Celestial Mechanics and Cosmic Dynamics usually contain more than three bodies. So it makes sense to add the third particle [P.sub.2] of mass [m.sup.2] to the system [P.sub.0][P.sub.1][P.sub.3] and to investigate an influence of its mass on the number and location of equilibrium positions of the particle [P.sub.3]. Obviously, the particle [P.sub.2] may be situated in any of the five equilibrium positions of the corresponding three-body problem which coincide with the libration points [L.sub.1], ..., [L.sub.5] in case of [m.sub.2] =0. In this way we obtain the restricted four-body problem that has been a subject of many papers (see, for example, [2, 6, 7, 8, 9,10,11,13,19]). It should be emphasized that in the framework of this model, motion of the massive particles [P.sub.0], [P.sub.1], [P.sub.2] is given and is determined by the corresponding solutions of the three-body problem, namely, by the Lagrange triangular solutions or by the Euler collinear solutions (see [17]).

The case when particle [P.sub.2] is situated in the vertex [L.sub.4] of the equilateral triangle [P.sub.0][P.sub.1][L.sub.4] was studied in detail in [3,6,9,19]. It was shown that four new equilibrium positions of particle [P.sub.3] arise from the point of libration [L.sub.4] if the second mass parameter [[mu].sub.2] = [m.sub.2]/[m.sub.0] becomes greater than zero. The rest four libration points [L.sub.1], [L.sub.2], [L.sub.3], [L.sub.5] only change their positions depending on the values of parameters [[mu].sub.1], [[mu].sub.2]. Besides, one or two new equilibrium positions may arise inside the triangle [P.sub.0][P.sub.1][P.sub.2]. Location and stability of all these equilibrium positions were investigated in [6,7,8,9,15].

To complete the study of the equilibrium positions in the restricted fourbody problem one needs to consider the case when particle [P.sub.2] is situated in one of the collinear libration points [L.sub.1], [L.sub.2], [L.sub.3] and, therefore, the massive particles [P.sub.0], [P.sub.1], [P.sub.2] are located on the same straight line. It should be emphasized that this problem differs essentially from the case of triangular configuration of the particles [P.sub.0], [P.sub.1], [P.sub.2]. Actually, in the latter case the particles [P.sub.0], [P.sub.1], [P.sub.2] are fixed in the vertices of the equilateral triangle for any values of parameters [[mu].sub.1], [[mu].sub.2], while in the collinear case mutual distances between particles depend on these parameters. Therefore, given the values of [[mu].sub.1], [[mu].sub.2] we have to look for both equilibrium configuration of the massive particles [P.sub.0], [P.sub.1], [P.sub.2] and equilibrium positions of the particle [P.sub.3] as solutions of the corresponding nonlinear algebraic equations. Only afterwards we can analyze the Hamiltonian function in the neighborhood of each equilibrium configuration and investigate the stability of equilibrium positions (see [14,15,16]).

In the present paper we focus on the search of the equilibrium positions of particle [P.sub.3] in the case when three massive particles [P.sub.0], [P.sub.1], [P.sub.2] form a collinear configuration. As the algebraic equations determining such positions are nonlinear we cannot find their general solutions in symbolic form. However, for small values of parameter [[mu].sub.2] we can obtain the corresponding solutions in the form of power series in terms of [[mu].sub.2]. Then we use these solutions to separate different equilibrium solutions and to investigate their dependence on parameter [[mu].sub.2]. using numerical methods. We also estimate possible perturbations of the equilibrium configurations in case of non-zero mass of particle [P.sub.3]. Realization of this task involves very advanced symbolic and numerical calculations which are performed here with the aid of the computer algebra system Wolfram Mathematica [24].

2 General analysis of equilibrium configurations

Let us consider the rotating coordinate system with an origin located at the point [P.sub.0]. In such a system the particles [P.sub.1], [P.sub.2] are immovable and located on the same straight line which may be considered as the Ox axis. Without loss of generality, one can consider also that the dimensionless x-coordinates of particles [P.sub.1], [P.sub.2] are equal to 1 and a, respectively, where the variable a is defined by the equation (see [10,17,22])

1+ [[mu].sub.1] + [[mu].sub.1] (a / [[absolute value of a].sup.3] a - 1 / [[absolute value of a - 1].sup.3] = [[mu].sub.1]/a (1 + a - 1/ [[absolute value of a - 1].sup.3]) + 1 + [[mu].sub.2] / [[absolute value of a].sup.3] = [Kappa], (2.1)

and parameter [kappa] > 0 determines an angular velocity of rotation of the coordinate system. We assume also that parameters [[mu].sub.1], [[mu].sub.2] belong to the intervals 0 < [[mu].sub.1] 1, 0 [less than or equal to] [[mu].sub.1] [less than or equal to] 1 which realize all physically different configurations of the massive particles.

In case of [[mu].sub.1] > 0, [[mu].sub.2] = 0, the roots of equation (2.1) correspond to the collinear libration points [L.sub.1], [L.sub.2], [L.sub.3] in the restricted three-body problem (see [12,22]). Increasing parameter [[mu].sub.2] results in changing the equilibrium position of the particle [P.sub.2] on the Ox axis but the number of real-valued roots of equation (2.1) remains the same and it has exactly one root on the each interval a < 0, 0 < a < 1, a > 1, for any values of parameters [[mu].sub.1] > 0 and [[mu].sub.2] [greater than or equal to] 0. Given [[mu].sub.1] > 0 and [[mu].sub.2] [greater than or equal to] 0, one can easily find each of these roots numerically with the aid of the built-in Mathematica function FindRoot, for example (see [24]).

Equilibrium position of the particle [P.sub.3] of negligible mass corresponding to some given configuration of the particles [P.sub.0], [P.sub.1], [P.sub.2] is determined by the system of two equations (see [16])

[mathematical expression not reproducible]. (22)

Besides of two parameters [[mu].sub.1], [[mu].sub.2], the system (2.2) contains the equilibrium position a of the particle [P.sub.2] which also depends on parameters [[mu].sub.1], [[mu].sub.2] (see (2.1)). Therefore, a general solution of the system (2.2) cannot be found in analytical form and numerical methods should be utilized for its investigation.

2.1 Collinear configurations

One can readily see that the second equation of system (2.2) is satisfied for any x if y = 0. In this case all the particles are located on the Ox axis and form a collinear configuration. The corresponding equilibrium positions of the particle [P.sub.3] are defined by the first equation of system (2.2) which takes the form

f (x) = [kappa]x - x / [[absolute value of x].sup.3] - [[mu].sub.1] + (1 + x - 1 / [[absolute value of x - 1].sup.3]) - [[mu].sub.2] (a / [[absolute value of a].sup.3] + x - a / [[absolute value of x - a].sup.3) = 0. (2.3)

The function f (x) is continuous and finite for x [member of] R, except for the points 0,1,a, where the particles [P.sub.0], [P.sub.1], [P.sub.2] are located. Note that for given values of parameters [[mu].sub.1], [[mu].sub.2] there exist three solutions of equation (2.1) and for the each solution the points 0,1,a divide the real axis Ox into four intervals, for example, (-[infinity], 0), (0, a), (a, 1), and (1, +[infinity]) in case of a [member of] (0,1). The function f(x) takes values of different signs at the ends of each interval and increases monotonically from -[infinity] to +[infinity] inside of the interval. Indeed, for a [member of] (0,1) and x [member of] (-[infinity], 0) the function f (x) takes the form

f (x) = [kappa]x + 1 / [x.sup.2] [[mu].sub.1] (1 - 1 / [(x - 1).sup.2] - [[mu].sub.2] (1 / [a.sup.2] - 1 / [(x - a).sup.2]. (2.4)

Obviously,

[mathematical expression not reproducible]

and the derivative of function (2.4) is

f'(x) = [kappa] - 2/[x.sup.3] - 2[[mu].sub.1] / [(x - 1).sup.3] - 2[[mu].sub.2] / [(x - a).sup.3] > 0.

For a [member of] (0,1) and x [member of] (0, a) we obtain

[mathematical expression not reproducible], (2.5)

and the derivative of function (2.5) is

f,(x) = [kappa] + 2/[x.sup.3] - 2[[mu].sub.1] / [(x - 1).sup.3] - 2[[mu].sub.2] / [(x - a).sup.2] > 0.

Similarly, for a [member of] (0, 1) and x [member of] (a, 1) we have

[mathematical expression not reproducible], (2.6)

and the derivative of function (2.6) is

f'(x) = [kappa] + 2/[x.sup.3] - 2[[mu].sub.1] / [(x - 1).sup.3] - 2[[mu].sub.2] / [(x - a).sup.3] > 0.

At last, for a [member of] (0,1) and x [member of] (1, +[infinity]) we obtain

[mathematical expression not reproducible], (2.7)

and the derivative of function (2.7) is

f'(x) = [kappa] + 2/[x.sup.3] + 2[[mu].sub.1] / [(x - 1).sup.3] - 2[[mu].sub.2] / [(x - a).sup.3] > 0.

Therefore, for any a [member of] (0,1) there exist only one root of equation (2.3) on the each interval (-[infinity], 0), (0, a), (a, 1), (1, +[infinity]) and totally there are four different equilibrium positions [S.sub.1], [S.sub.2], [S.sub.3], [S.sub.4] of the particle [P.sub.3] on the Ox axis (see Figure 1).

Doing similar calculations, we show that in the cases a < 0 and a > 1 equation (2.3) also has four different roots on the intervals (-[infinity], a), (a, 0), (0,1), (1, +[infinity]) and (-[infinity], 0), (0,1), (1,a), (a, +[infinity]), respectively. This proves the following theorem.

Theorem 1. Equations (2.1), (2.2) determine exactly twelve different collinear equilibrium configurations in the circular restricted four-body problem for any values of parameters 0 < [[mu].sub.1] [less than or equal to] 1, 0 < [[mu].sub.2] [[mu].sub.2] 12.2 Triangular configurations

Three massive particles [P.sub.0], [P.sub.1], [P.sub.2] are located on the Ox axis and form a sequence [P.sub.0], [P.sub.1], [P.sub.2] in case of a > 1, [P.sub.0], [P.sub.2], [P.sub.1] for 0 < a < 1, and [P.sub.2], [P.sub.0], [P.sub.1] in case of a < 0. Connecting the outer particles [P.sub.0] and [P.sub.2] for a > 1, for example, with [P.sub.3] in case of y [not equal to] 0 gives a triangle [P.sub.0][P.sub.2][P.sub.3] with particle [P.sub.1] located on its side [P.sub.0][P.sub.2]. In the other two cases we obtain the triangles [P.sub.0][P.sub.1][P.sub.3] and [P.sub.2][P.sub.1][P.sub.3] with particles [P.sub.2] and [P.sub.0] located on the sides [P.sub.0][P.sub.1] and [P.sub.1][P.sub.2], respectively. Therefore, the corresponding configurations of the particles may be called the triangular ones.

Given parameters [[mu].sub.1], [[mu].sub.2], one can find equilibrium positions of particle [P.sub.3] as solutions of the system (2.2), where a multiplier y [not equal to] 0 in the second equation is eliminated and parameter a is determined by equation (2.1). The system does not change if we subtract the second equation of (2.8) multiplied by x from the first one. As a result we rewrite the system (2.2) in the form

[mathematical expression not reproducible], (2.8)

[mathematical expression not reproducible]. (2.9)

In case of [[mu].sub.2] = 0 which corresponds to the restricted three-body problem it follows immediately from equation (2.8) that [(x - 1).sup.2] + [y.sup.2] = 1. As [kappa] = 1 + [[mu].sub.1] in this case (see (2.1)), equation (2.9) gives [x.sup.2] + [y.sup.2] = 1. As a result, we obtain the well-known Lagrange solutions x = 1/2, y = [+ or -][square root of (3/2)] to the restricted three-body problem. For [[mu].sub.2] > 0 these equilibrium points can change their positions in the plane Oxy or some other equilibrium positions can arise. To answer a question on the number of triangular equilibrium positions of particle [P.sub.3] for [[mu].sub.2] > 0 we need to analyze the functions [f.sub.1](x,y), [f.sub.2](x,y) more carefully.

Note that the functions [f.sub.1](x,y), [f.sub.2](x, y) are continuous and finite in the plane Oxy, except for the points (0,0), (1, 0), (a, 0), where the particles [P.sub.0], [P.sub.1], [P.sub.2] are located. As these functions depend on the square of variable y then solutions to the system (2.8), (2.9) form pairs (x,y) and (x, -y) situated symmetrically with respect to the axis Ox in the plane Oxy. Therefore, it is sufficient to analyze the case y > 0 to conclude on the number of triangular equilibrium configurations.

One can readily see that a multiplier of [[mu].sub.1] in (2.8) takes negative and positive values inside and outside the circle [(x - 1).sup.2] + [y.sup.2] = 1, respectively. Similarly, in case of a > 0 a multiplier of [[mu].sub.2] in (2.8) takes negative and positive values inside and outside the circle [(x - a).sup.2] + [y.sup.2] = [a.sup.2]. These two circles have a joint point x = y = 0, where [f.sub.1] (0, 0) = 0, and their centers are located on the Ox axis. The function [f.sub.1] (x, y) is negative in the domain [(x - 1).sup.2] + [y.sup.2] [less than or equal to] 1 and is positive for a > 1 in the domain [(x - a).sup.2] + [y.sup.2] [greater than or equal to] [a.sup.2]. Therefore, the condition [f.sub.1] (x, y) = 0 may be fulfilled for a > 1 only in the domain between the circles [(x - 1).sup.2] + [y.sup.2] = 1 and [(x - a).sup.2] + [y.sup.2] = [a.sup.2]. Moreover, the condition [f.sub.1] (x, y) = 0 determines only one continuous curve in this domain that is symmetrical with respect to the axis Ox, it is shown as a bold solid curve in Figure 2.

Actually, partial derivatives of the function f1(x, y) are given by

[partial derivative][f.sub.1] / [partial derivative]x = 3[[mu].sub.1] (x - 1) / [([(x - 1).sup.2] + [y.sup.2]).sup.5/2] + 3[[mu].sub.2]a (x - 1) / [([(x - a).sup.2] + [y.sup.2]).sup.5/2, (2.10)

[partial derivative][f.sub.1] / [partial derivative]y = 3[[mu].sub.1] y / [([(x - 1).sup.2] + [y.sup.2]).sup.5/2] + 3[[mu].sub.2]ay / [([(x - 1).sup.2] + [y.sup.2]).sup.5/2] (2.11)

Equation (2.11) shows that [partial derivative][f.sub.1] / [partial derivative]y > 0 for y > 0, a > 0 and so for any fixed x the function [f.sub.1](x, y) increases monotonically if variable y grows up. As [f.sub.1](x, y) takes values of different signs at the circles [(x - 1).sup.2] + [y.sup.2] = 1 and [(x - a).sup.2] + [y.sup.2] = [a.sup.2] for any x [member of] [0, 2] there exists only one value of y inside of the interval bounded by the circles for which [f.sub.1] (x,y) = 0. For the points on the axis Ox, where y = 0 and x increases from x = 2 to x = 2a, equation (2.10) gives [partial derivative][f.sub.1] / [partial derivative]x > 0. So function [f.sub.1](x, 0) increases monotonically and takes values of different signs at the boundary points x = 2 and x = 2a. Therefore, there exists only one point [x.sub.0] [member of] (2, 2a) for which [f.sub.1] ([x.sub.0], 0) = 0, and [f.sub.1] (x, 0) > 0 if x [member of] ([x.sub.0], 2a]. As [f.sub.1] (x, 0) < 0 for x [member of] [2, [x.sub.0]), [f.sub.1] (x,y) > 0 on the circle [(x - a).sup.2] + [y.sup.2] = [a.sup.2] and [partial derivative][f.sub.1] / [partial derivative]y > 0 for each value of x [member of] [2, [x.sub.0]] there exists only one value y > 0 in the domain between the Ox axis and the circle [(x - a).sup.2] + [y.sup.2] = [a.sup.2] for which [f.sub.1] (x, y) = 0. Thus, the condition [f.sub.1] (x, y) = 0 defines a unique continuous curve in the domain between the circles [(x - 1).sup.2] + [y.sup.2] = 1 and [(x - a).sup.2] + [y.sup.2] = [a.sup.2]. Although we consider above only positive values of y the function [f.sub.1] (x,y) depends on [y.sup.2] and so the second branch of the curve in the domain y < 0 is obtained by means of mirror reflection of the upper branch with respect to the axis Ox (see Figure 2).

Note that condition [f.sub.2] (x,y) = 0 also defines a unique continuous curve in the plane Oxy (bold dashed curve in Figure 2) and there are only two points of intersection of the curves for a > 1 which correspond to solutions of the system (2.8)-(2.9). To prove this let us use expressions (2.1) for parameter [kappa] and rewrite the function [f.sub.2] (x, y) in the two forms

[mathematical expression not reproducible]. (2.12)

Partial derivatives of the function [f.sub.2](x,y) are given by

[partial derivative][f.sub.2] / [partial derivative]x = 3x / [([x.sup.2] + [y.sup.2]).sup.5/2] + 3[[mu].sub.1] / [([(x - 1).sup.2] + [y.sup.2]).sup.5/2] + [[mu].sub.2] (x - a) / [([(x - a).sup.2] + [y.sup.2]).sup.5/2], (2.13)

[partial derivative][f.sub.2] / [partial derivative]y = 3y (1 / [([x.sup.2] + [y.sup.2]).sup.5/2] + [[mu].sub.1] / [([(x - 1).sup.2] + [y.sup.2]).sup.5/2] + [[mu].sub.2] / [([(x - a).sup.2] + [y.sup.2]).sup.5/2]). (2.14)

Analyzing equation (2.12), we conclude that [f.sub.2] (x, 0) < 0 for x [member of] [0, a] and a > 1. As partial derivative (2.13) is positive for x > a the function [f.sub.2] (x, 0) increases monotonically and becomes positive for some large value of x. It means that there exists only one value [x.sub.1] > a such that [f.sub.2] ([x.sub.1],0) = 0. To prove that [x.sub.1] < [x.sub.0] let us compute the function [f.sub.2] (x, y) at the points belonging to the curve [f.sub.1] (x, y) = 0. Rewriting equation (2.8) in the form

[[mu].sub.1] / [([(x - a).sup.2] + [y.sup.2]).sup.3/2] = [[mu].sub.1] / (1 - 1 / [([(x - 1).sup.2] + [y.sup.2]).sup.3/2] + [[mu].sub.1] / [[absolute value of a].sup.3/2] (2.15)

and substituting (2.15) into (2.12), we obtain

[f.sub.2] (x,y) = 1 / [[absolute value of a].sup.3] - 1 / [([(x - 1).sup.2] + [y.sup.2]).sup.3/2] + [[mu].sub.1]((a - 1) / a (1 / [[absolute value of a- 1].sup.3] - 1 / [([(x - 1).sup.2] + [y.sup.2]).sup.3/2]. (2.16)

Now it is clear that

[f.sub.2]([x.sub.0],0) = 1 / [a.sup.3] - 1 / [x.sup.3.sub.0] + [[mu].sub.1] (a - 1) / a (1 / [[absolute value of a- 1].sup.3] - 1 / [[absolute value of [x.sub.0] - 1].sup.3] > 0

and, therefore, [x.sub.1] < [x.sub.0]. Moreover, it follows from (2.16) that the condition [f.sub.2] (x, y) > 0 is fulfilled for any point (x, y) belonging to the curve [f.sub.1] (x,y) = 0 and x [member of] ([x.sub.2], [x.sub.1]), where coordinate [x.sub.2] corresponds to the point of intersection of the curve [f.sub.1] (x,y) = 0 and the circle [x.sup.2] + [y.sup.2] = [a.sup.2]. Besides, [f.sub.2] (x,y) > 0 for the points (x,y) belonging to the circle [x.sup.2] + [y.sup.2] = [a.sup.2] for x [member of] [0,[x.sub.2]]. As [f.sub.2] (x, 0) < 0 for x [member of] [0,[x.sub.1]) and partial derivative (2.14) is positive for y > 0, the function [f.sub.2] (x, y) increases monotonically if variable y grows up and takes positive values on the circle [x.sup.2] + [y.sup.2] = [a.sup.2] for x [member of] [0, [x.sub.2]) and on the curve [f.sub.1] (x, y) = 0 for x [member of] [[x.sub.2], [x.sub.1]). Therefore, for any x [member of] [0, [x.sub.1]] there is only one value of y for which [f.sub.2] (x, y) = 0 and so the condition [f.sub.2] (x, y) = 0 defines a unique continuous curve in the domain x [greater than or equal to] 0. Note that existence of a unique continuous curve [f.sub.2] (x, y) = 0 for x < 0 is quite obvious because function [f.sub.2] (x, y) is negative at the origin and takes positive values for [x.sup.2] + [y.sup.2] [right arrow] [infinity] while its derivatives (2.13), (2.14) are negative and positive, respectively, for x < 0, y > 0. Thus, the condition [f.sup.2] (x, y) = 0 defines a unique continuous curve in the plane Oxy.

Note that coefficient of [[mu].sub.1] in (2.16) is positive in the domain [(x - 1).sup.2] + [y.sup.2] > 1 while a sum of the first two terms is negative for [x.sup.2] + [y.sup.2] < [a.sup.2] and it tends to -[infinity] at the origin. And if variable x decreases from x = [x.sub.2] to x = 0 the function [f.sub.2] (x,y) on the curve [f.sub.1] (x,y) = 0 (see (2.16)) decreases to zero at the point [S.sub.5] and becomes negative for smaller values of x. It means that the curve [f.sub.2] (x, y) = 0 crosses the curve [f.sub.1] (x, y) = 0 only at the point [S.sub.5] and for smaller values of x it is situated above the curve [f.sub.1] (x, y) = 0. Thus, we can conclude that the curves [f.sub.1] (x, y) = 0 and [f.sub.2] (x, y) = 0 have only one point of intersection [S.sub.5] for y > 0, a > 1 (another point [S.sub.6] is situated in the domain y < 0) and there are only two triangular equilibrium configurations of the particles for a > 1.

Similar analysis is easily repeated in the case 0 < a < 1 when particle [P.sub.2] is located between particles [P.sub.0], [P.sub.1] and the circle [(x - a).sup.2] + [y.sup.2] = [a.sup.2] is situated inside of the circle [(x - 1).sup.2] + [y.sup.2] = 1. Again the functions [f.sub.1] (x,y) = 0 and [f.sub.2] (x, y) = 0 define two different continuous curves which have only two points of intersection [S.sub.5] and [S.sub.6] (see Figure 3).

In case of a < 0, the curve [f.sub.1] (x, y) = 0 is situated in the domain x [greater than or equal to] 0 if [[mu].sub.1] > [[mu].sub.2], or in the domain x [less than or equal to] 0 if [[mu].sub.1] < [[mu].sub.2]. Actually, analysis of equation (2.1) shows that equilibrium position of particle [P.sub.2] is located inside of the interval - 1 < a < 0 if [[mu].sub.1] > [[mu].sub.2] and a < -1 if [[mu].sub.1] < [[mu].sub.2]. In case of [[mu].sub.1] = [[mu].sub.2] particle [P.sub.2] is located at the point a = -1 and the curve [f.sub.1] (x, y) = 0 degenerates into the Oy axis (see (2.8)). The circles [(x - a).sup.2] + [y.sup.2] = [a.sup.2] and [(x - 1).sup.2] + [y.sup.2] = 1 are situated in the domains x < 0 and x > 0, respectively, and have a joint point (0,0), where [f.sub.1] (0,0) = 0. Function [f.sub.1] (x,y) takes positive values on the circle [(x - a).sup.2] + [y.sup.2] = [a.sup.2] and is negative on the circle [(x - 1).sup.2] + [y.sup.2] = 1. Asymptotic value of [f.sub.1] (x,y) for [x.sub.2] + [y.sub.2] [right arrow] [infinity] is obtained from (2.1), (2.8) and in case of [[mu].sub.1] > [[mu].sub.2] it is positive:

[[mu].sub.1] + [[mu].sub.1] 1 / [[absolute value of a].sup.3] = a / a - 1 (-1 +1 / [[absolute value of a].sup.3]) + 1 / [[absolute value of a - 1].sup.3] ([[mu].sub.2] + [[mu].sub.1] / a) > 0 (2.17)

As partial derivative (2.11) is positive in the domain x > 0, y > 0, function [f.sub.1] (x,y) increases monotonically if variable y grows up and for any x [member of] (0, 2] there exists only one value of y > 0 for which [f.sub.1] (x,y) = 0. For x > 2, y = 0, partial derivative (2.10) is also positive and so there exists only one value of x = [x.sub.0] > 2 for which [f.sub.1] ([x.sub.0],0) = 0 and [f.sub.1] (x, 0) < 0 for x [member of] [2,[x.sub.0]). As partial derivative (2.11) is positive in the domain x > 2, y > 0, there exists only one value of y > 0 for which [f.sub.1] (x,y) = 0 for any x [member of] (2,[x.sub.0]]. Thus, condition [f.sub.1] (x,y) = 0 defines a unique continuous curve situated outside the circle [(x - 1).sup.2] + [y.sup.2] = 1 in the domain 0 [less than or equal to] x [less than or equal to] [x.sub.0] in case of [[mu].sub.1] > [[mu].sub.2] (bold solid curve in Figure 4) while in the domain x < 0 function [f.sub.1] (x, y) takes positive values.

In case of [[mu].sub.1] < [[mu].sub.2], asymptotic value (2.17) of [f.sub.1] (x,y) is negative while it takes positive values on the circle [(x - a).sup.2] + [y.sup.2] = [a.sup.2]. Repeating the analysis above we show that in this case condition [f.sub.1] (x, y) = 0 defines a unique continuous curve situated outside the circle [(x - a).sup.2] + [y.sup.2] = [a.sup.2] in the domain x [less than or equal to] 0.

Analysis of the function [f.sub.2] (x, y) for a < 0 is done similarly to the case a > 1 and shows that the condition [f.sub.2] (x, y) = 0 defines a unique continuous curve in the plane Oxy (bold dashed curve in Figure 4) which has only two points of intersection with the curve [f.sub.1] (x, y) = 0. In case of [[mu].sub.1] > [[mu].sub.2], these points [S.sub.5], [S.sub.6] are located in the domain x > 0 while for [[mu].sub.1] < [[mu].sub.2] they are located in the domain x < 0. In case of [[mu].sub.1] = [[mu].sub.2], when the curve [f.sub.1] (x, y) = 0 degenerates into the axis Oy the points [S.sub.5], [S.sub.6] are located at this axis and equilibrium configuration of the particles becomes symmetrical with respect to the Ox and Oy axes.

Summarizing analysis of the functions [f.sub.1] (x,y), [f.sub.2] (x,y), we can formulate the following theorem.

Theorem 2. Equations (2.1), (2.8), (2.9) determine exactly six different triangular equilibrium configurations in the circular restricted four-body problem for any values of parameters 0 < [[mu].sub.1] [less than or equal to] 1, 0 < [[mu].sub.2] [less than or equal to] 1

3 Computation of equilibrium positions

Equilibrium positions of particle [P.sub.3] are defined by the system (2.2), where parameter a is a root of equation (2.1). So computing the equilibria for given values of parameters [[mu].sub.1], [[mu].sub.2] includes two steps. First, one has to solve equation (2.1) and to find equilibrium position of particle [P.sub.2]. Then the found value of a is substituted into system (2.2) and some numerical procedure is applied for searching its solution. Note that numerical algorithms for solving non-linear equations are well known and are implemented in modern software such as the computer algebra system Mathematica [24], for example. However, utilizing the corresponding procedure requires to specify a starting point for the search and in general it is difficult to predict which solution will be found if equation has several roots. Computation shows only that correct result is usually obtained if the starting point is located sufficiently close to the solution.

Recall that equation (2.1) has three roots for any reasonable values of parameters [[mu].sub.1], [[mu].sub.2] but exactly one root is located within each of the fixed intervals a < 0, 0 < a < 1, a > 1. Choosing the starting point within any of these intervals and utilizing the built-in Mathematica function FindRoot (see [24]), we obtain correct equilibrium point located inside of the same interval. These three intervals exist also in case of m2 =0 when both equations (2.1) and (2.3) take the same form given by

1 + [[mu].sub.1] - 1 / [[absolute value of [a.sub.0]].sup.3] - [[mu].sub.1] / [a.sub.0] (1 + [a.sub.0] - 1 / [[absolute value of [a.sub.0] - 1].sup.3] = 0 (3.1)

Note that equation (3.1) defines three libration points [L.sub.1], [L.sub.2], [L.sub.3] in the restricted three-body problem (see [22]) which are easily computed for any [[mu].sub.1]. In the problem under consideration for [[mu].sub.2] = 0 two particles [P.sub.2] and [P.sub.3] are situated either at the same or at different libration points. If locations of these particles are different and mass of particle [P.sub.2] becomes greater than zero ([[mu].sub.2] > 0) then both particles [P.sub.2] and [P.sub.3] change their positions continuously according to equations (2.1) and (2.3), respectively. But if they are located at the same point and mass of particle P2 increases then two new equilibrium positions of particle [P.sub.3] arise for [[mu].sub.2] > 0 while particle [P.sub.2] only changes its position according to equation (2.1). In both cases for [[mu].sub.2] > 0, three points 0,1,a divide the axis Ox into four intervals, where different collinear equilibrium positions [S.sub.1], [S.sub.2], [S.sub.3], [S.sub.4] are located (see subsection 2.1), and the boundaries of these intervals are not fixed but are determined by parameter a. It complicates a choice of the starting point when equilibrium positions of particle [P.sub.3] are computed for different values of parameters [[mu].sub.1], [[mu].sub.2]. To select and investigate some equilibrium point it is very helpful to find first an approximate solution to equation (2.3) which may be used for specifying the starting point.

For small values of parameter [[mu].sub.2] one can look for solutions of equations (2.1) and (2.3) in the form of power series

a = [a.sub.20] + [a.sub.1] [[mu].sub.2] + [a.sup.2] [[mu].sup.2.sub.2] + ..., x = [a.sub.30] + [x.sub.1] [[mu].sub.2] + [x.sub.2] [[mu].sub.2] + ..., (3.2)

where [a.sub.20] and [a.sub.30] are different roots of equation (3.1). To find unknown coefficients [a.sub.k], k =1, 2, ..., we substitute expression (3.2) for a into equation (2.1) and expand both its sides into power series in terms of [[mu].sub.2]. Equating coefficients of [[mu].sup.k.sub.2] in the left- and the right-hand sides of equation we obtain a system of linear equations determining coefficients [a.sub.k], k = 1, 2, .... Solution of this system can be found in a general symbolic form but the corresponding expressions for coefficients [a.sub.k] are quite bulky and so we write out only coefficient [a.sup.1] as an example

[mathematical expression not reproducible]. (3.3)

This algorithm can be easily realized for any given value of parameter [[mu].sub.1] and particle [P.sub.2] situated for [[mu].sub.2] = 0 in one of the points [L.sub.1], [L.sub.2] or [L.sub.3]. For example, if [[mu].sub.1] = 0.1 and for [[mu].sub.2] = 0, particle [P.sub.2] is located at the point of libration [L.sub.1] the solution of equation (3.1) is [a.sub.0] = 0.717513 and expression (3.3) gives [a.sub.1] = -0.548713. Doing calculations up to the third order, we obtain

a = 0.717513 - 0.548713 [[mu].sub.2] + 2.14053 [[mu].sub.2] - 11.1868 [[mu].sup.3.sub.2] +.... (3.4)

To find collinear equilibrium solutions arising from the points of libration [L.sub.2] and [L.sub.3] we substitute (3.4) and expression (3.2) for x into equation (2.3) and expand both its sides into power series in terms of [[mu].sub.2]. Equating coefficients of [[mu].sup.k.sub.2] in the left- and the right-hand sides of equation we obtain a system of linear equations determining coefficients [x.sub.k], k = 1, 2, .... Computing the coefficients [x.sub.k], we obtain two solutions of equation (2.3) in the form

[mathematical expression not reproducible], (3.5)

[mathematical expression not reproducible]. (3.6)

If particle P2 is located at the point of libration [L.sub.1] for [[mu].sub.2] = 0, we cannot use expression (3.2) for representing the collinear solutions arising from this point, because the function f (x) has singularity at the point x = [a.sub.20] = [a.sub.30] (see (2.3)). However, we can use the power series of the form

X = [a.sub.30] + [x.sub.1] [[mu].sup.1/3.sub.2] + [x.sub.2] [[mu].sup.2/3.sub.2] + [x.sub.3] [[mu].sub.2] + [x.sub.4] [[mu].sup.4/3.sub.2] + ... (3.7)

On substituting (3.4) and (3.7) into equation (2.3) and expanding both its sides into power series in terms of [[mu].sub.2], we equate coefficients of [[mu].sup.k/3.sub.2], k = 1, 2, ... in the left- and the right-hand sides of equation. Solving the obtained system of linear equations, we find the coefficients [x.sub.k] and the corresponding solutions of equation (2.3)

[mathematical expression not reproducible]. (3.8)

Choosing some small value of [[mu].sub.2] and expressions (3.5), (3.6), (3.8), we compute collinear equilibrium positions of particle [P.sub.3] which are used then as starting points by the Mathematica function FindRoot applied for numerical solving equation (2.3). The obtained results are again used as the starting points by the function FindRoot, which solves equation (2.3) for larger value of parameter [[mu].sub.2]. Increasing [[mu].sub.2] with sufficiently small step and repeating the calculations, we obtain four numerical solutions of equation (2.3) determining the collinear equilibrium positions [S.sub.1], [S.sub.2], [S.sub.3], [S.sub.4] of particle [P.sub.3] for different 0 [less than or equal to] [[mu].sub.2] [less than or equal to] 1 and [[mu].sub.1] = 0.1, they are depicted in Figure 5 by bold solid curves. The corresponding approximate solutions (3.5), (3.6), (3.8) are depicted in Figure 5 by bold dashed curves together with numerical solution of equation (2.1) shown by longer dashes. One can readily see that for 0 [less than or equal to] [[mu].sub.2] < 0.04 the approximate solutions agree quite good with the numerical ones. Repeating the calculations in the case when particle [P.sub.2] is located near the point of libration [L.sub.2], we obtain numerical solutions of equation (2.3) shown in Figure 6. Obviously, all the calculations can be easily performed for other values of parameter [[mu].sub.1].

Using the method described above we can also compute approximate solutions of the system (2.2) in case of the triangular equilibrium configurations. If particle [P.sub.2] is located near the point [L.sub.1] and [[mu].sub.1] = 0.1, for example, the corresponding solution of (2.2) is given by

[mathematical expression not reproducible]. (3.9)

If particle [P.sub.2] is located near the points [L.sub.2] and [L.sub.3], the corresponding triangular solutions are

[mathematical expression not reproducible], (3.10)

[mathematical expression not reproducible]. (3.11)

Computing approximate solutions (3.9)-(3.11) for some small value of [[mu].sub.2] and using the obtained coordinates x, y as starting points for the function FindRoot, we can find solution of the system (2.2) with any required precision. Then we repeat the calculation for larger value of [[mu].sub.2], using the results from the previous step as the starting point for the function FindRoot and so on. Finally, we obtain numerical solutions of the system (2.2) showing triangular equilibrium positions for different values of 0 < [[mu].sub.2] < 1 and [[mu].sub.1] = 0.1 and three possible locations of particle [P.sub.2] near the libration points [L.sub.1], [L.sub.2], [L.sub.3] (see Figure 7).

Figure 7 demonstrates how triangular configuration of the particles varies when parameter [[mu].sub.2] grows up from 0 to 1. In case of [[mu].sub.2] = [[mu].sub.1] and particle [P.sub.2] situated near the point [L.sub.3], it is an isosceles triangle [P.sub.1][P.sub.2][P.sub.3] with the axis of symmetry Oy. An isosceles triangle [P.sub.0][P.sub.2][P.sub.3] is obtained also in case of [[mu].sub.2] = 1 and particle [P.sub.2] situated near the point [L.sub.2]. In the last case the straight line [P.sub.1][P.sub.3] is an axis of symmetry for any value of [[mu].sub.1]. These results correspond to our intuitive expectation and confirm correctness of the calculations.

4 The case of non-zero mass of [P.sub.3]

In the framework of the restricted three-body problem it is assumed that negligibly small mass of particle [P.sub.3] ([[mu].sub.3] = [m.sub.3]/[m.sub.0] = 0) does not affect the motion of massive particles [P.sub.0], [P.sub.1], [P.sub.2], and so equations (2.1) and (2.2) determining equilibrium positions of particles [P.sub.2] and [P.sub.3] are solved separately. However, if mass of particle [P.sub.3] is taken into account ([[mu].sub.3] > 0) the system of equations determining equilibrium configurations of the particles becomes more complicated and may be written in the form (see [10])

-[kappa] + 1 + [[mu].sub.1] = [[mu].sub.2] (a - 1 / [r.sup.3.sub.12] - 1 / [r.sup.3.sub.2]) + [[mu].sub.3] (x - 1 / [r.sup.3.sub.13] - x / [r.sup.3.sub.3]), (4.1)

0 = [[mu].sub.2] (b/[r.sup.3.sub.12] - b/[r.sup.3.sub.2]) + [[mu].sub.3] (y/[r.sup.3.sub.13] - y/[r.sup.3.sub.3]), (4.2)

- [kappa]a + (1 + [[mu].sub.2])a / [r.sup.3.sub.2] = - [[mu].sub.1] b / [r.sup.3.sub.12] + [[mu].sub.3] (x - a / [r.sup.3.sub.23] - x / [r.sup.3.sub.3]), (4.3)

- [kappa]b + (1 + [[mu].sub.2])b / [r.sup.3.sub.2] = [[mu].sub.1] b / [r.sup.3.sub.12] + [[mu].sub.3] (y - b / [r.sup.3.sub.23] - y / [r.sup.), (4.4)

- [kappa]x + (1 + [[mu].sub.3) x / [r.sup.3.sub.3] = - [[mu].sub.1] (x - 1 / [r.sup.3.sub.13] + 1) - [[mu].sub.2] (x - a / [r.sup.3.sub.23] + a / [r.sup.3.sub.2]), (4.5)

- [kappa]y + (1 + [[mu].sub.3])y / [r.sup.3.sub.3] = -[[mu].sub.1] y / [r.sup.3.sub.13] - [[mu].sub.2] (y - b / [r.sup.3.sub.23] + b / [r.sup.3.sub.2]), (4.6)

where

[r.sub.2] = [square root of ([a.sup.2] + [b.sup.2])], [r.sub.3] = [square root of ([x.sup.2] + [y.sup.2])], [r.sub.12] = [square root of ([(a - 1).sup.2] + [b.sup.2])],

[r.sub.13] = [square root of ([(x - 1).sup.2] + [y.sup.2])], [r.sub.23] = [square root of ([(x - a).sup.2] + [(y - b).sup.2])].

Here two pairs of variables (a, b) and (x, y) correspond to equilibrium coordinates of particles [P.sub.2] and [P.sub.3], respectively. It is assumed also that particles [P.sub.0] and [P.sub.1] are located at the origin and at the point (1,0) on the axis Ox, respectively, and parameter [kappa] takes only positive values.

One can readily see that in case of [[mu].sub.3] = 0, the system (4.1)-(4.6) splits into two independent sub-systems. The first sub-system (4.1)-(4.4) determines three well-known collinear (b = 0) and two triangular (a = 1/2, b = [+ or -][square root of (3/2)]) equilibrium configurations in the three-body problem. The second sub-system (4.5)-(4.6) reduces to the system (2.2) if b = 0 and the particles [P.sub.0], [P.sub.1], [P.sub.2] form collinear equilibrium configuration; this case is analyzed in detail in Section 2. Recall that for each of the three collinear equilibrium positions of particle [P.sub.2] there are four equilibrium positions [S.sub.1], [S.sub.2], [S.sub.3], [S.sub.4] of particle [P.sub.3] at the axis Ox if [[mu].sub.3] = 0 (see Theorem 1). These collinear equilibrium configurations exist in case of [[mu].sub.3] > 0, as well. Indeed, for b = y = 0, equations (4.2), (4.4), (4.6) are satisfied identically. Solving the rest three equations (4.1), (4.3), (4.5) numerically for given parameters [[mu].sub.1], [[mu].sub.2], [[mu].sub.3], we obtain equilibrium positions a, x of particles [P.sub.2], [P.sub.3] and parameter [kappa], the corresponding solutions are shown in Figure 8. Each of the four pairs of curves depicted in the same style in Figure 8 shows the equilibrium coordinates a, x for different values of [[mu].sub.3] when particle [P.sub.2] is located in the neighborhood of the point of libration [L.sub.1] and particle [P.sub.3] of non-zero mass is located in one of the collinear equilibrium points [S.sub.1], [S.sub.2], [S.sub.3], [S.sub.4]. As [P.sub.3] is assumed to be the lightest particle of the system, we consider only small enough values of parameter 0 [less than or equal to] [[mu].sub.3] [less than or equal to] 0.01. This simulation shows that maximum deviations [absolute value of [DELTA]a] [approximately equal to] 0.0537 and [absolute value of [DELTA]a] [approximately equal to] 0.0325 of particle [P.sub.2] from its initial position corresponding to the case [[mu].sub.3] = 0 take place for [[mu].sub.3] = 0.01 if massive particle [P.sub.3] is located at the points [S.sub.2] and [S.sub.1], respectively. Position of particle [P.sub.3] also changes if parameter [[mu].sub.3] grows up, but its deviation from the initial position does not exceed [absolute value of [DELTA]x] [approximately equal to] 0.03. If massive particle [P.sub.3] is located at the point S3 or S4, the deviations [absolute value of [DELTA]a] and [absolute value of [DELTA]x] do not exceed 0.01. Similar results are obtained if particle [P.sup.2] is located in the neighborhood of the points [L.sub.2] and [L.sub.3]. Thus, in case of small mass of the particle [P.sub.3], we observe only small deviations of particles [P.sub.2], [P.sub.3] of the order of [[mu].sub.3] from their positions in case of [[mu].sub.3] = 0; the number and structure of the collinear configurations of four massive particles remain the same as in case of the restricted four-body problem.

Recall that in case of [[mu].sub.3] = 0 there exist two triangular equilibrium configurations when three massive particles are located at the axis Ox and particle [P.sub.3] is located in the neighborhood of the point of libration [L.sub.4] or [L.sub.5] (see Figure 7).

Direct calculation shows that in case of [[mu].sub.3] > 0 and y [not equal to] 0, equations (4.2), (4.4), (4.6) are not satisfied if b = 0. It means that three massive particles [P.sub.0], [P.sub.1], [P.sub.2] cannot be located at the straight line even for small values of mass of particle [P.sub.3] (or [[mu].sub.3] > 0) and triangle configuration is transformed into 4-gon. For given [[mu].sub.1], [[mu].sub.2] and small values of [[mu].sub.3], one can find solution of the system (4.1)-(4.6) in the form of power series in terms of [[mu].sub.3]. For example, in case of [[mu].sub.1] = 0.1, [[mu].sub.2] = 0.02 and particle [P.sub.2] located near the point [L.sub.1], the corresponding solution in linear approximation is

a = 0.707315 + 0.008724[[mu].sub.3], b = 0.074969[[mu].sub.3],

x = 0.482648 - 0.012279[[mu].sub.3], y = 0.800825 + 0.014733[[mu].sub.3]. (4.7)

Using solution (4.7) as starting point for the function FindRoot, we have computed equilibrium coordinates of particles [P.sub.2], [P.sub.3] for different 0 [less than or equal to] [[mu].sub.3] [less than or equal to] 0.01 and found that their deviations from the initial triangular configuration satisfy the following inequalities

[absolute value of [DELTA]a] < 0.000089, [absolute value of [DELTA]b] < 0.00075, [absolute value of [DELTA]x] < 0.00013, [absolute value of [DELTA]y] < 0.00015.

Similar calculation in the case when particle [P.sub.2] is located near the point [L.sub.3], gives

[absolute value of [DELTA]a] < 0.0060, [absolute value of [DELTA]b] < 0.0756, [absolute value of [DELTA]x] < 0.00127, [absolute value of [DELTA]y] < 0.00081.

The simulations above show that the mass of particle [P.sub.3] disturbs equilibrium configuration of the particles computed in the framework of the restricted four-body problem. However, deviations of the particles from their positions corresponding to the case of [[mu].sub.3] = 0 are of the order of [[mu].sub.3] or less. Therefore, the restricted four-body problem may be applicable to modelling such systems as the Sun, the Jupiter ([[mu].sub.1] [approximately equal to] 0.001) and binary asteroids or the artificial satellites of comets and asteroids ([[mu].sub.2,3] [much less than] [[mu].sub.1]).

5 Conclusions

In the present paper we study the equilibrium configurations in the planar circular restricted four-body problem formulated on the basis of the Euler collinear solutions of the three-body problem. Equilibrium positions of particle [P.sub.3] are defined by the system of two non-linear equations (2.2) which contains two mass parameters [[mu].sub.1], [[mu].sub.2] and parameter a describing equilibrium configuration of the massive particles. We have performed a general analysis of the system (2.2) and proved that it defines 12 different collinear configurations and 6 triangular configurations for any reasonable values of the system parameters [[mu].sub.1], [[mu].sub.2]. For small values of parameter [[mu].sub.2] we find the equilibrium solutions in the form of power series in terms of [[mu].sub.2] (see (3.5), (3.6), (3.8)-(3.11)). These solutions are used later for computing the starting point in the numerical procedure for computing the equilibria with required precision. Such combination of symbolic and numerical computations enables to separate different equilibrium solutions and investigate their dependence on parameters, the corresponding results are demonstrated in Figures 5-7.

Note that all relevant numerical and symbolic calculations and visualization of the obtained results are performed with the aid of the computer algebra system Wolfram Mathematica.

Received January 4, 2018; revised June 27, 2018; accepted June 27, 2018

https://doi.org/10.3846/mma.2018.030

References

[1] A. Albouy, H.E. Cabral and A.A. Santos. Some problems of the classical n-body problem. Celestial Mechanics and Dynamical Astronomy, 113(4):369-375, 2012. https://doi.org/10.1007/s10569-012-9431-1.

[2] M. Alvares-Ramirez, J.E.F. Skea and T.J. Stuchi. Nonlinear stability analysis in a equilateral restricted four-body problem. Astrophysics and Space Science, 358(1):3, 2015. https://doi.org/10.1007/s10509-015-2333-4.

[3] R.E. Arenstrof. Central configurations of four bodies with one inferior mass. Celestial Mechanics, 28(1-2):9-15, 1982. https://doi.org/10.1007/BF01230655.

[4] D. Boccaletti and G. Pucacco. Theory of orbits. Volume 1: Integrable systems and non-perturbative methods. 3rd edn. Astronomy and Astrophysics Library. Springer-Verlag, Berlin Heidelberg, 2004. https://doi.org/10.1007/978-3-66203319-7.

[5] V.A. Brumberg. Celestial mechanics: past, present, future. Solar System Research, 47(5):347-358, 2013. https://doi.org/10.1134/S0038094613040011.

[6] D.A. Budzko and A.N. Prokopenya. Symbolic-numerical analysis of equilibrium solutions in a restricted four-body problem. Programming and Computer Software, 36(2):68-74, 2010. https://doi.org/10.1134/S0361768810020039.

[7] D.A. Budzko and A.N. Prokopenya. On the stability of equilibrium positions in the circular restricted four-body problem. In V.P. Gerdt, W. Koepf, E.W. Mayr and E.V. Vorozhtsov(Eds.), Computer Algebra in Scientific Computing, volume 6885 of Lecture Notes in Computer Science, pp. 88-100, Berlin, Heidelberg, 2011. Springer-Verlag. https://doi.org/10.1007/978-3-642-23568-9_8.

[8] D.A. Budzko and A.N. Prokopenya. Stability of equilibrium positions in the spatial circular restricted four-body problem. In V.P. Gerdt, W. Koepf, E.W. Mayr and E.V. Vorozhtsov(Eds.), Computer Algebra in Scientific Computing, volume 7442 of Lecture Notes in Computer Science, pp. 72-83, Berlin, Heidelberg, 2012. Springer-Verlag. https://doi.org/10.1007/978-3-642-32973-9_7.

[9] D.A. Budzko and A.N. Prokopenya. Symbolic-numerical methods for searching equilibrium states in a restricted four-body problem. Programming and Computer Software, 39(2):74-80, 2013. https://doi.org/10.1134/S0361768813020035.

[10] E.A. Grebenikov, E.V. Ikhsanov and A.N. Prokopenya. Numeric-symbolic computations in the study of central configurations in the planar Newtonian four-body problem. In V.G. Ganzha, E.W. Mayr and E.V. Vorozhtsov(Eds.), Computer Algebra in Scientific Computing, volume 4194 of Lecture Notes in Computer Science, pp. 192-204, Berlin, Heidelberg, 2006. Springer, Berlin, Heidelberg. https://doi.org/10.1007/11870814_16.

[11] R. Kozera. Asymptotics for length and trajectory from cumulative chord piecewise-quarties. Fundamenta Informaticae, 61(3):267-283, 2004.

[12] A.P. Markeev. Libration points in celestial mechanics and cosmic dynamics. Nauka, Moscow, 1978. (in Russian)

[13] J.I. Palmore. Collinear relative equilibria of the planar n-body problem. Celestial Mechanics, 28(1):17-24, 1982. https://doi.org/10.1007/BF01230656.

[14] A.N. Prokopenya. Computing the stability boundaries for the Lagrange triangular solutions in the elliptic restricted three-body problem. Mathematical Modelling and Analysis, 11(1):95-104, 2006. https://doi.org/10.1080/13926292.2006.9637305.

[15] A.N. Prokopenya. Hamiltonian normalization in the restricted many-body problem by computer algebra methods. Programming and Computer Software, 38(3):156-166, 2012. https://doi.org/10.1134/S0361768812030048.

[16] A.N. Prokopenya. Symbolic-numerical analysis of the relative equilibria stability in the planar circular restricted four-body problem. In V.P. Gerdt, W. Koepf, W.M. Seiler and E.V. Vorozhtsov(Eds.), Computer Algebra in Scientific Computing, volume 10490 of Lecture Notes in Computer Science, pp. 329-345, Cham, 2017. Springer International Publishing. https://doi.org/10.1007/9783-319-66320-3_24.

[17] A.E. Roy. Orbital motion. 4th edn. Institute of Physics Publishing, Bristol and Philadephia, 2005.

[18] D.G. Saari. On the role and the properties of n-body central configurations. Celestial Mechanics, 21:9-20, 1980. https://doi.org/10.1007/BF01230241.

[19] C. Simo. Relative equilibrium solutions in the four-body problem. Celestial Mechanics, 18(2):165-184, 1978. https://doi.org/10.1007/BF01228714.

[20] S. Smale. Mathematical problems for the next century. The mathematical Intelligencer, 20(2):7-15, 1998. https://doi.org/10.1007/BF03025291.

[21] K.F. Sundman. Memoire sur le probleme des trios corps. Acta Mathematica, 36(1):105-179, 1912.

[22] V.G. Szebehely. Theory of orbits. The restricted problem of three bodies. Academic Press, New York, London, 1967. Translated from Astronomicheskii Zhurnal, vol. 46, no. 2, pp. 459-460

[23] A. Wintner. The analytical foundations of Celestial Mechanics. Princeton University Press, Princeton, New York, 1941.

[24] S. Wolfram. An elementary introduction to the Wolfram Language. 2nd edition. Wolfram Media, Champaign, IL, USA, 2017.

Alexander N. Prokopenya

Warsaw University of Life Sciences--SGGW

Nowoursynowska str. 159, 02-776 Warsaw, Poland

E-mail(corresp.): alexander_prokopenya@sggw.pl

Caption: Figure 1. Collinear equilibrium positions [S.sub.1], [S.sub.2], [S.sub.3], [S.sub.4], [[mu].sub.1] = 0.1, [[mu].sub.2] = 0.04, 0 < a < 1.

Caption: Figure 2. Points of intersections of the curves [f.sub.1] (x, y) = 0 and [f.sub.2] (x, y) = 0 define the triangular equilibrium positions [S.sub.5] and [S.sub.6], [[mu].sub.1] = 0.1, [[mu].sub.2] = 0.04, a > 1.

Caption: Figure 3. Triangular equilibrium positions [S.sub.5] and [S.sub.6], [[mu].sub.1] = 0.1, [[mu].sub.2] = 0.04, 0 < a < 1.

Caption: Figure 4. Triangular equilibrium positions [S.sub.5] and [S.sub.6], [[mu].sub.1] = 0.1, [[mu].sub.2] = 0.04, a < 0.

Caption: Figure 5. Approximate (dashed curves) and numerical (solid curves) solutions for [S.sub.1], [S.sub.2], [S.sub.3], [S.sub.4], [[mu].sub.1] = 0.1, 0 < a < 1.

Caption: Figure 6. Collinear equilibrium positions for 0 [less than or equal to] [[mu].sub.2] [less than or equal to] 1, particle [P.sub.2] is located near the point [L.sub.2], a > 1, [[mu].sub.1] = 0.1.

Caption: Figure 7. Triangular equilibrium configurations for three possible positions of particle [P.sub.2], [[mu].sub.1] = 0.1, 0 [less than or equal to] [[mu].sub.2] [less than or equal to] 1.

Caption: Figure 8. Collinear equilibrium configurations for 0 < [[mu].sub.3] < 0.01, particle [P.sub.2] is located near the point [L.sub.1], 0 < a < 1, [[mu].sub.1] = 0.1, [[mu].sub.2] = 0.02.