Simulation of the algorithms and their visualization forthe solutions to the restricted problems of the cosmic dynamics of the fourteen bodies with three rings

– The restricted problem of the fourteen bodies with three rings is considered. The results of the visualization and dynamic investigations with the Mathematica system are given. The equilibrium positions are found with the use of analytical, numerical and graphical possibilities of the system. The stability of the equilibrium positions is then considered. Visualization and animation techniques are used for the observations of the motion processes.


Introduction
The computer algebra methods are the effective means for the search of problems of mathematical and theoretical physics, celestial mechanics, astrodynamics and other natural sciences. With the systems of computer mathematics put into practice, the solutions of many classical natural science problems have been improved. In such a way H. Poincaré [1] formulated the classical restricted three-body problem and provided the method for its solution. The generalization of this problem for n-bodies is given in the papers by Grebenikov and Elmabsout [2,3,4]. In the models like these the gravitational field is created by the bodies which make up regular polygons revolving around the central body. In such a gravitational field the motion of passively gravitating mass is studied. The systems of differential equations which describe such models are nonlinear and that is why it is impossible to integrate them by quadratures. In such cases H. Poincaré recommended to discover the problem of existence of equilibrium points [1]. This problem, as the authors have shown [5,6,7], can be reduced to the solution of rather complicated systems of the nonlinear algebraic (non-differential) equations. It's impossible to find the exact solution to these algebraic systems in the analytical form because of their essential nonlinearity. But when using modern computer mathematical systems (for example,the Mathematica system by Wolfram Research [8] the coordinates of equilibrium points can be found with the arbitrary calculation accuracy. In this paper the restricted problem of the fourteen bodies with three rings is studied. That means that in the arbitrary Euclidean space there are 13 bodies (one central and 12 bodies situated at the apexes of three squares, which have the general center and revolve around the central body with the equal angular velocity) and one body m with a negligibly small mass which does not gravitate with other bodies. It is necessary to study the motion of body m with a negligibly small mass. While solving this problem we use the analytical, numerical and graphical possibilities of the Mathematica system.
We should point out that such restricted problems with the number of bodies from 4 to 12 and the number of rings equal to one or two are studied in papers [5, 6, ?, 9, 10, 11, 12]. The restricted problem with the number equal fourteen and the number of rings equal one or two is studied in papers [13,14,15].

Algorithms for solutions to the restricted problem of the cosmic dynamics
Let us define the next algorithm of solving the considered three-ring restricted problem: The program modules realized in the Mathematica system codes, and the graphs built with the use of the visualization capabilities of the system are given below.
Application of the picture [n_, α_, β_] module: using, for example, the following function picture[4,0.5,0.15] we will obtain the visualization of the three-ring model (   Fig.1 [16], we obtain the analytical condition for the angular velocity of the rotating squares

Geometric visualization of the partial solution represented in
and two analytical conditions for the masses m i (= 1, 2, 3) and distances α, β relationship Conditions (1) 3 By means of the Mathematica we can solve system (1)-(3) as relates to omega, m 2 , m 3 using the function Solve. Using the animational capabilities of the Mathematica system [7] makes it easy to prove that there are the areas of the positive dynamic m 0 .m 1 , and geometric parameters α, β, in which the values ω, m 2 , m 3 will be positive as well. At that we can accept -without the restriction of generality -that m 0 = 1.
In Figs 2 and 3 such areas built with the use of the Plot3D and Manipulate functions [8] are depicted for the defined parameters. Below we can see the graphics of the dependence between the masses, distances and angular velocity.  The general method for finding the equilibrium positions in the restricted problem of the fourteen bodies, the problem visualization, as well as finding the coordinates of the equilibrium position for the cases of one or two rings are presented in [13]. We will use the method and the module sys [n]. In this problem n = 4. As a result, function sys [4] gives the analytical form for the g(x,y) and h(x,y) functions. The intersection of these functions results in the localization of the equilibrium positions on the plane. In Fig.4 all the equilibrium positions are shown for the parameters α = 0.8, β = 0.9, m 0 = 1, m 1 = 0.0001. Bold points are the bodies P k (vertexes of the three squares). Fig.5 presents the enlarged view of the portion of Fig.4 (0.8 ≤ x ≤ 1.2, −0.1 ≤ y ≤ 0.3), where the two of the equilibrium positions are located. The other the equilibrium positions are located symmetrically.    N 1 is a radial equilibrium position as it is located on the line going through the two bodies, whereas S 1 is a non-radial equilibrium position. Likewise using the function FindRoot we can easily find the coordinates of all the equilibrium positions with the arbitrary required accuracy.

2.5
Let us do numerical research of the non-radial equilibrium position S 1 . In papers [11,12] it was proved that in the case of one-and two-ring problems with some accurately defined values of masses, the similar equilibrium positions are stable in the Lyapunov sense. We will study the point S 1 with the following coordinates Let us set the initial conditions and choose for instance the time t (0 < t < 10000). Using the numerical capabilities of the Mathematica system we can solve the Cauchy problem [12,15]: where the functions g and h are defined by function sys [4], with the angular velocity ω obtained from (1). We obtained the results in the form of interpolating functions. We build graphics of these interpolating functions for different values of time t and different initial conditions x 1 , y 1 , x ′ 1 , y ′ 1 . Let us build the graphics of these interpolating functions [12,15]    With the changes in the initial conditions and the defined values of masses, the trajectory does not deviate from the position of equilibrium S 1 . We can propose that the point S 1 is stable or -to be more precise -does not contradict to the stability property. Exact proof for the stability of the point S 1 requires the use of the results of the Kolmogorov-Arnold-Moser theory [17]. Having done the numerical study for the point N 1 similar to that which was done for the point S 1 (with such initial conditions the trajectory moves away from the equilibrium position N 1 ), let us make sure that this point has the non-stability characteristics.