Numerical solution for solving procedure for 3D motions near libration points in the Circular Restricted Three Body Problem (CR3BP)

In a recent paper in Astrophysics and Space Science Vol. 364 no. 11 (2019), S. Ershkov & D. Leschenko presented a new solving procedure for Euler-Poisson equations for solving momentum equations of the CR3BP near libration points for uniformly rotating planets having inclined orbits in the solar system with respect to the orbit of the Earth. The system of equations of the CR3BP has been explored with regard to the existence of an analytic way of presentation of the approximated solution in the vicinity of libration points. A new and elegant ansatz has been suggested in their publication whereby, in solving, the momentum equation is reduced to a system two coupled Riccati ODEs. In this paper, we presented a numerical solution of such coupled Riccati ODEs using Mathematica software package.


Introduction.
The equations of motion of the R3BP describe the dynamics of infinitesimal mass m under the action of gravitational forces affected by two celestial bodies of giant masses (in this problem M Sun and m planet, m planet < M Sun), which are rotating around their common centre of masses on Kepler's trajectories.The aforesaid infinitesimal mass m is supposed to be moving (as first approximation) inside the restricted region of space around the planet of mass m planet or inside the so-called Hill sphere (where a p is semimajor axis of the planet).[1] It is worth to note that a large amount of previous and recent works concerning analytical development with respect to these equations exist.We can mention a few of them here.
Motivated by trajectory design applications, analytical work from the 1970s includes Farquhar and Kamel's third-order expansion of orbits in the vicinity of the Earth-Moon L2 libration point using a Poincar´e-Lindstedt method.5During the same time period, Richardson and Cary developed a third-order approximation of quasi-periodic motion near the Sun-Earth L1 and L2 libration points via the method of multiple time scales.More recent semi-analytical work leverages computer algebra tools.G´omez et al. use a Poincar´eLindstedt method to generate high-order quasihalo orbit expansions.Jorba and Masdemont investigate the dynamics around the libration points using center manifold reduction.These local methods offer a thorough view of the dynamics in the libration point vicinity, but they are limited by their regions of convergence.Specialized algebraic manipulators are often required as well, which can be an implementation obstacle.Still, the solutions have been successfully exploited for generating heteroclinic connections.[5] In the 1980s Howell and Pernicka presented a numerical shooting approach for correcting approximate quasi-periodic orbits, though lacking control of orbit parameters.Gomez and Mondelo later developed a scheme for computing twodimensional quasi-periodic tori by using an invariant circle parameterized by Fourier coefficients of a stroboscopic map.Kolemen et al. use a similar approach except with multiple Poincare maps and directly parameterizing the invariant circle using states.1 The current approach uses similar concepts to these methods but with a more general formulation.We use a stroboscopic map, which requires less knowledge of the solution structure than a sequence of Poincare maps, and consider tori of arbitrary dimension.We use a grid of states on the invariant torus allowing us to avoid a convolution that is necessary when using a Fourier coefficient parameterization.In addition, we incorporate a different set of solution and continuation constraints that are broadly applicable.This allows the method to be used unchanged between families of solutions.An alternative approach presented by Schilder et al. computes a torus of a flow directly, and some of the constraints used have been adapted from this approach.The flow approach has been applied to the CR3BP; however, this requires dealing with a torus of dimension one larger than the torus of an associated map.[5] We should especially emphasize the theory of orbits, which was developed in profound work [4] for the case of the Circular Restricted Problem of Three Bodies (CR3BP) (primaries M Sun and m planet are rotating around their common centre of masses on circular orbits).According to [4], equations of motion of the infinitesimal mass m should be presented in the co-rotating frame of a Cartesian coordinate system r  ={x, y, z} in case of CR3BP (at given initial conditions).
The suggested ansatz could be useful from a practical point of view in celestial mechanics (for example, to obtain a regular data of astrometric observations).Indeed, sometimes it is convenient for an observer to consider such the celestial dynamical system at a fixed angle to the plane of mutual circular rotation of the chosen primaries e.g. in the dynamics of planetary rotation in the Solar system.[1] The most common example is Pluto's angle of orbit's inclination with respect to Earth's orbit (which is 17°9'); namely, Pluto's orbit is inclined, or tilted, circa 17.1 degrees from the ecliptic of Earth's orbit.Except Mercury's inclination of 7 degrees, all the other orbits of planets are closer to the orbit of Earth.[1]

Coupled Riccati ODEs.
According to the ansatz for solving Poisson equations, the system of Equations for three problem has the analytical way to present the general solution in regard to the time t: [1] here σ = σ(x, y, z) is some arbitrary (real) function, given by the initial conditions; the real-valued coefficients a(x,y,z, t), b(x,y,z, t) in ( 1) are solutions of the mutual system of two Riccati ordinary differential equations, ODEs in regard to the time t System (2) describes the evolution of function a in dependence on the function b in regard to the time t (and vice versa); we should especially note that the Riccati ODE has no analytical solution in general case.

Numerical solution of equation (2)
As an alternative approach to aforementioned solution as described in [1], we would like to prove that numerical solutions exist for equation (2).
We will obtain the numerical solution using Mathematica software package.
Case 1:     The above shematic diagrams show numerical results of equation ( 2) in three different sets of parameters.

Discussion and conclusion.
As we can see in [1], the system of equations for the Circular Restricted Problem of Three Bodies (CR3BP) (which governs the dynamics of infinitesimal mass m under the action of gravitational forces affected by two mutually rotating celestial bodies of giant masses) has been proved to be very complicated to solve analytically in 3D case.
Indeed, numerical solutions of the mutual system (2) of two Riccati ordinary differential equations in regard to the time t, can be found using Mathematica.
We also note that due to the special character of the solutions of Riccati-type ODEs, there is the possibility for sudden jumping in the magnitude of the solution at some time t₀ (see e.g. the solutions of Riccati-type for Abel ODE of 1-st order).In the physical sense, such jumping of Riccati-type solutions could be associated with the effect of a sudden acceleration/deceleration of the celestial body's velocity at a