Comparison on Sufficient Conditions for the Stability of Hill Equation: An Arnold’s Tongues Approach

Carlos A. Franco, Joaquin Collado

Department of Automatic Control, Centro de Investigacion y de Estudios Avanzados del Instituto Politecnico Nacional, Mexico City, Mexico.

DOI: 10.4236/am.2017.810109                

Abstract

It is known that the solutions of a second order linear differential equation with periodic coefficients are almost always analytically impossible to obtain and in order to study its properties we often require a computational approach. In this paper we compare graphically, using the Arnold Tongues, some sufficient criteria for the stability of periodic differential equations. We also present a brief explanation on how the authors, of each criterion, obtained them. And a comparison between four sufficient stability criteria and the stability zones found by perturbation methods is presented.

Keywords

Periodic Differential Equations, Stability Sufficient Criteria, Hill Equation, Hamiltonian, Lyapunov Characteristic Constant

1. Introduction

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 [1] to the suspension bridges [2] or the flow of the light in a 1D photonic crystal [3] . 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; $t\in \left[\mathrm{0,}T\right]$ , gives us sufficient information to know the solution in any $t\in ℝ$ . 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 [4] gave one of the most known sufficient condition for stability, namely; let the differential equation $\ddot{x}+p\left(t\right)x=0$ with $p\left(t+T\right)=p\left(t\right)$ , if

$\frac{n^2{\text{π}}^2}{T^2}\le p\left(t\right)\le \frac{{\left(n+1\right)}^2{\text{π}}^2}{T^2}$ then the solutions of the periodic system are stable.

¹Let $x_1$ and $x_2$ be two linearly independent solutions of a periodic differential equation subject to the initial conditions $x_1\left(t_0\right)=1$ , ${\stackrel{\dot{}}{x}}_1\left(t_0\right)=0$ , $x_2\left(t_0\right)=0$ and ${\stackrel{\dot{}}{x}}_1\left(t_0\right)=1$ , then the discriminant is equal to $x_1\left(T\right)+{\stackrel{\dot{}}{x}}_2\left(T\right)$ , $T$ is the minimum period of the differential equation.

Lyapunov in his celebrated work “The general problem of the stability of motion” [5] developed an approximation of the discriminant¹ 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 [4] Yakubovich [6] [7] and Tang [8] had taken advantage of the properties of the Hamiltonian structure in order to obtain sufficient conditions for stability. Hochstadt [9] [10] and Xu [11] 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 $\ddot{x}+\left(\alpha +\beta p\left(t\right)\right)x=0$ , $p\left(t+T\right)=p\left(t\right)$ , 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 [12] . Each stability criterion will be used in order to find zones in the $\alpha -\beta$ 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

$\begin{array}{l}\text{Mathieu}\ddot{x}+\left(\alpha +\beta \sqrt{2}cos\left(t\right)\right)x=0\\ \text{Meissner}\ddot{x}+\left(\alpha +\beta sign\left(cos\left(t\right)\right)\right)x=0\\ \text{Lyapunov}\ddot{x}+\left(\alpha +\beta \sqrt{\frac{32}{25}}\left(cos\left(t\right)+\frac{3}{4}cos\left(2t\right)\right)\right)x=0\end{array}$

the scaling factor $\sqrt{2}$ and $\sqrt{\frac{32}{25}}$ of the excitation functions related to the Mathieu and Lyapunov equations were added so its ${\mathcal{L}}_2\left[\mathrm{0,2}\text{π}\right]$ norm be equal to $2\text{π}$ that is, ${\displaystyle {\int }_0^{\text{2π}}}{\left|p\left(t\right)\right|}^2\text{d}t=2\text{π}$ , for each excitation function. The excitation function of the classical Mathieu equation is $p\left(t\right)=\cos \left(t\right)$ , if we want the excitation function to have ${\mathcal{L}}_2\left[\mathrm{0,2}\text{π}\right]=2\text{π}$ we must multiply the function $\cos \left(t\right)$ by a constant a, since ${\displaystyle {\int }_0^{\text{2π}}}{\left|a\cos \left(t\right)\right|}^2\text{d}t=a^2\text{π}$ then, $a=\sqrt{2}$ ; by doing a similar procedure for the classical Meissner ( $p\left(t\right)=sign\left(\cos \left(t\right)\right)$ ) and Lyapunov ( $p\left(t\right)=\cos \left(t\right)+\frac{3}{4}\cos \left(2t\right)$ ) 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 [13] ; 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.

2. Preliminaries

Consider the linear differential equation with periodic coefficients

$\ddot{y}+\left(\alpha +\beta p\left(t\right)\right)y=0$ (1)

where $p\left(t+T\right)=p\left(t\right)$ is a periodic piece-wise continuous real function with zero average, i.e. $\frac{1}{T}{\displaystyle {\int }_0^T}p\left(t\right)\text{d}t=0$ , and $\alpha ,\beta$ are real constants; Equation (1) is known as Hill equation. One can prove that (1) is equivalent to the system

$\ddot{y}+\left(\lambda +q\left(t\right)\right)y=0$ (2)

where $q\left(t\right)$ is a T periodic real function and $\lambda$ 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 $\alpha$ and $\beta$ (or $\lambda$ ) the Equation (1) has bounded solutions and for some others the solutions grow without bounds. The plane of parameters $\alpha -\beta$ 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 [14] . Unstable regions are called Arnold tongues after remarkable works of V. A. Arnold on Mathieu equation [15] . It can be proved that, for $\beta ={\beta }_0$ constant, the stable regions alternate with unstable zones, see Theorem 2.1 below.

By the usual change of variables, $x_1=y$ and $x_2=\stackrel{\dot{}}{y}$ , Equation (1) may be rewritten as

$\stackrel{\dot{}}{x}=\underset{\triangleq B\left(t\right)}{\underset{︸}{\left[\begin{array}{cc}0& 1\\ -\alpha -\beta p\left(t\right)& 0\end{array}\right]}}x$ (3)

where $x=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]$ ; and $B\left(t\right)=B\left(t+T\right)\in {ℝ}^{2\times 2}$ . Notice that (3) may also be written as the Hamiltonian system

$\stackrel{\dot{}}{x}=JH\left(t\right)x$ (4)

where $H\left(t\right)$ is the symmetric and T periodic matrix $H\left(t\right)=\left[\begin{array}{cc}\alpha +\beta p\left(t\right)& 0\\ 0& 1\end{array}\right]$ and $J=\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]$ . The Equation (4) may be taken as a definition of Hamiltonian system and it will be used in section 4.

Let $\Phi \left(t\mathrm{,}{t}_0\right)$ be the state transition matrix of the system (3), from Floquet Theorem it is known that the matrix $\Phi \left(t\mathrm{,}{t}_0\right)$ 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 [16] or [14] . Such factorization is

$\Phi \left(t,t_0\right)=F^{-1}\left(t\right)e^{R\left(t-t_0\right)}F\left(t_0\right)$ (5)

where $F\left(t\right)=F\left(t+T\right)$ is a real bounded matrix, and R is a $2\times 2$ constant matrix not necessarily real. Factorization (5) is due to Floquet [17] . Notice that if we assume $t_0=0$ then, $F\left(0\right)=I_2$ , and if we define the Monodromy Matrix $M$ as $M\triangleq \Phi \left(T\mathrm{,0}\right)$ then, from Equation (5) one can infer

$\Phi \left(T,0\right)=M=e^{RT}$ (6)

so, the growth of the system solutions are related with monodromy matrix $M$ . For better understanding notice that for all $t\ge 0$ , $t=kT+\tau$ , $\tau \in \mathrm{0,}T\mathrm{<mo\; stretchy="false">\; )\; </mo>}$ and for some positive integer k any solution of (3) can be written as

$\begin{array}{c}x\left(t\right)=\Phi \left(t,0\right)x\left(0\right)\\ =\Phi \left(kT+\tau ,kT\right)\Phi \left(kT,\left(k-1\right)T\right)\cdots \Phi \left(T,0\right)x\left(0\right)\\ =\Phi \left(kT+\tau ,kT\right){\Phi }^k\left(T,0\right)x\left(0\right)\\ =\Phi \left(\tau ,0\right)M^k\left(T\right)x\left(0\right)\end{array}$ (7)

thus one can conclude the following

Theorem 2.1. Let ${\mu }_i$ be the eigenvalues of the monodromy matrix $M$ . The system (3) is:

a) Asymptotically stable if and only if all the eigenvalues ${\mu }_i$ of the monodromy matrix are $\left|{\mu }_i\right|<1$ ;

b) Stable if and only if all $\left|{\mu }_i\right|\le 1$ and if any ${\mu }_i$ has modulo one then it must be a simple root of the minimal polynomial of $M$ ; and

c) Unstable if and only if there is a ${\mu }_i$ such that $\left|{\mu }_i\right|>1$ or if all $\left|{\mu }_i\right|\le 1$ and there exists a ${\mu }_j$ such that $\left|{\mu }_i\right|=1$ and it is a multiple root of the minimal polynomial of $M$ .

Since Equation (3) can be written as a Hamiltonian system, the matrix $\Phi \left(t\mathrm{,}{t}_0\right)$ is a symplectic matrix² [ [14] , Chapter 2 § 3.4], in other words the characteristic multipliers of the monodromy matrix $\Phi \left(T\mathrm{,0}\right)$ 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 $f\left(t\right)$ and $g\left(t\right)$ be real solutions of (2) subject to the initial conditions

$\begin{array}{cc}f\left(0\right)=1& g\left(0\right)=0\\ \stackrel{\dot{}}{f}\left(0\right)=0& \stackrel{\dot{}}{g}\left(0\right)=1\end{array}$ (8)

Then the state transition matrix associated to (2) is

$\Phi \left(t,0\right)=\left[\begin{array}{cc}f\left(t\right)& g\left(t\right)\\ \stackrel{\dot{}}{f}\left(t\right)& \stackrel{\dot{}}{g}\left(t\right)\end{array}\right]$ (9)

and the characteristic multipliers³ are solutions of the characteristic equation

$\det \left[M-{\mu }_i{I}_2\right]={\mu }_i^2-2A\left(\lambda \right){\mu }_i+1=0$ (10)

²We say that a matrix $B$ is symplectic if the condition $B^{\prime }JB=J$ is fulfilled. Symplectic matrices are of even order and they form a group [18]

³The monodromy matrix eigenvalues ${\mu }_i$ are usually called characteristic multipliers or just multipliers.

⁴The 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 theorem⁴ [19] and the function

$A\left(\lambda \right)=\frac{1}{2}\left(f\left(T\right)+\stackrel{\dot{}}{g}\left(T\right)\right)=\frac{1}{2}Trace\left(M\right)$ (11)

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, $f\left(t\right)$ and $g\left(t\right)$ , to be stable may be restated as: If $A^2\left(\lambda \right)<1$ , both solutions are stable; if $A^2\left(\lambda \right)>1$ one solution is stable and one unstable; if $A^2\left(\lambda \right)=1$ then one solution is periodic when $A\left(\lambda \right)=1$ , (or anti-periodic when $A\left(\lambda \right)=-1$ ); and the second solution may or may not be periodic (anti-periodic).

The Haupt oscillation Theorem asserts that the $\lambda$ real line can be split into alternating intervals known as stability and instability intervals, the former are characterized by $\left|A\left(\lambda \right)\right|<1$ and the latter by $\left|A\left(\lambda \right)\right|>1$ , the endpoints of the intervals $\left|A\left(\lambda \right)\right|=1$ are characterized by values of $\lambda$ for which the system (2) has at least one periodic solution. The Haupt oscillation Theorem [20] is as follows

Theorem 2.2. For the differential equation $\ddot{x}+\left(\lambda +q\left(t\right)\right)x=0$ . There exists an infinite sequence

${\lambda }_0<{\lambda }_1\le {\lambda }_2<{\lambda }_3\le {\lambda }_4<\cdots$ (12)

such that $A\left({\lambda }_i\right)=1$ . There exists a second infinite sequence

${{\lambda }^{\prime }}_1\le {{\lambda }^{\prime }}_2<{{\lambda }^{\prime }}_3\le {{\lambda }^{\prime }}_4<\cdots$ (13)

such that $A\left({{\lambda }^{\prime }}_i\right)=-1$ . Both sequences do not have accumulation points, ${\lambda }_i\underset{i\to \infty }{\to }\infty$ and ${{\lambda }^{\prime }}_i\underset{i\to \infty }{\to }\infty$ . These two sequences interlace such that

${\lambda }_0<{{\lambda }^{\prime }}_1\le {{\lambda }^{\prime }}_2<{\lambda }_1\le {\lambda }_2<{{\lambda }^{\prime }}_3\le {{\lambda }^{\prime }}_4<{\lambda }_3\le {\lambda }_4<\cdots$ (14)

whenever $\lambda$ lies in one of the intervals

$\left({\lambda }_0,{{\lambda }^{\prime }}_1\right),\left({{\lambda }^{\prime }}_2,{\lambda }_1\right),\left({\lambda }_2,{{\lambda }^{\prime }}_3\right),\cdots ,\left|A\left(\lambda \right)\right|<1$ (15)

if $\lambda$ lies in

$\left(-\infty \mathrm{,}{\lambda }_0\right)\mathrm{,}\left({{\lambda }^{\prime }}_1\mathrm{,}{{\lambda }^{\prime }}_2\right),\left({\lambda }_1\mathrm{,}{\lambda }_2\right)\mathrm{,}\cdots ,\left|A\left(\lambda \right)\right|\mathrm{<mo>\; >\; </mo>\; 1}$ (16)

The proof of the Theorem is based on the fact that the functions $A\left(\lambda \right)-1=0$ and $A\left(\lambda \right)+1=0$ are entire functions of the real variable $\lambda$ and its order of

growth for $\left|\lambda \right|\to \infty$ is exactly $\frac{1}{2}$ , see [20] [21] [22] . Then, $A\left(\lambda \right)-1=0$ and $A\left(\lambda \right)+1=0$ have infinitely many zeros.

Notice that the Theorem 2.2 establishes conditions, in terms of $\lambda$ , for the solutions of (2) to be stable; Theorem 2.2 is easily restated so that the stability conditions depend on the parameters $\alpha$ and $\beta$ 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 $\alpha -\beta$ , unstable zones are called Arnold tongues. Figure 1 shows the stable and unstable zones, in the $\alpha -\beta$ plane, of the Mathieu, Meissner and Lyapunov equations.

3. Stability Criteria Based on Lyapunov Approximation

In [5] Lyapunov proposed a method to approximate the discriminant $A\left(\lambda \right)$ associated to the equation $\ddot{x}+\bar{p}\left(t\right)x=0$ when $\bar{p}\left(t\right)\ge 0$ , and in [23] he delved into the study of the approximation properties. In order to be consistent with Theorem 2.2 we will consider $\bar{p}\left(t\right)=\lambda +q\left(t\right)$ . The estimation consist in the alternating series

$A\left(\lambda \right)=A_0-A_1+A_2+\cdots +{\left(-1\right)}^n{A}_n+\cdots$ (17)

where each coefficient is defined as a definite n-multiple integral, that is

$A_0=1,A_1=\frac{T}{2}{\displaystyle {\int }_0^T}\bar{p}\left(t\right)\text{d}t$

$A_2=\frac{1}{2}{\displaystyle {\int }_0^T}\text{d}{t}_1{\displaystyle {\int }_0^{t_1}}\left(T-t_1+t_2\right)\left(t_1-t_2\right)\bar{p}\left(t_1\right)\bar{p}\left(t_2\right)\text{d}{t}_2$ (18)

Figure 1. Stability chart for (a) Mathieu, (b) Meissner and (c) Lyapunov equations, stable areas in white and unstable areas in gray.

$⋮$

$A_n=\frac{1}{2}{\displaystyle {\int }_0^T}\text{d}{t}_1{\displaystyle {\int }_0^{t_1}}d{t}_2\cdots {\displaystyle {\int }_0^{t_{n-1}}}\left(T-t_1+t_n\right)\left(t_1-t_2\right)\cdots \left(t_{n-1}-t_n\right)\cdot \bar{p}\left(t_1\right)\cdots \bar{p}\left(t_n\right)\text{d}{t}_n$ (19)

notice that the sub-index of each coefficient $A_n$ is equal to the order of the n-multiple integral, that is, the coefficient $A_3$ requires a triple definite integral to be calculated, $A_4$ requires a forth order definite integral to be calculated, and so on.

In [23] Lyapunov studied in depth the properties of the approximation for $A\left(\lambda \right)$ such as the convergence of the series and the monotonic decrease of the series coefficients. One can prove that for $n>1$ the coefficients $A_n,n>1$ satisfy the inequality

$A_{n-1}{A}_1-n{A}_n>0$ (20)

For details of the inequality (20) see appendix A. Using (20) Lyapunov got the following

Criterion 3.1. (Lyapunov 1). If $\bar{p}\left(t\right)\ge 0$ , and if it satisfies the condition

$T{\displaystyle {\int }_0^T}\bar{p}\left(t\right)\text{d}t\le 4$ (21)

then the solutions of $\ddot{x}+\bar{p}\left(t\right)x=0$ are stable.

Remark 3.2. As far as knowledge of the authors, the only reference in which this criterion is graphically shown is [ [24] , pp. 205]

Before proceeding, we shall call the interval $\left(-\infty \mathrm{,}{\lambda }_0\right)$ the zeroth instability zone, the intervals $\left({{\lambda }^{\prime }}_1\mathrm{,}{{\lambda }^{\prime }}_2\right)$ and $\left({\lambda }_1\mathrm{,}{\lambda }_2\right)$ the first and the second instability zone and so on. Similarly the intervals $\left({\lambda }_0\mathrm{,}{{\lambda }^{\prime }}_1\right)\mathrm{,}\left({{\lambda }^{\prime }}_2\mathrm{,}{\lambda }_1\right)\mathrm{,}\left({\lambda }_2\mathrm{,}{{\lambda }^{\prime }}_3\right)$ are called the zeroth, first and second stability zones respectively.

Next criterion can be proved by using the above described Lyapunov method [14] .

Criterion 3.3 (Lyapunov 2). If

$\bar{p}\left(t\right)\ge -a^2$ (22)

and

${\displaystyle {\int }_0^T}\bar{p}\left(t\right)\text{d}t\ge \mathrm{0,\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}{\displaystyle {\int }_0^T}\left(\bar{p}\left(t\right)+a^2\right)\le 2acoth\left(\frac{aT}{2}\right)$ (23)

where $0\le a^2\le \frac{{\text{π}}^2}{T^2}$ . Then, $\bar{p}\left(t\right)$ 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 $0<a<\frac{1}{2}$ . 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 $\left[\mathrm{0,0.5}\right[$ , 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 $cos\left(t\right)$ , $sign\left(cos\left(t\right)\right)$ and

$cos\left(t\right)+\frac{3}{4}cos\left(2t\right)$ are all equal to zero, so both criteria depends only on the real constants $\alpha$ and $a$ .

For more criterion based on Lyapunov method see [14] where a whole chapter is dedicated to study stability of the characteristic constant $A\left(\lambda \right)$ .

4. Stability Criteria Based on Properties of Canonical Forms

Consider a second order system in canonical (Hamiltonian) form

$\stackrel{\dot{}}{x}=JH\left(t\right)x$ (24)

where $x=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]$ , $H\left(t\right)$ is a symmetric real periodic matrix $H\left(t+T\right)=H\left(t\right)=H^{\text{T}}\left(t\right)$ and $J$ is the skew symmetric matrix defined in section 2.

The state transition matrix of (24) may be expressed as

$\begin{array}{l}\Phi \left(t,0\right)=F\left(t\right)e^{tK}\\ \Phi \left(0,0\right)=I_2\end{array}$ (25)

⁵In [25] , Theorem 1.4.2, gives necessary and sufficient conditions for a non-singular real matrix to have a real logarithm.

where $F\left(t+T\right)=F\left(t\right)$ or $F\left(t+T\right)=-F\left(t\right)$ , $\det \left(F\left(t\right)\right)=1$ , $F\left(t\right)$ is a real matrix function, and K is a real matrix with $Tr\left(K\right)=0$ , see appendix B. It can be proved that matrix K could be defined as⁵

$K=\frac{1}{T}\ln \left(\pm M\right)$ (26)

where the sign $\pm$ 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 $e^{KT}$ 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 [14] , we will say that: a) $K\in \mathcal{H}$ if $\det \left(K\right)<0$ ; b) $K\in \Pi$ if $\det \left(K\right)=0$ ; and c) $K\in \mathcal{O}$ if $\det \left(K\right)>0$ . So, if $K\in \mathcal{H}$ then (24) has one stable solution and one unstable solution (the characteristic multipliers are real and distinct from $\pm 1$ ); if $K\in \Pi$

Table 1. Relation between the subsets $\mathcal{H}$ , $\Pi$ and $\mathcal{O}$ and the stability of canonical system solutions (24), $ϵ$ , $\lambda$ and $\varphi$ are real positive constants.

⁶Coexistence refers to the simultaneous existence of two linearly independent solutions of period T or 2T of (1) or (24).

with $ϵ\ne 0$ , see Table 1, then (24) has one periodic solution $x\left(t\right)$ and one unstable solution $tx\left(t\right)$ , with $ϵ=0$ both solutions are stable and periodic (the characteristic multipliers are $+1$ or $-1$ ), these solutions corresponds to coexistence⁶; and, if $K\in \mathcal{O}$ then both solutions are bounded (the characteristic multipliers are complex, lie on the unit circle and are distinct from $\pm 1$ ). All of these properties are summarized in the Table 1.

Remark 4.1. The subset $\Pi$ defines the transition boundaries, i.e., $\Pi$ defines lines on the $\alpha -\beta$ plane separating the stable zones from the unstable ones. Moreover if $K\in \Pi$ and $ϵ\ne 0$ then (24) has one periodic stable solution $x\left(t\right)$ and one unstable solution of the form $y\left(t\right)=y_1\left(t\right)+tx\left(t\right)$ .

Let $\Omega$ be the set of real continuous function matrices $F\left(t\right)$ with $F\left(t+T\right)=\pm F\left(t\right)$ and $\det \left(F\left(t\right)\right)=1$ . And, let $x=F\left(t\right)a$ be a solution of Hamiltonian system (24), where a is a non-zero arbitrary vector and let ${\phi }_x$ denote the rotation of the solution x in time T, since $F\left(T\right)=\pm I_2$ it follows that ${\phi }_x=n\text{π}$ . Let ${\Omega }_n$ denote the set of matrices $F\left(t\right)$ such that the rotation over a period T is ${\phi }_x=n\text{π}$ , and $\Omega ={\displaystyle {\cup }_{n=-\infty }^{\infty }{\Omega }_n}$ , each of the subsets ${\Omega }_n$ are disjoint sets [14] .

By the above mentioned properties, of the subsets $\Pi$ , $\mathcal{H}$ , $\mathcal{O}$ and ${\Omega }_n$ , we can say that the symmetric matrix $H\left(t\right)$ in (24) belongs to one of the subsets ${\mathcal{H}}_n\triangleq {\Omega }_n\times \mathcal{H}$ , ${\Pi }_n\triangleq {\Omega }_n\times \Pi$ or ${\mathcal{O}}_n\triangleq {\Omega }_n\times \mathcal{O}$ if and only if the pair $\left(F\left(t\right)\mathrm{,}K\right)$ defined by the state transition matrix of (24) are: $F\left(t\right)\in {\Omega }_n$ , and K belongs to $\mathcal{H}$ , $\Pi$ or $\mathcal{O}$ respectively, that is, $H\left(t\right)\in {\mathcal{H}}_n$ if and only if $\Phi \left(t\mathrm{,0}\right)\in {\mathcal{H}}_n$ and so on.

Let $H\left(t\right)\in {\mathcal{O}}_n$ , i.e. the system (24) is stable and $F\left(t\right)\in {\Omega }_n$ and $K\in \mathcal{O}$ , then the rotation $\phi$ of all the solutions of (24) will be $n\text{π}<\phi <\left(n+1\right)\text{π}$ . The sketch of the proof is as follows, by assumption $F\left(t\right)\in {\Omega }_n$ and $K\in \mathcal{O}$ , from (25) we know that any solution of (24) can be written as $x=F\left(t\right)y$ where $y\left(t\right)=e^{tK}y\left(0\right)$ , from form (c) it follows

$y\left(t\right)=S\left[\begin{array}{cc}\cos \left(\varphi t\right)& -\sin \left(\varphi t\right)\\ \sin \left(\varphi t\right)& \cos \left(\varphi t\right)\end{array}\right]S^{-1}y\left(0\right),0<\varphi <\frac{\text{π}}{T}$ (27)

with $b=S^{-1}y\left(0\right)$ and $\det \left(S\right)=1$ . Thus $S^{-1}y\left(t\right)$ moves uniformly in a circle, describing and angle $\varphi T$ in time T. Since $0<\varphi T<\text{π}$ then $0<{\phi }_y<\text{π}$ and finally $n\text{π}<\phi <\left(n+1\right)\text{π}$ .

Remark 4.2. Let $x=F\left(t\right)a$ , $y=F\left(t\right)b$ , be two linear independent solutions of (24), and $Z=\left[x,y\right]$ then $\det \left(Z\right)=\det F\left(t\right)\det \left(\left[a,b\right]\right)=\det \left([a,b]\right)=\Vert a\Vert \Vert b\Vert \sin \left(\theta \right)$ so the area of the parallelogram defined by $x$ and $y$ , do not change do to the influence of $F\left(t\right)$ and vectors $x$ , $y$ cannot overlap each other neither $\theta \ge \text{π}$ . Then the rotation number of $x$ and $y$ must coincide.

Remark 4.3. The definition of the subsets ${\Pi }_n$ , ${\mathcal{H}}_n$ and ${\mathcal{O}}_n$ 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

$\begin{array}{l}\stackrel{\dot{}}{p}=-\frac{\partial \bar{H}\left(p,q\right)}{\partial q}\\ \stackrel{\dot{}}{q}=\frac{\partial \bar{H}\left(p,q\right)}{\partial p}\end{array}$ (28)

where $\bar{H}\left(p\mathrm{,}q\right)$ is the quadratic form

$\bar{H}\left(p,q\right)=\frac{1}{2}{x}^{\text{T}}Hx$ (29)

the matrix $H$ is the symmetric matrix associated to (24) and $x\left(t\right)=\left[\begin{array}{c}q\left(t\right)\\ p\left(t\right)\end{array}\right]$ is

a solution of the same equation. Defining ${\omega }_x\left(t\right)$ as the argument of $x\left(t\right)$ and deriving

$\begin{array}{l}{\omega }_x\left(t\right)=\arctan \left(\frac{q}{p}\right)\\ \frac{\text{d}{\omega }_x\left(t\right)}{\text{d}t}=\frac{\stackrel{\dot{}}{q}p-\stackrel{\dot{}}{p}q}{p^2+q^2}=\frac{2\bar{H}\left(p,q\right)}{p^2+q^2}\end{array}$ (30)

it follows from the multiplication $\left[\begin{array}{cc}p& -q\end{array}\right]\left[\begin{array}{c}\stackrel{\dot{}}{q}\\ \stackrel{\dot{}}{p}\end{array}\right]=\left[\begin{array}{cc}q& p\end{array}\right]H\left[\begin{array}{c}q\\ p\end{array}\right]$ . Integrating (30) we get

${\omega }_x\left(t\right)={\omega }_x\left(0\right)+{\displaystyle {\int }_0^t}\frac{2\bar{H}\left(p\left(t_1\right),q\left(t_1\right)\right)}{p^2\left(t_1\right)+q^2\left(t_1\right)}\text{d}{t}_1$ (31)

and then

$\varphi ={\displaystyle {\int }_0^T}\frac{2\bar{H}\left(p\left(t_1\right),q\left(t_1\right)\right)}{p^2\left(t_1\right)+q^2\left(t_1\right)}\text{d}{t}_1$ (32)

and then ${\displaystyle {\int }_0^T}{h}_{\min }\left(t\right)\text{d}t\le {\phi }_x\le {\displaystyle {\int }_0^T}{h}_{\max }\left(t\right)\text{d}t$ , where $h_{min}\left(t\right)$ and $h_{max}\left(t\right)$ be the smallest and largest eigenvalues $H\left(t\right)$ . The latter procedure was obtained from [6] . Following the aforementioned, V. A. Yakubovich proved the following

Criterion 4.4. (Yakubovich 1). Let $h_{min}\left(t\right)$ and $h_{max}\left(t\right)$ be the smallest and largest characteristic values of the matrix $H\left(t\right)$ in (24). And let

$n\text{π}<{\displaystyle {\int }_0^T}{h}_{\min }\left(t\right)\le {\displaystyle {\int }_0^T}{h}_{\max }\left(t\right)<\left(n+1\right)\text{π}$ (33)

then the Hamiltonian system (24) belongs to the n-th stability region ${\mathcal{O}}_n$ ( $n=0,1,2,3,\cdots$ ).

One must notice that Hill’s Equation (1) could be written as in (24), setting $y_1=x$ , $y_2=\stackrel{\dot{}}{x}$ and $q\left(t\right)\triangleq \alpha +\beta p\left(t\right)$ , one gets

$J^{\text{'}}\left[\begin{array}{c}{\stackrel{\dot{}}{y}}_1\\ {\stackrel{\dot{}}{y}}_2\end{array}\right]=\underset{H\left(t\right)}{\underset{︸}{\left[\begin{array}{cc}q\left(t\right)& 0\\ 0& 1\end{array}\right]}}\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right]$ (34)

thus Hill’s equation could be seen as a Hamiltonian system.

It is easy to verify that (34) may be written as

$J^{\text{'}}\left[\begin{array}{c}c{\stackrel{\dot{}}{y}}_1\\ {\stackrel{\dot{}}{y}}_2\end{array}\right]=\left[\begin{array}{cc}\frac{1}{c}q\left(t\right)& 0\\ 0& c\end{array}\right]\left[\begin{array}{c}c{y}_1\\ y_2\end{array}\right]$ (35)

where c is a positive real constant. Then the smallest and largest eigenvalues of $H\left(t\right)$ are

$\begin{array}{l}{h}_{\min }\left(t\right)=\{\begin{array}{l}\frac{1}{c}q\left(t\right)\text{if}q\left(t\right)<c^2\\ c\text{if}{c}^2<q\left(t\right)\end{array}\\ h_{\max }\left(t\right)=\{\begin{array}{l}\frac{1}{c}q\left(t\right)\text{if}q\left(t\right)>c^2\\ c\text{if}{c}^2>q\left(t\right)\end{array}\end{array}$ (36)

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 $h_{min}\left(t\right)$ and $h_{max}\left(t\right)$ are

${\displaystyle {\int }_0^{\text{2π}}}{h}_{\min }\left(t\right)\text{d}t=\{\begin{array}{l}\frac{\text{2π}}{c}\alpha \text{if}\alpha +\beta sign\left(\cos \left(t\right)\right)<c^2\\ 2\text{π}c\text{if}{c}^2<\alpha +\beta sign\left(\cos \left(t\right)\right)\end{array}$ (37)

${\displaystyle {\int }_0^{2\pi }}{h}_{\max }\left(t\right)\text{d}t=\{\begin{array}{l}\frac{\text{2π}}{c}\alpha \text{if}\alpha +\beta sign\left(\cos \left(t\right)\right)>c^2\\ 2\text{π}c\text{if}{c}^2>\alpha +\beta sign\left(\cos \left(t\right)\right)\end{array}$ (38)

obviously the inequalities $\alpha +\beta sign\left(\cos \left(t\right)\right)<c^2$ and $\alpha +\beta sign\left(\cos \left(t\right)\right)>c^2$ are fulfilled if $\alpha +\beta <c^2$ and $\alpha -\beta >c^2$ respectively. So, for the $\alpha +\beta <c^2$ case we get that the sufficient stability condition reduces to

$n<\frac{2}{c}\alpha \le 2c<2\left(n+1\right)$ (39)

and for the $\alpha -\beta >c^2$ case we get

$n<2c\le \frac{2}{c}\alpha <2\left(n+1\right)$ (40)

In fact the inequalities (39) and (40) are invariant as long as the periodic excitation function $p\left(t\right)$ has zero mean, $\frac{1}{T}{\displaystyle {\int }_0^T}p\left(t\right)\text{d}t=0$ . Then the conditions for Mathieu, Meissner and Lyapunov equations to have stable solutions are:

$\begin{array}{l}\ln <\frac{2}{c}\alpha \le 2c<2\left(n+1\right)\text{if}\alpha +\beta \underset{t\in \left[0,T\right]}{\max }p\left(t\right)<c^2\\ n<2c\le \frac{2}{c}\alpha <2\left(n+1\right)\text{if}\alpha +\beta \underset{t\in \left[0,T\right]}{\min }p\left(t\right)>c^2\end{array}$ (41)

where the function $p\left(t\right)$ 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 $h_{min}\left(t\right)$ and $h_{max}\left(t\right)$ 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 $0<c\le 3$ and then we merge the resulting areas.

By Yakubovich 1 criterion, a sufficient condition for (35) to belong to ${\mathcal{O}}_n$ with $q\left(t\right)\equiv c^2$ is

$\frac{n\text{π}}{T}<c<\frac{\left(n+1\right)\text{π}}{T}$ (42)

Set

$\begin{array}{l}{{\rho }^{\prime }}_n\left(c\right)=\inf {\displaystyle {\int }_0^T}\left|q\left(t\right)-c^2\right|\text{d}t\text{over}\text{all}q\left(t\right)\in {\Pi }_n\\ {{\rho }^{″}}_{n+1}\left(c\right)=\inf {\displaystyle {\int }_0^T}\left|q\left(t\right)-c^2\right|\text{d}t\text{over}\text{all}q\left(t\right)\in {\Pi }_{n+1}\end{array}$ (43)

where ${\Pi }_n$ are the boundaries of instability areas ${\mathcal{H}}_n$ , the boundary of ${\mathcal{O}}_n$ are the sets ${\Pi }_n$ and ${\Pi }_{n+1}$ . Suppose that for some function $q(t)$

$\inf {\displaystyle {\int }_0^T}\left|q\left(t\right)-c^2\right|\text{d}t<\min \left[{{\rho }^{\prime }}_n\left(c\right),{{\rho }^{″}}_{n+1}\left(c\right)\right]$ (44)

in words, this condition establishes that the distance from $c^2\in {\mathcal{O}}_n$ to the function $q\left(t\right)$ is less than the distance from $c^2$ to the boundary of ${\mathcal{O}}_n$ , so the function $q\left(t\right)\in {\mathcal{O}}_n$ . Thus inequality (44) is a test for the stability of (24). This test requires the explicit expressions for ${{\rho }^{\prime }}_n\left(c\right)$ and ${{\rho }^{″}}_n\left(c\right)$ . The next three criteria use inequality (44) and the expressions of ${{\rho }^{\prime }}_n\left(c\right)$ and ${{\rho }^{″}}_n\left(c\right)$ to found sufficient conditions for the stability of the Hamiltonian system (24).

Now consider the differential equation

$\ddot{x}+q\left(t\right)x=0$ (45)

where $q\left(t\right)$ 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 [7] .

Remark 4.5. Notice that in these criterion, the regions which are guaranteed stable are not convex.

Criterion 4.6 (Yakubovich 2). Let

$q\left(t\right)\ge \frac{n^2{\text{π}}^2}{T^2},\frac{n\text{π}}{T}\le c<\frac{\left(n+1\right)\text{π}}{T}$ (46)

if the following inequality holds,

${\displaystyle {\int }_0^T}\left|q\left(t\right)-c^2\right|\text{d}t\le 2c\left(n+1\right)\cot \frac{Tc}{2\left(n+1\right)}$ (47)

then the solution of equation (45) is stable and $q\left(t\right)$ belongs to the n-th zone of stability, $n=0,1,2,\cdots$ .

Criterion 4.7 (Yakubovich 3). Let

$q\left(t\right)\le \frac{{\left(n+1\right)}^2{\text{π}}^2}{T^2},\frac{n\text{π}}{T}\le c<\frac{\left(n+1\right)\text{π}}{T}$ (48)

if for some $n=0,1,2,\cdots$ , we have the inequality

${\displaystyle {\int }_0^T}\left|q\left(t\right)-c^2\right|\text{d}t\le c\left(Tc-n\text{π}\right)$ (49)

then the trivial solution of Equation (45) is stable and $q\left(t\right)$ 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 ( $n=0,1,\cdots$ ) depends on the parameter c. The Yakubovich 2 criterion assumptions, $\frac{n}{2}\le c<\frac{n+1}{2}$ and $\zeta \le 2c\left(n+1\right)cot\left(\frac{\text{π}c}{\left(n+1\right)}\right)$ , define parallelograms whose vertexes depends on the current value of c, for example, for the stability zone ${\mathcal{O}}_1$ the vertexes are: $\left(\alpha ,\beta \right)=\left(c^2+\frac{2}{\text{π}}c\cot \left(\frac{\text{π}c}{2}\right),0\right)$ ; $\left(\alpha ,\beta \right)=\left(c^2-\frac{2}{\text{π}}c\cot \left(\frac{\text{π}c}{2}\right),0\right)$ ; $\left(\alpha ,\beta \right)=\left(c^2+\frac{2}{\text{π}}c\cot \left(\frac{\text{π}c}{2}\right),\frac{2}{\text{π}}c\cot \left(\frac{\pi c}{2}\right)\right)$ ; and $\left(\alpha ,\beta \right)=\left(c^2-\frac{2}{\text{π}}c\cot \left(\frac{\text{π}c}{2}\right),\frac{2}{\text{π}}c\cot \left(\frac{\text{π}c}{2}\right)\right)$ . And the inequality $\alpha -\beta \ge \frac{n^2}{4}$ define triangular sectors. If we denote each parallelogram as $\mathcal{P}\left(c\right)$ and the triangular sector as $\mathcal{T}\left(n\right)$ then, the n-th stability zone given by Yakubovich 2 criterion is

${\displaystyle \underset{\frac{n}{2}\le c<\frac{n+1}{2}}{\cup }\mathcal{P}\left(c\right)\cap \mathcal{T}\left(n\right)}\subset {\mathcal{O}}_n$

Figure 4 Shows the parallelograms defined by different values of $c$ and the triangular sector defined by the inequality $\alpha -\beta \ge \frac{n^2}{4}$ for $n=1$ stability zone.

From the Figure 4 we notice that the stability zone obtained by the Yakubovich 2 criterion belongs to the cone $V=\left\{\left(\alpha ,\beta \right)|\alpha -\beta \ge \frac{1}{4}\text{and}\alpha +\beta \le 1\right\}$ , moreover, by following the same procedure we can prove that each subset of ${\mathcal{O}}_n$ , defined by the Yakubovich 2 criterion, is content in the cone

Table 2. Yakubovich 2 and Yakubovich 3 criteria for the Meissner equation, where $\zeta =\left(\left|\alpha +\beta -c^2\right|+\left|\alpha -\beta -c^2\right|\right)\text{π}$ .

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 $\frac{1}{2}\le c<1$ and $\zeta \le 4ccot\left(\frac{\text{π}c}{2}\right)$ .

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 $\frac{1}{2}\le c<1$ and $\zeta \le c\left(2c-1\right)$ .

$V_n=\left\{\left(\alpha ,\beta \right)|\alpha -\beta \ge \frac{n^2}{4}\text{and}\alpha +\beta \le \frac{{\left(n+1\right)}^2}{4}\right\}$ (50)

By doing a similar procedure for Yakubovich 3 criterion we obtain the Figure 5.

Notice that the cone $V=\left\{\left(\alpha ,\beta \right)|\alpha -\beta \ge \frac{1}{4}\text{and}\alpha +\beta \le 1\right\}$ is inside the

stable zone obtained with the Yakubovich 3 criterion. It can be proved that the cones $V_n$ defined in (50) are inside of the stable zones (subsets of ${\mathcal{O}}_n$ ) 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 [4] . An alternative proof was made in [14] .

Criterion 4.8 (Borg). Consider the system (45) and let

$q_{av}=\frac{1}{T}{\displaystyle {\int }_0^T}q\left(t\right)\text{d}t$ (51)

Suppose that for some integer n

$\frac{n^2{\text{π}}^2}{T^2}<q_{av}<\frac{{\left(n+1\right)}^2{\text{π}}^2}{T^2}$ (52)

${\displaystyle {\int }_0^T}\left|q\left(t\right)-q_{av}\right|\text{d}t<2\sqrt{q_{av}}\left(T\sqrt{q_{av}}-n\text{π}\right)\mathrm{,}\text{if}n\ge \mathrm{1,}$ (53)

${\displaystyle {\int }_0^T}\left|q\left(t\right)-q_{av}\right|\text{d}t\le 4\sqrt{q_{av}}\left(n+1\right)\cot \left(\frac{T\sqrt{q_{av}}}{2\left(n+1\right)}\right)$ (54)

Then all solutions of Equation (45) are bounded on $\left(-\infty \mathrm{,}+\infty \right)$ and the corresponding Hamiltonian in (24) is in the stability domain ${\mathcal{O}}_n$ .

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 $q_{av}$ is written instead of $c^2$ and the distance between the function $q\left(t\right)$ and the constant c is divided by 2. Moreover for Mathieu, Meissner and Lyapunov equations the constant $q_{av}$ is equal to $\alpha$ , so Borg criterion may be rewritten as:

The solutions of Mathieu, Meissner and Lyapunov equations belongs to the stability domain ${\mathcal{O}}_n$ if for some integer n the inequalities

$\frac{n^2}{4}<\alpha <\frac{{\left(n+1\right)}^2}{4}$

$\beta {\displaystyle {\int }_0^{\text{2π}}}\left|p\left(t\right)\right|\text{d}t<2\sqrt{\alpha }\left(2\sqrt{\alpha }-n\right)\text{π},\text{if}n\ge 1,$

$\beta {\displaystyle {\int }_0^{\text{2π}}}\left|p\left(t\right)\right|\text{d}t\le 4\sqrt{\alpha }\left(n+1\right)\cot \left(\frac{\text{π}\sqrt{\alpha }}{\left(n+1\right)}\right)$

are fulfilled. The integral ${\displaystyle {\int }_0^{\text{2π}}}\left|p\left(t\right)\right|\text{d}t$ for Mathieu, Meissner and Lyapunov equations are equal to $4\sqrt{2}$ , $2\text{π}$ and 5.40537 respectively.

5. Stability Criterion Based on Properties of the Sturm Liouville Equation

⁷These class of systems are called reversible by V. I. Arnold, see [26] .

In [13] Hochstadt proved that given the differential equation⁷

$\ddot{x}+\left(\lambda +p\left(t\right)\right)y=0,\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}p\left(t+T\right)=p\left(t\right),p\left(t\right)=p\left(-t\right)$ (55)

may, under suitable conditions, be solved by assuming solutions of the form

$\begin{array}{l}{x}_1=A_1\left(t\right)\cos {\phi }_1\left(t\right)\hfill \\ x_2=A_2\left(t\right)\sin {\phi }_2\left(t\right)\hfill \\ {\stackrel{\dot{}}{x}}_1=-A_1\left(t\right)\sqrt{\lambda +p\left(t\right)}\sin {\phi }_1\left(t\right)\hfill \\ {\stackrel{\dot{}}{x}}_2=A_2\left(t\right)\sqrt{\lambda +p\left(t\right)}\cos {\phi }_2\left(t\right)\hfill \end{array}\}$ (56)

where

${\stackrel{\dot{}}{\phi }}_n\left(t\right)=\sqrt{\lambda +p\left(t\right)}+{\left(-1\right)}^n\frac{1}{4}\frac{\stackrel{\dot{}}{p}\left(t\right)}{\lambda +p\left(t\right)}\sin 2{\phi }_n\left(t\right),n=1,2$

${\stackrel{\dot{}}{A}}_1\left(t\right)=-A_1\left(t\right)\frac{\stackrel{\dot{}}{p}\left(t\right)}{2\lambda +2p\left(t\right)}{\sin }^2{\phi }_1(t)$

${\stackrel{\dot{}}{A}}_2\left(t\right)=-A_2\left(t\right)\frac{\stackrel{\dot{}}{p}\left(t\right)}{2\lambda +2p\left(t\right)}{\cos }^2{\phi }_2(t)$

$A_n\left(0\right)=1,{\phi }_n\left(0\right)=0,n=1,2$

notice that $\phi$ 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 $\left|x_1\left(T\right)+{\stackrel{\dot{}}{x}}_2\left(T\right)\right|<2$ most be fulfilled. If $p\left(t\right)=p\left(-t\right)$ it is possible to establish relations between $x_1$ , ${\stackrel{\dot{}}{x}}_1$ , $x_2$ and ${\stackrel{\dot{}}{x}}_2$ at $t=T$ and $x_1$ , ${\stackrel{\dot{}}{x}}_1$ , $x_2$ , ${\stackrel{\dot{}}{x}}_2$ at

$t=\frac{T}{2}$ as follows

$\begin{array}{l}{x}_1\left(T\right)=2x_1\left(\frac{T}{2}\right){\stackrel{\dot{}}{x}}_2\left(\frac{T}{2}\right)-1=1+2{\stackrel{\dot{}}{x}}_1\left(\frac{T}{2}\right)x_2\left(\frac{T}{2}\right)\hfill \\ {\stackrel{\dot{}}{x}}_1\left(T\right)=2x_1\left(\frac{T}{2}\right){\stackrel{\dot{}}{x}}_1\left(\frac{T}{2}\right)\hfill \\ x_2\left(T\right)=2x_2\left(\frac{T}{2}\right){\stackrel{\dot{}}{x}}_2\left(\frac{T}{2}\right)\hfill \\ {\stackrel{\dot{}}{x}}_2\left(T\right)=x_1\left(T\right)\hfill \end{array}\}$ (57)

since if $x_1\left(t\right)$ and $x_2\left(t\right)$ are linearly independent solutions of (55) then $x_1\left(t+T\right)$ and $x_2\left(t+T\right)$ are also solutions and can be written as

$\begin{array}{l}{x}_1\left(t+T\right)=x_1\left(t\right)x_1\left(T\right)+x_2\left(t\right){\stackrel{\dot{}}{x}}_1\left(T\right)\\ x_2\left(t+T\right)=x_1\left(t\right)x_2\left(T\right)+x_2\left(t\right){\stackrel{\dot{}}{x}}_2\left(T\right)\end{array}$ (58)

differentiating both sides

$\begin{array}{l}{\stackrel{\dot{}}{x}}_1\left(t+T\right)={\stackrel{\dot{}}{x}}_1\left(t\right)x_1\left(T\right)+{\stackrel{\dot{}}{x}}_2\left(t\right){\stackrel{\dot{}}{x}}_1\left(T\right)\\ {\stackrel{\dot{}}{x}}_2\left(t+T\right)={\stackrel{\dot{}}{x}}_1\left(t\right)x_2\left(T\right)+{\stackrel{\dot{}}{x}}_2\left(t\right){\stackrel{\dot{}}{x}}_2\left(T\right)\end{array}$ (59)

setting $t=-\frac{T}{2}$ and noticing that the solutions $x_1\left(t\right)$ and $x_2\left(t\right)$ are even

and odd respectively, that is, $x_1\left(t\right)=x_1\left(-t\right)$ and $x_2\left(t\right)=-x_2\left(-t\right)$ . Solving the equations for $x_1\left(T\right)$ , ${\stackrel{\dot{}}{x}}_1\left(T\right)$ , $x_2\left(T\right)$ and ${\stackrel{\dot{}}{x}}_2\left(T\right)$ we arrive to (57), see [20] .

From (57) one can easily rewrite the stability inequality $\left|x_1\left(T\right)+{\stackrel{\dot{}}{x}}_2\left(T\right)\right|<2$ as

$\left|4x_1\left(\frac{T}{2}\right){\stackrel{\dot{}}{x}}_2\left(\frac{T}{2}\right)-2\right|<2$ (60)

or

$\left|2+4{\stackrel{\dot{}}{x}}_1\left(\frac{T}{2}\right)x_2\left(\frac{T}{2}\right)\right|<2$ (61)

The stability of the solution is then determined by the examination of the signs of $x_1\left(\frac{T}{2}\right)$ , $x_2\left(\frac{T}{2}\right)$ , ${\stackrel{\dot{}}{x}}_1\left(\frac{T}{2}\right)$ and ${\stackrel{\dot{}}{x}}_2\left(\frac{T}{2}\right)$ , and the number of zeros of $x_1\left(t\right)$ , $x_2\left(t\right)$ in the open interval $\left(\mathrm{0,}\frac{T}{2}\right)$ . This follows from the Sturm oscillation theorem. Moreover, if $x_1\left(\frac{T}{2}\right)$ , $x_2\left(\frac{T}{2}\right)$ and ${\stackrel{\dot{}}{x}}_2\left(\frac{T}{2}\right)$ are positive and ${\stackrel{\dot{}}{x}}_1\left(\frac{T}{2}\right)$ is negative and $x_1\left(t\right)$ , $x_2\left(t\right)$ have no zeros in $t\in \left(\mathrm{0,}\frac{T}{2}\right)$ then $x_1$ and $x_2$ belongs to the first stability zone [27] . So, if $x_1\left(T/2\right)\ge 0$ and ${\stackrel{\dot{}}{x}}_2\left(T/2\right)\ge 0$ then the solutions of (55) are bounded. On the other hand, by Equation (56), if ${\phi }_1\left(t\right)<\frac{\text{π}}{2}$ and ${\phi }_2\left(t\right)\mathrm{<mo>\; <\; </mo>}\frac{\text{π}}{2}$ at $0\le t\le T/2$ then $x_1\left(t\right)\ge 0$ and ${\stackrel{\dot{}}{x}}_2\left(t\right)\ge 0$ at $0\le t\le T/2$ . In [9] it is proved that the solutions “rotation” is bounded by

$\begin{array}{c}\left|{\phi }_n\left(\frac{T}{2}\right)\right|=\left|{\displaystyle {\int }_0^{\frac{T}{2}}}\left[\sqrt{\lambda +p\left(t\right)}+{\left(-1\right)}^n\frac{1}{4}\frac{\stackrel{\dot{}}{p}\left(t\right)}{\lambda +p\left(t\right)}\sin 2\phi \left(t\right)\right]\text{d}t\right|\\ \le {\displaystyle {\int }_0^{\frac{T}{2}}}\sqrt{\lambda +p\left(t\right)}\text{d}t+\frac{1}{4}{\displaystyle {\int }_0^{\frac{T}{2}}}\left|\frac{\stackrel{\dot{}}{p}\left(t\right)}{\lambda +p\left(t\right)}\right|\text{d}t\le \frac{\text{π}}{2}\end{array}$ (62)

and we get the following

Criterion 5.1 (Hochstadt 1). A sufficient condition for the boundedness of all solutions of the periodic differential equation $\ddot{x}+q\left(t\right)x=0$ is

${\displaystyle {\int }_0^{\frac{T}{2}}}{\left(q\left(t\right)\right)}^{\frac{1}{2}}\text{d}t+\frac{1}{4}{\displaystyle {\int }_0^{\frac{T}{2}}}\left|\frac{\stackrel{\dot{}}{q}\left(t\right)}{q\left(t\right)}\right|\text{d}t\le \frac{\text{π}}{2}$ (63)

It is well known that, for $q\left(t\right)=\alpha +\beta p\left(t\right)$ , the transition curves are defined by points in the $\alpha -\beta$ plane for which there is at least one periodic (anti-periodic) solution of (55). In [13] Hochstadt shows that for periodic solutions of (55) the condition

${\phi }_n\left(T\right)=2n\text{π}$ (64)

for some positive integer n, must be satisfied. And for anti-periodic solutions the condition

${\phi }_n\left(T\right)=\left(2n+1\right)\text{π}$ (65)

must be satisfied.

If the solutions of (55) are stable then, the condition

$n\text{π}\le {\phi }_n\left(T\right)\le \left(n+1\right)\text{π}$ (66)

must be satisfied. By noticing the above Hochstadt generalized criterion 13 as follows [10]

Criterion 5.2 (Hochstadt 2). A sufficient condition for the boundedness of all solutions of the periodic differential equation $\ddot{x}+q\left(t\right)x=0$ is

$n\text{π}\le {\displaystyle {\int }_0^T}{\left(q\left(t\right)\right)}^{\frac{1}{2}}\text{d}t-\frac{1}{4}{\displaystyle {\int }_0^T}\left|\frac{\stackrel{\dot{}}{q}\left(t\right)}{q\left(t\right)}\right|\text{d}t\le {\displaystyle {\int }_0^T}{\left(q\left(t\right)\right)}^{\frac{1}{2}}\text{d}t+\frac{1}{4}{\displaystyle {\int }_0^T}\left|\frac{\stackrel{\dot{}}{q}\left(t\right)}{q\left(t\right)}\right|\text{d}t\le \left(n+1\right)\text{π}$

Both criteria, Hochstadt 1 and Hochstadt 2, require the derivative of the function $p\left(t\right)$ , obviously the functions $cos\left(t\right)$ and $cos\left(t\right)+\frac{3}{4}cos\left(2t\right)$ do not

have any trouble, but the function $sign\left(cos\left(t\right)\right)$ 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

$sign\left(sin\left(t\right)\right)={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}\frac{4}{\left(2n+1\right)\text{π}}sin\left(\left(2n+1\right)t\right)$ (67)

instead of $sign\left(cos\left(t\right)\right)$ . The change of the excitation function is possible since: if $\Phi \left(t\mathrm{,0}\right)$ is the transition matrix of $\ddot{x}+\left(\alpha +\beta sign\left(\cos \left(t\right)\right)\right)x=0$ and $\bar{\Phi }\left(t\mathrm{,0}\right)$ is the transition matrix of $\ddot{x}+\left(\alpha +\beta sign\left(\cos \left(t\right)\right)\right)x=0$ then, $Trace\left(\Phi \left(T,0\right)\right)=Trace\left(\bar{\Phi }\left(T,0\right)\right)$ . 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 [12] the authors gave an approximation for the discriminant associated to the Equation (2). This approximation was obtained by rewriting Equation (2) as:

$\left[\begin{array}{c}{\stackrel{\dot{}}{x}}_1\\ {\stackrel{\dot{}}{x}}_2\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -\lambda q\left(t\right)& \frac{\stackrel{\dot{}}{q}\left(t\right)}{2q\left(t\right)}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]+\left[\begin{array}{c}0\\ -\frac{\stackrel{\dot{}}{q}\left(t\right)}{2q\left(t\right)}{x}_2\end{array}\right]$ (68)

where $q\left(t\right)$ 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

$\left[\begin{array}{c}{\stackrel{\dot{}}{y}}_1\\ {\stackrel{\dot{}}{y}}_2\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -\lambda q\left(t\right)& \frac{\stackrel{\dot{}}{q}\left(t\right)}{2q\left(t\right)}\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right]$ (69)

can be solved and its state transition matrix is

$X\left(t\right)=\left[\begin{array}{cc}\cos \left(w\Phi \left(t\right)\right)& \frac{1}{w\varphi \left(0\right)}\sin \left(w\Phi \left(t\right)\right)\\ -w\varphi \left(t\right)\sin \left(w\Phi \left(t\right)\right)& \frac{\varphi \left(t\right)}{\varphi \left(0\right)}\cos \left(w\Phi \left(t\right)\right)\end{array}\right]$ (70)

where $w=\sqrt{\lambda }$ , $\varphi =\sqrt{q\left(t\right)}$ and $\Phi \left(t\right)={\displaystyle {\int }_0^t}\varphi \left(\tau \right)\text{d}\tau$ . Then the solution of (68), with ${\left[\begin{array}{cc}{x}_1\left(0\right)& x_2\left(0\right)\end{array}\right]}^{\prime }={\left[\begin{array}{cc}{x}_{10}& x_{20}\end{array}\right]}^{\prime }$ , has the form

$\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=X\left(t\right)\left[\begin{array}{c}{x}_{10}\\ x_{20}\end{array}\right]+{\displaystyle {\int }_0^t}X\left(t\right)X^{-1}\left(\tau \right)\left[\begin{array}{c}0\\ -\frac{\stackrel{\dot{}}{q}\left(\tau \right)}{2q\left(\tau \right)}{x}_2\left(\tau \right)\end{array}\right]\text{d}\tau$ (71)

suppose that $\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]$ and $\left[\begin{array}{c}{\bar{x}}_1\\ {\bar{x}}_2\end{array}\right]$ are solutions of (68) and they are subject to the initial conditions $\left[\begin{array}{c}{x}_1\left(0\right)\\ x_2\left(0\right)\end{array}\right]=\left[\begin{array}{c}1\\ 0\end{array}\right]$ , $\left[\begin{array}{c}{\bar{x}}_1\left(0\right)\\ {\bar{x}}_2\left(0\right)\end{array}\right]=\left[\begin{array}{c}0\\ 1\end{array}\right]$ then the discriminant (2) is

$A\left(\lambda \right)=x_1\left(T\right)+{\bar{x}}_2\left(T\right)$ and one can obtain $2A\left(\lambda \right)$ by means of successive approximations. For details see [12] .

The approximation given in [12] is as follows

$\begin{array}{c}A\left(\lambda \right)=2\cos {\displaystyle {\int }_0^T}\sqrt{\lambda q\left(t\right)}\text{d}t\\ +{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}}\frac{1}{2^{4n-1}}{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^{t_1}}\cdots {\displaystyle {\int }_0^{t_{2n-1}}}\cos \psi \left(t_1{t}_2\cdots t_{2n}\right){\displaystyle {\prod }_{i=1}^{2n}}\frac{\stackrel{\dot{}}{q}\left(t_i\right)}{q\left(t_i\right)}\text{d}{t}_{2n}\cdots \text{d}{t}_1\\ =2\cos {\displaystyle {\int }_0^T}\sqrt{\lambda q\left(t\right)}\text{d}t+{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}}{\Delta }_n\end{array}$ (72)

where

$\begin{array}{c}\psi \left(t_1{t}_2\cdots t_{2n}\right)={\displaystyle {\int }_0^T}\sqrt{\lambda q\left(s\right)}\text{d}s-{\displaystyle {\int }_{t_2}^{t_1}}\sqrt{\lambda q\left(s\right)}\text{d}s-{\displaystyle {\int }_{t_4}^{t_3}}\sqrt{\lambda q\left(s\right)}\text{d}s\\ -\cdots -{\displaystyle {\int }_{t_{2n}}^{t_{2n-1}}}\sqrt{\lambda q\left(s\right)}\text{d}s\end{array}$

Following [11] and noticing that

${\displaystyle {\int }_0^T}{\displaystyle {\int }_0^{t_1}}\cdots {\displaystyle {\int }_0^{t_{k-1}}}f\left(t_1\right)f\left(t_2\right)\cdots f\left(t_k\right)\text{d}{t}_k\cdots \text{d}{t}_1=\frac{1}{k!}{\left({\displaystyle {\int }_0^T}f\left(t\right)\text{d}t\right)}^k$ (73)

and

$2{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}}\frac{x^{2n}}{\left(2n\right)!}=e^x+e^{-x}-2$ (74)

the second term of the right hand side of (72) is

$\left|{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}}{\Delta }_n\right|\le e^{\frac{v}{4}}+e^{-\frac{v}{4}}-2$ (75)

where

$\nu \triangleq {\displaystyle {\int }_0^T}\frac{\left|\stackrel{\dot{}}{q}\left(t\right)\right|}{q\left(t\right)}\text{d}t$ (76)

then

$\left|A\left(\lambda \right)\right|\le 2\left|cos{\displaystyle {\int }_0^T}\sqrt{\lambda q\left(t\right)}\text{d}t\right|+e^{\frac{v}{4}}+e^{-\frac{v}{4}}-2$ (77)

for the solutions of (2) to be stable the inequality

$2\left|cos{\displaystyle {\int }_0^T}\sqrt{\lambda q\left(t\right)}\text{d}t\right|+e^{\frac{v}{4}}+e^{-\frac{v}{4}}-2<2$ (78)

must be fulfilled.

From (78) it is not so hard to prove the next

Criterion 6.1 (Xu). Suppose $v={\displaystyle {\int }_0^{\text{π}}}\frac{\left|\stackrel{\dot{}}{q}\left(\frac{T}{\text{π}}\tau \right)\right|}{q\left(\frac{T}{\text{π}}\tau \right)}\frac{T}{\text{π}}\text{d}\tau$ and

$\begin{array}{l}\lambda \in \left({\left(\frac{n\text{π}+{\cos }^{-1}\left({\phi }_0\right)}{{\displaystyle {\int }_0^T}\sqrt{q\left(t\right)}\text{d}t}\right)}^2,{\left(\frac{\left(n+1\right)\text{π}-{\cos }^{-1}\left({\phi }_0\right)}{{\displaystyle {\int }_0^T}\sqrt{q\left(t\right)}\text{d}t}\right)}^2\right)\\ n=0,1,2,\cdots \end{array}$ (79)

where ${\phi }_0=\frac{1}{2}\left(4-e^{v/4}-e^{v/4}\right)$ . Then (2) is stable.

Xu criterion, as the Hochstadt 1 and Hochstadt 2 criteria, needs the derivative of the function $q\left(t\right)$ . 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 $sign\left(sin\left(t\right)\right)$ and then substitute them instead of $sign\left(cos\left(t\right)\right)$ .

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: [4] 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 $\ddot{y}+\mu P\left(t\right)y=0$ ; In [ [28] , 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 [ [14] , 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.

(Hamiltonian) systems.

7. Conclusions

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 $\alpha \in \left[\mathrm{0,3}\right]$ and $\beta \in \left[\mathrm{0,1.5}\right]$ , for each of the equations: Mathieu, $\ddot{x}+\left(\alpha +\beta \sqrt{2}\cos \left(t\right)\right)x=0$ ; Meissner, $\ddot{x}+\left(\alpha +\beta sign\left(\cos \left(t\right)\right)\right)x=0$ ; and

Lyapunov, $\ddot{x}+\left(\alpha +\beta \sqrt{\frac{32}{25}}\left(\cos \left(t\right)+\frac{3}{4}\cos \left(2t\right)\right)\right)x=0$ . 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 $\alpha \in \left[\mathrm{0,3}\right]$ and $\beta \in \left[\mathrm{0,1.5}\right]$ .

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 [ [28] , pag. 61].

The best two criteria for stability of the solution for the whole $\alpha -\beta$ plane ( $\alpha \in \left[\mathrm{0,3}\right]$ and $\beta \in \left[\mathrm{0,1.5}\right]$ ) 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 $\beta \ll 1$ , and then seeking the value of $\alpha$ such that the solutions be T or 2T periodic. The general solution and the coefficient $\alpha$ are written in terms of powers of $\beta$ (perturbation expansion)

Table 3. Percentage of stable area obtained by numerical calculation in the plane $\alpha -\beta$ , $\alpha \in \left[\mathrm{0,3}\right]\mathrm{,}\beta \in \left[\mathrm{0,1.5}\right].$

Table 4. Percentage of stability zone obtained by applying each criterion on the Mathieu, Meissner and Lyapunov equation.

$x\left(t,\beta \right)=x_0\left(t\right)+\beta x_1\left(t\right)+{\beta }^2{x}_2\left(t\right)+\cdots$ (80)

$\alpha ={\alpha }_0+\beta {\alpha }_1+{\beta }^2{\alpha }_2+\cdots$ (81)

then, the series (80) and (81) are substituted into the Hill equation and by grouping terms of like powers of $\beta$ , one obtains a set of recursive differential equations. The initial condition of $\alpha$ , i.e. ${\alpha }_0$ , is the value of $\alpha$ where the

Arnold tongues rise, and depends on the excitation function period, ${\alpha }_0={\left(\frac{\text{π}n}{T}\right)}^2$ ,

$n=0,1,2,\cdots$ . 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 [29] .

In [30] the calculation of the first transition curves of the classical form of the Mathieu equation

$\ddot{x}+\left(\alpha +\beta \cos \left(t\right)\right)=0$ (82)

are obtained as

$\alpha =-\frac{1}{2}{\beta }^2+O\left({\beta }^3\right)$ Transition curve associated to the 0^th tongue

$\alpha =\frac{1}{4}\pm \frac{1}{2}\beta -\frac{1}{8}{\beta }^2+O\left({\beta }^3\right)$ Transition curves associated to the 1^st tongue

$\alpha =1+\frac{5}{12}{\beta }^2+O\left({\beta }^3\right)$ Transition curve associated to the 2^th tongue

$\alpha =1-\frac{1}{12}{\beta }^2+O\left({\beta }^3\right)$ Transition curve associated to the 2^th 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 $\alpha \in \left[0,1.5\right]$ and $\beta \in \left[0,0.05\right]$ .

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 $\alpha \in \left[\mathrm{0,1.5}\right]$ and $\beta \in \left[\mathrm{0,0.05}\right]$ .

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 $\alpha \in \left[0,1.5\right],$ $\beta \in \left[0,0.05\right]$ and $\alpha \in \left[0,1.5\right],\beta \in \left[0,0.1\right]$ .

work for small values of $\beta$ , the sufficient stability criteria, here presented, do not have that limitation.

Acknowledgements

The first author, Carlos A. Franco acknowledges the financial support of CONACyT and CINVESTAV.

Appendix

A.

Let $f\left(t\right)$ and $g\left(t\right)$ be two linear independent solutions of

$\ddot{x}+\lambda q\left(t\right)x=0,q\left(t+T\right)=q\left(t\right),q\left(t\right)>0,\forall t$ (83)

subject to the initial condition

$\begin{array}{l}f\left(0\right)=1g\left(0\right)=0\\ f^{\prime }\left(0\right)=0g^{\prime }\left(0\right)=1\end{array}$ (84)

The approximation of the characteristic constant $A\left(\lambda \right)=\frac{1}{2}\left[f\left(T\right)+g^{\prime }\left(T\right)\right]$

associated to the periodic differential Equation (83), consists in write $A\left(\lambda \right)$ and the linear independent solutions $f\left(t\right)$ and $g\left(t\right)$ as alternating series of powers of $\lambda$

$A\left(\lambda \right)=A_0-A_1\lambda +A_2{\lambda }^2+\cdots +{\left(-1\right)}^n{A}_n{\lambda }^n+\cdots$

$f\left(t\right)=f_0\left(t\right)-\lambda f_1\left(t\right)+{\lambda }^2{f}_2\left(t\right)-\cdots$ (85)

$g\left(t\right)=g_0\left(t\right)-\lambda g_1\left(t\right)+{\lambda }^2{g}_2\left(t\right)-\cdots$

where

$A_k=\frac{1}{2}\left(f_k\left(T\right)+{\stackrel{\dot{}}{g}}_k\left(T\right)\right),k=0,1,2,\cdots$

$f_k\left(t\right)={\displaystyle {\int }_0^t}\text{d}{t}_1{\displaystyle {\int }_0^{t_1}}q\left(t_2\right)f_{k-1}\left(t_2\right)\text{d}{t}_2$ (86)

$g_k\left(t\right)={\displaystyle {\int }_0^t}\text{d}{t}_1{\displaystyle {\int }_0^{t_1}}q\left(t_2\right)g_{k-1}\left(t_2\right)\text{d}{t}_2$

with $f_0\left(t\right)=1$ and $g_0\left(t\right)=t$ . And then solve the recurrence system of equations.

$A_0=1,A_1=\frac{T}{2}{\displaystyle {\int }_0^T}q\left(t\right)\text{d}t$

$A_2=\frac{1}{2}{\displaystyle {\int }_0^T}\text{d}{t}_1{\displaystyle {\int }_0^{t_1}}\left(T-t_1+t_2\right)\left(t_1-t_2\right)q\left(t_1\right)q\left(t_2\right)\text{d}{t}_2$

$⋮$ (87)

$A_n=\frac{1}{2}{\displaystyle {\int }_0^T}\text{d}{t}_1{\displaystyle {\int }_0^{t_1}}\text{d}{t}_2\cdots {\displaystyle {\int }_0^{t_{n-1}}}\left(T-t_1+t_n\right)\left(t_1-t_2\right)\cdots \left(t_{n-1}-t_n\right)\cdot q\left(t_1\right)\cdots q\left(t_n\right)\text{d}{t}_n$

In [23] Lyapunov studied in depth the properties of the approximation for $A\left(\lambda \right)$ such as the convergence of the series and the monotonic decrease of the series coefficients. One can prove that for $n>1$ the coefficients $A_n,n>1$ satisfy the inequality

$A_{n-1}{A}_1-n{A}_n>0$ (88)

this can be demonstrated following [5] . By the definitions of the coefficients $A_n$ we can restate (88) as

$\left(f_{n-1}\left(t\right)+{\stackrel{\dot{}}{g}}_{n-1}\left(t\right)\right)t{\displaystyle {\int }_0^t}q\left(t\right)\text{d}t-2n\left(f_n\left(t\right)+\stackrel{\dot{}}{g}\left(t\right)\right)>0$ (89)

for $n>2$ . Defining $S_n\triangleq \left(f_{n-1}\left(t\right)+{\stackrel{\dot{}}{g}}_{n-1}\left(t\right)\right)t{\displaystyle {\int }_0^t}q\left(t\right)\text{d}t-2n\left(f_n\left(t\right)+\stackrel{\dot{}}{g}\left(t\right)\right)$ and noticing that it can be rewritten as

$S_n={\displaystyle {\int }_0^t}\left(F_n+q\left(t\right){\Phi }_n\right)\text{d}t$ (90)

where

$\begin{array}{l}{F}_n=t{\stackrel{\dot{}}{f}}_{n-1}\left(t\right){\displaystyle {\int }_0^t}q\left(t\right)\text{d}t+\left(f_{n-1}\left(t\right)+{\stackrel{\dot{}}{g}}_{n-1}\left(t\right)\right){\displaystyle {\int }_0^t}q\left(t\right)\text{d}t-2n{\stackrel{\dot{}}{f}}_n\left(t\right)\\ {\Phi }_n=t{g}_{n-2}\left(t\right){\displaystyle {\int }_0^t}q\left(t\right)\text{d}t+\left(f_{n-1}\left(t\right)+{\stackrel{\dot{}}{g}}_{n-1}\right)t-2n{g}_{n-1}\left(t\right)\end{array}$ (91)

then (89) is fulfilled if all the coefficients $F_n$ and ${\Phi }_n$ are positive, $F_n>0$ and ${\Phi }_n>0$ . Rewriting $F_n$ and ${\Phi }_n$ we have

$\begin{array}{l}{F}_n={\displaystyle {\int }_0^t}\left(2{\stackrel{\dot{}}{f}}_{n-1}\left(t\right){\displaystyle {\int }_0^t}q\left(t\right)\text{d}t+q\left(t\right)u_n\right)\text{d}t\\ {\Phi }_n={\displaystyle {\int }_0^t}\left(2q\left(t\right)t{\stackrel{\dot{}}{g}}_{n-2}\left(t\right)+{\nu }_n\right)\text{d}t\end{array}$ (92)

where

$\begin{array}{l}{u}_n={\displaystyle {\int }_0^t}\left(2q\left(t\right)\left(g_{n-2}\left(t\right)+t{f}_{n-2}\left(t\right)\right)+F_{n-1}\right)\text{d}t\\ {\nu }_n={\displaystyle {\int }_0^t}\left(2{\stackrel{\dot{}}{f}}_{n-1}\left(t\right)+2{\stackrel{\dot{}}{g}}_{n-2}\left(t\right){\displaystyle {\int }_0^t}q\left(t\right)\text{d}t+q\left(t\right){\Phi }_{n-1}\left(t\right)\right)\text{d}t\end{array}$ (93)

since the excitation function $q\left(t\right)$ is positive for all t, one can notice that $f_n\left(t\right)>0$ , $g_n\left(t\right)>0$ , ${\stackrel{\dot{}}{f}}_n\left(t\right)>0$ and ${\stackrel{\dot{}}{g}}_n\left(t\right)>0$ $\forall t>0$ , these follow from definition. So, $u_n$ and $v_n$ are positive which implies that $F_n$ and ${\Phi }_n$ are positive and the inequality (89) is fulfilled. For details of the proof see [5] .

B. A Brief Introduction to Hamiltonian System Solutions

Consider a second order differential equation in canonical form

$\stackrel{\dot{}}{x}=JH\left(t\right)x$ (94)

where $x=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]$ , $H\left(t\right)$ is a symmetric real periodic matrix $H\left(t+T\right)=H\left(t\right)$ $=H^{\text{T}}\left(t\right)$ and $J$ is the skew symmetric matrix $J=\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]$ .

The state transition matrix of (94) may be expressed as

$\begin{array}{l}\Phi \left(t,0\right)=F\left(t\right)e^{tK}\\ \Phi \left(0,0\right)=I_2\end{array}$ (95)

where $F\left(t+T\right)=F\left(t\right)$ or $F\left(t+T\right)=-F\left(t\right)$ , $\det \left(F\left(t\right)\right)=1$ , $F\left(t\right)$ is a real matrix function, and K is a real matrix with $Tr\left(K\right)=0$ [14] . 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 $ln\left(M\right)$ , in the second a real matrix $ln\left(-M\right)$ and in the third case both $ln\left(M\right)$ and $ln\left(-M\right)$ are real matrices. So K could be defined as

$K=\ln \left(\pm M\right)$ (96)

where the sign $\pm$ is chosen so that K be real.

Remembering that, for periodic systems, $\Phi \left(t+T,0\right)=\Phi \left(t+T,T\right)\Phi \left(T,0\right)=$ $\beta \in \left[\mathrm{0,1.5}\right]$ and substituting (95) into $\Phi \left(t+T\mathrm{,0}\right)$ we have

$F\left(t+T\right)e^{\left(t+T\right)K}=F\left(t\right)e^{tK}F\left(T\right)e^{TK}=F\left(t\right)e^{\left(t+T\right)K}$ (97)

since $F\left(T\right)=M{e}^{-TK}=I_2$ . Then $F\left(t+T\right)=\pm F\left(t\right)$ where the sign is the same as the one in the definition of K. From (95) one can easily see that $\det \left(F\left(t\right)\right)=$ $\alpha \in \left[\mathrm{0,3}\right]\mathrm{,}\beta \in \left[\mathrm{0,1.5}\right].$ . Since $\det \left(F\left(t+T\right)\right)=\det \left(F\left(t\right)\right)$ and $Tr\left(K\right)$ is a real numbers, we can infer that that $Tr\left(K\right)=0$ and $\det \left(F\left(t\right)\right)=1$ .

It follows from matrix similarity that the K matrix may be brought by a similarity transformation to one of following forms:

$\text{(a)}K=S\left[\begin{array}{cc}\lambda & 0\\ 0& -\lambda \end{array}\right]S^{-1}$

$\text{(b)}K=S\left[\begin{array}{cc}0& \epsilon \\ 0& 0\end{array}\right]S^{-1}$

$\text{(c)}K=S\left[\begin{array}{cc}0& -\varphi \\ \varphi & 0\end{array}\right]S^{-1}$

where $\lambda \mathrm{,}\epsilon$ and $\varphi$ are real numbers and S is a real non-singular matrix. Following the nomenclature of [14] , let $\mathcal{H}$ , $\Pi$ and $\mathcal{O}$ be subsets of all admissible matrix K. We will say that the matrix $K\in \mathcal{H}$ if it has the form (a), $K\in \Pi$ if it has the form (b), and $K\in \mathcal{O}$ 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 $e^{KT}$ , and therefore on K i. e. if $K\in \mathcal{H}$ then (94) has one stable solution and one unstable solution (the characteristic multipliers are real and distinct from $\pm 1$ ). If $K\in \Pi$ with $ϵ\ne 0$ then (94) has one periodic solution $x_1\left(t\right)$ and one unstable solution $t{x}_2\left(t\right)$ , if $ϵ=0$ both solutions are stable and periodic (the characteristic multipliers are $\pm 1$ ). And, if $K\in \mathcal{O}$ then both solutions, of (94), are bounded (the characteristic multipliers are complex, lie on the unit circle and are distinct from $\pm 1$ ). All the latter properties are summarized in the Table S1.

Let $\Omega \mathrm{<mi>\; </mi>}$ be the set of real continuous function matrices $F\left(t\right)$ with $F\left(t+T\right)=\pm F\left(t\right)$ and $\det \left(F\left(t\right)\right)=1$ . Now, let $x=F\left(t\right)a$ be a solution of Hamiltonian system (94), where a is a non-zero arbitrary vector and let ${\phi }_x$ denote the rotation of the solution x in time T, since $F\left(T\right)=\pm I_2$ it follows that ${\phi }_x=n\text{π}$ , where n is even if $F\left(T+t\right)=F\left(t\right)$ and n is odd if $F\left(T+t\right)=-F\left(t\right)$ . Let ${\Omega }_n$ denote the set of matrices $F\left(t\right)$ such that the rotation over a period T is $\phi =n\text{π}$ , $\Omega ={\displaystyle {\cup }_{n=-\infty }^{\infty }{\Omega }_n}$ , each of these ${\Omega }_n$ are disjoint sets [14] .

Let $\mathfrak{L}$ be the set of all symmetric matrices $H\left(t\right)$ in (94). As we know each matrix $H\left(t\right)\in \mathfrak{L}$ determines a unique state transition matrix $\Phi \left(t\mathrm{,0}\right)$ . By (95)

Table S1. Relation between the subsets $\mathcal{H}$ , $\Pi$ and $\mathcal{O}$ and the stability of canonical system solutions (24), $ϵ$ , $\lambda$ and $\varphi$ are real positive constants.

the matrix $\Phi \left(t\mathrm{,0}\right)$ determines the pair $F\left(t\right)$ , $K$ where $F\left(t\right)\in \Omega$ and $K\in \mathcal{H}\cup \Pi \cup \mathcal{O}$ . So one can say that $\mathfrak{L}=\Omega \times \left(\mathcal{H}\cup \Pi \cup \mathcal{O}\right)$ . As we have seen $\Omega ={\displaystyle {\cup }_{n=-\infty }^{\infty }{\Omega }_n}$ then

$\mathfrak{L}={\displaystyle \underset{\forall n}{\cup }\left({\mathcal{H}}_n\cup {\Pi }_n\cup {\mathcal{O}}_n\right)}$ (98)

where

${\mathcal{H}}_n={\Omega }_n\times \mathcal{H}$

${\Pi }_n={\Omega }_n\times \Pi$ (99)

${\mathcal{O}}_n={\Omega }_n\times \mathcal{O}$

In words, the set $\mathfrak{L}$ is divided into the subsets ${\mathcal{H}}_n$ , ${\Pi }_n$ and ${\mathcal{O}}_n$ where the matrices K and $F\left(t\right)$ , which are defined by the state transition matrix of the Hamiltonian system (94), belong to one of the sets $\mathcal{H}$ , $\Pi$ or $\mathcal{O}$ and $F\left(t\right)\in {\Omega }_n$ . We say that $H\left(t\right)\in {\mathcal{H}}_n$ if and only if $\Phi \left(t\mathrm{,0}\right)\in {\mathcal{H}}_n$ and so on.

Let $H\left(t\right)\in {\mathcal{H}}_n$ then $F\left(t\right)\in {\Omega }_n$ and $K\in \mathcal{H}$ , let $v_{+}$ and $v_{-}$ be eigenvectors of K such that

$K{v}_{+}=\lambda v_{+},\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}K{v}_{-}=-\lambda v_{-}$

$e^{tK}{v}_{+}=e^{t\lambda }{v}_{+},\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}{e}^{tK}{v}_{-}=e^{-t\lambda }{v}_{-}$

we can notice that the rotation of the two linear independent solutions of (94), $x_1=e^{t\lambda }F\left(t\right)v_{+}$ and $x_2=e^{-t\lambda }F\left(t\right)v_{-}$ , is ${\phi }_{x_1}={\phi }_{x_2}=n\text{π}$ , these follows from the fact that the rotation of a solution doesn’t depend on the eigenvectors $v_{+}$ and $v_{-}$ but on the matrix $F\left(t\right)$ . For any other solution the rotation is either $n\text{π}<\phi <\left(n+1\right)\text{π}$ or $\left(n-1\right)\text{π}<\phi <n\text{π}$ (See [14] § VIII. 1.9, lemma I and II).

Similarly, one can prove that if $H\left(t\right)\in {\mathcal{O}}_n$ then the rotation $\phi$ of all the solutions of (94) will be $n\text{π}<\phi <\left(n+1\right)\text{π}$ . The sketch of the proof is as follows, by assumption $F\left(t\right)\in {\Omega }_n$ and $K\in \mathcal{O}$ , from (95) we know that any solution of (94) can be written as $x=F\left(t\right)y$ where $y\left(t\right)=e^{tK}y\left(0\right)$ , from form (c) it follows

$y\left(t\right)=S\left[\begin{array}{cc}\cos \left(\varphi t\right)& -\sin \left(\varphi t\right)\\ \sin \left(\varphi t\right)& \cos \left(\varphi t\right)\end{array}\right]S^{-1}y\left(0\right),0<\varphi <\frac{\text{π}}{T}$ (100)

with $b=S^{-1}y\left(0\right)$ and $\det \left(S\right)=1$ . Thus $S^{-1}y\left(t\right)$ moves uniformly in a circle, describing and angle $\varphi T$ in time $T$ . Since $0<\varphi T<\text{π}$ then $0<{\phi }_y<\text{π}$ and finally $n\text{π}<\phi <\left(n+1\right)\text{π}$ .

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Richards, J.A. (2012) Analysis of Periodically Time-Varying Systems. Springer Science & Business Media.
[2]	Berchio, E. and Gazzola, F. (2015) A Qualitative Explanation of the Origin of Torsional Instability in Suspension Bridges. Nonlinear Analysis: Theory, Methods & Applications, 121, 54-72.
[3]	Nusinsky, I. and Hardy, A.A. (2006) Band-Gap Analysis of One-Dimensional Photonic Crystals and Conditions for Gap Closing. Physical Review B, 73, Article ID: 125104. https://doi.org/10.1103/PhysRevB.73.125104
[4]	Starzinskii, V. (1955) Survey of Works on Conditions of Stability of the Trivial Solution of a System of Linear Differential Equations with Periodic Coefficients. American Mathematical Society.
[5]	Lyapunov, A.M. (1992) The General Problem of the Stability of Motion. International Journal of Control, 55, 531-534. https://doi.org/10.1080/00207179208934253
[6]	Yakubovich, V. (1951) Stability Criteria for a System of Two Canonical Equations with Periodic Coefficients. Doklady Akademii nauk SSSR, 28, 221-224.
[7]	Yakubovich, V. (1950) On the Boundedness of Solutions of the Equation: y + p(t)y = 0, p (t + w) = p(t). Doklady Akademii Nauk SSSR, No. 5, 901-903.
[8]	Tang, X.-H. and Zhang, M. (2012) Lyapunov Inequalities and Stability for Linear Hamiltonian Systems. Journal of Differential Equations, 252, 358-381.
[9]	Hochstadt, H. (1962) A Stability Criterion for Hill’s Equation. Proceedings of the American Mathematical Society, 13, 601-603.
[10]	Hochstadt, H. (1962b) A Class of Stability Criteria for Hills Equation. Quarterly of Applied Mathematics, 20, 92-93. https://doi.org/10.1090/qam/132883
[11]	Xu, R. (2007) Eigenvalue Estimates for Hills Equation with Periodic Coefficients1. Applied Mathematical Sciences, 1, 2601-2608.
[12]	Shi, J. (1999) A New Form of Discriminant for Hill Equation. Annals of Differential Equations, 15, 191-210.
[13]	Hochstadt, H. (1961) Asymptotic Estimates for the Sturm-Liouville Spectrum. Communications on Pure and Applied Mathematics, 14, 749-764. https://doi.org/10.1002/cpa.3160140408
[14]	Yakubovich, V. and Starzhinskii, V. (1975) Linear Differential Equations with Periodic Coefficients 2 vols. Wiley, New York.
[15]	Arnol’d, V.I. (1983) Remarks on the Perturbation Theory for Problems of Mathieu Type. Russian Mathematical Surveys, 38, 215-233. https://doi.org/10.1070/RM1983v038n04ABEH004210
[16]	Gelfand, I. and Lidskii, V. (1958) On the Structure of the Regions of Stability of Linear Canonical Systems of Differential Equations with Periodic Coefficients. American Mathematical Society Translations, 8, 143-181.
[17]	Brockett, R. (1970) Finite Dimensional Linear Systems. Wiley.
[18]	Meyer, K., Hall, G. and Offin, D. (2008) Introduction to Hamiltonian Dynamical Systems and the N-Body Problem. Vol. 90, Springer Science & Business Media.
[19]	Lancaster, P. and Tismenetsky, M. (1985) The Theory of Matrices: With Applications. Elsevier.
[20]	Magnus, W. and Winkler, S. (2013) Hill’s Equation. Dover Publications.
[21]	Ungar, P. (1961) Stable Hill Equations. Communications on Pure and Applied Mathematics, 14, 707-710. https://doi.org/10.1002/cpa.3160140403
[22]	Hochstadt, H. (1963) Functiontheoretic Properties of the Discriminant of Hill’s Equation. Mathematische Zeitschrift, 82, 237-242. https://doi.org/10.1007/BF01111426
[23]	Lyapunov, A.M. (1902) On the Series, Encountered in the Theory of Linear Second Order Differential Equations with Periodic Coefficients. Academy of Science Notes on the Physics-Mathematical Department, VIII, 1-70. (In Russian)
[24]	Demidovich, B.P. (1967) Lectures on Stability Theory. Nauka, Moscow. (In Russian)
[25]	Adrianova, L.Y. (1995) Introduction to Linear Systems of Differential Equations. American Mathematical Soc.
[26]	Sevryuk, M. (1986) Reversible Systems. Volume 1211 of Lecture Notes in Mathematics. Springer, Berlin. https://doi.org/10.1007/BFb0075877
[27]	Loud, W. (1975) Stability Regions for Hill’s Equation. Journal of Differential Equations, 19, 226-241.
[28]	Cesari, L. (2012) Asymptotic Behavior and Stability Problems in Ordinary Differential Equations. Vol. 16, Springer Science & Business Media.
[29]	Nayfeh, A.H. and Mook, D.T. (2008) Nonlinear Oscillations. John Wiley & Sons.
[30]	Kevorkian, J. and Cole, J.D. (2013) Perturbation Methods in Applied Mathematics. Vol. 34, Springer Science & Business Media.