Survey of Eight Modern Methods of Hamiltonian Mechanics

: Here we describe eight new methods, arisen in the last 60 years, to study solutions of a Hamiltonian system with n degrees of freedom. The ﬁrst six of them are intended for systems with small parameters or without them. The methods allow to ﬁnd families of periodic solutions and families of invariant n -dimensional tori by means of analytic computation near a stationary solution, near a periodic solution and near an invariant torus, using the corresponding normal form of a Hamiltonian. Then we can continue the founded families by means of numerical computation. In a Hamiltonian system without parameters, only periodic solutions and invariant n -dimensional tori form one-parameter families. The last two methods are intended for systems with not small parameters, which do not depend on time. They allow computing sets of parameters, which guarantee the stability of some solutions for linear (method seven) and nonlinear (method eight) systems. We do not consider chaotic behaviors, but only regular ones.


Introduction
The following methods, arisen in the last 60 years, are considered here.

1.
A normal form method that allows one to study regular perturbations near a stationary solution [1] (Ch. I), near a periodic solution [1] (Ch. II) [2,3], near the invariant torus [1] (Ch. II) and near families of such solutions [1] (Ch. VII, VIII), as well as bifurcations of periodic solutions and invariant tori and their stability. 2.
The method of truncated systems obtained with the help of Newton polyhedra, which allows the study of singular perturbations. For the theory and three applications, see [4] (Ch. IV). Other applications: Beletskiy's equation on satellite oscillations [5,6], problems of periodic flyby of the Moon and planets. 3.
The method of generating families of periodic solutions (regular and singular). Generating families are the limits of families of periodic solutions as the perturbing parameters tend to zero. The solutions of the generating families consist of certain parts of the solutions to the limit problem. If the limit problem is integrable, then the generating families are found analytically. Applications: the restricted three-body problem, where the limit problem is the two-body problem and the generating families are one-parameter [1] (Ch. III-V) [7][8][9]; Hill's problem [10], where the limit problem is an intermediate Henon problem and each generating family consists of one solution [11,12]. This approach can be applied to families of invariant tori as well. 4.
Methods of numerical computation of families of periodic solutions and of families of invariant tori. 5.
The method of generalized problems admitting in Celestial Mechanics bodies with negative masses [13]. Such problems have unified complete families of periodic solutions, which facilitates their calculation. Example: generalized Hill's problem [13].

6.
Calculation of the network of families of periodic solutions and families of invariant tori as a "skeleton" of a part of the phase space. Poincaré [14] wrote about the benefits of such "skeletons". Examples: the Hill problem [13] and partly the restricted three-body problem [9,15]. 7.
Method of computation of a set of stability of stationary solutions of a linear multiparameter Hamiltonian system combines modern techniques of elimination theory and power geometry [4]. It can be applied in the case when the Hamiltonian function depends on parameters in a polynomial way and gives the description of the boundary of the set of stability in the parameter space [16], and breaking it into cells in which nonlinear terms cannot impact the type of stability [17]. Examples: one gyroscopic problem with three-dimensional space of parameters [16], a double pendulum with a following force. 8.
Application of the q-analog of the classical subresultant for the characteristic polynomial of the matrix of a linear multi-parameter Hamiltonian system allows one to find resonant manifolds [18] and invariant coordinate subspace of the normal form of a Hamiltonian system. Resonant manifolds together with normal form in the vicinity of a stationary solution provide the method of dividing the set of stability into cells where formal stability takes place. The invariant coordinate subspaces allow reducing the phase flow of the initial system to a subspace of essentially less dimension.
There are many more works on these methods. The authors and their collaborators contributed to the development and application of these eight methods. These methods are discussed below in the order shown in Sections 2-9.
Theorem 1 ([1] (Ch. I)). There is a canonical formal transformation (4), which reduces the Hamiltonian (3) to the normal form g(x, y) = ∑ g pq x p y q , where the series g(x, y) contains only resonant terms with p − q, λ = 0 and the quadratic part g 2 (x, y) of g(x, y) has its normal form (so that the matrix of the linear part of the system is the Hamiltonian analog of the Jordan normal form).
For the real original system (1), the coefficients g pq of the complex normal form (5) satisfy special realness relations, and under the standard canonical linear change of coordinates x, y → X, Y, the system with Hamiltonian (5) goes over into the real system. There are several ways to calculate the coefficients g pq of the normal form (5). The simplest is described in the book [19] by Zhuravlev, Petrov, Shunderyuk. The resonant normal form of the autonomous Hamiltonian system near the stationary solution, which takes into account only the eigenvalues of the matrix A of its linear part and without restrictions on this matrix A, was introduced in 1972. Previous variants of normal form were proposed by Birkhoff [20], Cherry [21], Gustavson [22] and others. They assumed some restrictions on the matrix of linear parts of the system. Later, a slightly simpler superresonant normal form was introduced, which took into account the Jordan cells of the normal form of the matrix A [23]. However, these additional simplifications did not allow an additional decrease in the number of degrees of freedom. For a more general approach to the normal form, see [24].

Theorem 2.
There is a formal canonical substitution (7) that reduces the Hamiltonian function (6) to the normal form where the power series h pq (µ) are nonzero only when p − q, λ = 0.
We point out some features of this theorem. First, the parameter does not change under a normalizing transformation. Second, the µ-dependence of the normal form lies in the fact that for resonant indices p, q the coefficient h pq depends on µ in an arbitrary fashion. Third, the eigenvalues λ j do not depend on µ and are computed for the value µ = 0.
The theory of a resonant normal form near a stationary solution is described in detail in the book [1] (Chapter I).
Usually the normalizing transformation diverges in the whole neighborhood of the stationary point. However, near the point there are such analytical families of stationary points, of periodic solutions and of n-dimensional invariant tori, which adjoin the point. Next we describe these families in coordinates of the normal form and small parameters.

Families of Periodic Solutions
All non-zero imaginary eigenvalues λ j , we divide into such blocks that all eigenvalues λ j from the same block are pairwise commensurable, and they are not commensurable with λ k from other blocks. A block can consist of one eigenvalue.
Let one such block correspond to a set of indices I = {i 1 , i 2 , . . . , i m }. Then I = {i 1 , ..., i m } is a set of increasing natural indices i : l < i n. Here 1 m n − l. Consider the coordinate subspace Then in the coordinate subspace K I , a family of periodic solutions satisfies the following system [25] ∂g ∂y j = λ j x j a, ∂g ∂x j = λ j y j a, j = 1, . . . , l, and j ∈ I.
From its solutions, we have to exclude the family of stationary points (8), (9). In (10) a is a free parameter. We can exclude it and obtain the system (8) and In total, we obtain 2(l + m) − 1 equations for 2(l + m) unknowns, which describes a one-parameter family of periodic solutions.
We consider the set A in the cartesian coordinates ρ = (ρ 1 , . . . , ρ n ). In the generic case, each coordinate subspace (with respect to ρ) contains one one-dimensional (with respect to ρ) component of the set A that does not lie in a smaller coordinate subspace. Consequently, the set A consists of 2 n − 1 such components; for each d n there are exactly n!/[d!(n − d)!] of these components situated in d-dimensional (with respect to ρ) coordinate subspaces. In particular, there is one component, situated outside the coordinate subspaces. It is the formal family of the tori. Let where T is a symmetric matrix. In the generic case, det T = 0 and the system of equations in (11) has a one-dimensional solution, System (11) can be written in the form ∂h λ j ∂ρ j = ∂h λ 1 ∂ρ 1 , j = 2, . . . , n.
The formal family A 0 will be analytical if eigenvalues λ satisfy the following condition on small divisors [25].

Condition ω
Let ω k = min| p, λ | over p, It is a very weak numerical restriction on the eigenvalues λ.

Stability
The stationary point ξ = η = 0 can be stable if all eigenvalues λ are pure imaginary.
has no integer solutions p with ||p|| k.
This condition means that there are no resonances up to and including the order k. If it is satisfied, then in the normal form (5) where g l (r) are homogeneous polynomials from r j = x j y j , j = 1, . . . , n, of degree l, andg (k) is a series from x, y starting with powers above k and [k/2] is an integer part of number k/2. Thus it is possible to obtain a Hamiltonian of the form (12) with partial normalization only up to order k, wheng (k) contains not only resonance terms.
In particular, under the condition A n 2 we have Under the condition A n 4 , we have µ jk r j r k +g (5) (x, y).
of the system (1), where f l (ξ, η) is the homogeneous form of degree l. In other words, and f l (ξ, η) does not return to zero at any ξ, η except the point ξ = η = 0.
Stability is possible only if λ = 0. If the condition A n 2 is satisfied, then all λ j are different and non-zero. In this case, the complex coordinates x, y are related to the real coordinates X, Y by the canonical substitution With complex conjugation x j = −iy j ,ȳ j = −ix j , j = 1, . . . , n, and the Hamiltonian function g(x, y) goes into itself, i.e., into (5): Then in real coordinates R j 0, α j is real, Theorem 3 ([26]). If the condition A n 2 is satisfied and the numbers α 1 , . . . , α n are of the same sign, then the fixed point ξ = η = 0 is stable according to Lyapunov.
Here, the role of the integral f is played by the Hamiltonian γ itself, for it is an integral, the notation (13) has the form (12) with k = 2 and the form γ By formal we will mean power series, about the convergence of which nothing is known.

Definition 2 ([27]).
A stationary point (2) of a real Hamiltonian system (1) is formally stable if there exists a formal real sign-defined integral (15) of the system (1), i.e., the formal identity (16) is satisfied and the homogeneous form f l is zero only at ξ = η = 0.
Formal stability means that the departure of solutions from the fixed point, if anything, is very slow: slower than any finite degree t. Since r j r k = − 1 4 R j R k , then under the condition A n 4 the entry (14) takes the form µ jk R j R k +g (5) .
Hence, all the coefficients of µ jk are real. Let K ⊂ R n be a linear hull of integers q satisfying the equation α, q = 0, and Q = {q 0, q = 0} ⊂ R n is a non-negative orthant without origin.

Theorem 4 ([29]
). If the condition A n 4 is met and in (19) n ∑ j,k=1 then the point ξ = η = 0 is formally stable in the sense of the Definition 2 Here, the normal form of the Hamiltonian (5) is used to construct the formal integral.
According to (17) in real coordinates, the normal form (12) is where the homogeneous polynomials h l = i 2 l g l (R) are real. The following generalization of Theorem 4 is proved verbatim as well.
Theorem 5. If the condition A n k is satisfied and in the normal form (21) then the point ξ = η = 0 is formally stable in the sense of the Definition 3.

Condition M2 ([28] (Ch. 8, § 3)). System of equations
µ jk q j q k = 0 has no solution q ∈ Q, i.e., q 0, q = 0. Under conditions A n 4 and M2, the conditions of Theorem 4 are fulfilled and there is formal stability. However, the condition M2 is easier to check than the (20) condition.
If n = 2, the M2 takes the form: system of two equations α 1 q 1 + α 2 q 2 = 0, µ 20 q 2 1 + 2µ 11 q 1 q 2 + µ 02 q 2 2 = 0 has no solution q 1 , q 2 0 with q 1 + q 2 = 0. However, the solutions of the first equation have the form q 1 = − α 2 α 1 q 2 . For them, q 1 , q 2 > 0 only when α 1 α 2 < 0, i.e., the first equation has no solutions with q 1 , q 2 > 0 under the condition of Theorem 3 α 1 α 2 > 0. Substituting them into the second equation and reducing by q 2 2 /α 2 1 , we obtain the condition which is called the Arnold-Moser condition. Under this condition there is not only formal stability, but also Lyapunov stability, because there are one-parameter families of two-dimensional invariant tori with similar sets of frequencies that lock the origin of coordinates. However, J. Moser [30] in 1968 and V.I. Arnold [31] in 1963 made mistakes in proving this fact. At the end of the article [25] there is a criticism of the first proof of Theorem 7 in [30].
Most works on stability use conditions such as condition M2, where the numbertheoretic character of frequencies α j is not taken into account. However, the structure of the normal form depends on them. For example, if the equation α, q = 0 has no solutions in integer q = 0, then the condition A n ∞ is satisfied and the normal form of the Hamiltonian (5) is g(r). Then any r j is a formal integral and the fixed point is formally stable. In particular, at n = 2 this is satisfied if the ratio α 1 /α 2 is an irrational number. Example 1. In the book [28] (Ch. 7), the stability of the libration points of the planar circular restricted three-body problem is studied. There n = 2, the frequencies ω 1 = α 1 , ω 2 = −α 2 with 1 ω 1 > ω 2 > 0 satisfy the equation where µ is the ratio of the masses of the two bodies and the only parameter of the problem (0 µ 1).
In this case, the stability is studied for It is shown in § 4 that according to (4.7) in the normal form Let us show that at these values the frequencies of ω 1 and ω 2 are incommensurable, i.e., there is formal stability.

However, on the interval
Therefore, it follows from equality (27) that Assume z = x/y, i.e., x = zy. Here z is the ratio of the squares of the frequencies. According to (26) we obtain y = 1 z + 1 Substituting this and x = zy in (28), we obtain Consequently, z satisfies the quadratic equation .
Given (28), we see that both values of z are irrational. Consequently, the ratio of frequencies √ z is also irrational.

Example 2.
In the book [28] (Ch. 8), the stability of libration points of a spatial circular restricted three-body problem is studied. There n = 3, the frequencies ω 1 and ω 2 are the same as in Example 1, and ω 3 = 1. In § 3 on p. 136, the formal stability theorem is formulated for all values of µ such that 0 < 27µ(1 − µ) < 1, except where there is double resonance. Let us show that in this problem, the double resonance is impossible. Indeed, in the case of double resonance, the frequencies ω 1 and ω 2 are commensurate with each other and commensurate with unity. Let where p, q, r, s-are integers, According to (26) That is, Let us put k = qr, l = pr, m = qs.
Then the Equation (30) takes the form As we know, all solutions to the Equation (32) in integer non-negative numbers have the form where κ is a non-negative integer. According to (29) and (31) l < k. Therefore, the Equations (33) will apply when κ > 2, and when κ = 0, κ = 1 and κ = 2 we put By direct verification, we make sure that when 0 κ < 3, the Equations (31) and (34) are impossible for integers. When κ > 3, the Equations (31) and (33) follow Therefore, The numbers κ − 1, κ, κ + 1 have no common factor, and the numbers κ 2 − 1 and κ 2 + 1 have no common factor other than 2. Therefore, the ratio (35) cannot be an integer.

Local Coordinates
Let a real Hamiltonian system with n + 1 degrees of freedom have a real 2π-periodic solution M and the Hamiltonian function is analytic in some neighborhood of it.
There exists a complex formal invertible 2π-periodic in ψ and ϕ canonical transformation of coordinates in the form of Poisson series which reduces the Hamiltonian γ into normal form where x, y ∈ C n , 0 p, q ∈ Z n , l 0 and m are integers, and all terms of the second sum are resonant, i.e., The normal form preserves small parameters and linear automorphisms of the initial system.

Theorem 7.
There exists a canonical transformation with rational β j , which reduces the normal form of the Hamiltonian (39), (40) to an autonomous power series where in the first sum all nonzero γ j are irrational numbers and in the second sum 0 p, q ∈ Z n , 0 l ∈ Z, h pql = const ∈ C and present only resonant terms with where γ = (γ 1 , . . . , γ n ) and γ 1 = . . . = γ l = 0.
A similar theorem is in [3]. Variable s is now a formal integral of the following systeṁ If the initial Hamiltonian γ is real for real coordinates ξ, ψ, η, ρ, then in Theorem 6 variables x, y are complex but variables ψ, ρ and ϕ, r are real. Here according to [1] (Chs. I, II) variables x, y are connected with real variables X = (X 1 , . . . , X n ), Y = (Y 1 , . . . , Y n ) by the linear standard transformation

Families of Periodic Solutions
Let all imaginary numbers among eigenvalues λ be λ 1 , . . . , λ m , i.e., λ j = iα j , j = 1, . . . , m. Let all rational numbers α j have eigenvalues λ 1 , . . . , λ l . Then the family of periodic solutions going through the solution M satisfies the following system [25] ∂g where a is a free parameter. Excluding it, we obtain the system After transformation (41) we obtain γ 1 = . . . = γ l = 0 and Equation (45) takes the form The subsystem of Equation (46) defines the set of all periodic solutions of the subsystem (43). Equation (44)φ = ∂h/∂s gives dependence of ϕ from t for each of these solutions.

Families of (n + 1)-Dimensional Irreducible Invariant Tori
Such a family goes through the periodic solution M only if all eigenvalues λ are pure imaginary, i.e., λ = 0, and equation p 0 + p, λ = 0 has no integral solutions p 0 ∈ Z, p ∈ Z n . Then in normal form g = g(ρ, r), and our family is defined by the system where a is a free parameter. Excluding it, we obtain the system That formal family is analytic if eigenvalues λ satisfy the following condition on small divisors [25].

Stability
The periodic solution M is stable only if all eigenvalues λ have λ = 0.
Here there is a notion of Lyapunov stability, but at n = 1 the conditions for its existence coincide with those for formal stability.

Definition 4. Periodic solution
orbitally formally stable if there exists such a power real series on ξ, η, ρ almost periodic on ψ which may diverge, but is a formal sign-defined integral of the system (49).
In other words, all the coefficients of a power series ∂γ ∂ψ must convert to zero and the homogeneous in ξ, η, Recall that a function f (ψ) is periodic if it has a single frequency, conditionally (or quasi) periodic if it has a finite number of frequencies, and almost periodic if it has a countable number of frequencies. In our case, there will be quasi-periodic functions F pql (ψ).
Definition 4 is similar to Definition 2, but one can also define formal orbital stability similar to Definition 3.
Condition B n k . For all integer p with ||p|| |p 1 | + . . . + |p n | k the scalar products p, α are not integers, i.e., the comparison p, α ≡ 0 (mod 1) has no solutions with such p. 3]). Under the condition B n 2 , there exists a complex formal reversible 2π-periodic on ψ and φ canonical coordinate transformation which brings the Hamiltonian γ to the normal form α j x j y j + ∑ g pqlm x p y q r l e imφ , (51) where x, y ∈ C n , 0 p, q ∈ Z n , l 0 and m are integers, all terms of the second sum have order in x, y, √ r above two and resonant, that is, Let us put r j = x j y j , j = 1, . . . , n; r = (r 1 , . . . , r n ).

Corollary 1.
If the condition B n 4 is satisfied, then the normal form (51), (52) has the form where δ = const ∈ C n , ε = const ∈ C.
). The canonical transformation leads the normal form of the Hamiltonian (51) to an autonomous power series h(u, v, s) = s + ∑ h pqlm u p v q s l , corresponding to the second sum in (51).
Note that the returns from the variables u, v, s to the original variables are given by formal power series on ξ, η, ρ with quasi-periodic coefficients on ψ. Let us call the Hamiltonian (56) a reduced normal form.
The variable s is now the formal integral of the systeṁ The orbital stability problem of the periodic solution M has now been reduced to the stability problem of the fixed point u = v = 0, s = 0 in the system (57).

Corollary 2.
If the condition B n 4 is satisfied, then according to (53) and (55) the reduced normal form (56) is

Real Case
If the original Hamiltonian γ is real under the real variables ξ, ψ, η, ρ, then in Theorem 8 the variables x, y are complex and the variables ψ, ρ and φ, r are real.
If the condition B n 2 is satisfied, then according to [1] (Chapters I and II) the complex variables x, y are related to the real variables X, Y by the formulae The complex variables x j , y j and their conjugate variablesx j ,ȳ j are related by the relations With complex conjugation, the Hamiltonian (51) is preserved: g(x, φ, y, r) = g(x, φ, y, r).
Indeed, iα j x j y j =īα jxjȳj = iα j x j y j , and we can show that Note that according to (60) ir j = −ix jȳj = (−i) 3 x j y j = ix j y j = ir j , j = 1, . . . , n.
According to (59) (5) . (62) All quantities here are real. All integer vectors q that satisfy the comparison α, q ≡ 0 (mod 1), form in R n the lattice L. Let M be its linear hull and Q = {q 0, q = 0} is a non-negative orthant in R n with no origin. µ jk q j q k + 2 α, q ∆, q − ε α, q 2 = 0 for all q ∈ M ∩ Q, then the periodic solution (48)

Now (58) takes the form
have an invariant torus T k of dimension k. Using the normal form, we can study the solutions of the system (63) in a neighborhood of the torus T k for any k n [25]. Here we confine ourselves to the most important case k = n. We call a torus T n regular if in its neighborhood there are local coordinates ρ = (ρ 1 , . . . , ρ n ), ψ = (ψ 1 , . . . , ψ n ) that have the following properties.

1.
The coordinates ρ and ψ are canonically conjugate and are analytic functions of ξ and η.

3.
The torus T n is specified by the equations ρ = 0. 4.
On T n the system (63) induces the systeṁ Then, in the neighborhood of the regular torus T n the system (63) takes the forṁ In the neighborhood under consideration, the Hamiltonian function g is 2π-periodic in each ψ j and is expanded in a convergent series where the integer vectors l are non-negative and the coefficients g l are analytic 2π-periodic functions. By hypothesis, on the torus T n , the system (65) takes the form (64), that is, where ω = (ω 1 , . . . , ω n ) is its frequency basis. Now we try to simplify the system (65) by means of a formal canonical local coordinate transformation ρ, ψ → r, ϕ, where ρ = r + · · · , ψ = ϕ + · · · , and we have not written out terms of higher degrees in r. As a result of such a transformation, let g(ρ, ψ) = h(r, ϕ). We expand h(r, ϕ) in a Poisson series We call the Hamiltonian function h a normal form if in (68) only those coefficients h lm are nonzero for which m, ω = 0. Theω-condition: where the limit is taken over the integers m such that m, ω = 0. i.e., h = h(r).The formal family of our tori satisfies the system [25] ∂h ∂r j = −ω j a, j = 1, . . . , n, where a is a free parameter, i.e., ∂h ω j ∂r j − ∂h ω 1 ∂r 1 = 0, j = 2, . . . , n.
Let the MacLaurin series for h have the form where T is a symmetric square matrix. In our case, ∂h/∂ϕ j = 0, and so the system (70) is Taking account of (71), we obtain ω − Tr + · · · = ωa that is, In the generic case, det T = 0, and so the system has a unique one-parametric solution r = r(a − 1) = T −1 ω(1 − a) + · · · that is, the set A is a one-parameter family of n-dimensional regular tori with frequency basis ωa. The formal family (58) is analytic, if frequencies ω satisfy the following condition on small divisors.
Condition ω n Let δ k = min| ω, p | for p ∈ Z n , ||p|| 2 k . Then Thus, we have proved the following result.

Corollary 3.
In a generic Hamiltonian system with n degrees of freedom, an invariant regular torus T n , whose frequency basis ω is non-resonant and satisfies the ω-condition, lies on an analytic one-parameter family of n-dimensional tori, filled by conditionally-periodic solutions with frequency basis aω. The Hamiltonian function changes monotonically in this family. Torus T n is not bifurcating, if det T = 0.

Stability
Non-resonant torus is formally stable.

The Truncated Systems Method
If an equation (or a system of equations) contains a linear part, then sometimes it can be reduced to normal form. However, if the equation does not have a linear part, then the question arises: what should be considered as the first approximation of the equation (or the system)? The answer to it is given by the method of truncated equations, which makes it possible to write out several first approximations and for each indicate the region in the space of variables and parameters where it dominates.

Restricted Three Bodies Problem (RTBP)
Let two bodies P 1 and P 2 with masses 1 − µ and µ, respectively, revolve around their common center of mass with a period of 2π. The plane circular restricted three-body problem studies the plane motion of a body P 3 of infinitesimal mass under the action of the Newtonian attraction of bodies P 1 and P 2 . In a rotating (synodic) coordinate system, the problem is described by an autonomous Hamiltonian system with two degrees of freedom and one parameter µ. It was introduced by Euler in 1772 [35]. The Hamilton function has the form [1] Here body P 1 = {x, y : x 1 = x 2 = 0} and body P 2 = {x, y : x 1 = 1, x 2 = 0}, where x = (x 1 , x 2 ), y = (y 1 , y 2 ). Consider small values of the mass ratio µ 0. For µ = 0 the problem becomes the problem of two bodies P 1 and P 3 . However, here it is necessary to remove from the phase space the points corresponding to the collisions of bodies P 2 and P 3 . Collision points split the solutions of the problem of two bodies P 1 and P 3 into parts. For small µ > 0 near the body P 2 there is a singular perturbation of the case µ = 0.
In order to find all the first approximations of the restricted three-body problem, it is necessary to introduce local coordinates ξ 1 = x 1 − 1, ξ 2 = x 2 , η 1 = y 1 , η 2 = y 2 − 1 near the body P 2 and expand the Hamiltonian function in a power series in these coordinates. After expanding 1/ (ξ 1 + 1) 2 + ξ 2 2 in a Maclaurin series, Hamilton's function (73) takes the form where f is a convergent power series that does not contain terms of order less than three. Let The set S of these points (p, q, r) consists of the points   Figure 1 depicts a polyhedron Γ for series (74) in p, q, r, coordinates, which is a semi-infinite trihedral prism with an oblique base. It has four faces and six edges. Let us consider them.
It describes Hill's problem [10], found in 1878, which is non-integrable. A power transformationξ reduces the corresponding Hamiltonian system to the Hamiltonian system with the Hamiltonian function of the form (75), where ξ i , η i , µ must be replaced byξ iηi , 1, respectively.
Face Γ 2 , which is obtained from the function h at µ = 0. It describes the problem of two bodies P 1 and P 3 , which is integrable.
Consider the edges. Of the six edges, one is improper. It passes through the point (0, 2, 0) parallel to the vector (1, 0, 0). On three edges, q = 0, that is, for them the truncated Hamiltonian functions do not depend on η 1 , η 2 , and the solutions of the corresponding truncated Hamiltonian systems have ξ 1 , ξ 2 = const, which is not interesting. Two edges remain.
Edge Γ It describes the problem of two bodies P 2 and P 3 . Power transformation (76) reduces it into a Hamilton system with a Hamilton function of the form (77), where ξ i , η i , µ is replaced byξ i ,η i , 1, respectively.
Therefore, very close to the body P 2 the first approximation of the original restricted problem with the Hamiltonian function (74) is the problem of two bodies P 2 and P 3 with Hamiltonian (77), just close is the Hill problem with Hamiltonian (75), further from the body P 2 is the intermediate problem, and far from the body P 2 is the problem of two bodies P 1 and P 3 . Near the body P 2 , the periodic solutions of the restricted problem are perturbations of both periodic solutions of all the above four first approximations and the results of gluing the hyperbolic orbits of the two-body problem P 2 , P 3 with arcsolutions of either the two-body problem P 1 , P 3 , or an intermediate problem. In [38][39][40][41][42], periodic solutions of the intermediate problem were used as generators for finding periodic quasi-satellite orbits of the restricted problem.
Therefore, the restricted three-body problem was stated by L. Euler in 1772, one of its first approximations was found by G. Hill in 1878 as a result of long work, another of its first approximations was found by M. Hénon in 1969 again as a result of non-trivial work. However, using the polyhedron, these approximations and the others can be found without difficulty.

Truncated Algebraic Systems
Consider now the set of polynomials x, y, µ), . . . , f m (x, y, µ). (78) Each f j (x, y, µ) has its own support S j ⊂ R 2n+s and all accompanying objects: Newton's polyhedron Γ j its generalized faces Γ (dj) . Moreover, for every non-empty intersection there corresponds to a set of truncationŝ which is the first approximation of the set (78), for (log |x|, log |y|, log |µ|) → ∞ near the intersection (79) and is called the truncation of the set (78). Consider now the system of equations corresponding to the set (78). System (81) corresponds to all the objects indicated for the set (78), as well as the truncated systems of equationŝ each of which corresponds to one set of truncations (80). Each truncated system (82) is the first approximation of the complete system (81).

Analytical Computation of Local Families
In Section 2, we obtained systems of equations of the form (78), (81), the roots of which are families of stationary points, of periodic solutions and of invariant tori, when coordinates tend to zero. The method, described in Section 3.3, allows to calculate such roots in the form of power series from some parameters. Examples of such computations are shown in [43]. The method of truncated equations and systems allows to calculate asymptotic solutions for ordinary differential equations: the Beletskii equation of satellite oscillations [5], the problem of periodic flyby of planets with a close approach to the Earth and for partial differential equations: boundary layer on a needle.

Generating Families of Periodic Solutions and Generating Families of Invariant Tori
As soon as electronic computers appeared, scientists began to calculate families of periodic solutions of the restricted three-body problem for different cases: Sun-Jupiter (µ ≈ 10 −3 ), Earth-Moon (µ ≈ 10 −2 ), etc. It turned out that these families are very similar, and their periodic solutions resemble solutions to the two-body problem. In 1968, M. Hénon [44] realized that it was necessary to consider the limits of these families for µ → 0. Generating families allow to understand and to prescribe by analytic computation singularities of families for small parameter µ.

Method
Let the Hamilton function H(µ) depend analytically on small parameters µ = (µ 1 , . . . , µ s ) and the corresponding Hamiltonian system has families of periodic solutions F j (µ). Some of these families may have limits F j (0) at µ → 0. Families F j (0) are called generating. Their solutions are formed by parts of solutions of the Hamilton limit system with µ = 0. If this limit system is integrable, then the generating families can be described analytically. This approach was proposed by Hénon [44]. It was used for the Hill problem and for the restricted three-body problem [1] (Ch. III-V), [7,8].

The Hill Problem
Its Hamilton function has the form The corresponding systeṁ describes the motion of the Moon (P 3 ) with zero mass under the influence of the attraction of the Sun (P 1 ), located at infinity, and the Earth (P 2 ) with mass 1, located at the origin of coordinates. The Hamilton function (83) is analytic in We apply the canonical coordinate transformation and we obtain the Hamiltonian systeṁ where Let ε = 2 | H | and H → −∞. Then, in the limit, we obtain system (4.2) with This is the Hénon problem [36]. For h 0 , system (84) is linear and, therefore, integrable. Since the Hamiltonian h 0 is homogeneous, it suffices to consider it for h 0 = 1/2. It has one regular periodic solution If the orbit (X 1 (t), X 2 (t)) of the solution to the Henon problem passes through the point then the body P 3 collides with the body P 2 and the solution cannot be continued through the collision. Therefore, point (4.3) divides the solution into independent parts. Hénon [36] found all the arc-solutions that start and end with such collisions. They form a countable set of two types. The arc-solutions of the first type are denoted by the symbols ±j, j ∈ N, and their orbits are epicycloids. For j = +1, +2, +3 they are shown in Figure 2. The orbits of arc-solutions with negative j values are symmetric to them about the X 2 axis. The arc-solutions of the second type are denoted by the letters i and e, their orbits are ellipses passing through the point (85). They are shown in Figure 3.

Theorem 12. ([45]) A sequence of arc-solutions that does not contain two successively identical arc-solutions of the second type is a generating solution to the Hill problem.
Here the generating family of periodic solutions consists of one solution. All known families of periodic solutions to the Hill problem include at least one generating solution.
In the restricted three-body problem, there is a countable set of one-parameter generating families of periodic solutions. Some of them are quite complex [7].
The same approach is applicable to families of n-dimensional invariant tori. The Hill problem has no generating families of them, but the RTBP has an infinite amount: indeed all solutions to the two-body problem with fixed irrational mean motion form such a family [1]. However, generating families of invariant tori are not studied yet.

Numerical Computation of Families of Periodic Solutions and of Invariant Tori
In Sections 2 and 3, we describe a way for computation of local such families near a stationary point, near a periodic solution and near n-dimensional invariant torus. For periodic solutions, there are many methods for numerical continuation of their families (see [46][47][48][49]) and a lot applications of these methods.
For invariant tori, there is one method for numerical computation of their family proposed by C. Simó [50] (see also [51]). Therefore, the method of computation of families of invariant tori appeared only since 1998.

Generalized Problems
Usually in celestial mechanics, bodies with non-negative masses are considered. However, Batkhin [13] proposed to consider problems where some masses are negative. In the Hill problem with body mass equal to −1 (called the anti-Hill problem), families of periodic solutions are extensions of families of periodic solutions to the classical Hill problem. Therefore, it is more convenient to calculate families of periodic solutions for both problems at once: for the Hill problem and for the anti-Hill ones. This approach provides new families of periodic solutions for the classical Hill problem. Figure 4 shows a diagram of the connections between these families of the Hill (left) and the anti-Hill (right) problems. The center column gives the generating solutions for these families.
The usefulness of negative masses in Celestial Mechanics was only realized after 2014.

Skeletons
In some parts of the phase space of the Hamilton system, there are many families of periodic solutions, invariant tori and they form the "skeleton" of this part of the phase space. Therefore, the calculation of such families is very useful for studying the structure of the phase space.

Stability in a Linear Multi-Parameter Hamiltonian System
The last two sections presented the description of methods providing an investigation of stability for the case of a Hamiltonian that depends in a polynomial way on the parameter vector P of the parameter space Π. Assuming a generic case, we describe in Section 8 the method which allows finding in the parameter space such domains that SP is Lyapunov stable either in the full problem or in the linear approximation of the full problem. This method essentially explores modern elimination theory and can be implemented in any computer algebra system.
The study of Lyapunov stability of the SP of a Hamiltonian system in the case in which the number of degrees of freedom is greater than two requires considerable effort. This is due, on the one hand, to the fact that the stability of the equilibrium position in the linear approximation can be destroyed by any arbitrarily small perturbation of higher order. On the other hand, the Arnold-Mozer theorem on the stability of the equilibrium position is inapplicable for large dimensions. However, for many applications, the formal stability proposed by J. Moser [27] is quite sufficient.
Let us consider the simplest invariant manifold of dimension 0, namely, the stationary points of a Hamiltonian function H(z) depending in a polynomial way on the parameter vector P ∈ Π ≡ R s , where s is the dimension of the parameter space Π. The proposed methods can be applied to the case of manifolds of large dimensions, i.e., in the case of studying Hamiltonian phase flow near a periodic solution or near an invariant torus.
In the generic case, an analytic time-independent Hamiltonian function H(z) in the vicinity of the SP, coinciding with the origin, is expanded into a convergent series of homogeneous polynomials H k of degree k of its phase variables z = (x, y) The well-known Lagrange-Dirichlet theorem [26] and [55] ( § 29) states that SP is stable if the quadratic form H 2 (z) is sign-determined (see Theorem 3).
If the number of degrees of freedom is not more than two • Stability is determined by the Arnold-Moser theorem in the absence of resonances of order four or less, which requires normalizing H to order four; • For resonances of order less than four, the stability conditions are derived in the works of A. P. Markeev and A. G. Sokolsky (see [28] and also Section 2.1.5).
The series (86) starts with the quadratic Hamiltonian H 2 (z; P) defining the local dynamics near the SP. The behavior of the phase flow in the first approximation is described by a linear Hamiltonian systeṁ Let us recall here the main properties of a linear Hamiltonian system.

1.
If λ j is an eigenvalue of the matrix B, then −λ j is also its eigenvalue. All eigenvalues λ j , j = 1, . . . , 2n, of the matrix B can be reordered in such a way that λ j+n = −λ j , j = 1, . . . , n.

2.
The characteristic polynomialf (λ) of the matrix B contains only even powers of λ, so it is a polynomial in µ = λ 2 . The following [16] polynomial is called semi-characteristic.

3.
If λ j = 0 for any j, i.e., the SP is hyperbolic, then it is structurally stable according to the Hartman-Grobman theorem.

4.
For an elliptic SP, the behavior of the phase flow in its vicinity can only be obtained by taking into account the nonlinear terms. Usually this is performed using KAM-theory, but here such study is performed using the Hamiltonian normal form described in Section 2.

Definition 5.
The stability set Σ of the system (87) is the set of all values of parameter P ∈ Π for which the SP z = 0 is Lyapunov stable.
In terms of roots of a semi-characteristic polynomial (88), the condition of stability of the SP is given by Let us first briefly recall the definition of a k-th order subdiscriminant of an arbitrary monic polynomial of n degree (for more details, see [59]).
Definition 6. Let f n (x) = x n + a 1 x n−1 + · · · + a n−1 x + a n (89) be some monic polynomial from the x variable. Then its k-th subdiscriminant D (k) ( f n ), where x j are roots of the polynomial (89), I is any non empty subset of the set {1, 2, . . . , n}, #(I) is its cardinality (number of elements in I). For k = n − 1 we put D (n−1) ( f n ) = n and for k = n we put D (n) ( f n ) = 1. For k = 0 we obtain D (0) ( f n ) = D( f n ) the classical discriminant.
The criteria of reality and negativity of roots is given by the following statement.

Proposition 1.
For all roots of a polynomial f (µ) of degree n to be real, negative, and distinct, it is necessary and sufficient that the conditions where D (k) ( f ) is the k-th sub-discriminant of the polynomial f (µ). Assuming that functions f n (P) and D( f ) are polynomials in P, we obtain the problem of description of the affine varieties F 0 and F 2 , which divide the space of parameters Π into cells. The appropriate cells are included in the stability set Σ according to Proposition 1. The boundaries of these cells are selected according to Proposition 2 with the help of multiple roots founded from (94).

Remark 1.
According to Proposition 1 the set of stability Σ is a solution of the system (90), i.e., it is a semi-algebraic set. It can be computed with the algorithm of cylindrical algebraic decomposition (CAD) [59] (Sect. 5). The CAD algorithm usually demonstrates rather high computational complexity.
For the system (87), let the stability set Σ be calculated. Then for each value of P ∈ Σ the eigenvalues of the matrix B(P) are purely imaginary λ j = −λ j+n = iω j , j = 1, . . . , n and elementary divisors are simple. In this case, the Hamiltonian H 2 (z) is reduced to the normal form with a set of invariants σ i (P), i = 1, . . . , n, The change of the signature s of the quadratic form H 2 can occur while crossing the hypersurface F 0 , which divides the stability set Σ into regions Σ j with the same signature value s. Let us distinguish those Σ j regions for which s = ±2n, i.e., for them the Lagrange-Dirichlet theorem is inapplicable. The formal stability of the SP should be studied in these regions Σ j of the stability set Σ.

Remark 2.
To determine the invariants σ i and consequently the signature s of the normal form (91), it is not necessary to perform a normalization procedure. It is sufficient to use the method described in [17].

Studying of Formal Stability of Stationary Point
Here, we give a schematic description of the method for studying the formal stability of the SP. This method is based on the following key results: normal form of the Hamiltonian at the SP, Bruno's Theorem 4 on formal stability and q-analog of the classical elimination theory.
In the absence of strong resonances between eigenvalues of a linearized Hamiltonian system in the neighborhood of the SP, the condition for its formal stability is formulated by the Bruno theorem [29]. In the paper [17], a scheme for investigating formal stability was proposed. This scheme assumes, firstly, that the set of stability Σ of the SP in linear approximation is computed, and secondly, that the so-called resonant sets R q of the characteristic polynomialf (λ) corresponding to strong resonances 2:1, 3:1 between eigenvalues are found. Here we consider a method for investigating the stability of the SP of a multi-parameter Hamiltonian system with more than two degrees of freedom, based on the description of the discriminant and resonant sets of the real polynomial proposed in [18].
One of the earliest results on the formal stability of the SP is Theorem 14 (Birkhoff). If for some value of the parameter vector P ∈ Σ all components of the vector λ of eigenvalues are rationally incommensurable, i.e., the equation λ, p = 0 has no integer solution p ∈ Z n , p = 0, then there exists a formal canonical transformation w of the variables (x, y) → (ρ,ϕ) such that the Hamiltonian H(w(x, y)) is a power series over ρ. Here ρ, ϕ are action-angle variables.
The Birkhoff theorem is difficult to apply because resonance sets are dense everywhere in the parameter space Π.
Theorem 15 (Bruno [29]). If on any pair of nonzero integer vectors K 1 and K 2 of ortant k i 0, i = 1, . . . , n, which are solutions to the equation bilinear form K T 1 BK 2 = 0 at λ = 0, then the SP z = 0 of the Hamilton system is formally stable.
Note that the condition (92) of the Bruno theorem is equivalent to the fact that the semi-algebraic system g 2 (ρ) = g 4 (ρ) = 0, ρ 0 is incompatible.
Thus, to apply the Bruno theorem on formal stability, it is necessary to find the boundaries of regions in the parameter space Π defined by resonance sets.
Here we consider a scheme for studying the formal stability of the equilibrium position of a Hamiltonian system under the following assumptions: • Number of degrees of freedom more than two; • The quadratic form H 2 (z) in the expansion (86) is nondegenerate and not sign-defined; • Hamiltonian function H(z) smoothly depends on the parameter vector P.
The proposed research scheme essentially uses the normal form of the Hamiltonian system near the SP, a method for computing the stability set Σ of a linear multi-parameter Hamiltonian system [16] and a method for computing the resonance set of a real polynomial [18]. These methods, in the case of polynomial dependence on the parameter P, essentially use so-called quantum calculus [60].
The application of q-subdiscriminants D (k) q ( f n ) of the characteristic polynomial f n allows not only to find out if it has commensurable roots, but also, under certain conditions, to find these roots without having to calculate all its eigenvalues. If the coefficients of polynomial (89) depend on parameters, both q-subdiscriminants are functions of these parameters, which makes it possible to determine at what values of q resonance takes place and to find its multiplicity and order.
To determine the rational comparability of the roots of the polynomial, we will use q-analogs of the classic derivative and subdiscriminant. Definition 7. Jackson derivative (q-derivative, q-differential Jackson operator): The Jackson derivative has all the properties of an ordinary derivative. In addition to the above properties, we note that the q-derivative function x n is equal to [n] q x n−1 , and applying q-derivative to q-binomial of degree n again gives q-binomial of degree n − 1 multiplied by [n] q .
With the q-derivative, the q-analog of the classic polynomial discriminant is now determined.
Definition 8. Define q-discriminant D q ( f n ) of polynomial f n (x) as the resultant of a polynomial pair f n (x) and (A q f n )(x): D q ( f n ) = (−1) n(n−1)/2 Res x ( f n (x), (A q f n )(x)).
The equality to zero of the q-discriminant of polynomial f n (x) with fixed q is a signature of the existence of at least one pair of q-commensurable roots; however, the detailed structure of all commensurable roots can be obtained using the sequence of qsubdiscriminants of various orders of the polynomial f n (x).
Theorem 16 ([18]). The f n (x) polynomial has exactly n − d different sequences of q-commensurable roots if and only if the first non-zero q-subdiscriminant in sequence (93) of k-th q-subdiscriminants D (k) q ( f n ), k = 0, . . . , n − 2, has the index d.
All commensurable roots of the polynomial f n (x) are the roots of the largest common polynomial divisor of f n (x) and its q-derivative (A q f n )(x): f q (x) = GCD( f n (x), (A q f n )(x)).
Theorem 16 states that the degree of the polynomialf q (x) is equal to the number d of the first non-zero q-subdiscriminant in the sequence S q ( f n ).
Note that q-subdiscriminants are calculated using any of the matrix methods to calculate the classical subresultants of the polynomial pair f n (x) and (A q f n )(x) [61].
Let, in the conditions of Theorem 16, the first non-zero q-discriminant have the index d which is the determinant of this inner. Then, as shown in [18], the following proposition takes place. Proposition 3. If in the sequence (93) the first one different from zero q-subdiscriminant D (d) q ( f n ) has the index d, thenf The roots of polynomial (94) are either q-commensurable roots for q = 1 or multiple roots for q = 1 of the initial polynomial (89).
Let us describe schematically the procedure of formal stability studying.
In the first step, the stability set Σ ⊂ Π of the linear Hamiltonian system (87) is computed. Then the open regions Σ j in which the signature s of the quadratic form H 2 (z) is not equal to ±2n, i.e., the Lagrange-Dirichlet theorem on Lyapunov stability is not applicable.
In the next step, the resonance sets R 2 ( f ) and R 3 ( f ) are computed for the semicharacteristic polynomial (88), dividing each region Σ j into subregions that have resonances of order four and higher. Then, in these regions, the procedure of reducing the Hamiltonian (86) to the fourth-order normal form g(ρ, ϕ) = g 2 (ρ) + g 4 (ρ) is applied. After normalization, we check whether the condition of the Bruno theorem is satisfied.
On sections of resonance sets R k ( f ), k = 2, 3, where there is a single resonance, the theorems of [28] (Ch. 4, § § 2, 3) are applied. For the case of multiple resonance, as well as for resonance of order two, the formal stability study can be carried out using a different approach.
On the one hand, the presence of resonance leads to additional first formal integrals, which allows decreasing the number of degrees of freedom of the normalized system [1] (Ch. I, §3), and in some cases even to integrate it. On the other hand, finding the invariant coordinate subspaces of the Hamiltonian normal form reduces the phase flow to subspaces of fewer dimensions, on which one can search the realization of conditions of analyticity of the normalizing transformation.