Coupled Mathieu Equations: γ-Hamiltonian and μ-Symplectic

Abstract

Several theoretical studies deal with the stability transition curves of coupled and damped Mathieu equations utilizing numerical and asymptotic methods. In this contribution, we exploit the fact that symplectic maps describe the dynamics of Hamiltonian systems. Starting with a Hamiltonian system, a particular dissipation is introduced, which allows the extension of Hamiltonian and symplectic matrices to more general γ -Hamiltonian and μ -symplectic matrices. A proof is given that the state transition matrix of any γ -Hamiltonian system is μ -symplectic. Combined with Floquet theory, the symmetry of the Floquet multipliers with respect to a μ -circle, which is different from the unit circle, is highlighted. An attempt is made for generalizing the particular dissipation to a more general form. The methodology is applied for calculation of the stability transition curves of an example system of two coupled and damped Mathieu equations.

1. Introduction

Dynamical systems represented by nonlinear or linear ordinary differential equations with periodic coefficients occur in many engineer problems (see for instance [1, 2]). The simplest example of such a system is the Mathieu equation. Most investigations in literature deal with the corresponding stability transition curves [3]. Some works analyze the stability of two coupled Mathieu equations [4, 5, 6]. In general, an asymptotic or a numerical analysis method is required for analyzing this class of systems. Perturbation techniques may lead to cumbersome expression, at least for second-order perturbation [7], and a numerical analysis may require considerable computation time. In this contribution, an extension of the theory developed in [8] is exposed in which coupled Mathieu equations are analyzed in the context of a Hamiltonian system.

The literature on Hamiltonian systems is vast. We focus on the two main references [9, 10] that are relevant for the present work. The latter focuses on linear periodic Hamiltonian systems. Although every periodic mechanical system possesses at least a small amount of dissipation, the main literature on linear Hamiltonian systems does not incorporate a dissipation. The dynamics of Hamiltonian systems can be described by symplectic maps [11]. A key fact here is that a symplectic transformation preserves the Hamiltonian structure of the underlying dynamic system. In this work we attempt to derive an appropriate formalism for linear Hamiltonian systems incorporating a very particular dissipation. For this purpose we redefine and develop the properties of the so-called $\gamma$ -Hamiltonian and $\mu$ -symplectic matrices. With the last definitions, we prove that the state transition matrix of any $\gamma$ -Hamiltonian system is $\mu$ -symplectic. The relevance of the symplectic matrices or symplectic maps lies on their symmetry which allows simplifying many computations and analysis [12]. The formalism is benchmarked for two coupled and damped Mathieu equations highlighting its advantages. Due to the symmetry of the symplectic matrices, the parametric resonance zones are characterized, which allows faster computations, and with higher accuracy, of the stability transition curves. This work is an extension of the contribution presented in [8, 13].

2. Preliminaries on matrices

2.1 Symplectic matrices

Definition 1 The matrix $A\in {\mathbb{R}}^{2n\times 2n}$ is called symplectic if it satisfies

$$A^T\mathit{JA}=J,$$ E1

with

$$J=\left[\begin{array}{cc}0& I_n\\ -I_n& 0\end{array}\right]$$ E2

and $I_n$ is the $n\times n$ identity matrix.

Note that for $J$ the following relations hold: $J^T=-J$ , $J^{-1}=J^T$ , $J^2=-I_{2n}$ , and $det\left(J\right)=1$ . The determinant of a symplectic matrix is 1 ([9]), and $I_{2n}$ and $J$ are symplectic matrices themselves. If $A$ and $B$ are of the same dimensions and symplectic, then $\mathit{AB}$ is also symplectic because ${\left(\mathit{AB}\right)}^{T}J\left(\mathit{AB}\right)=B^T{A}^T\mathit{JAB}=B^T\mathit{JB}=J$ . Finally and importantly, the inverse of a symplectic matrix always exists and is also symplectic:

$$A^{-1}=J^{-1}{A}^{T}J:{\left(J^{-1}{A}^{T}J\right)}^{T}J\left(J^{-1}{A}^{T}J\right)=J^T{\mathit{AJA}}^{T}J=J.$$ E3

The set of the symplectic matrices of dimension $2n\times 2n$ forms a group. The corresponding characteristic polynomial of a symplectic matrix $A\in {\mathbb{R}}^{2n\times 2n}$

$$P_A\left(\lambda \right)=det\left(\lambda I_{2n}-A\right)={\lambda }^{2n}+a_{2n-1}{\lambda }^{2n-1}+\dots +a_1\lambda +1$$

is a reciprocal polynomial:

$$P_A\left(\lambda \right)={\lambda }^{2n}{P}_A\left(\frac{1}{\lambda }\right)$$ E4

This is equivalent to stating that the coefficients of $P_A\left(\lambda \right)$ satisfy the relation $a_k=a_{2n-k}$ or rewriting as a matrix product

$$\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{2n-1}\\ a_{2n}\end{array}\right]=\left[\begin{array}{ccccccc}0& 0& \cdots & 0& \cdots & 0& 1\\ 0& 0& \cdots & 0& \cdots & 1& 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 1& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮& ⋮\\ 0& 1& \cdots & 0& \cdots & 0& 0\\ 1& 0& \cdots & 0& \cdots & 0& 0\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{2n-1}\\ a_{2n}\end{array}\right]$$ E5

Since $A$ is real, if $\lambda$ is an eigenvalue of $A$ , then so are ${\lambda }^{-1}$ , $\overline{\lambda }$ , and ${\overline{\lambda }}^{-1}$ , where the bar indicates the complex conjugate. Equivalently, the eigenvalues of a symplectic matrix are reciprocal pairs. This property is called reflexivity [11]. Consequently, the eigenvalues are symmetric with respect to the unit circle, namely, if there is an eigenvalue inside of the unit circle, then there must be a corresponding eigenvalue outside of the unit circle. As a result of the coefficient symmetry of a symplectic matrix $A$ , the following transformation is proposed in [12]:

$$\delta =\lambda +\frac{1}{\lambda },$$ E6

where $\lambda \in \sigma \left(A\right)$ . This transforms the characteristic polynomial $P_A\left(\lambda \right)$ of degree $2n$ to an auxiliary polynomial $Q_A\left(\delta \right)$ of degree $n$ , while keeping all pertinent information of the original polynomial [12].

2.2 Hamiltonian matrices

Definition 2 The matrix $A\in {\mathbb{R}}^{2n\times 2n}$ ( $A\in {\mathbb{C}}^{2n\times 2n}$ ) is said to be Hamiltonian if and only if

$$A^{T}J+\mathit{JA}=\mathbf{0}.$$ E7

Let $P_A\left(s\right)$ be the characteristic polynomial of $A$ , then $P_A\left(s\right)$ is an even polynomial, and it only has even powers. Thus, the eigenvalues of $A$ are symmetric with respect to the imaginary axis, i.e., if $s$ is an eigenvalue of $A$ , then $-s$ is an eigenvalue, too. Furthermore, if the matrix $A$ is real, $\overline{s}$ and $-\overline{s}$ are eigenvalues as well. Then the eigenvalues of the Hamiltonian matrix are located symmetrically with respect to both real and imaginary axis. The eigenvalues appear in real pairs, purely imaginary pairs, or complex quadruples [9, 14].

2.3 $\mathit{\mu }$ -symplectic matrices

The next definitions and properties attempt to generalize the classical definitions above.

Definition 3 $M\in {\mathbb{R}}^{2n\times 2n}$ is called $\mu$ -symplectic matrix if

$$M^T\mathit{JM}=\mathit{\mu J}$$ E8

is satisfied for $\mu \in \left(0,1\right]$ .

Lemma 4 The determinant of a $\mu$ -symplectic matrix $M\in {\mathbb{R}}^{2n\times 2n}$ is ${\mu }^n$ .

To see the proof of the last lemma, see Appendix A. If $M$ is a $\mu$ -symplectic matrix, $M^2$ is a ${\mu }^2$ -symplectic matrix, and the set of $\mu$ -symplectic matrix matrices does not form a group.

Lemma 5 The characteristic polynomial of a $\mu$ -symplectic matrix $M\in {\mathbb{R}}^{2n\times 2n}$ satisfies

$$P_M\left(\frac{\mu }{\lambda }\right)=\frac{{\mu }^n}{{\lambda }^{2n}}P\left(\lambda \right).$$ E9

Proof 6 $P_M\left(\lambda \right)=det\left(\lambda I_{2n}-M^T\right)=det\left(\lambda I_{2n}-\mu {\mathit{JM}}^{-1}{J}^{-1}\right)$

$$\begin{array}{l}=det\left(J\right)det\left(\lambda I_{2n}-\mu M^{-1}\right)det\left(J^{-1}\right)=det\left(\frac{\lambda }{\mu }M-I_{2n}\right)det\left(\mu M^{-1}\right)\\ ={\mu }^{n}det\left(\left(-\frac{\lambda }{\mu }\right)\left(\frac{\mu }{\lambda }{I}_{2n}-M\right)\right)={\mu }^n{\left(-\frac{\lambda }{\mu }\right)}^{2n}det\left(\frac{\lambda }{\mu }{I}_{2n}-M\right)=\frac{{\lambda }^{2n}}{{\mu }^n}{P}_M\left(\frac{\mu }{\lambda }\right)\end{array}$$

Corollary 7 The eigenvalues of a $\mu$ -symplectic matrix $M$ satisfy the symmetry

$$\lambda \in \sigma \left(M\right)\Rightarrow \left(\frac{\mu }{\lambda }\right)\in \sigma \left(M\right).$$ E10

The product of each pair of eigenvalues contributes with $\mu$ to $det\left(M\right)$ , and there are $n$ of these pairs; therefore, $det\left(M\right)={\mu }^n$ . If all eigenvalues have the same magnitude, i.e., ${\lambda }_i=rexp\left({\theta }_i\right)$ , then ${\prod }_{i=1}^{2n}\left|\lambda \right|={\prod }_{i=1}^{2n}\left|{\mathit{re}}^{{\theta }_i}\right|=r^{2n}=det\left(M\right)={\mu }^n$ . From this we find that $r=\sqrt{\mu }$ , independent of $n$ . This may be interpreted as a “symmetry” with respect to a circle of radius $r=\sqrt{\mu }$ . Since $M$ is real if $\lambda$ is an eigenvalue of $M$ , then $\overline{\lambda }$ , $\frac{\mu }{\lambda }$ , and $\frac{\mu }{\overline{\lambda }}$ are also eigenvalues of $M$ . Moreover, the eigenvalues are symmetric with respect to the $\mu$ -circle: if there is an eigenvalue inside of the $\mu$ -circle, then there must be another eigenvalue outside (see Figure 1a for a visualization).

Figure 1.

Symmetries in the spectra in the complex plane: (a) $\mu$ -symplectic matrix and (b) $\gamma$ -Hamiltonian matrix.

Remark 8 Due to Eq. (9) , the characteristic polynomial $P_M\left(\lambda \right)=m_{2n}{\lambda }^{2n}+\dots m_1\lambda +m_0$ of the $\mu$ -symplectic matrix $M$ satisfies the following relations:

$$\begin{array}{c}{m}_0=m_{2n}{\mu }^n\\ m_1=m_{2n-1}{\mu }^{n-1}\\ ⋮\\ m_n=m_n\\ ⋮\\ m_{2n-1}=m_1{\mu }^{1-n}\\ m_{2n}=1=m_0{\mu }^{-n}\end{array}$$ E11

rewritten as a product of matrices yields

$$\left[\begin{array}{c}{m}_0\\ m_1\\ ⋮\\ m_{2n-1}\\ m_{2n}\end{array}\right]=\left[\begin{array}{ccccccc}0& 0& \cdots & 0& \cdots & 0& {\mu }^{-n}\\ 0& 0& \cdots & 0& \cdots & {\mu }^{-n+1}& 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 1& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮& ⋮\\ 0& {\mu }^{n-1}& \cdots & 0& \cdots & 0& 0\\ {\mu }^n& 0& \cdots & 0& \cdots & 0& 0\end{array}\right]\left[\begin{array}{c}{m}_0\\ m_1\\ ⋮\\ m_{2n-1}\\ m_{2n}\end{array}\right]$$ E12

For $\mu =1$ , the relations in Eq. (12) reduce to Eq. (5) .

Remark 9 By applying the transformation

$$\delta =\lambda +\frac{\mu }{\lambda }$$ E13

the characteristic polynomial $P_M\left(\lambda \right)$ of degree $2n$ , associated to a $\mu$ -symplectic matrix, is reduced to an auxiliary polynomial $Q_M\left(\delta \right)$ of degree $n$ . For instance,

$$\begin{array}{l}n=2\left\{\begin{array}{l}{P}_M\left(\lambda \right)={\lambda }^4+m_3{\lambda }^3+m_2{\lambda }^2+m_3\mathit{\mu \lambda }+{\mu }^2\\ Q_M\left(\delta \right)={\delta }^2+m_3\delta +m_2-2\mu \end{array}\right)\\ n=3\left\{\begin{array}{l}{P}_M\left(\lambda \right)={\lambda }^6+m_5{\lambda }^5+m_4{\lambda }^4+m_3{\lambda }^3+\mu m_4{\lambda }^2+{\mu }^2{m}_5\lambda +{\mu }^3\\ Q_M\left(\delta \right)={\delta }^3+m_5^2\delta +\left(m_4-3\mu \right)\delta +m_3-2m_5\mu \end{array}\right)\\ n=4\left\{\begin{array}{l}{P}_M\left(\lambda \right)={\lambda }^8+m_7{\lambda }^7+m_6{\lambda }^6+m_5{\lambda }^5+m_4{\lambda }^4+\mu m_5{\lambda }^3\\ +{\mu }^2{m}_6{\lambda }^2+{\mu }^3{m}_7\lambda +{\mu }^4\\ Q_M\left(\delta \right)={\delta }^4+m_7^3\delta +\left(m_6-4\mu \right){\delta }^2+\left(m_5-3m_7\mu \right)\delta +m_4\\ -2m_6\mu +2{\mu }^2\end{array}\right)\\ n=5\left\{\begin{array}{l}{P}_M\left(\lambda \right)={\lambda }^{10}+m_9{\lambda }^9+m_8{\lambda }^8+m_7{\lambda }^7+m_6{\lambda }^6+m_5{\lambda }^5\\ +\mu m_6{\lambda }^4+{\mu }^2{m}_7{\lambda }^3+{\mu }^3{m}_8{\lambda }^2+{\mu }^4{m}_9\lambda +{\mu }^5\\ Q_M\left(\delta \right)={\delta }^5+m_9{\delta }^4+\left(m_8-5\mu \right){\delta }^3+\left(m_7-4m_9\mu \right){\delta }^2\\ +\left(m_6-3m_8\mu +5{\mu }^2\right)\delta +m_5-2\mu m_7+2m_9{\mu }^2\end{array}\right)\end{array}$$

Note that the property of the characteristic polynomial of a $\mu$ -symplectic matrix in Eq. (9) reduces to Eq. (4) at $\mu =1$ . Then Eq. (12) represents the “symmetry” of the characteristic polynomial for all $\mu \in \left(0,1\right]$ . Although the definition of $\mu$ -symplectic matrices appears in [9], no further properties were developed within this reference. In the next section, we reveal their relationship as a generalized definition of Hamiltonian matrices, the so-called $\gamma$ -Hamiltonian matrices.

2.4 $\mathit{\gamma }$ -Hamiltonian matrices

Definition 10 A matrix $A\in {\mathbb{R}}^{2n\times 2n}$ ( $A\in {\mathbb{C}}^{2n\times 2n}$ ) is called $\gamma$ -Hamiltonian matrix if for some $\gamma \ge 0$ ,

$$A^{T}J+\mathit{JA}=-2\mathit{\gamma J}$$ E14

Lemma 11 $A$ is $\gamma$ -Hamiltonian if and only if $A+\gamma I_{2n}$ is Hamiltonian.

Proof 12 If $A$ is $\gamma$ -Hamiltonian, then $A^{T}J+\mathit{JA}=-2\mathit{\gamma J}$ which can be rewritten as ${\left[A+\gamma I_{2n}\right]}^{T}J+J\left[A+\gamma I_{2n}\right]=0$ . Hence, $\left[A+\gamma I_{2n}\right]$ is Hamiltonian.

Lemma 13 If $A$ is $\gamma$ -Hamiltonian and if $s+\gamma \in \sigma \left(A\right)$ , then $-s+\gamma \in \sigma \left(A\right)$ .

Proof 14 Recall that if $\sigma \left(R\right)=\left\{r_1,\dots r_{2n}\right\}$ , then $\sigma \left(R+\gamma I_{2n}\right)=\left\{r_1+\gamma ,\dots ,r_{2n}+\gamma \right\}$ . Then if $s+\gamma \in \sigma \left(A\right)$ , then $s\in \sigma \left(A+\gamma I_{2n}\right)$ , since $\left[A+\gamma I_{2n}\right]$ is Hamiltonian and consequently $-s\in \sigma \left(A+\gamma I_{2n}\right)$ which is equivalent to $-s+\gamma \in \sigma \left(A\right)$ .

Remark 15 If in the last lemma all the eigenvalues of the Hamiltonian matrix $A+\gamma I_{2n}$ have zero real parts, then the real parts of the eigenvalues of the $\gamma$ -Hamiltonian matrix $A$ are identical to $-\gamma$ . Thus, the eigenvalues of the $\gamma$ -Hamiltonian matrix $A$ are symmetric with respect to the vertical line $-\gamma$ in the complex plane (see Figure 1b for a visualization).

Figure 2.

All configurations of multiplier positions with respect to the unit (in solid line) and $\mu$ -circle (dashed line).

Notice that real Hamiltonian matrices have their spectrum symmetric with respect to the real and imaginary axes, whereas the spectrum of real $\gamma$ -Hamiltonian matrices is symmetric with respect to the real axis and a vertical line at $\mathbf{Re}\left(s\right)=-\gamma$ . Then the eigenvalues of a real $\gamma$ -Hamiltonian matrix are placed: (i) in quadruples symmetrically with respect the real axis and the line $\mathbf{Re}\left(s\right)=-\gamma$ , (ii) pairs on the line $\mathbf{Re}\left(s\right)=-\gamma$ and symmetric with the real axis, and (iii) real pairs symmetric with the line $\mathbf{Re}\left(s\right)=-\gamma$ . All cases are shown in Figure 1b .

By the last lemma, the characteristic polynomial of the $\gamma$ -Hamiltonian $A$ satisfies

$$P_A\left(s+\gamma \right)=P_A\left(\gamma -s\right)$$

with

$$\begin{array}{l}{P}_A\left(\gamma -s\right)={\left(\gamma -s\right)}^{2n}+a_{2n-1}{\left(\gamma -s\right)}^{2n-1}+\dots +a_1\left(\gamma -s\right)+a_0\\ P_A\left(s+\gamma \right)={\left(s+\gamma \right)}^{2n}+a_{2n-1}{\left(s+\gamma \right)}^{2n-1}+\dots +a_1\left(s+\gamma \right)+a_0\end{array}$$

Thus, $P_A\left(s\right)$ depends only on $n$ coefficients. For instance, for $n=1$ ,

${\left(s+\gamma \right)}^2+a_1\left(s+\gamma \right)+a_0={\left(\gamma -s\right)}^2+a_1\left(\gamma -s\right)+a_0$ . Equating the coefficients leads to $a_1=-2\gamma$ , $a_0=a_0$ , and finally to

$$P_A\left(s\right)=s^2-2\mathit{\gamma s}+a_0.$$

Similarly, the polynomials for the lowest values of $n$ read

$$n=2:$$

$$P_A\left(s\right)=s^4-4\gamma s^3+a_2{s}^2+\left(8{\gamma }^3-2\gamma a_2\right)+a_0$$

$$n=3:$$

$$P_A\left(s\right)=s^6-6\gamma s^5+a_4{s}^4+\left(40{\gamma }^3-4\gamma a_4\right)s^3+a_2{s}^2+\left(-96{\gamma }^5+8{\gamma }^3{a}_4-2\gamma a_2\right)+a_0$$

$$n=4:$$

$$P_A\left(s\right)=s^8-8\gamma s^7+a_6{s}^6+\left(112{\gamma }^3-6\gamma a_6\right)s^5+a_4{s}^4+\left(-896{\gamma }^5+40{\gamma }^3{a}_6-4\gamma a_4\right)s^3$$

$$+a_2{s}^2+\left(2176{\gamma }^7-96{\gamma }^5{a}_6+8{\gamma }^3{a}_4-2\gamma a_2\right)+a_0.$$

Furthermore, by applying the transformation

$$\varphi =s-\gamma ,$$ E15

the polynomial $P_A\left(s\right)$ can be written as an auxiliary polynomial $Q_A\left(\varphi \right)$ which only has even coefficients, namely,

$$\begin{array}{ll}{P}_A\left(s\right)& =s^{2n}+a_{2n-1}{s}^{2n-1}+a_{2n-2}{s}^{2n-2}+\dots +a_2{s}^2+a_{1}s+a_0\\ Q_A\left(\varphi \right)& ={\varphi }^{2n}+q_{2n-2}{\varphi }^{2n-2}+\dots +q_2{\varphi }^2+q_0\end{array}$$

For instance,

$$n=1:$$

$$Q_A\left(\varphi \right)={\varphi }^2+a_0-{\gamma }^2$$

$$n=2:$$

$$Q_A\left(\varphi \right)={\varphi }^4+\left(a_2-6{\gamma }^2\right){\varphi }^2+5{\gamma }^4-a_2{\gamma }^2+a_0$$

$$n=3:$$

$$Q_A\left(\varphi \right)={\varphi }^6+\left(a_4-15{\gamma }^2\right){\varphi }^4+\left(a_2-6a_4{\gamma }^2+75{\gamma }^4\right){\varphi }^2-61{\gamma }^6+5a_4{\gamma }^4-a_2{\gamma }^2+a_0$$

$$n=4:$$

$$Q_A\left(\varphi \right)={\varphi }^8+\left(a_6-28{\gamma }^2\right){\varphi }^6+\left(a_4-15a_6{\gamma }^2+350{\gamma }^4\right){\varphi }^4$$

$$+\left(a_2-6a_4{\gamma }^2+75a_6{\gamma }^4-1708{\gamma }^6\right){\varphi }^2+1385{\gamma }^8-61a_6{\gamma }^6+5a_4{\gamma }^4-a_2{\gamma }^2+a_0.$$

3. Linear $\mathit{\gamma }$ -Hamiltonian systems

Definition 16 If there is a differentiable function called Hamiltonian function (energy) $\mathcal{H}\left(t,x,y\right)$ , $\mathcal{H}:\mathbb{R}\times {\mathbb{R}}^n\times {\mathbb{R}}^n\mapsto \mathbb{R}$ , which satisfies

$$\dot{x}={\left(\frac{\partial \mathcal{H}}{\partial y}\right)}^T\mathit{\text{and}}\dot{y}=-\left(\frac{\partial \mathcal{H}}{\partial y}\right),$$

then it is called a Hamiltonian system. If $\mathcal{H}\left(t,x,y\right)$ is a quadratic function with respect to $x$ and $y$ , then the system is a linear Hamiltonian system.

It is easy to prove that if $\mathcal{H}$ does not depend on $t$ , $\mathcal{H}\left(x,y\right)$ is a first integral. However, this is no longer true in the time-periodic case. In the time-periodic case, even for $n=1$ , the integration of the equations is not possible. Any linear Hamiltonian system can be written as

$$\dot{z}=\mathit{JH}\left(t\right)z$$ E16

where $H^T\left(t\right)=H\left(t\right)$ is a symmetric matrix (Hermitian in the complex case). Herein, the variables used in the definition satisfy $z={\left[x^T,y^T\right]}^T$ . Therefore, the dimension of real Hamiltonian systems is always even. Finally, note that the product $\mathit{JH}$ satisfies the condition for a Hamiltonian matrix. The fundamental property of any linear Hamiltonian system is that the state transition matrix of the system in Eq. (16) is a symplectic matrix (see [9] for more details).

If $A$ is $\gamma$ -Hamiltonian matrix, or equivalently, $A+\gamma I_{2n}$ is a Hamiltonian matrix for some $\gamma >0$ ; then it follows from Eq. (16) that

$$\dot{x}=\left[A+I_{2n}\right]x=\mathit{JHx}$$

for some matrix $H=H^T$ . From the last equation $\left[A+I_{2n}\right]=\mathit{JH}$ , we obtain

$$A=J\left[H+\gamma I_{2n}\right].$$ E17

Any $\gamma$ -Hamiltonian matrix $A$ may be written as in Eq. (17), which motivates the next definition.

Definition 17 Any linear system that can be written as

$$\dot{x}=A\left(t\right)x=J\left[H\left(t\right)+\mathit{\gamma J}\right]x$$ E18

with $x\in {\mathbb{R}}^{2n}$ , $H^T\left(t\right)=H\left(t\right)$ , and $\gamma \ge 0$ is called a linear $\mathit{\gamma }$ -Hamiltonian system.

Lemma 18 The state transition matrix of a linear $\gamma$ -Hamiltonian system in Eq. (18) is $\mu$ -symplectic with $\mu =e^{-2\mathit{\gamma t}}$ .

Proof 19 Let be $N\left(t\right)=\mathrm{\Phi }\left(t,0\right)$ be the state transition of Eq. (17) , and then

$$\dot{N}\left(t\right)=A\left(t\right)N\left(t\right).$$

Differentiating the product $N^T\mathit{JN}$ gives

$$\begin{array}{c}\frac{d}{\mathit{dt}}{N}^T\mathit{JN}={\dot{N}}^T\mathit{JN}+N^{T}J\dot{N}={\left(\mathit{\text{AN}}\right)}^T\mathit{JN}+N^{T}J\left(\mathit{\text{AN}}\right)\\ =N^T\left(A^{T}J+\mathit{JA}\right)N=N^T\left({\left(J\left(H+\mathit{\gamma J}\right)\right)}^{T}J+J\left(J\left(H+\mathit{\gamma J}\right)\right)\right)N\\ =-2\gamma N^T\mathit{JN}\end{array}$$ E19

Since $N^T\left(0\right)\mathit{JN}\left(0\right)=J$ , we get ¹

$$N^T\left(t\right)\mathit{JN}\left(t\right)=e^{-2\mathit{\gamma t}}J=\mathit{\mu J}.$$

Therefore, $N$ is $\mu$ -symplectic.

Lemma 20 Consider the transformation

$$x=S\left(t\right)z$$ E20

with $S\left(t\right)$ a symplectic matrix for all $t$ . Then the transformation in Eq. (20) preserves the $\gamma$ -Hamiltonian form of the system, Eq.(18).

Proof 21 From the definition $S^T\mathit{JS}=0\to {\dot{S}}^T\mathit{JS}+S^{T}J\dot{S}=0$ , thus ${\dot{S}}^T\mathit{JS}=-S^{T}J\dot{S}$ , and from Eq.(20)

$$\dot{x}=S\dot{z}+\dot{S}z\to S^{-1}\dot{x}=\dot{z}+S^{-1}\dot{S}z$$

then applying the transformation Eq.(20) into Eq.(18) it is obtained as $\dot{z}+S^{-1}\dot{S}z=S^{-1}J\left(H+\mathit{\gamma J}\right)\mathit{Sz}$ ; then from the symplectic definition matrix $S^{-1}=J^{-1}{S}^{T}J$ ,

$$\begin{array}{cc}\dot{z}& =S^{-1}J\left(H+\mathit{\gamma J}\right)\mathit{Sz}-S^{-1}\dot{S}z=\left(J^{-1}{S}^{T}J\right)\mathit{\text{JHSz}}-\mathit{\gamma Iz}-\left(J^{-1}{S}^{T}J\right)\dot{S}z\\ & ={\mathit{JS}}^T\mathit{HSz}-\mathit{\gamma Iz}+{\mathit{JS}}^{T}J\dot{S}z=J\left(S^T\mathit{HS}+S^{T}J\dot{S}+\mathit{\gamma J}\right)z=J\left(\tilde{H}+\mathit{\gamma J}\right)z\end{array}$$

where $\tilde{H}=S^T\mathit{HS}+S^{T}J\dot{S}$ , but ${\left(S^{T}J\dot{S}\right)}^T={\dot{S}}^T{J}^{T}S=-{\dot{S}}^T\mathit{JS}=-\left(-S^{T}J\dot{S}\right)=S^{T}J\dot{S}$ ; therefore $\tilde{H}={\tilde{H}}^T■$ .

3.1 Mechanical, linear $\mathit{\gamma }$ -Hamiltonian system

Consider any mechanical system described by the equation

$$\tilde{M}\ddot{y}+\tilde{D}\dot{y}+\tilde{K}\left(t\right)y=0$$ E21

where $y\left(t\right)\in {\mathbb{R}}^n$ , $\tilde{K}\left(t\right)={\tilde{K}}^T\left(t\right)\in {\mathbb{R}}^{n\times n}$ , and the constant matrices $\tilde{M}$ and $\tilde{D}\in {\mathbb{R}}^{n\times n}$ such that $\tilde{M}={\tilde{M}}^T>0$ and $\tilde{D}={\tilde{D}}^T$ . Then there always exists a linear transformation $T$ such that

$$\begin{array}{c}{T}^T\tilde{M}T=I_n\\ T^T\tilde{D}T=D=diag\left\{d_1,d_2,\dots ,d_n\right\}\\ \sigma \left({\tilde{M}}^{-1}\tilde{D}\right)=\left\{d_1,d_2,\dots ,d_n\right\}\end{array}$$

(e.g., see [15]). Therefore, applying the transformation $y=\mathit{Tz}$ yields

$$\ddot{z}+D\dot{z}+K\left(t\right)z=0,$$ E22

where $K\left(t\right)=T^T\tilde{K}\left(t\right)T$ . Eq. (22) can be rewritten as a first-order system by introducing the state vector $x={\left[z^T{\dot{z}}^T\right]}^T$ :

$$\dot{x}=\left[\begin{array}{cc}{0}_{n\times n}& I_n\\ -K\left(t\right)& -D\end{array}\right]x$$ E23

where $x\in {\mathbb{R}}^{2n\times 2n}$ . Let

$$Q=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_n& I_n\\ -I_n& I_n\end{array}\right]$$ E24

be an orthogonal matrix satisfying ${\mathit{QQ}}^T=Q^{T}Q=I_{2n}$ , and also $\mathit{JQ}=\mathit{QJ}$ , one can introduce the transformation $w=Q^{T}x$ , and Eq. (23) gives

$$\dot{w}=Q^T\left[\begin{array}{cc}{0}_{n\times n}& I_n\\ -K\left(t\right)& -D\end{array}\right]\mathit{Qw}=\frac{1}{2}\left[\begin{array}{cc}K\left(t\right)-I_n-D& K\left(t\right)+I_n+D\\ -K\left(t\right)+D-I_n& -K\left(t\right)+I_n-D\end{array}\right]w,$$

or equivalently,

$$\dot{w}=J\left(\frac{1}{2}\left[\begin{array}{ll}K\left(t\right)+I_n-D& K\left(t\right)-I_n\\ K\left(t\right)-I_n& K\left(t\right)+I_n+D\end{array}\right]+\frac{1}{2}\left[\begin{array}{cc}D& 0_{n\times n}\\ 0_{n\times n}& D\end{array}\right]J\right)w.$$ E25

Since $D=diag\left\{d_1,d_2,\dots ,d_n\right\}$ and $K=K^T$ , the matrix

$$H\left(t\right)=\frac{1}{2}\left[\begin{array}{cc}K\left(t\right)+I_n-D& K\left(t\right)-I_n\\ K\left(t\right)-I_n& K\left(t\right)+I_n+D\end{array}\right]$$

is also symmetric $H\left(t\right)=H{\left(t\right)}^T$ . Therefore, Eq. (25) can be cast into the $\gamma$ -Hamiltonian linear system form $\dot{w}=J\left(H+\mathit{\gamma J}\right)w$ if $\gamma$ is approximated as $\gamma \approx \frac{1}{2n}{\sum }_{i=1}^n{d}_i$ . In the special case $d=d_1=d_2=\cdots =d_n$ , $\gamma$ is given exactly given by $\gamma =\frac{d}{2}$ .

3.2 Periodic linear systems

This section summarizes the main results on periodic linear systems. The proofs and details are omitted and can be found in [16, 17]. Consider the linear periodic system:

$$\dot{x}=B\left(t\right)x\text{with}B\left(t\right)=B\left(t+\mathrm{\Omega }\right)$$ E26

where $x\in {\mathbb{R}}^n$ , $B\in {\mathbb{R}}^{n\times n}$ , and $\mathrm{\Omega }$ are the fundamental periods.

Theorem 22 (Floquet) The state transition matrix $\mathrm{\Phi }\left(t,t_0\right)$ of the system in Eq. (26) may be factorized as

$$\mathrm{\Phi }\left(t,t_0\right)=P^{-1}\left(t\right)e^{R\left(t-t_0\right)}P\left(t_0\right)$$ E27

where

$$P^{-1}\left(t\right)=\mathrm{\Phi }\left(t,0\right)e^{-\mathit{Rt}}.$$ E28

In addition $P^{-1}\left(t\right)=P^{-1}\left(t+\mathrm{\Omega }\right)$ is a periodic matrix of the same period $\mathrm{\Omega }$ , and $R$ is in general a complex constant matrix [ 18 ].

Definition 23 We define the monodromy matrix $M$ associated to the Eq. (26) as

$$M=\mathrm{\Phi }\left(\mathrm{\Omega },0\right).$$ E29

The monodromy matrix may be defined as $M_{t_0}=\mathrm{\Phi }\left(\mathrm{\Omega },\mathrm{\Omega }+t_0\right)$ , but we use only the spectrum of the monodromy matrix, $\sigma \left(M\right)$ . From.

$\mathrm{\Phi }\left(t,t_0\right)={\left.P^{-1}\left(t\right)e^{R\left(t-t_0\right)}P\left(t_0\right)\right|}_{t=t_0+\mathrm{\Omega }}=\mathrm{\Phi }\left(\mathrm{\Omega },\mathrm{\Omega }+t_0\right)=P^{-1}\left(t_0+\mathrm{\Omega }\right)e^{R\mathtt{W}}P\left(t_0\right)=P^{-1}\left(t_0\right)e^{R\mathtt{W}}P\left(t_0\right)$ , because $P$ and $P^{-1}$ are $\mathrm{\Omega }$ -periodic. This last relation shows that $M$ and $M_{t_0}$ are similar matrices and possess the same spectrum. Moreover, if $t_0=0$ in the Floquet theorem, then $\mathrm{\Phi }\left(t,0\right)=Q\left(t\right)e^{\mathit{Rt}}$ based on $Q\left(t\right)=Q\left(t+\mathrm{\Omega }\right)$ and $Q\left(0\right)=I_n$ ; we have

$$M=\mathrm{\Phi }\left(\mathrm{\Omega },0\right)=Q\left(\mathrm{\Omega }\right)e^{R\mathrm{\Omega }}=Q\left(0\right)e^{R\mathrm{\Omega }}=e^{R\mathrm{\Omega }}.$$ E30

Definition 24 The eigenvalues ${\lambda }_i$ of the monodromy matrix are called characteristic multipliers or multipliers. The numbers ${\rho }_i$ , not unique, defined as ${\lambda }_i=e^{{\rho }_i\mathrm{\Omega }}$ , are called characteristic exponents or Floquet exponents.

Corollary 25 (Lyapunov-Floquet Transformation) If we define the change of coordinates

$$z\left(t\right)=P\left(t\right)x\left(t\right)$$ E31

where $P$ fulfills Eq. (28), then the periodic linear system in Eq. (26) can be transformed into a linear time-invariant system

$$\dot{z}\left(t\right)=\mathit{Rz}\left(t\right)$$ E32

where $R$ is a constant matrix as introduced in the Floquet theorem.

The transformation in Eq. (31) is a Lyapunov transformation which means that the stability properties of the linear system in Eq. (26) are preserved. Therefore any periodic system as in Eq. (26) is reducible to a system in Eq. (32) with constant coefficients² ([16]). However, the matrix $R$ is not always real (e.g., see [10, 20]). In the present discussion, we only use its spectrum $\sigma \left(R\right)$ .

For analyzing $x\left(t\right)$ as $t\to \infty$ , we assume that the initial conditions are given at $t_0=0$ . Then for any $t>0$ , $t$ may be expressed as $t=k\mathrm{\Omega }+\tau$ , where $k\in {\mathbb{Z}}_{+}$ and $\tau \in [0,\mathrm{\Omega })$ . Applying the well-known properties of the state transition matrix, the solution can be written as

$$\begin{array}{c}x\left(t\right)=\mathrm{\Phi }\left(t,0\right)x_0=\mathrm{\Phi }\left(k\mathrm{\Omega }+\tau ,0\right)x_0\\ =\mathrm{\Phi }\left(k\mathrm{\Omega }+\tau ,k\mathrm{\Omega }\right)\underset{k-\text{terms}}{\underbrace{\mathrm{\Phi }\left(k\mathrm{\Omega },\left(k-1\right)\mathrm{\Omega }\right)\mathrm{\Phi }\left(\left(k-1\right)\mathrm{\Omega },\left(k-2\right)\mathrm{\Omega }\right)\dots \mathrm{\Phi }\left(\mathrm{\Omega },0\right)}}{x}_0\\ =\mathrm{\Phi }\left(\tau ,0\right){\left[\mathrm{\Phi }\left(\mathrm{\Omega },0\right)\right]}^k{x}_0=\mathrm{\Phi }\left(\tau ,0\right)M^k{x}_0\end{array}$$

Analyzing the last expression, the terms $\mathrm{\Phi }\left(\tau ,0\right)$ and $x_0$ are bounded; the following three cases can be distinguished:

$$x\left(t\right)\to 0\iff \underset{k\to \infty }{lim}{M}^k=0\iff \sigma \left(M\right)\subset D=\left\{z\in \mathbb{C}:\left|z\right|<1\right\}.$$

1. the solution $x\left(t\right)$ is bounded $\iff \underset{k\to \infty }{lim}{M}^k=0$ is bounded $\iff \sigma \left(M\right)\subset \overline{D}=\left\{z\in \mathbb{C}:\left|z\right|\le 1\right\}$ , and if $\lambda \in \sigma \left(M\right)$ and $\left|\lambda \right|=1$ , $\lambda$ is a simple root of the minimal polynomial of $M$ .

2. $‖x\left(t\right)‖\to \infty \iff \exists \lambda \in \sigma \left(M\right):\left|\lambda \right|>1$ or $\exists \lambda \in \sigma \left(M\right):\left|\lambda \right|=1$ and $\lambda$ is a multiple root of the minimum polynomial of $M$ .

Theorem 26 (Lyapunov-Floquet) Considering the linear periodic system in Eq. (26) , then the system is (a) asymptotic stable if and only if Eq. (1) is satisfied, (b) stable if and only if Eq. (2) is satisfied, and (c) Unstable if and only if Eq. (3) is satisfied.

Due to the Lyapunov-Floquet transformation in Eq. (31), the stability of the periodic linear system in Eq. (26) can be determined by analyzing the system in Eq. (32).

Corollary 27 The system in Eq. (26) is:

1. Asymptotically stable $\iff \sigma \left(R\right)\subset Z=\left\{z\in \mathbb{C}:\mathbf{Re}\left(z\right)<0\right\}$ .

2. Stable $\iff \sigma \left(R\right)\subset \overline{Z}=\left\{z\in \mathbb{C}:Re\left(z\right)\le 0\right\}$ , if $\mathbf{Re}\left(z_i\right)=0$ are simple roots of the minimum polynomial of $R$ .

3. Unstable $\iff \exists$ ${\rho }_i\in \sigma \left(R\right):Re\left(z\right)>0$ or $\sigma \left(R\right)\subset \overline{Z}$ $\&$ $\exists$ $\mathbf{Re}\left(z_i\right)=0$ which is a multiple root of minimum polynomial of $R$ .

4. Periodic $\mathit{\gamma }$ -Hamiltonian systems

Once the linear Hamiltonian systems become periodic, i.e., the matrix $H\left(t\right)$ of the system in Eq. (18) possesses a periodically time-varying $H\left(t\right)=H\left(t+\mathrm{\Omega }\right)$ , the underlying monodromy matrix becomes $\mu$ -symplectic and $\gamma$ -Hamiltonian.

Definition 28 Any linear periodic system that can be written as

$$\dot{x}=J\left[H\left(t\right)+\mathit{\gamma J}\right]x$$ E33

with $H\left(t\right)=H\left(t+\mathrm{\Omega }\right)$ will be named linear periodic $\gamma$ -Hamiltonian system, where $x\in {\mathbb{R}}^{2n}$ and $H^T\left(t\right)=H\left(t\right)$ are a $2n\times 2n$ matrix and $\gamma \ge 0$ .

Remark 29 According to Lemma 18, the state transition matrix $\mathrm{\Phi }\left(t,t_0\right)$ of Eq. (33) is $\mu$ -symplectic, in particular, the state transition matrix evaluated over one period $\mathrm{\Omega }$ .

Corollary 30 The monodromy matrix $M=e^{R\mathrm{\Omega }}$ and the matrix $R$ of the periodic system in Eq. (33) are $\mu$ -symplectic and $\gamma$ -Hamiltonian matrices, respectively, with $\mu =e^{-2\gamma \mathrm{\Omega }}$ .

Proof 31 From the definition of a $\mu$ -symplectic matrix $M^T\mathit{JM}={\left(e^{R\mathrm{\Omega }}\right)}^{T}J\left(e^{R\mathrm{\Omega }}\right)=\mathit{\mu J}$ , we obtain $e^{R^T\mathrm{\Omega }}=\mu {\mathit{Je}}^{-R\mathrm{\Omega }}{J}^{-1}=\mathit{\mu J}\left\{I_{2n}-R\mathrm{\Omega }+\frac{\mathit{RR}{\mathrm{\Omega }}^2}{2}-\frac{\mathit{RRR}{\mathrm{\Omega }}^3}{3!}+\dots +\frac{R^k{\mathrm{\Omega }}^4}{k!}+\dots \right\}{J}^{-1}=\mu e^{-{\mathit{JRJ}}^{-1}\mathrm{\Omega }}=e^{-2\gamma \mathrm{\Omega }}{e}^{-{\mathit{JRJ}}^{-1}\mathrm{\Omega }}$ thus $e^{R^T\mathrm{\Omega }}=e^{-2\gamma \mathrm{\Omega }-{\mathit{JRJ}}^{-1}\mathrm{\Omega }}\Rightarrow R^{T}J+\mathit{JR}=-2\mathit{\gamma J}$ .

This corollary states the main relation in our analysis. The symmetry of the $\mu$ -symplectic matrix will be utilized for obtaining the stability conditions of the system in Eq. (33). Furthermore, by applying the Lyapunov transformation

$$z\left(t\right)=P\left(t\right)x\left(t\right)$$ E34

we conclude that any linear periodic $\gamma$ -Hamiltonian system can be reduced to a linear time-invariant $\gamma$ -Hamiltonian system

$$\dot{z}\left(t\right)=\mathit{Rz}\left(t\right).$$ E35

The next two subsections are based on [12] and are adapted for characteristic polynomials of $\mu$ -symplectic matrices.

4.1 Stability of a system with one degree of freedom

For $n=1$ , the characteristic polynomial of the monodromy matrix $M$ associated with the system in Eq. (33) becomes $P_M\left(\lambda \right)={\lambda }^2+\mathit{a\lambda }+\mu$ with $a=-\mathrm{tr}\left(M\right)$ . According to the Lemma 18, $M$ is $\mu$ -symplectic. Then, there are two multipliers symmetric to the circle of radius $r$ and the real axis. Therefore, the multipliers only can leave the unit circle at the coordinates $\left(1,0\right)$ or $\left(-1,0\right)$ (see Figure 2 ). Note that the term $-a$ is equal to the transformation in Eq. (13):

$$\delta =\lambda +\frac{\mu }{\lambda }=\mathrm{tr}\left(M\right)=-a$$

Theorem 32 For $n=1$ , the system in Eq. (33) is asymptotically stable if and only if the inequality

$$\left|a\right|<\left(1+\mu \right)$$

is satisfied.

Proof 33 Since the multipliers only leave the unit circle on the points $\lambda =1$ or $\lambda =-1$ , the stability boundaries are given by

$$\begin{array}{c}{P}_M\left(1\right)={\left(1\right)}^2+a\left(1\right)+\mu =a+\left(\mu +1\right)\\ P_M\left(-1\right)={\left(-1\right)}^2+a\left(-1\right)+\mu =-a+\left(\mu +1\right)\end{array}.$$

This means that $a+\left(\mu +1\right)>0$ and $-a+\left(\mu +1\right)>0$ must be fulfilled; thus, $\left|a\right|<\left(1+\mu \right)$ .

4.2 Stability of a system with two degrees of freedom

For $n=2$ , the characteristic polynomial of the monodromy matrix $M$ associated with the system in Eq. (33) reads

$$P_M\left(\lambda \right)={\lambda }^4+a{\lambda }^3+b{\lambda }^2+\mathit{a\mu \lambda }+{\mu }^2$$ E36

where $a=-\mathrm{tr}\left(M\right)$ and $2b={\left(\mathrm{tr}\left(M\right)\right)}^2-\mathrm{tr}\left(M^2\right)$ . There are four multipliers, and due to the symmetry with respect to the $\mu$ -circle, they can be categorized in the position configurations depicted in Figure 2 .

Respecting that the characteristic polynomial is associated with a $\mu$ -symplectic matrix, we can use the transformation

$$\delta =\lambda +\frac{\mu }{\lambda }$$ E37

to obtain the auxiliary polynomial

$$Q_M\left(\delta \right)={\delta }^2+\mathit{a\delta }+b-2\mu .$$ E38

The symmetry of the eigenvalues yield

$$a=-\mathrm{tr}\left(M\right)={\lambda }_1+\frac{\mu }{{\lambda }_1}+{\lambda }_3+\frac{\mu }{{\lambda }_3}={\delta }_1+{\delta }_2.$$

The transition boundaries are characterized by having at least one eigenvalue at $\left|\lambda \right|=1$ . The simplest cases are if $\lambda =1$ ( $\delta =1+\mu$ ) or $\lambda =-1$ ( $\delta =-1-\mu$ ). These points overlay if a real-valued multiplier leaves the unit circle at the point $\left(1,0\right)$ or $\left(0,-1\right)$ (see the cases c, d, e, f, or g in Figure 2 ). Substituting these two values into Eq. (36) gives

$$b=-a\left(1+\mu \right)-\left(1+{\mu }^2\right)$$ E39

and

$$b=a\left(1+\mu \right)-\left(1+{\mu }^2\right).$$ E40

Considering the case $\lambda \in \mathbb{C}$ , we search the transition boundary line when two complex multipliers leave the unit circle at points different to $\left(1,0\right)$ and $\left(0,-1\right)$ (see cases a or b in Figure 2 ). Then the transition boundary line can be obtained by considering the symmetry of the multipliers with respect to the real axis and the circle of the radius $r=\sqrt{\mu }$ . Here, ${\lambda }_1=x+\mathit{iy}$ , ${\lambda }_2=\frac{\mu }{{\lambda }_1}$ , ${\lambda }_3=x-\mathit{iy}$ , and ${\lambda }_4=\frac{\mu }{{\lambda }_3}$ . At $\left|\lambda \right|=1$ $\left(\sqrt{x^2+y^2}=1\right)$ , it follows that

$${\lambda }_1=x+\mathit{iy},{\lambda }_2=\mu \left(x-\mathit{iy}\right),{\lambda }_3=y-\mathit{iy},{\lambda }_4=\mu \left(x+\mathit{iy}\right).$$

Hence, the transformation in Eq. (13) follows:

$$\begin{array}{l}{\delta }_1={\lambda }_1+\frac{\mu }{{\lambda }_1}=x\left(1+\mu \right)+\mathit{iy}\left(1-\mu \right)\\ {\delta }_2={\lambda }_3+\frac{\mu }{{\lambda }_3}=x\left(1+\mu \right)-\mathit{iy}\left(1-\mu \right)\end{array}$$

Adding ${\delta }_1$ and ${\delta }_2$ gives

$${\delta }_1+{\delta }_2=2x\left(1+\mu \right).$$ E41

From Eq. (38) we obtain

$${\delta }_{1,2}=\frac{-a}{2}\pm \frac{\sqrt{a^2+8\mu -4b}}{2}.$$ E42

Note that for ${\delta }_1$ and ${\delta }_2$ to become complex, the inequality

$$4b>a^2+8\mu$$

must be fulfilled. Adding ${\delta }_1$ and ${\delta }_2$ , one obtains

$${\delta }_1+{\delta }_2=-a$$ E43

Equating Eqs. (41) and (43) yields

$$2x\left(1+\mu \right)=-a.$$ E44

The real part $x$ of the eigenvalues results from Eq. (37)

$${\lambda }^2-\mathit{\lambda \delta }+\mu =0$$

and

$$\lambda =\frac{\delta \pm \sqrt{{\delta }^2-4\mu }}{2}.$$ E45

Substituting Eq. (42) into Eq. (45) and choosing only the positive signs gives

$${\lambda }_1=\frac{1}{4}\left(-a+\sqrt{w-4\mu -2b+a^2}\right)+\frac{i}{4}\left(\sqrt{w+4\mu +2b-a^2}+\sqrt{-\left(a^2+8\mu \right)+4b}\right)$$

with the abbreviation $w=2\sqrt{-4a^2\mu +{\left(b+2\mu \right)}^2}$ . Consequently, the real part of $\lambda$ is

$$x=\frac{1}{4}\left(-a+\sqrt{w-4\mu -2b+a^2}\right),$$

and substituting into Eq. (44) results in

$$\frac{1}{2}\left(-a+\sqrt{w-4\mu -2b+a^2}\right)\left(1+\mu \right)=-a$$

which can be solved for $b$ to obtain the transition boundary curve

$$b=\frac{{\mu }^4+2{\mu }^3+2{\mu }^2+2\mu +a^2\mu +1}{{\left(1+\mu \right)}^2}.$$ E46

Two intersection points exist on each line in Eqs. (39) and (40) with the curve defined by Eq. (46). These points are

$$b=\frac{1}{\mu }\left({\mu }^4+{\mu }^3+2{\mu }^2+\mu +1\right)$$ E47

$$b={\mu }^2+4\mu +1$$ E48

and are highlighted in Figure 3 . The line in Eq. (47), dashed line in the figure, is a necessary condition for stability.

Figure 3.

Multiplier map in the case of $n=2$ : Horizontal and vertical axes are the coefficients a and $b$ of the characteristic polynomial of the monodromy matrix M in Eq. (36). The solid lines represent the borders of the inequalities in Theorem 34, Eqs. (49), (50), and (51). The dots indicate the position of multipliers and the unit circle, in solid line, associated with the system in Eq. (33) in the case of $n=2$ . The dashed circle depicts the $\mu$ -circle.

Theorem 34 The Eq. (33) when $n=2$ is asymptotically stable if and only if the inequalities are fulfilled:

$$b\ge -a\left(1+\mu \right)-\left(1+{\mu }^2\right),$$ E49

$$b\ge a\left(1+\mu \right)-\left(1+{\mu }^2\right),$$ E50

$$b\le \frac{{\mu }^4+2{\mu }^3+2{\mu }^2+a^2\mu +2\mu +1}{{\left(1+\mu \right)}^2}.$$ E51

From this analysis, the multipliers position in relation to the unit circle and $\mu$ -circle are defined by inequalities. These split the complex plane into four regions as it is shown in the Figure 3 .

5. Coupled Mathieu equations

Figure 4.

Multiplier map and stability charts, for example, systems in Eq. (52). Multiplier map corresponds to Figure 3 but now with colored regions for the different unstable multiplier configurations. Stability charts are given for different values of damping. (a) Multiplier map: $a$ and $b$ are the coefficients of the corresponding characteristic polynomial of the monodromy matrix. All colored zones correspond to unstable positions multipliers configurations. (b) Stability chart of coupled Mathieu equations, Eq. (52), with small damping: ${\mathit{\Theta }}_{12}={\mathit{\Theta }}_{21}=\mathit{0}$ and ${\mathit{\Theta }}_{11}={\mathit{\Theta }}_{22}=\mathit{0.1}.$ Each color code is according to the position of the multipliers as in Figure 4a . (c) Stability chart of coupled Mathieu equations, Eq. (52), with damping: ${\mathit{\Theta }}_{12}={\mathit{\Theta }}_{21}=\mathit{0}$ and ${\mathit{\Theta }}_{11}={\mathit{\Theta }}_{22}=\mathit{0.3}$ .

Consider two coupled and damped Mathieu equations of the following form:

$$\left[\begin{array}{c}{\ddot{z}}_1\\ {\ddot{z}}_2\end{array}\right]+\left[\begin{array}{cc}{\mathrm{\Theta }}_{11}& {\mathrm{\Theta }}_{12}\\ {\mathrm{\Theta }}_{21}& {\mathrm{\Theta }}_{22}\end{array}\right]\left[\begin{array}{c}{\dot{z}}_1\\ {\dot{z}}_2\end{array}\right]+\left(\left[\begin{array}{cc}{\omega }_1^2& 0\\ 0& {\omega }_2^2\end{array}\right]+\beta \left[\begin{array}{cc}{Q}_{11}& Q_{12}\\ Q_{21}& Q_{22}\end{array}\right]\left[\begin{array}{c}{z}_1\\ z_2\end{array}\right]cos\left(\mathit{\nu t}\right)\right)\left[\begin{array}{c}{z}_1\\ z_2\end{array}\right]=0.$$ E52

Following the procedure presented in Section 3.1, the system in Eq. (52) can be cast into the $\gamma$ -Hamiltonian form in Eq. (33) if ${\mathrm{\Theta }}_{12}={\mathrm{\Theta }}_{21}$ and $Q_{12}=Q_{21}$ , i.e., the coefficient matrices $\mathrm{\Theta }$ and $Q$ are symmetric. In this case, the coupled Mathieu equations present all the properties of the periodic $\gamma$ -Hamiltonian system defined in Eq. (33) for $n=2$ and $\mathrm{\Omega }=2\pi /\nu$ . Hence, all the above analysis on Hamiltonian systems can be applied. The monodromy matrix is computed by numerical methods, and the stability chart is obtained by applying the Theorem 34.

The following numerical values are chosen for the analysis of a specific system ${\omega }_1^2=8$ , ${\omega }_2^2=26$ , $Q_{11}=Q_{22}=2$ , $Q_{12}=Q_{21}=-2$ . Figure 4a depicts the multiplier chart similar to Figure 3 . The unstable regions are colored and the stable regions are kept white. Each color depicts a specific configuration of the multiplier positions within the unit circle and the $\mu$ -circle according to the inequalities stated in Theorem 34 and visualized in Figures 3 and 4a . The description of each color is relevant because each color describes the parametric resonance phenomenon. Thus, yellow, magenta, and cyan colors refer to the configuration of four real-valued multipliers, two of them inside and two outside of the unit circle. These multipliers are either all negative (magenta region), all positive (yellow region), or two positive and two negative (cyan region). The blue and red regions indicate two complex conjugate multipliers on the $\mu$ -circle, while the other two are real with $\left|\lambda \right|>1$ . The two real multipliers are either positive (blue) or negative (red). Then, all four multipliers are complex conjugate within the green region. In this case, two multipliers lie inside and two outside of the unit circle.

Additionally, parametric primary resonances occur at parametric excitation frequencies $\nu =2{\omega }_i/k$ , with $k\in {\mathrm{\mathbb{N}}}^{+}$ , and parametric combination resonances of summation type occur at $\nu =\left({\omega }_1+{\omega }_2\right)/k$ [7, 10]. These frequencies are also observed for the example system in Figure 4 . The green regions mark parametric combination resonances. The blue and red regions correspond to parametric primary resonances. The presented calculation technique can be categorized as a semi-analytical method. After rewriting the original system into the form in Eq. (33), the monodromy matrix is constructed by integrating the equations of motion using numerical methods.

Subsequently, the coefficients of the characteristic polynomial of the monodromy matrix can be computed as $a=-\mathrm{tr}\left(M\right)$ and $2b={\left(\mathrm{tr}\left(M\right)\right)}^2-\mathrm{tr}\left(M^2\right)$ . This technique avoids the computation of the eigenvalues itself. This has the main advantage that numerical problems on the computation of the eigenvalues are avoided, e.g., numerical sensitivity of multipliers [21].

The definitions of $\mu$ -symplectic and $\gamma$ -Hamiltonian matrices allow the analysis of a linear periodic Hamiltonian system with a particular dissipation. The main result of the proposed theory lies in Corollary 30 which states that the state transition matrix of any $\gamma$ -Hamiltonian system is $\mu$ -symplectic. The symmetry properties of the eigenvalues of $\mu$ -symplectic matrices lead to an efficient calculation of the stability boundaries of this type of system. The general framework is applied for the example analysis of two damped and coupled Mathieu equations confirming the faster and robust computation of the stability chart. The procedure can be extended to a higher number of coupled Mathieu equations as outlined above.

Acknowledgments

M.Ramírez is thankful to Conacyt for the provided support during his Ph.D. studies in México and to the Austrian EFRE/LEADER project LaZu-CLLD IWB TIROL OSTT 005 for his research stay at Campus Technik Lienz.

Proof 35 In [22], Rim proposed an elementary proof that real symplectic matrices have determinant one; following the same procedure, we prove that the symplectic matrices have determinant ${\mu }^n$ . From the definition $det\left(M^T\mathit{JM}\right)=det\left(M^T\right)det\left(J\right)det\left(M\right)=det\left(\mathit{\mu J}\right)={\mu }^{2n}\to det\left(M\right)=\pm {\mu }^n$ , therefore it is necessary to prove that $det\left[M\right]=-{\mu }^{2n}$ is false. Considering the matrix $S=M^{T}M+\mu I_{2n}$ since $M^{T}M\ge 0$ and $\mu \in \left(0,1\right]$ , the matrix $S$ has real and greater than $\mu$ eigenvalues:

$$det\left(S\right)=det\left(M^{T}M+\mu I_{2n}\right)>\mu .$$ E53

Now from the definition $M^{-T}={\mu }^{-1}{\mathit{JMJ}}^{-1}$ and rewriting $S$ ,

$$S=M^{T}M+\mu I_{2n}=M^T\left(M+\mu M^{-T}\right)=M^T\left(M+{\mathit{JMJ}}^{-1}\right)$$

denotes the subblocks of $M$ as follows:

$$M=\left[\begin{array}{cc}{M}_{11}& M_{12}\\ M_{21}& M_{21}\end{array}\right]$$

with $M_{11},M_{12},M_{21},M_{21}\in {\mathbb{R}}^{n\times n}$ ; thus $M+{\mathit{JMJ}}^{-1}=\left[\begin{array}{cc}{M}_{11}+M_{22}& M_{12}-M_{21}\\ M_{21}-M_{12}& M_{11}+M_{22}\end{array}\right]$ if $A=M_{11}+M_{22}$ and $B=M_{12}-M_{21}$ ; then

$$M+{\mathit{JMJ}}^{-1}=\left[\begin{array}{cc}A& B\\ -B& A\end{array}\right].$$ E54

We rewrite (54) with the unitary transformation $T$ :

$$T=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_n& I_n\\ {\mathit{iI}}_n& -{\mathit{iI}}_n\end{array}\right]\iff T^{-1}=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_n& -{\mathit{iI}}_n\\ I_n& {\mathit{iI}}_n\end{array}\right]$$

$$M+{\mathit{JMJ}}^{-1}=\left[\begin{array}{cc}A& B\\ -B& A\end{array}\right]=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_n& I_n\\ {\mathit{iI}}_n& -{\mathit{iI}}_n\end{array}\right]\left[\begin{array}{cc}A+\mathit{iB}& 0\\ 0& A-\mathit{iB}\end{array}\right]\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_n& -{\mathit{iI}}_n\\ I_n& {\mathit{iI}}_n\end{array}\right]$$

therefore

$$\begin{array}{cc}det\left(S\right)& =det\left(M^T\left(M+{\mathit{JMJ}}^{-1}\right)\right)=det\left(M^T\right)det\left(M+{\mathit{JMJ}}^{-1}\right)\\ & =det\left(M\right)det\left(A+\mathit{iB}\right)det\left(A-\mathit{iB}\right)=det\left(M\right)det\left(A+\mathit{iB}\right)det\left(\overline{A+\mathit{iB}}\right)\end{array}$$

since $A,B$ are real, the complex conjugation commute with the determinant, and then

$$det\left(S\right)=det\left(M\right)det\left(A+\mathit{iB}\right)\overline{det\left(A-\mathit{iB}\right)}=det\left(M\right){\left|det\left(A+\mathit{iB}\right)\right|}^2$$

and from Eq. (53)

$$0<\mu <det\left(S\right)=det\left(M\right){\left|det\left(A+\mathit{iB}\right)\right|}^2$$

then ${\left|det\left(A+\mathit{iB}\right)\right|}^2>0$ and $det\left(M\right)>0$ ; therefore $det\left(M\right)={\mu }^n$ . $■$