The second order differential equations are often encountered in engineering and physical problems, and they have been studied for more than a hundred years. The Hill equation is a particular and representative equation among the linear periodic equations, and it receives its name after the work of W. Hill on the lunar perigee. The Hill equation can be used to describe from the simples dynamical systems to the more complex systems; from a child playing on a swing or a spring mass system  to the suspension bridges  or the flow of the light in a 1D photonic crystal  . The complexities of studying the properties of periodic systems often demand a computational approach. However we are not so much interested in the exact shape of the solutions but in the question whether the solution is stable or unstable. By stable solution we mean one that is bounded in the whole real line.
For any linear T-periodic ODE, the knowledge of the state transition matrix in one period; , gives us sufficient information to know the solution in any . The main problem with the determination of the stability of a periodic system solutions lies in the calculation of the state transition matrix, which is almost always impossible to obtain analytically. Hence it is desirable to find some criteria for determining the stability of solutions without the need of obtain them. Zhukovskii  gave one of the most known sufficient condition for stability, namely; let the differential equation with , if
then the solutions of the periodic system are stable.
1Let and be two linearly independent solutions of a periodic differential equation subject to the initial conditions , , and , then the discriminant is equal to , is the minimum period of the differential equation.
Lyapunov in his celebrated work “The general problem of the stability of motion”  developed an approximation of the discriminant1 of the periodic system (see section 2) and obtained a variety of sufficient stability conditions, being the most known the one here presented. Authors such as Borg  Yakubovich   and Tang  had taken advantage of the properties of the Hamiltonian structure in order to obtain sufficient conditions for stability. Hochstadt   and Xu  had made use of the properties of the Sturm-Liouville equation and by means of successive approximations, respectively, in order to develop a new discriminant approximation and then finding sufficient criteria for stability.
The aim of the present work is to collect and graphically display some of the most known and efficient sufficient criteria for the stability of the periodic differential equation , , known as Hill’s equation. There is a vast number of stability criteria, in the literature, they are based on different approaches, here we present an easy explanation of four of those approaches starting with: a) the approximation of the discriminant due to Lyapunov followed by; b) canonical forms (Hamiltonian structure); c) Sturm-Liouville equation properties and ending with d) another discriminant approximation due to Shi  . Each stability criterion will be used in order to find zones in the plane, where the solutions are stable, for the best known forms of Hill’s equation, Mathieu, Meissner and Lyapunov equations i.e., for equations of the form
the scaling factor and of the excitation functions related to the Mathieu and Lyapunov equations were added so its norm be equal to that is, , for each excitation function. The excitation function of the classical Mathieu equation is , if we want the excitation function to have we must multiply the function by a constant a, since then, ; by doing a similar procedure for the classical Meissner ( ) and Lyapunov ( ) we get the equations given above.
This work is structured as follows: Section 2 is dedicated to introduce some basic concepts on periodic differential equations and its solutions; in section 3 we give the discriminant approximation made by Lyapunov and the first two criteria are presented; section 4 is devoted to the study of canonical (Hamiltonian) systems, some properties of such systems solutions are described and four criteria are presented; In section 5 two criteria based on Sturm-Liouville equation properties are exhibit, both criteria follows from solutions proposed by Hochstadt in  ; in section 6 we briefly present a new approximation of the discriminant of the Hill equation obtained by Shi and a criterion do to Xu. In section 7 some conclusions, about the efficiency of the criteria, are presented and a comparison between four stability criteria and the stable zones found by perturbation methods (Linsdtedt-Poincare method) is done. A brief introduction to periodic Hamiltonian systems solutions is given in the appendix.
Consider the linear differential equation with periodic coefficients
where is a periodic piece-wise continuous real function with zero average, i.e. , and are real constants; Equation (1) is known as Hill equation. One can prove that (1) is equivalent to the system
where is a T periodic real function and is a real constant. Throughout the document we will use system (1) or (2) indistinctly.
We say that a system is stable if and only if all its solutions are bounded in the whole real line. It is known that for some values of the parameters and (or ) the Equation (1) has bounded solutions and for some others the solutions grow without bounds. The plane of parameters can be splited into stable zones (where all solutions are stable) and unstable zones (where at least one solution is unstable). Stable and unstable regions are separated by some curves known as transition curves, these transition curves are defined by having at least one periodic or anti-periodic solution  . Unstable regions are called Arnold tongues after remarkable works of V. A. Arnold on Mathieu equation  . It can be proved that, for constant, the stable regions alternate with unstable zones, see Theorem 2.1 below.
By the usual change of variables, and , Equation (1) may be rewritten as
where ; and . Notice that (3) may also be written as the Hamiltonian system
where is the symmetric and T periodic matrix and . The Equation (4) may be taken as a definition of Hamiltonian system and it will be used in section 4.
Let be the state transition matrix of the system (3), from Floquet Theorem it is known that the matrix may be written as a multiplication of three matrices, two of them are time dependent matrices and one constant matrix; one of the time independent matrices is bounded and periodic, and the other is an exponential one, the former gives us information about the phase of the solutions and the latter contains information about the growth of the solutions, see  or  . Such factorization is
where is a real bounded matrix, and R is a constant matrix not necessarily real. Factorization (5) is due to Floquet  . Notice that if we assume then, , and if we define the Monodromy Matrix as then, from Equation (5) one can infer
so, the growth of the system solutions are related with monodromy matrix . For better understanding notice that for all , , and for some positive integer k any solution of (3) can be written as
thus one can conclude the following
Theorem 2.1. Let be the eigenvalues of the monodromy matrix . The system (3) is:
a) Asymptotically stable if and only if all the eigenvalues of the monodromy matrix are ;
b) Stable if and only if all and if any has modulo one then it must be a simple root of the minimal polynomial of ; and
c) Unstable if and only if there is a such that or if all and there exists a such that and it is a multiple root of the minimal polynomial of .
Since Equation (3) can be written as a Hamiltonian system, the matrix is a symplectic matrix2 [  , Chapter 2 § 3.4], in other words the characteristic multipliers of the monodromy matrix are symmetric with respect to the unitary circle, thus the solutions of (3) can not be asymptotically stable. And the only chance for the solutions of (3) are to be bounded or unstable.
Let and be real solutions of (2) subject to the initial conditions
Then the state transition matrix associated to (2) is
and the characteristic multipliers3 are solutions of the characteristic equation
2We say that a matrix is symplectic if the condition is fulfilled. Symplectic matrices are of even order and they form a group 
3The monodromy matrix eigenvalues are usually called characteristic multipliers or just multipliers.
4The fact that the independent term be equal to one follows from the symplectic matrix property that states that the determinant of a symplectic matrix is equal one.
where the independent term is equal to one because of the Liouville theorem4  and the function
is known as the characteristic constant or the discriminant associated to (2) and it plays a fundamental roll on the determination of the stability of the system, the Trace operator in (11) is the sum of the main diagonal entries of a square matrix. By Floquet theory, the condition for the solutions, and , to be stable may be restated as: If , both solutions are stable; if one solution is stable and one unstable; if then one solution is periodic when , (or anti-periodic when ); and the second solution may or may not be periodic (anti-periodic).
The Haupt oscillation Theorem asserts that the real line can be split into alternating intervals known as stability and instability intervals, the former are characterized by and the latter by , the endpoints of the intervals are characterized by values of for which the system (2) has at least one periodic solution. The Haupt oscillation Theorem  is as follows
Theorem 2.2. For the differential equation . There exists an infinite sequence
such that . There exists a second infinite sequence
such that . Both sequences do not have accumulation points, and . These two sequences interlace such that
whenever lies in one of the intervals
if lies in
The proof of the Theorem is based on the fact that the functions and are entire functions of the real variable and its order of
growth for is exactly , see    . Then, and have infinitely many zeros.
Notice that the Theorem 2.2 establishes conditions, in terms of , for the solutions of (2) to be stable; Theorem 2.2 is easily restated so that the stability conditions depend on the parameters and of the system (1). It is well known that the stability of solutions of any Hill equation, of the form (1), can be represented as a stability chart in the plane of parameters , unstable zones are called Arnold tongues. Figure 1 shows the stable and unstable zones, in the plane, of the Mathieu, Meissner and Lyapunov equations.
3. Stability Criteria Based on Lyapunov Approximation
In  Lyapunov proposed a method to approximate the discriminant associated to the equation when , and in  he delved into the study of the approximation properties. In order to be consistent with Theorem 2.2 we will consider . The estimation consist in the alternating series
where each coefficient is defined as a definite n-multiple integral, that is
Figure 1. Stability chart for (a) Mathieu, (b) Meissner and (c) Lyapunov equations, stable areas in white and unstable areas in gray.
notice that the sub-index of each coefficient is equal to the order of the n-multiple integral, that is, the coefficient requires a triple definite integral to be calculated, requires a forth order definite integral to be calculated, and so on.
In  Lyapunov studied in depth the properties of the approximation for such as the convergence of the series and the monotonic decrease of the series coefficients. One can prove that for the coefficients satisfy the inequality
For details of the inequality (20) see appendix A. Using (20) Lyapunov got the following
Criterion 3.1. (Lyapunov 1). If , and if it satisfies the condition
then the solutions of are stable.
Remark 3.2. As far as knowledge of the authors, the only reference in which this criterion is graphically shown is [  , pp. 205]
Before proceeding, we shall call the interval the zeroth instability zone, the intervals and the first and the second instability zone and so on. Similarly the intervals are called the zeroth, first and second stability zones respectively.
Next criterion can be proved by using the above described Lyapunov method  .
Criterion 3.3 (Lyapunov 2). If
where . Then, belongs to the zero stability domain.
Figure 2 shows the stability zones obtained by applying the criteria Lyapunov
Figure 2. Sufficient stability zones of (a) Mathieu, (b) Meissner and (c) Lyapunov equations by using Lyapunov 1 (blue) and Lyapunov 2 (red) criteria.
1 and Lyapunov 2 on Mathieu, Meissner and Lyapunov equations.
The blue and red areas are the stability zones found by criterion 0.3 and criterion 0.5 respectively. From the Lyapunov 2 criterion, one can see that the
parameter a must fulfill the inequality . The stability area found by the
Criterion 0.5 is obtained by calculating the stability area defined by the equation (23) and 500 different values of a, equispaced in the interval , and then, merging the 500 resulting areas.
From Figure 2 one can notice that the blue and red regions are the same for the three equations. This is because the integrals of , and
are all equal to zero, so both criteria depends only on the real constants and .
For more criterion based on Lyapunov method see  where a whole chapter is dedicated to study stability of the characteristic constant .
4. Stability Criteria Based on Properties of Canonical Forms
Consider a second order system in canonical (Hamiltonian) form
where , is a symmetric real periodic matrix and is the skew symmetric matrix defined in section 2.
The state transition matrix of (24) may be expressed as
5In  , Theorem 1.4.2, gives necessary and sufficient conditions for a non-singular real matrix to have a real logarithm.
where or , , is a real matrix function, and K is a real matrix with , see appendix B. It can be proved that matrix K could be defined as5
where the sign is chosen so that K be real. From the factorization (25) we can notice that the stability of the canonical system (24) depends on the exponential matrix and therefore on the matrix K; it can be proved that the matrix K is similar to one of the three matrices shown in Table 1. The set of all possible matrices K could be divided into subsets depending on the determinant sign of each matrix K, following the nomenclature of  , we will say that: a) if ; b) if ; and c) if . So, if then (24) has one stable solution and one unstable solution (the characteristic multipliers are real and distinct from ); if
Table 1. Relation between the subsets , and and the stability of canonical system solutions (24), , and are real positive constants.
6Coexistence refers to the simultaneous existence of two linearly independent solutions of period T or 2T of (1) or (24).
with , see Table 1, then (24) has one periodic solution and one unstable solution , with both solutions are stable and periodic (the characteristic multipliers are or ), these solutions corresponds to coexistence6; and, if then both solutions are bounded (the characteristic multipliers are complex, lie on the unit circle and are distinct from ). All of these properties are summarized in the Table 1.
Remark 4.1. The subset defines the transition boundaries, i.e., defines lines on the plane separating the stable zones from the unstable ones. Moreover if and then (24) has one periodic stable solution and one unstable solution of the form .
Let be the set of real continuous function matrices with and . And, let be a solution of Hamiltonian system (24), where a is a non-zero arbitrary vector and let denote the rotation of the solution x in time T, since it follows that . Let denote the set of matrices such that the rotation over a period T is , and , each of the subsets are disjoint sets  .
By the above mentioned properties, of the subsets , , and , we can say that the symmetric matrix in (24) belongs to one of the subsets , or if and only if the pair defined by the state transition matrix of (24) are: , and K belongs to , or respectively, that is, if and only if and so on.
Let , i.e. the system (24) is stable and and , then the rotation of all the solutions of (24) will be . The sketch of the proof is as follows, by assumption and , from (25) we know that any solution of (24) can be written as where , from form (c) it follows
with and . Thus moves uniformly in a circle, describing and angle in time T. Since then and finally .
Remark 4.2. Let , , be two linear independent solutions of (24), and then so the area of the parallelogram defined by and , do not change do to the influence of and vectors , cannot overlap each other neither . Then the rotation number of and must coincide.
Remark 4.3. The definition of the subsets , and allows us to discriminate between stable (unstable) zones where the solutions rotation have similar properties, see appendix B.
Remember that every linear Hamiltonian system (24) may be defined as
where is the quadratic form
the matrix is the symmetric matrix associated to (24) and is
a solution of the same equation. Defining as the argument of and deriving
it follows from the multiplication . Integrating (30) we get
and then , where and be the smallest and largest eigenvalues . The latter procedure was obtained from  . Following the aforementioned, V. A. Yakubovich proved the following
Criterion 4.4. (Yakubovich 1). Let and be the smallest and largest characteristic values of the matrix in (24). And let
then the Hamiltonian system (24) belongs to the n-th stability region ( ).
One must notice that Hill’s Equation (1) could be written as in (24), setting , and , one gets
thus Hill’s equation could be seen as a Hamiltonian system.
It is easy to verify that (34) may be written as
where c is a positive real constant. Then the smallest and largest eigenvalues of are
Notice that the Yakubovich 1 criterion inequalities could be reformulated as follow. For simplicity consider the Meissner equation, the integral over one period of the functions and are
obviously the inequalities and are fulfilled if and respectively. So, for the case we get that the sufficient stability condition reduces to
and for the case we get
In fact the inequalities (39) and (40) are invariant as long as the periodic excitation function has zero mean, . Then the conditions for Mathieu, Meissner and Lyapunov equations to have stable solutions are:
where the function is the excitation function of each periodic differential equation, see section 1.
Figure 3 shows the stability zones obtained by applying Yakubovich 1
Figure 3. Sufficient stability zones of (a) Mathieu (b) Meissner and (c) Lyapunov equations by using Yakubovich 1 (blue) criterion.
criterion and the definitions of and given in (36) to the Mathieu, Meissner and Lyapunov equation.
For blue stable areas we have calculated the Yakubovich 1 stability criterion for 200 different values of and then we merge the resulting areas.
By Yakubovich 1 criterion, a sufficient condition for (35) to belong to with is
where are the boundaries of instability areas , the boundary of are the sets and . Suppose that for some function
in words, this condition establishes that the distance from to the function is less than the distance from to the boundary of , so the function . Thus inequality (44) is a test for the stability of (24). This test requires the explicit expressions for and . The next three criteria use inequality (44) and the expressions of and to found sufficient conditions for the stability of the Hamiltonian system (24).
Now consider the differential equation
where is a non-negative, non-identically equal to zero, piecewise continuous periodic function with period T. This and the following results were obtained by V. A. Yakubovich  .
Remark 4.5. Notice that in these criterion, the regions which are guaranteed stable are not convex.
Criterion 4.6 (Yakubovich 2). Let
if the following inequality holds,
then the solution of equation (45) is stable and belongs to the n-th zone of stability, .
Criterion 4.7 (Yakubovich 3). Let
if for some , we have the inequality
then the trivial solution of Equation (45) is stable and belongs to the n-th zone of stability.
For Meissner equation, the criteria 0.11 and 0.12 may be rewritten just as we did in the criterion 0.9. Table 2 shows the sufficient conditions for stability of Meissner equation solutions
From Table 2 one can notice that the n-th stability zone ( ) depends on the parameter c. The Yakubovich 2 criterion assumptions, and , define parallelograms whose vertexes depends on the current value of c, for example, for the stability zone the vertexes are: ; ; ; and . And the inequality define triangular sectors. If we denote each parallelogram as and the triangular sector as then, the n-th stability zone given by Yakubovich 2 criterion is
Figure 4 Shows the parallelograms defined by different values of and the triangular sector defined by the inequality for stability zone.
From the Figure 4 we notice that the stability zone obtained by the Yakubovich 2 criterion belongs to the cone , moreover, by following the same procedure we can prove that each subset of , defined by the Yakubovich 2 criterion, is content in the cone
Table 2. Yakubovich 2 and Yakubovich 3 criteria for the Meissner equation, where .
Figure 4. Red continuous line represents the boundary of the stability zone obtained by Yakubovich 2 criterion. Parallelograms in blue discontinuous lines are defined by the inequalities and .
Figure 5. Red continuous line represents the boundary of the stability zone obtained by Yakubovich 3 criterion. Parallelograms in blue discontinuous lines are defined by the inequalities and .
By doing a similar procedure for Yakubovich 3 criterion we obtain the Figure 5.
Notice that the cone is inside the
stable zone obtained with the Yakubovich 3 criterion. It can be proved that the cones defined in (50) are inside of the stable zones (subsets of ) obtained by Yakubovich 3 criterion. Then, for the case of Meissner equation, we can say that the stable zones obtained by Yakubovich 2 criterion are contained in the zones defined by the Yakubovich 3 criterion. This is not the same for the cases of Mathieu and Lyapunov equations, see Figure 6.
Next criterion was developed by Borg  . An alternative proof was made in  .
Criterion 4.8 (Borg). Consider the system (45) and let
Suppose that for some integer n
Then all solutions of Equation (45) are bounded on and the corresponding Hamiltonian in (24) is in the stability domain .
Figure 6 shows the stable zones obtained by the Yakubovich 2 (red), Yakubovich
Figure 6. Stability zones, for (a) Mathieu, (b) Meissner and (c) Lyapunov equations, obtained by Yakunovich 2 (red), Yakubovich 3 (green) and Borg (blue) criteria.
3 (green) and Borg (blue) criteria for the Mathieu, Meissner and Lyapunov equations.
It is worth to notice that Borg criterion uses almost the same statements of the criteria 0.11 and 0.12 but is written instead of and the distance between the function and the constant c is divided by 2. Moreover for Mathieu, Meissner and Lyapunov equations the constant is equal to , so Borg criterion may be rewritten as:
The solutions of Mathieu, Meissner and Lyapunov equations belongs to the stability domain if for some integer n the inequalities
are fulfilled. The integral for Mathieu, Meissner and Lyapunov equations are equal to , and 5.40537 respectively.
5. Stability Criterion Based on Properties of the Sturm Liouville Equation
7These class of systems are called reversible by V. I. Arnold, see  .
In  Hochstadt proved that given the differential equation7
may, under suitable conditions, be solved by assuming solutions of the form
notice that could be seen as a simile of the rotation of the solutions, see section 4.
We know that for the solutions of (55) to be stable the inequality most be fulfilled. If it is possible to establish relations between , , and at and , , , at
since if and are linearly independent solutions of (55) then and are also solutions and can be written as
differentiating both sides
setting and noticing that the solutions and are even
and odd respectively, that is, and . Solving the equations for , , and we arrive to (57), see  .
From (57) one can easily rewrite the stability inequality as
The stability of the solution is then determined by the examination of the signs of , , and , and the number of zeros of , in the open interval . This follows from the Sturm oscillation theorem. Moreover, if , and are positive and is negative and , have no zeros in then and belongs to the first stability zone  . So, if and then the solutions of (55) are bounded. On the other hand, by Equation (56), if and at then and at . In  it is proved that the solutions “rotation” is bounded by
and we get the following
Criterion 5.1 (Hochstadt 1). A sufficient condition for the boundedness of all solutions of the periodic differential equation is
It is well known that, for , the transition curves are defined by points in the plane for which there is at least one periodic (anti-periodic) solution of (55). In  Hochstadt shows that for periodic solutions of (55) the condition
for some positive integer n, must be satisfied. And for anti-periodic solutions the condition
must be satisfied.
If the solutions of (55) are stable then, the condition
must be satisfied. By noticing the above Hochstadt generalized criterion 13 as follows 
Criterion 5.2 (Hochstadt 2). A sufficient condition for the boundedness of all solutions of the periodic differential equation is
Both criteria, Hochstadt 1 and Hochstadt 2, require the derivative of the function , obviously the functions and do not
have any trouble, but the function does, i.e. both criteria can be used to obtain stability zones for Mathieu and Lyapunov equations but, we can not directly apply the criteria for Meissner equation.
For the Meissner case we must use the expansion on Fourier series of
instead of . The change of the excitation function is possible since: if is the transition matrix of and is the transition matrix of then, . We take just the first 20 terms of the expansion.
Figure 7 shows in red and blue the stability areas obtained applying Hochstadt 1 and Hochstadt 2 criteria respectively to Mathieu, Meissner and Lyapunov equations
6. Stability Criterion Based on Shi Approximation
In  the authors gave an approximation for the discriminant associated to the Equation (2). This approximation was obtained by rewriting Equation (2) as:
where is assumed to be positive and differentiable for all t, then it is
Figure 7. Stability zones obtained, for (a) Mathieu, (b) Meissner and (c) Lyapunov equations, by Hochstadt 1 (green) and Hochstadt 2 (blue) criteria.
proved that the system
can be solved and its state transition matrix is
where , and . Then the solution of (68), with , has the form
suppose that and are solutions of (68) and they are subject to the initial conditions , then the discriminant (2) is
and one can obtain by means of successive approximations. For details see  .
The approximation given in  is as follows
Following  and noticing that
the second term of the right hand side of (72) is
for the solutions of (2) to be stable the inequality
must be fulfilled.
From (78) it is not so hard to prove the next
Criterion 6.1 (Xu). Suppose and
where . Then (2) is stable.
Xu criterion, as the Hochstadt 1 and Hochstadt 2 criteria, needs the derivative of the function . In order to avoid troubles we will do the same as in the previous section, i.e. take the first 20 elements of the expansion of and then substitute them instead of .
Applying Xu criterion to Mathieu, Meissner and Lyapunov equations we obtain the Figure 8.
There is a vast number of sufficient conditions for the stability of periodic differential equations, we have just mentioned and explained some of them. For more stability criteria see for example:  where Starzhinskii not only collect sufficient stability criteria for Equation (1) but he consider second order periodic differential equations with dissipation, n-th order systems and some particular cases of the vector equation ; In [  , Chapter 2 § 4.3] some conditions for stability of solutions of (1) that belongs to the first stability zone are presented and some general stability criteria are described; And in [  , Chapters 7 and 8] Yakubovich and Starshinskii, in their outstanding monograph on Linear differential equations with periodic coefficients, they presented and proved some stability criteria starting with the Lyapunov approach (approximation of the characteristic constant), and doing a deep analysis on canonical
Figure 8. Stability zones obtained by Xu criterion for a) Mathieu b) Meissner and c) Lyapunov equations.
We have reviewed some of the most known stability criteria for second order differential equations; these criteria were obtained by four different approaches; by an approximation of the discriminant made by Lyapunov; by properties of canonical (Hamiltonian) systems; by Sturm-Liouville equation properties and by a discriminant approximation due to Shi. We have given an easy explanation of each approach.
From the figures of the present work we can see, by simple inspection, that the best criterion among the ones presented is Hochstadt 2 (criterion 5.2) which is based on some properties of the solutions of the Sturm-Liouville equation and the rotation of the solutions, the second best is Xu’s criterion which is based on the approximation of the discriminant of the Hill’s equation made by Shi, and the discriminant approximation was obtained by means of successive approximations.
We have numerically calculated the percentage of the stability zones, in the ranges and , for each of the equations: Mathieu, ; Meissner, ; and
Lyapunov, . See Table 3.
Table 4 shows us the percentage of the stability zone obtained by using each sufficient condition for stability. Being 100% the whole stability zone of each stability chart in the ranges and .
From Table 4 we can see that the best criterion among the criteria for obtaining the first stability zone of Mathieu, Meissner and Lyapunov equations is Lyapunov 2. Some others conditions that assure the solutions of (1) are stable and belong to the first zone of stability can be found in [  , pag. 61].
The best two criteria for stability of the solution for the whole plane ( and ) are Hochstadt 2 and Xu, both of them are based on an approximation of the discriminant of the Hill equation, the former approximation is based on 2 proposed linear independent solutions of the Sturm-Liouville problem, see section 5, and the later is obtained by means of successive approximations, see section 6.
Perturbation methods, such as strained parameters (Lindstedt-Poincare method) and multiple scales methods, are frequently used for analysing the stability of the periodic differential equations. They are based on the assumption that the variable-coefficient terms are small in some sense. The stability boundaries associated to a Hill equation may be determined by the strained parameters method, that is, assuming that , and then seeking the value of such that the solutions be T or 2T periodic. The general solution and the coefficient are written in terms of powers of (perturbation expansion)
Table 3. Percentage of stable area obtained by numerical calculation in the plane ,
Table 4. Percentage of stability zone obtained by applying each criterion on the Mathieu, Meissner and Lyapunov equation.
then, the series (80) and (81) are substituted into the Hill equation and by grouping terms of like powers of , one obtains a set of recursive differential equations. The initial condition of , i.e. , is the value of where the
Arnold tongues rise, and depends on the excitation function period, ,
. The accuracy of the method depends on how many elements of the series (80) and (81) are obtained. Notice that the same procedure must be done for each transition curve, see  .
In  the calculation of the first transition curves of the classical form of the Mathieu equation
are obtained as
Transition curve associated to the 0th tongue
Transition curves associated to the 1st tongue
Transition curve associated to the 2th tongue
Transition curve associated to the 2th tongu
Figure 9(a) shows the transition curves, associated to the first tongue, obtained by the strained parameters method (Lindstedt-Poincare method) and the actual one. Figures 9(b)-(e) show the stability zones obtained by Yakubovich 1, Hochstadt 2, Borg and Xu criteria respectively, in the ranges and .
Notice that the stability areas found with the strained parameters method are larger than the ones obtained with the sufficient stability criteria (Hochstadt 2, Xu, Yakubovich 1 and Borg). Table 5 shows the percentage of the stability zones obtained by using the four stability criteria and the strained parameters method on the classic Mathieu equation, being 100% the whole stability area in and .
Even though the stability areas found with strained parameters method are larger than the ones obtained with the stability criteria, for the Mathieu equation case (82), the complexities of the method and the set of differential equations that has to be solved, for obtaining the stable zones, make the strained parameters method less suitable than the simplicity of the criteria statements for the stability analysis of periodic differential equations. Besides, we must remember that the strained parameters method has the limitation that it just
Figure 9. Comparison between Strained parameters method and sufficient stability criteria. (a) Unstable zone (UZ) vs. Strained parameters stability zone (SPSZ), (b) stable zone obtained by Yakubovich 1 criterion, (c) stable zone obtained by Hochstadt 2 criterion, (d) stable zone obtained by Borg criterion, (e) stable zone obtained by Xu criterion. TCSP = Transition curves obtained by strained parameters, CSZ = Criterion stable zone.
Table 5. Percentage of stability zones obtained by four stability criteria and the strained parameters method in and .
work for small values of , the sufficient stability criteria, here presented, do not have that limitation.
The first author, Carlos A. Franco acknowledges the financial support of CONACyT and CINVESTAV.
Let and be two linear independent solutions of
subject to the initial condition
The approximation of the characteristic constant
associated to the periodic differential Equation (83), consists in write and the linear independent solutions and as alternating series of powers of
with and . And then solve the recurrence system of equations.
In  Lyapunov studied in depth the properties of the approximation for such as the convergence of the series and the monotonic decrease of the series coefficients. One can prove that for the coefficients satisfy the inequality
this can be demonstrated following  . By the definitions of the coefficients we can restate (88) as
for . Defining and noticing that it can be rewritten as
then (89) is fulfilled if all the coefficients and are positive, and . Rewriting and we have
since the excitation function is positive for all t, one can notice that , , and , these follow from definition. So, and are positive which implies that and are positive and the inequality (89) is fulfilled. For details of the proof see  .
B. A Brief Introduction to Hamiltonian System Solutions
Consider a second order differential equation in canonical form
where , is a symmetric real periodic matrix and is the skew symmetric matrix .
The state transition matrix of (94) may be expressed as
where or , , is a real matrix function, and K is a real matrix with  . The fact that K is a real matrix follows from the eigenvalues properties of the monodromy matrix M, specifically on the multipliers position with respect to the unit circle. For the multipliers of M there exist only three possibilities; a) both eigenvalues are real, positive and reciprocal; b) both eigenvalues real negative and reciprocal; or c) the eigenvalues are complex conjugated numbers on the unit circle. In the first case, there exists a real matrix , in the second a real matrix and in the third case both and are real matrices. So K could be defined as
where the sign is chosen so that K be real.
Remembering that, for periodic systems, and substituting (95) into we have
since . Then where the sign is the same as the one in the definition of K. From (95) one can easily see that . Since and is a real numbers, we can infer that that and .
It follows from matrix similarity that the K matrix may be brought by a similarity transformation to one of following forms:
where and are real numbers and S is a real non-singular matrix. Following the nomenclature of  , let , and be subsets of all admissible matrix K. We will say that the matrix if it has the form (a), if it has the form (b), and if it has the form (c).
From the Floquet factorization (95) one can say that the stability of the solutions of (94) depends on the exponential matrix , and therefore on K i. e. if then (94) has one stable solution and one unstable solution (the characteristic multipliers are real and distinct from ). If with then (94) has one periodic solution and one unstable solution , if both solutions are stable and periodic (the characteristic multipliers are ). And, if then both solutions, of (94), are bounded (the characteristic multipliers are complex, lie on the unit circle and are distinct from ). All the latter properties are summarized in the Table S1.
Let be the set of real continuous function matrices with and . Now, let be a solution of Hamiltonian system (94), where a is a non-zero arbitrary vector and let denote the rotation of the solution x in time T, since it follows that , where n is even if and n is odd if . Let denote the set of matrices such that the rotation over a period T is , , each of these are disjoint sets  .
Let be the set of all symmetric matrices in (94). As we know each matrix determines a unique state transition matrix . By (95)
Table S1. Relation between the subsets , and and the stability of canonical system solutions (24), , and are real positive constants.
the matrix determines the pair , where and . So one can say that . As we have seen then
In words, the set is divided into the subsets , and where the matrices K and , which are defined by the state transition matrix of the Hamiltonian system (94), belong to one of the sets , or and . We say that if and only if and so on.
Let then and , let and be eigenvectors of K such that
we can notice that the rotation of the two linear independent solutions of (94), and , is , these follows from the fact that the rotation of a solution doesn’t depend on the eigenvectors and but on the matrix . For any other solution the rotation is either or (See  § VIII. 1.9, lemma I and II).
Similarly, one can prove that if then the rotation of all the solutions of (94) will be . The sketch of the proof is as follows, by assumption and , from (95) we know that any solution of (94) can be written as where , from form (c) it follows
with and . Thus moves uniformly in a circle, describing and angle in time . Since then and finally .