Determination of the stability boundaries for the hamiltonian systems with periodic coefficients

Abstract. We consider the hamiltonian system of linear differential equations with periodic coefficients. Using the infinite determinant method based on the existence of periodic solutions on the boundaries between the domains of stability and instability in the parameter space we have developed the algorithm for analytical computation of the stability boundaries. The algorithm has been realized for the second and the fourth order hamiltonian systems arising in the restricted many-body problems. The stability boundaries have been found in the form of powers series, accurate to the sixth order in a small parameter. All the computations are done with the computer algebra system Mathematica.


Introduction
Let us consider the linear hamiltonian system of differential equations dx dt = JH(t, ε)x, (1.1) where x T = (x 1 , x 2 , . . . , x 2n ) ia a 2n-dimensional vector whose components x k and x n+k are the canonically conjugated variables, J = 0 E n −E n 0 and E n is the n × n identity matrix, H(t, ε) is the real-valued 2n × 2n matrix function which can be represented in the form of the converging series H(t, ε) = H 0 + εH 1 (t) + ε 2 H 2 (t) + . . . , (1.2) where ε is a small parameter. The matrix functions H k (t) (k = 1, 2, . . .) in (1.2) are continuous and periodic with a period T , while H 0 is a constant matrix. Besides, H 0 and H k (t) can depend on some parameters. Equations of the form (1.1) describe dynamical systems with intrinsic periodicity and appear in many branches of science and engineering (see, for example, [11]). Our interest to the system (1.1) arises because such a system occurs in studying the stability of equilibrium solutions in the elliptic restricted many-body problems [3,9]. We are interested also in determination of the boundaries between the domains of stability and instability in the parameter space for the system (1.1). According to the general theory of differential equations with periodic coefficients [11], the behaviour of solutions of the system (1.1) is determined by its characteristic exponents which are continuous functions of ε. And the system may be stable only if none of the eigenvalues of the matrix JH 0 has a positive real part. However, the system being stable for ε = 0 may become unstable even for very small values of ε > 0. So, in order to analyse the stability of system (1.1) we have to calculate its characteristic exponents for ε > 0. Using the method of a small parameter, we can find them in the form of power series in ε as it was done in [3,8], for example. But if we are looking for the stability boundaries the method of infinite determinant turns out to be more effective [4,7]. It was developed first by Bolotin for a restricted class of differential equations, namely, for a system of uncoupled canonical equations [1]. Then Lindh and Likins extended this method for completely damped mechanical systems [6]. But in both cases the stability boundaries were determined numerically. Now there are modern computer algebra systems such as Mathematica [10], for example, that essentially increases our ability in doing symbolic calculations. The main aim of the present paper is to develop the algorithm for analytical calculation of the stability boundaries and to use it for the hamiltonian systems of the second and the fourth order. All calculations are done with the computer algebra system Mathematica.

Properties of the Characteristic Multipliers
The hamiltonian systems of linear differential equations with periodic coefficients and their general properties have been studied quite well. It is known that their characteristic multipliers obey the Liapunov-Poincare reciprocal root Theorem (see [5,11]). It indicates that a hamiltonian system with characteristic multiplier ρ of multiplicity m must also have characteristic multiplier ρ −1 of the same multiplicity m. Besides, coefficients in the equation (1.1) are real-valued functions. Therefore, we can formulate the following theorem, which restricts substantially possible values of the characteristic multipliers for the hamiltonian systems. Theorem 1. If ρ is a characteristic multiplier for the system (1.1) then ρ −1 , ρ, ρ −1 are its characteristic multipliers as well, where ρ is a complex-conjugate value for ρ.
Hence, characteristic multipliers of the hamiltonian system are divided into groups each of which in general case consists of four elements, namely, ρ, ρ −1 , ρ, ρ −1 . If |ρ| = 1 they are situated in the complex plane symmetrically both with respect to the unit circle and with respect to the real axis (see Fig. 1). Of course, theorem 1 will be satisfied if ρ = ρ and two characteristic multipliers are situated on the real axis symmetrically with respect to the unit circle or |ρ| = 1 and they are on the unit circle symmetrically with respect to the real axis. But in any case, if |ρ| = 1 then there exists at least one characteristic multiplier of modulus exceeding unity which means instability of the system. Thus, the hamiltonian system may be stable only if all of its characteristic multipliers are situated on the unit circle in the complex plane. Let us suppose that for ε = 0 all characteristic multipliers of the system (1.1) are different complex numbers of unit modulus, i.e., the system is stable. According to Theorem 1, there exist n pairs of complex-conjugate characteristic multipliers situated on the unit circle in the complex plane symmetrically with respect to the real axis. If for ε > 0 modulus of one characteristic multiplier becomes greater than 1, for example, and it leaves the circle then, according to Theorem 1, three additional characteristic multipliers shown on Fig. 1 must arise. But in this case the number of characteristic multipliers will become greater than 2n what is impossible. Hence, all characteristic multipliers must stay on the circle and the system will be stable for sufficiently small values of ε > 0. However, if for ε = 0 there exist multiple characteristic multipliers, for example, one pair of complex-conjugate characteristic multipliers has a multiplicity 2, then they can leave the circle without disturbing Theorem 1. Indeed, they can move along radial directions and form a configuration shown on Fig. 1. Thus, the existence of multiple characteristic multipliers is the necessary condition of instability of the hamiltonian system (1.1) for ε > 0. It should be noted also that if there exists at least one characteristic multiplier of the system (1.1) whose modulus is not equal to unity for ε = 0 then the system will be unstable for sufficiently small values of ε > 0 as well because its characteristic multipliers are continuous functions of ε.

Computing the Stability Boundaries for the Second Order Hamiltonian System
Let us consider the linear second order differential equation of the form where a, ε are some positive parameters. It arises in studying the stability of equilibrium solutions in the elliptic restricted many-body problems, where the motion of a particle of infinitesimal mass in the gravitational field generated by (N + 1) point particles P 0 , P 1 , . . . , P N is investigated [3,9]. The particles P 1 , . . . , P N have equal masses and move in elliptic orbits about their common center of mass being at any instant of time in the vertices of a regular polygon with N sides. The polygon rotates about its center where the particle P 0 is resting. It is supposed that orbits of the particles P 1 , . . . , P N are situated in the xOy plane of the barycentric inertial frame of reference and its origin is at the point P 0 . Then equation (3.1) describes the disturbed motion of the particle along Oz axis. Obviously, equation (3.1) can be written in the form of the second order hamiltonian system (1.1) with the matrix function The matrix function (3.2) is periodic with the period T = 2π and can be represented and the corresponding series converges in the domain |ε| < 1 for any t. The second order hamiltonian system has two characteristic multipliers ρ 1 and ρ 2 which must satisfy Theorem 1. Hence, they must be situated in the complex plane either on the real axis or on the unit circle (see Fig. 2). The first case corresponds to unstable behaviour of the system because one of its characteristic multipliers has modulus exceeding unity. In the second case |ρ 1 | = |ρ 2 | = 1 and ρ 2 = ρ 1 and the system is stable. Changing parameters of the system, we can force characteristic multipliers to move in the complex plane. But its transition from stable to unstable behaviour and back is possible only via the points ρ = ±1 where the system has multiple characteristic multipliers. Thus, the cases ρ 1 = ρ 2 = 1 and ρ 1 = ρ 2 = −1 correspond to the boundaries between stable and unstable behaviour of the system.
The eigenvalues of the matrix JH 0 are easily found and can be written as λ 1,2 = ±i √ a. The corresponding characteristic multipliers are complex-conjugate numbers of unit modulus. Obviously, the conditions ρ 1 = Although this is a Fourier series for the function z = z(t) of period 4π, it can also be used to obtain the solution with period 2π by setting to zero the Fourier coefficients corresponding to k being an odd integer. By substituting (3.4) into equation (3.1) and equating coefficients of cos( k 2 t) and sin( k 2 t) to zero we obtain the following infinite sequence of equations determining coefficients of the Fourier series (3.4): Of course, it's impossible to calculate a determinant of the infinite matrix. So, in order to find the stability boundaries a = a(ε) we should truncate the infinite subsequences of equations (3.5) -(3.8) after the k-th term, where k is a suitably large number. The corresponding determinant for the system (3.5), for instance, can be written as Equating determinant (3.9) to zero we obtain an algebraic equation giving an approximation for the stability boundary a = a(ε). An exact expression for the boundary is obtained when k → ∞. This approach giving the equation for the stability boundary is known as the infinite determinant method [6].
Determinant (3.9) is most efficiently evaluated from the following recurrence relation which is readily established from (3.9). To start the iterative process we observe that A similar procedure can be followed for the other systems (3.6) -(3.8). For instance, determinant of the system (3.7) is just the same as (3.9) with the first row and column deleted. The recurrence relation is again (3.10) for k ≥ 3, but the starting values are now given by 1)(a − 4).

(3.13)
It is easy to show from (3.10), (3.12) that in order to find the curves (3.13) in the vicinity of the point a = 1 4 k 2 with accuracy o(ε 2n ), it is sufficient to calculate the determinant D k+n . Then we should substitute (3.13) into the expression for D k+n and expand it in powers of ε. Afterwards, equating coefficients of ε k (k = 1, 2, 3, . . . ) to zero, we obtain a system of algebraic equations determining the coefficients a k in the expansion (3.13). As a result we have found the following curves  It should be emphasized that, increasing the order n of determinant D n of the systems (3.5) -(3.7), we'll be able to find only small correction terms in the equations of the stability boundaries (3.13). The terms we have already found in (3.14), (3.15) will be the same.

Characteristic multipliers for the fourth order system
The fourth order hamiltonian system has four characteristic multipliers which must obey Theorem 1 as well. Hence, the system may be stable only if all its characteristic multipliers are complex-valued with unit magnitude a b Figure 3. Motion of the characteristic multipliers in the complex plane.
Geometrically these characteristic multipliers are situated on the unit circle in the complex plane, symmetrically in pairs with respect to the real axis (see Fig. 3a). System (1.1) becomes unstable if at least one characteristic multiplier leaves the circle. But Theorem 1 imposes restrictions on possible motion of the characteristic multipliers in the complex plane. One possibility is shown in Fig. 3b, when two characteristic multipliers, being in the same semi-plane, move toward each other on the circle until their coincidence and then start to move away of the circle along radial directions. It means that the system becomes unstable because modulus of two characteristic multipliers becomes greater than 1. If such a case is realized then system (1.1) has no periodic solutions and we have to calculate its characteristic multipliers explicitly in order to find the stability boundaries. There is also another possibility when two characteristic multipliers, moving on the circle toward each other, coincide in the point ρ = −1 and then continue their motion along the real axis (see Fig. 4). Similar situation can occur if two characteristic multipliers coincide in the point ρ = 1. In both cases the other two characteristic multipliers remain on the unit circle, being symmetrical with respect to the real axis. Again the cases ρ = ±1 correspond to the boundary between stable and unstable behaviour of the system (1.1), similarly as it is in the case of the second order hamiltonian system. The interesting and significant result from this analysis is that for ρ = ±1 the system (1.1) has periodic solutions of period T and 2T respectively and this property may be used for determination of the boundaries between the domains of stability and instability in the parameter space.

Computing the stability boundaries
Let us consider now the hamiltonian system (1.1) of the fourth order with the matrix function where b and ε are some positive parameters. Such a system describes the disturbed motion of the particle in the xOy plane in the elliptic restricted problem of four bodies [2]. The matrix function (4.1) is periodic with the period T = 2π and may be represented in the form (1.2), where and the corresponding series converges in the domain |ε| < 1 for any t. Hence, characteristic exponents for the system are continuous functions of ε. In the case of ε = 0 they are just the eigenvalues of the matrix JH 0 and can be represented as They are distinct pure imaginary numbers if parameter b satisfies the following inequalities Thus, the considered hamiltonian system may be stable for ε > 0 only if parameter b belongs to the intervals (4.2). In other cases there exists at least one characteristic exponent with a positive real part and the system will be unstable for sufficiently small ε. Let us consider the first interval in (4.2). The corresponding intervals for σ 1,2 can be easily found and are given by Obviously, there is only one possibility for the system to have multiple characteristic multipliers. It is just the case σ 2 = 1 2 when two characteristic multipliers ρ 3,4 = exp(±2πσ 2 i) = −1.
The corresponding geometrical configuration is shown in Fig. 4. From the analysis above it follows that the domain of instability can arise only in the vicinity of the point for which σ 2 = 1 2 . Hence, the boundaries between the domains of stability and instability in the b − ε plane are some curves b = b(ε) which are characterized by the presence of periodic solutions with the period 2T = 4π and cross the b-axis at the point (4.3).
In order to find the stability boundaries let us rewrite the system (1.1) with matrix (4.1) in the form of two linear second order differential equations where a dot means the derivative d dt . Now we can attempt to seek a solution of the system (4.4) in the form of Fourier series By substituting (4.5) into equations (4.4) and equating coefficients of cos( k 2 t) and sin( k 2 t) to zero we obtain two infinite sequences of linear homogeneous equations. The first system is for the odd coefficients p 1 , p 3 , ... , p 2k−1 and β 1 , ... , β 2k−1 and is given by The second system of equations determines the coefficients q 1 , q 3 , ... , q 2k−1 and α 1 , ... , α 2k−1 and it can be written as Extracting coefficients of cos(k t) and sin(k t) we can easily obtain two similar sequences of equations for the even coefficients p 2k , q 2k , α 2k , β 2k . For a solution to exist, the corresponding determinants of infinite systems (4.6), (4.7) must vanish, thus determining the stability boundaries in the b − ε plane. These boundaries obviously must reduce to the point b = (6 − √ 33)/4 when ε → 0.
In order to find the stability boundaries b = b(ε) we should truncate the infinite sequences of equations (4.6) -(4.7) after the k-th term, where k is a suitably large number. For k = 3, for example, the corresponding determinants may be written as . Substituting (4.9) into (4.8) and equating coefficients of ε k to zero we get the system of algebraic equations determining the coefficients b k (k = 1, 2, . . .). Solving this system we obtain the stability boundaries in the form where the error term is O(ε 7 ). The second interval in (4.2) may be analysed similarly. As a result we obtain the boundaries of the second domain of instability as Thus, we can formulate the following theorem. It should be emphasized that increasing the order n of the determinants D n of the systems (4.6) -(4.7) we'll be able to find only the higher order coefficients b 7 , b 8 , . . . in the equations of the stability boundaries (4.9). The coefficients b 1 , b 2 , . . . , b 6 are the same as in (4.10), (4.11).

Conclusions
Using the infinite determinant method based on the existence of periodic solutions on the boundaries between the domains of stability and instability in the parameter space we have developed the algorithm for analytical computation of the stability boundaries for the hamiltonian systems of linear differential equations with periodic coefficients. The algorithm has been implemented in the case of the second and the fourth order hamiltonian systems arising in the elliptic restricted many-body problems. The obtained results are in a good agreement with similar results of [2,3], where the calculations are done with smaller accuracy and another method is used.