ON MOTIONS OF A CONSERVATIVE SYSTEM ON INVARIANT MANIFOLDS

In [Irtegov and Burlakova, 2017], the algorithms for the qualitative analysis of conservative systems have been presented. These are based on the RouthLyapunov method [Lyapunov, 1954] and some its modifications [Irtegov and Titorenko, 2009] as well as computer algebra methods [Cox, Little, and O’Shea, 1997]. In the paper proposed, the application of the algorithms is demonstrated by analysing a conservative system, the study of which is also of interest. We conduct qualitative analysis for the differential equations describing the rotational motion of a rigid body with a fixed point in two constant force fields. Similar problems arise, e.g., in space dynamics [Sarychev and Gutnik, 2015], quantum mechanics [Adler, Marikhin, and Shabat, 2012], [Smirnov, 2008]. In the phase space of the problem, we isolate the invariant manifolds of maximal dimension and study the equations of motion on them. For these equations, solutions (and their families) corresponding in the original phase space of the problem to permanent rotations and pendulum-like oscillations of the body as well as the invariant manifolds of 2ndlevel, which these solutions belong to, have been found and their Lyapunov’s stability has been investigated. The possibility of stabilization for the motions of conservative systems, whose stability conditions have the form of some constraints on the constants of first integrals, was discussed.


Introduction
Let us consider the differential equations 2ṗ = q r + bδ 3 ,γ 3 = γ 1 q − γ 2 p, 2q = x 0 γ 3 − p r,δ 1 = δ 2 r − δ 3 q, describing the rotational motion of a rigid body around a fixed point in two constant force fields. Here p, q, r are the projections of the angular velocity vector onto the coordinate axes rigidly attached to the body; γ i , δ i (i = 1, 2, 3) are the components of the direction vectors of the 1st and 2nd force fields, respectively; x 0 , b are the components of the radius vectors of the 1st and 2nd force centers. The distribution of mass in the body corresponds to the Kowalewski integrable case [Kowalewski, 1889]. Eqs. (1) admit the following first integrals: 2H = 2(p 2 + q 2 ) + r 2 + 2(x 0 γ 1 − b δ 2 ) = 2h, V 1 = (p 2 − q 2 − x 0 γ 1 − b δ 2 ) 2 + (2p q − x 0 γ 2 + b δ 1 ) 2 = c 1 , V 2 = γ 2 1 + γ 2 2 + γ 2 3 = 1, V 3 = δ 2 1 + δ 2 2 + δ 2 3 = 1, V 4 = γ 1 δ 1 + γ 2 δ 2 + γ 3 δ 3 = c 2 , V 5 = x 2 0 (pγ 1 + qγ 2 + r 2 γ 3 ) 2 + b 2 (pδ 1 + qδ 2 + r 2 δ 3 ) 2 − x 0 b r[(γ 2 δ 3 − γ 3 δ 2 )p +(γ 3 δ 1 − γ 1 δ 3 )q + r 2 (γ 1 δ 2 − γ 2 δ 1 )] +x 0 b 2 γ 1 (δ 2 1 + δ 2 2 + δ 2 3 ) − x 2 0 b δ 2 (γ 2 1 + γ 2 2 + γ 2 3 ) − bx 0 (bδ 1 − γ 2 x 0 )(δ 1 γ 1 + δ 2 γ 2 +δ 3 γ 3 ) = c 3 , where V 5 is the additional first integral obtained in [Reiman and Semenov-Tyan-Shanskii, 1988], [Bobenko, Reyman, and Semenov-Tian-Shansky, 1989]. Thus, the system under consideration is completely integrable. When b = 0, Eqs. (1) correspond to the Kowalewski integrable case of motion of the top in gravitational force field. There exists a series of works devoted to the topological analysis of system (1) (see, e.g., [Ryabov P.E. et al., 2012], [Kharlamov and Ryabov, 2017] and the bibliography therein). The analysis was made on the basis of the methods and approaches described in [Bolsinov and Fomenko, 1999]. In the work presented, we used the Routh-Lyapunov method and some its modifications in combination with computer algebra methods to analyse Eqs. (1) that enabled us to obtain new results in the problem under consideration. As is well-known, the problem of qualitative analysis of differential equations is to find special solutions (equilibria, periodic motions) of these equations and to investigate their stability and bifurcations. In the case of conservative systems, the variety of the special solutions is expanded through stationary sets. By these sets, we mean sets of any finite dimension on which the problem's first integrals (or their combinations) take a stationary value. Zero-dimensional sets having this property are traditionally called stationary solutions. By analogy, we call positive dimensional sets the stationary invariant manifolds (IMs) [Irtegov, 1986]. In the phase space of system (1), such sets were studied in [Irtegov and Titorenko, 2016]. In order to obtain the complete phase portrait of this system, it is necessary to construct phase portraits on its special IMs, i.e. to find the special solutions of the equations of motion on the IMs and to analyse their qualitative properties. The present paper solves this problem. The paper is organized as follows. In section 2, we find the special IMs of maximal dimension for system (1) on the basis of the Routh-Lyapunov method. In addition to previously known IMs, new IM has been obtained. Section 3 is devoted to finding the special solutions on these IMs. In section 4, the Lyapunov stability of the found solutions is investigated. In section 5, the possibility of stabilization for the motions of conservative systems is discussed. Finally, we give a conclusion in section 6.
2 On the Special IMs of Equations of Motion in the Original Phase Space of the Problem We shall seek the IMs of maximal dimension for Eqs.
(1). Since the problem's first integrals are IMs of codimension 1, let us begin with IMs of codimension 2. According to the Routh-Lyapunov method, the stationary sets of the differential equations under study can be obtained from the necessary extremum conditions for the elements of the algebra of the problem's first integrals. Following this method, we construct the complete linear combination from the first integrals (the family of the integrals): Other combinations of the integrals also are acceptable, including nonlinear ones. Here λ j (j = 0, . . . , 5) are the parameters of the family of the integrals K.
Then, we write the necessary conditions for the family K to have an extremum with respect to the phase variables: (2) The solutions of system (2) allow one to define the IMs (or their families) for differential equations (1) which correspond to the family of the first integrals K. So, the problem of finding the IMs of Eqs. (1) by the Routh-Lyapunov method is reduced to solving the system of algebraic equations (2). It is the system of 9 cubic equations with the parameters λ j , b, x 0 . To simplify further computations, it is convenient to use the variables of the kind [Ryabov P.E. et al., 2012]: where i = √ −1.
Their use allowed us to avoid cumbersome computations in the problem under consideration as well as to represent the obtained results in more compact form than when the original variables were employed (see [Irtegov and Titorenko, 2016]).
In variables (3), Eqs. (1) and first integrals (2) can be written as and respectively. The stationary conditions for the integral K in new variables take the form: We shall seek the IMs of codimension 2 for differential equations (4) as solutions of system (6). For this purpose, the Gröbner basis method [Cox, Little, and O'Shea, 1997] will be applied. Software codes of this method are included in many widespread computer algebra systems, e.g. Mathematica [Wolfram, 2003], Maple [Char, Geddes et al., 1992]. We shall employ the Mathematica program GroebnerBasis based on the Buchberger algorithm [Buchberger, 1976]. In order to compute a Gröbner basis for a polynomial system the program uses the system and its variables as input data. A resulting system is equivalent to the initial one, but has other form.
In our work, to obtain the desired IMs, we construct a lexicographical Gröbner basis for the polynomials of system (6) with respect to part of the phase variables and parameters λ i . In this case, a resulting system is not, generally speaking, equivalent to the initial one. It has the special structure: one part of its equations define the sought IM, and the other part enables to obtain the first integrals for the equations of motion on this IM. Let us take, e.g., λ 0 , λ 1 , λ 2 , λ 3 , λ 5 , y 1 , y 2 as unknowns. Here the number of the phase variables determines the codimension of the desired IM. Next, we compute a lexicographical basis for system (6) with respect to the above variables. A resulting system can be split up into two subsystems. The first subsystem has the form: x 2 z 1 − y 2 z 2 = 0, x 1 z 2 − y 1 z 1 = 0.
It is easy to verify by IM definition that Eqs. (8) determine the IM of codimension 2 of differential equations (4): the derivative of expressions (8) calculated by virtue of Eqs. (4) vanishes on the given expressions. The motions on IM (8) are described by the equations: They are derived by elimination of the variables y 1 , y 2 from Eqs. (4) with the use of (8).
Each of the latter two Eqs. (7) defines the family of IMs of codimension 1 for differential equations (9). It is easy verify by IM definition. Let us resolve these equations with respect to λ 2 , λ 3 : Having computed the derivatives of the right-hand sides of relations (10) in virtue of Eqs. (9), we find that they are identically equal to zero. So, one can conclude that these expressions are the first integrals of differential equations (9).
In the original variables, Eqs. (8) can be written as: It is easy to verify by IM definition that the equations define the IM of codimension 2 of system (1). After analogous transformations, Eqs. (12) become and determine the IM of codimension 2 of differential equations (1).
Having constructed a lexicographical basis for the polynomials of system (6) with respect to, e.g., λ 0 , λ 1 , λ 2 , λ 3 , λ 4 , x 1 , x 2 , we have found two IMs of codimension 2 of differential equations (4). One of them coincides with IM (12). The 2nd differs from the above IMs and is given by the equations In the original variables, Eqs. (15) can be written as and determine the IM of codimension 2 for the equations of motion (1). Here For the polynomials of system (6) there were constructed lexicographical bases with respect to other similar combinations of the parameters λ i and phase variables, but they did not give new results.
Another IM of codimension 2 in this problem has been found in [Bogoyavlenskii, 1984]. When the constant c 1 of integral V 1 (2) is zero, the equations define the IM of codimension 2 of system (1).
3 Finding the Special Solutions of the Equations of Motion on the IMs 3.1 On IM (13) The equations of motion on IM (13) can be written as: They have the following first integrals: Eqs. (18) and integrals (19) are derived by the elimination of the variables δ 1 , γ 2 from Eqs. (1) and integrals (2), respectively, with the use of (13).

Finding the IMs of 2nd-level
Consider the problem of finding the special IMs for Eqs. (18). Such IMs we call the IMs of 2nd-level. In order to obtain them, the technique of section 2 is applied. First, we solve the stated problem for differential equations (9) on IM (8). The first integrals of these equations are derived by the elimination of the variables y 1 , y 2 from expressions (5) with the use of (8): We take independent integrals from the above ones (e.g.,H,Ṽ 1 ,V 2 ,V 4 ,V 5 ), and construct the linear combination (21) from them and write the stationary conditions for the integralK with respect to the variables w 1 , w 2 , w 3 , x 1 , x 2 , z 1 , z 2 : Here µ i (i = 0, . . . , 5) are the parameters of the family of the integralsK. Then, we compute a lexicographical basis for the left-handed part of system (22) with respect to µ 0 , µ 1 , µ 2 , µ 4 , x 1 , x 2 and, as a result, obtain the IM of codimension 2 of differential equations (9). Its equations are given by In the original variables, these equations take the form and define the IM of codimension 2 of differential equations (18).
Finding the stationary solutions. Now, we shall seek the stationary solutions of Eqs. (18). As mentioned before, stationary solutions are usually found by the Routh-Lyapunov method from the necessary conditions of stationarity for a family of problem's first integrals. In the problem under consideration, this approach leads to solving a system of 7 polynomial equations of 5-8 degrees. To simplify computations, we apply the technique described in [Irtegov and Titorenko, 2015]. First, we find solutions of differential equations (18) under zero derivatives, and then, obtain the families of the integrals which take stationary values on these solutions. Equate the right-hand sides of differential equations (18) to zero and add relationsV 2 = 1,V 3 = 1 (19) to them: For the subsystem (24), construct a lexicographical basis with respect to part of the phase variables, e.g., q, r, γ 1 , γ 3 , δ 2 , δ 3 . A resulting system is split up into several subsystems. Let us consider two of them. The 1st subsystem: The 2nd subsystem: bp It is easy to verify by IM definition that Eqs. (25), (26) define two one-dimensional IMs of differential equations (18). The vector field on each of the IMs is described by the equationṗ = 0 which has the following family of solutions: Eqs. (25), (26) together with (27) determine 4 families of solutions of differential equations (18): Here p 0 is the parameter of the family, α = ( x 2 0 p −4 0 − 1. All the above solutions can be "lifted up" into the phase space of system (1). To this end, it is sufficient to add the equations of IM (13) to that of IMs (23), (25), (26). In the case of the families of solutions, the values of the variables γ 2 , δ 1 , obtained from Eqs. (13) under the corresponding values of the rest of the variables, it is necessary to add to relations (28), (29).
In order to find out the direction of the force fields corresponding to solutions (28), (29), we substitute the solutions into integralV 4 (19). It becomes identically " − 1 ′′ and " + 1 ′′ , respectively. So, with mechanical viewpoint, the elements of the families of solutions under study in the original phase space correspond to the permanent rotations of the top in the opposite in direction (or parallel) force fields. Resolve Eqs. (25) and (26) with respect to q, r, γ 1 , γ 3 , δ 2 , δ 3 and, then, substitute the resulting expressions into (23). The latter equations turn into identities. Whence, one can conclude that the IMs defined by Eqs.
From Eqs. (31), we find the constraints on µ i under which solutions (28) satisfy these equations: From now and further, we use the following denotations: On substituting the latter expressions into (30), one obtains: As one can see from (32), the family of the integralŝ K 1 is divided into 3 subfamilies corresponding to the coefficients of µ 0 , µ 1 , µ 4 . Both the family itself and each of its subfamilies takes a stationary value on the elements of the family of solutions (28). It is easy to verify by direct computations. Similarly, we have the family of integrals 2K 2 = K 0 + µ 4 (V 2 +V 3 − 2V 4 ). Both this family and its subfamilies (the coefficients of µ 0 , µ 1 , µ 4 ) takes a stationary value on the elements of families (29).

On IM (14)
The equations of motion on IM (14) have the form: They are derived by the elimination of the variables δ 1 , γ 2 from Eqs.
(1) with the use of (14). Analogously, the first integrals of Eqs. (33) are derived by the elimination of the variables δ 1 , γ 2 from expressions (2) with the use of (14): Finding the stationary solutions For Eqs. (33), by the technique of section 3.1, we have found 2 one-dimensional IMs and 4 families of solutions, which belong to them. The equations of the 1st IM are given by The family of solutions belonging to IM (35): The equations of the 2nd IM are written as: bp The family of solutions belonging to IM (37): Here p 0 is the parameter of the families, α = (x 4 0 − (b 2 + x 2 0 ) p 4 0 ) 1/2 . With the use of the technique of section 3.1, the above solutions can be "lifted up" into the phase space of system (1).
On substituting solutions (36), (38) into integralV 4 (34) it becomes identically " − 1" and " + 1", respectively. So, with mechanical viewpoint, the elements of the given families of solutions in the original phase space correspond to the permanent rotations of the top in the opposite in direction (or parallel) force fields. As one can see from (35)-(38), the solutions formally coincide with solutions (25), (26), (28), (29) of section 3.1. It is not difficult to show that in the original phase space of the problem the same IM corresponds to, e.g., IM (25) and IM (35).
They determine the one-dimensional IM of differential equations (1). Similarly, one can establish that in the original phase space the following families of solutions correspond to the families of solutions (28) on IM (13) and (36) on IM (14): They belong to IM (39). Resolve Eqs. (39) with respect to q, r, γ i , δ i (i = 1, 2, 3) and, then, substitute the resulting expressions into Eqs. (13), (14). They become identities. Whence, one can conclude that IM (39) belongs to both IM (13) and IM (14), i.e. it lies at their intersection. Analogously, we find that in the phase space of system (1) the same IM corresponds to IM (26) on IM (13) and IM (37) on IM (14) as well as the same families of solutions correspond to the families of solutions (29) on IM (13) and (38) on IM (14). Using the technique of section 3.1, we have obtained the families of integrals the elements of which take stationary values on the elements of families (36) and (38), respectively: ] . (40) The integralsV 1 ,V 2 +V 3 ± 2V 4 take stationary values on IMs (35), (37), i.e. these IMs are stationary.
On special solutions corresponding to pendulum-like motions In the problem under study, the solutions, for which p = r = 0 or q = r = 0, correspond to the oscillations and rotations of the top in a plane around one of its principal axes. Consider the following solutions of Eqs. (33): and It is easy to verify by IM definition that Eqs. (41), (42) define four IMs of codimension 4 of Eqs. (33). The equations of motion on IM (41) have the form: In the original phase space of the problem, with mechanical viewpoint, IMs (41) correspond to pendulumlike oscillations of the body around its immobile principal axis Oy. The equations of motion on IM (42) can be written as: In the original phase space, IMs (42) correspond to pendulum-like oscillations of the body around its immobile principal axis Ox. On substituting IMs (41), (42) into integralV 4 (34) it is identically equal to zero. So, with mechanical viewpoint, in the original phase space of the problem the IMs correspond to the pendulum-like oscillations of the top in orthogonal force fields. By direct computations, one can verify that the inte- ] take stationary values on corresponding IMs (41). Analogously, the in- take stationary values on IMs (42). Thus, these IMs are stationary.
Finding the IMs of 2nd-level Now, for differential equations (33), we shall seek the IMs of maximal dimension, which the above solutions belong to. First, we solve this problem for the equations of motion on IM (12): They are derived by the elimination of the variables y 1 , y 2 from Eqs. (4) with the use of (12). Consider some of explicitly written relations (11), e.g., They define the IMs of codimension 1 for Eqs. (43).
In the original variables, expressions (44) are:

respectively. Each of the equations
Having substituted the latter expressions into (45), we have the equations for the families of IMs of codimension 1, which solutions (36) belong to: Here p 0 is the parameter of the families. Thus, solutions (36) lie at intersection of the families of IMs defined by Eqs. (46). Let us consider the expressions in the square brackets (46) as the polynomials of p 0 . To eliminate p 0 from one of them, we compute the resultant of these polynomials. It has the form: Res = P 1 P 2 P 3 , where P i (i = 1, 2, 3) are the polynomials of p, q, r, γ 1 , γ 3 , δ 2 , δ 3 .
It is easy to verify by IM definition that each of the equations P i = 0 (i = 1, 2, 3) determines the IM of codimension 1 for differential equations (33). Both solutions (36) and one-dimensional IM (35) belong to one of these IMs. Its equation is written as: Similarly, we find the IM of codimension 1, which solutions (38) and one-dimensional IM (37) belong to. This IM is given by the equation: The equation of the IM, which solutions (41), (42) belong to, is derived from Eqs. (45) when λ 4 = 0. It is given by G 1 = 0.

On IM (16)
The equations of motion on IM (16) are rather bulky. It is more convenient to analyse these equations in variables (3), i.e. the equations of motion on IM (15). They are derived by the elimination of the variables x 1 , x 2 from Eqs. (4) with the use of (15): where σ = w 2 z 2 (w 2 1 y 2 + z 2 1 ) + w 1 z 1 (w 2 2 y 1 + 2w 2 w 3 z 2 + z 2 2 ). As one can see from Eqs. (47), they have no solutions corresponding to equilibria of the body, i.e. for which w 1 = w 2 = w 3 = 0. In this case, the equations have a singularity. For the same reason, we could not obtain the solutions corresponding to the permanent rotations and pendulum-like motions of the top. However, Eqs. (47) can be integrated on one of their IMs. Consider the first integral of these equations:
Here λ 1 , λ 5 are some constants. This integral has been obtained from stationary equations (6) together with IM (15) in section 2.
It is easy to verify by IM definition, the equations w 2 = 0, w 3 = 0, y 1 = 1, y 2 = 1, z 1 ± w 1 = 0. (49) define two one-dimensional IMs of differential equations (48). The motions on these IMs are described by the equations respectively. They can be integrated in the elementary functions: w 1 (t) = C 1 exp ∓it/2 , where C 1 is a constant of integration.
4 On Stability of the Stationary Solutions 4.1 On IM (13) Let us investigate the stability of the families of solutions (28) on IM (13) on the basis of Lyapunov's linear stability theorems [Lyapunov, 1992]. As is wellknown, they provide necessary stability conditions. We introduce the deviations from the elements of the family under study, where p 0 , q 0 , r 0 , δ 0 2 , δ 0 3 , γ 0 1 , γ 0 3 are the values of the variables in unperturbed solution (28). Next, we write the equations of first approximation in the vicinity of the elements of the family: From now and further, the denotations of section 3 are used: . Taking into account the conditions for solutions (28) to be real, the characteristic equation of system (50) has only zero and pure imaginary roots under the following constraints on the parameters b, x 0 , p 0 : The above-represented inequalities provide the necessary conditions of the stability for solutions (28). We derive the sufficient stability conditions for these solutions by the Routh-Lyapunov method. This method uses first integrals and their families taking stationary values on solutions under study. The problem is to verify the sign-definiteness conditions for the 2nd variation of the corresponding family of integrals obtained in the vicinity of the studied solution on the linear manifold defined by the first variations of "conditional" integrals.
We shall use the family of integralsK 1 (32) taking stationary values on the elements of the families of solutions (28) under the condition µ 4 = 0. In this case, it has the form: The 2nd variation of the integralK 0 in the vicinity of the elements of the studied family of solutions in the deviations y i (i = 1, . . . , 7) on the linear manifold δH = p 0 x 0 (x 0 y 3 − by 1 ) + 2p 2 0 (x 0 y 5 + by 6 ) ± αy 7 = 0, is written as: 2δ 2K 0 = α 11 y 2 1 + α 12 y 1 y 2 + α 15 y 1 y 5 + α 17 y 1 y 7 + α 22 y 2 2 + α 25 y 2 y 5 + α 27 y 2 y 7 + α 55 y 2 5 + α 57 y 5 y 7 + α 77 y 2 7 , where The conditions for the quadratic form δ 2K 0 to be positive definite are sufficient for the stability of the elements of the family of solutions (28). In the form of Sylvester's inequalities, they are: Taking into account the conditions for solutions (28) to be real, the above inequalities are consistent under the following constraints on the parameters µ 0 , µ 1 , b, x 0 , p 0 : and As one can see from (54), inequalities (54) are divided into 2 groups. The first (the constraints on the parameters p 0 , b, x 0 ) provides the sufficient conditions for the stability of the elements of families (28). The 2nd separates some subfamily from the family of the integralŝ K 0 , the elements of which give us a possibility to derive these conditions. On comparing inequalities (54) and (51) we conclude that the conditions (51) are necessary and sufficient for the stability of the elements of families (28) up to the boundary of the conditions. We have derived the above stability conditions for the families of solutions (29). The family of integralsK 0 (52) was used to obtain the sufficient conditions.

On IM (14)
The stability of the stationary solutions on IM (14) is analyzed also as above. To obtain the sufficient stability conditions for the elements of families (36) we use the family of integralsǨ (40) when λ 2 = (b 2 − x 2 0 )(λ 1 + 1/(b 2 + x 2 0 )), λ 3 = 0. Under these constraints the family of the integralsǨ becomes: The 2nd variation of the integralǨ 1 obtained in the vicinity of the elements of the studied family of solutions on the linear manifold δȞ = p 0 x 0 (x 0 ξ 3 − bξ 1 ) + 2p 2 0 (x 0 ξ 5 + bξ 6 ) ± αξ 7 = 0, δV 2 = 2αp 3 0 x 2 0 (x 0 ξ 3 − bξ 1 ) + α 2 p 0 x 2 0 ξ 4 +2p 5 0 x 2 0 (b 2 ξ 2 − x 2 0 ξ 4 ) +D 2 (αξ 5 + p 2 0 x 0 ξ 7 ) = 0, is written as: δ 2Ǩ 1 = ξ 2 1 + α 15 ξ 1 ξ 5 + α 16 ξ 1 ξ 6 + ξ 2 3 + α 35 ξ 3 ξ 5 + α 36 ξ 3 ξ 6 + α 55 ξ 2 5 + α 56 ξ 5 ξ 6 + α 66 ξ 2 6 , where Here ξ i (i = 1, . . . , 7) are the deviations from the elements of the studied family of solutions (36). The conditions of positive definiteness for the quadratic form δ 2Ǩ 1 in the form of Sylvester's inequalities are given by Taking into account the conditions for solutions (36) to be real, the inequalities are consistent under the following constraints on the parameters λ 1 , p 0 , b, x 0 : λ 1 + 1 D < 0 and p 0 ̸ = 0 and (p 2 0 + x 0 < 0 or x 0 > p 2 0 ) and where σ = x 0 √ x 2 0 p −4 0 − 1. Whence, we can conclude on the stability of the elements of the family under study. Likewise as above, the obtained inequalities are divided into 2 groups. The first (the constraints on the parameters p 0 , b, x 0 ) provides the sufficient stability conditions for the elements of families (36). The second separates some subfamily from the family of the inte-gralsǨ 1 , the elements of which enable us to derive these conditions. The necessary stability conditions for the studied family of solutions (36) were obtained from the analysis of differential equations (33) linearized in the vicinity of the elements of this family. They are given by p 0 ̸ = 0 and (p 2 0 + x 0 < 0 or x 0 > p 2 0 ) and Their comparison with conditions (56) shows that they are necessary and sufficient for the stability of the elements of the families of solutions (36) up to the boundary of the conditions. The above stability conditions have also been derived for the families of solutions (38). To obtain the sufficient conditions, we used the family of integrals which is similar to (55).

On Stabilization of Motions in
Conservative Systems Let us consider the following equations: These define two IMs of codimension 6 for differential Eqs.
The above presented problem of stability can be considered as a stabilization problem. When stability conditions for the solution of conservative system have the form of some constraints on the constants of problem's first integrals, as in the above case, one can speak about the stabilization of this solution by means of these constrains. Moreover, in the case discussed, we can speak about an optimal stabilization, because the corresponding first integral takes an extremum value on this solution.

Conclusion
In this work, with the use of the algorithms presented in the previous paper, the equations of motion for the generalized Kowalewski top were studied. We have made qualitative analysis for the differential equations obtained as a result of reduction of the initial equations of motion to their IMs. Some of these IMs are previously known, and one of them has been obtained in this work. For the equations of motion on the IMs, the special solutions and their families have been found. The linear and nonlinear combinations of the first integrals taking stationary values on the found solutions have been constructed. These combinations were used to analyse the stability of the solutions. For a series of the solutions, sufficient stability conditions have been compared with necessary ones. The possibility of stabilization of motions in conservative systems was discussed. The obtained results prove the efficiency of the approach and algorithms which were applied. The algorithms and the results can be used in the qualitative analysis of conservative systems.