ABSTRACT

In this work we study two types of Discrete Hill’s equation. The first comes from the discretization process of a Continuous-time Hill’s equation, we called Discretized Hill’s equation. The Second is a naturally obtained in Discrete-Time and will be called Discrete-time Hill’s equation. The objective of discretization is preserving the continuous-time behavior and we show this property. On the contrary a completely different dynamic property was found for the Discrete-Time Hill’s equation. At the end of the paper is shown that both types share the nonoscillatory behavior of solutions in the 0-th Arnold Tongue.

Keywords:

Arnold Tongues, Discrete Hill’s Equation, Monodromy Matrix, Discretized Hamiltonian

1. Introduction

Hill equation $\ddot{z}+\left(\alpha +\beta p\left(t\right)\right)z\left(t\right)=0$ where $p\left(t+T\right)=p\left(t\right)$ was introduced in 1877 by George W. Hill on his work on the part of the motion of the lunar perigee [1] , this model equation arises in many areas of applied mathematics where the stability of periodic motions is an issue. The Floquet theory gives us a factorization of the solution of linear systems with periodic coefficients [2] [3] even more gives us tools to determine the stability of the system in terms of the Monodromy matrix, i.e. the stability of the system is given in terms of the state transition matrix evaluated after one period T.

In 1912 Hamel found that as far as $\beta$ is sufficiently small there is an ${\alpha }_0$ such that for $\alpha \le {\alpha }_0$ all the solutions of the Hill equation are nonoscillatory and oscillatory if $\alpha >{\alpha }_0$ [4] [5] .

In late 1940s until 1960s, two Russian scientists Krein and Yakubovich, established the foundation of linear Hamiltonian with periodic coefficients [6] and [7] . Other important contributions were made by Gelfand-Lidskii [8] , Starzhinskii [9] , Bolotin [10] , Atkinson [11] , Eastham [12] and first of all Lyapunov [13] , more recent works on Hill’s equation and Arnold Tongues can be found on [14] [15] [16] , however this results are mostly for continuous time case, results on the discrete time case are frequently found in solid-state physics topics [17] [18] [19] .

Through the Arnold tongues, we compare both types of equations: Discretized and Discrete Hill’s equations. We will plot for the first time the Arnold tongues for the Discrete Hill’s equation and give a formula of where the Arnold tongues touches the $\beta =0$ axis. We will also proof the discrete version of the nonoscillatory criteria for discrete Hill’s equation.

This paper is organized as follows: the first section is an introduction and historical overview, in Section 2 we give the preliminary necessary for this work as well as the characterization of the continuous time Hill’s equation, Section 3 is focused on the Discretization of linear Hamiltonian systems and the methods that preserve such structure, in Section 4 we discretize the Continuous time Hill’s equation and analyze its stability conditions, in Section 5 we present the discrete Hill’s equation and develop its properties, finally in Section 6 we present some conclusions.

2. Preliminaries

In this section we will give a brief introduction to Floquet theory [20] . We define the Monodromy matrix and give stability conditions for the Hill’s equation in terms of the Monodromy matrix. We will mention the Hamiltonian systems and prove that the Hill’s equation is a Hamiltonian system. We will define some discretization methods and give the form of the solutions of a linear difference equation with constant coefficients.

Consider the following system

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

known as the Hill’s equation, where $\alpha$ , $\beta$ are real parameters, $p\left(t\right)=p\left(t+T\right)$ is a periodic function with period T. The parameter $\alpha$ represents the square of the natural frequency for $\beta =0$ ; the parameter $\beta$ is the amplitude of the parametric excitation, and the periodic function $p\left(t\right)$ is called the excitation function.

If we define the two-dimensional vector $y=\left[\begin{array}{cc}z& \dot{z}\end{array}\right]$ , the Equation (1) can be rewritten as

$\dot{y}\left(t\right)=A\left(t\right)y\left(t\right)$ (2)

where $A\left(t\right)=\left[\begin{array}{cc}0& 1\\ -\left(\alpha +\beta p\left(t\right)\right)& 0\end{array}\right]$ , $A\left(t+T\right)=A\left(t\right)$ .

2.1. Floquet Theory

The solution of any linear system $\overline{y}\left(t\right)=\overline{A}\left(t\right)\overline{y}\left(t\right)$ can be given in terms of its state transition matrix as $\overline{y}\left(t\right)=\Phi \left(t,t_0\right)\overline{y}\left(t_0\right)$ . For linear systems with periodic coefficients as (2), Floquet asserts that the state transition matrix¹ may be factorized as:

Theorem 1. Floquet [20] Consider a linear periodic system $\dot{y}\left(t\right)=A\left(t\right)y\left(t\right)$ as in (2), if $\Phi \left(t\mathrm{,}{t}_0\right)$ is the state transition matrix solution of (2) then $\Phi \left(t+T\mathrm{,}{t}_0\right)$ is also a solution of (2). Even more, exists an invertible matrix $Q\left(t\right)=Q\left(t+T\right)$ of the same period as the system (2) with initial value $Q\left(0\right)=I$ such that its state transition matrix satisfies:

$\Phi \left(t,t_0\right)=Q^{-1}\left(t\right)e^{B\left(t-t_0\right)}Q\left(t_0\right)$ (3)

where $B$ is a constant matrix, not necessarily real².

If we make $t_0=0$ in (3), we get the most well-known version:

Corollary 1. Floquet theorem Consider a linear periodic system $\dot{y}\left(t\right)=A\left(t\right)y\left(t\right)$ as in (2), exists an invertible matrix $Q\left(t\right)=Q\left(t+T\right)$ of the same period as the system (2) with initial value $Q\left(0\right)=I$ such that its state transition matrix satisfies:

$\Phi \left(t,0\right)=Q^{-1}\left(t\right)e^{Bt}$ (4)

$B$ is a constant matrix, again not necessarily real.

Now if we evaluate (4) at $t=T$ , we get the following definition.

Definition 1. The matrix

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

is a matrix particularly important and is known as the Monodromy Matrix.

Remark 1. The Monodromy matrix defined by (5) is dependent of the initial time $t_0$ ; but not its spectrum. Let us designate $M_{t_0}=\Phi \left(T+t_0,t_0\right)$ then using (3) for $t=T+t_0$ ,

$\begin{array}{c}\Phi \left(T+t_0,t_0\right)=Q^{-1}\left(T+t_0\right)e^{BT}Q\left(t_0\right)\\ =Q^{-1}\left(t_0\right)e^{BT}Q\left(t_0\right)=Q^{-1}\left(t_0\right)MQ(\; t\; 0\; )\end{array}$

This shows that $M_{t_0}$ and $M$ are similar. Therefore as long as our use of the Monodromy matrix is reduced to its spectrum, we can use either $M$ or $M_{t_0}$ .

Stability

Recall the stability definition in the sense of Lyapunov [13] :

Definition 2. The zero solution of $\dot{y}=A\left(t\right)y$ is

1) Stable, if $\forall ϵ>0$ , $\exists \delta >0$ such that $\Vert y\left(t_0\right)\Vert <\delta \Rightarrow \Vert y\left(t_0\right)\Vert <ϵ$ , $\forall t\ge t_0$

2) Unstable if it is not stable.

¹To review the properties of the state transition matrix you can see [20] or [21] .

²See Adrianova pp.17 [2] , The necessary and sufficient condition for which the matrix S, has a real logarithm. Even more S² always has a real logarithm.

3) Asymptotically stable if the zero solution is stable and ${\lim }_{t\to \infty }y\left(t\right)=0$ .

In our system (2) $\dot{y}=A\left(t\right)y$ , for $t\ge 0$ , $t$ may be expressed as: $t=kT+\tau$ , where $k$ is a non-negative integer and $\tau \in \left[\mathrm{0,}T\right)$ ; then the solution satisfies (2) for $t_0=0$ and $y\left(0\right)=y_0$ :

$y\left(t\right)=\Phi \left(kT+\tau ,0\right)y_0=\Phi \left(\tau ,0\right)M^k{y}_0$

from which, we can conclude that, system (2) is:

1) Stable: $y\left(t\right)$ remains bounded $\forall t\ge 0$ if and only if $\sigma \left(M\right)\subset {\overline{D}}_1\triangleq \left\{z\in ℂ\mathrm{<mo>\; :\; </mo>}\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) Unstable if $\exists \mu \in \sigma \left(M\right):\left|\mu \right|>1$ or if $\sigma \left(M\right)\subset {\overline{D}}_1$ , $\exists \mu :\left|\mu \right|=1$ , $\mu$ is not a simple root of the minimal polynomial of $M$ .

2.2. Hamiltonian Systems

Hamiltonian systems consist in a pair of $2n$ ordinary differential equations of the form:

$\begin{array}{l}\dot{q}={\left(\frac{\partial \mathcal{H}\left(t\mathrm{,}q\mathrm{,}p\right)}{\partial p}\right)}^{\text{T}}\\ \dot{p}=-{\left(\frac{\partial \mathcal{H}\left(t\mathrm{,}q\mathrm{,}p\right)}{\partial q}\right)}^{\text{T}}\end{array}$ (6)

where $\mathcal{H}\left(t\mathrm{,}q\mathrm{,}p\right)$ is called a Hamiltonian function which represents the energy of the system, when the Hamiltonian function is being preserved constant along the solutions of (6) then the Hamiltonian system (6) is called a Conservative system, and this happens if $\mathcal{H}\left(t\mathrm{,}q\mathrm{,}p\right)$ is independent of time t.

In this work we will deal only with linear Hamiltonian systems, so the Hamiltonian function can be rewritten as a quadratic homogeneous form, i.e.

$\mathcal{H}\left(t,q,p\right)={\left[\begin{array}{c}q\\ p\end{array}\right]}^{\text{T}}H\left(t\right)\left[\begin{array}{c}q\\ p\end{array}\right]$ (7)

where $H\left(t\right)$ is a $2n\times 2n$ symmetric matrix, in this case the Hamiltonian System (7) may be expressed as:

$\frac{\text{d}}{\text{d}t}\left[\begin{array}{c}q\\ p\end{array}\right]=JH\left(t\right)\left[\begin{array}{c}q\\ p\end{array}\right]$ (8)

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

Notice that $J$ is orthogonal and skew-symmetric; i.e., $J^{-1}=J^{\text{T}}=-J$ and that $J^2=-I_{2n}$ .

Definition 3. Hamiltonian matrix [22]

An even-order matrix $A\in {ℝ}^{2n\times 2n}$ is called Hamiltonian Matrix, if

$A^{\text{T}}J+JA=0$ (9)

From $A^{\text{T}}J+JA=0$ , we get $A=J^{-1}\left(-A^{\text{T}}\right)J$ , i.e., $A$ is similar to $-A^{\text{T}}$ therefore they have the same spectrum:

$\sigma \left(A\right)=\sigma \left(-A^{\text{T}}\right)=\sigma (\; -\; A\; )$

In other words, if $\lambda \in \sigma \left(A\right)$ , then $-\lambda \in \sigma \left(A\right)$ .

Theorem 2. Let $A\in {ℝ}^{2n\times 2n}$ be a Hamiltonian matrix, then if $\lambda \in \sigma \left(A\right)\Rightarrow -\lambda \in \sigma \left(A\right)$ . Equivalently, the characteristic polynomial of a Hamiltonian matrix has only even powers or it is an even polynomial.

Definition 4. Symplectic Matrix [6]

An even order real matrix $S\in {ℝ}^{2n\times 2n}$ is called a Symplectic Matrix, if

$S^{\text{T}}JS=J$ (10)

The key property of constant symplectic matrices is that its spectrum is symmetric with respect to the unit circle, this can be proven from the definition and the fact that a symplectic matrix is always invertible, then $S^{\text{T}}=J{S}^{-1}{J}^{-1}$ , i.e. $\sigma \left(S^{\text{T}}\right)=\sigma \left(S^{-1}\right)=\sigma \left(S\right)$ $\Rightarrow$ if $\lambda \in \sigma \left(S\right)$ then ${\lambda }^{-1}\in \sigma \left(S\right)$ .

Theorem 3. Let $S\in {ℝ}^{2n\times 2n}$ be a symplectic matrix, then if $\lambda \in \sigma \left(A\right)\Rightarrow {\lambda }^{-1}\in \sigma \left(S\right)$ . Equivalently, the characteristic polynomial of a Symplectic matrix is self-reciprocal, i.e. $p_S\left(\lambda \right)={\lambda }^{2n}{p}_S\left({\lambda }^{-1}\right)$

The Figure 1 shows graphically the properties of Theorems 2 and 3.

Definition 5. A continuous Linear time invariant system $\dot{x}\left(t\right)=Ax\left(t\right)$ is said to be a Hamiltonian system if and only if $A$ is a Hamiltonian matrix [22] .

Definition 6. Similarly $x_{k+1}=B{x}_k$ is a discrete Hamiltonian system if and only if $B$ is a Symplectic matrix [23] .

Theorem 4. [22] The state transition matrix $\Phi \left(t\mathrm{,}{t}_0\right)$ is symplectic, for linear Hamiltonian system $\dot{x}\left(t\right)=A\left(t\right)x\left(t\right)$ .

2.3. Hill’s Equation: Continuous Time Case

Here we will see that the Hill’s equation is a Hamiltonian system, we will talk about the discriminant of the Monodromy matrix $\varphi \left(\alpha \mathrm{,}\beta \right)$ and its relation to the roots of the Monodromy matrix, we will enunciate the Hochstadt theorem, that gives intervals of stability and unstability, this can be easily seen on the Arnold tongues. We give a formula of where the Arnold tongues begin.

Lemma 5. The Hill’s equation

$\dot{y}\left(t\right)=A\left(t\right)y(\; t\; )$

where $A\left(t\right)=\left[\begin{array}{cc}0& 1\\ -\left(\alpha +\beta p\left(t\right)\right)& 0\end{array}\right]$ , $A\left(t+T\right)=A\left(t\right)$ is a Hamiltonian system.

Proof. according with the definition 5 we only need to prove that the matrix $A\left(t\right)$ is a Hamiltonian matrix

$A^{\text{T}}\left(t\right)J={\left[\begin{array}{cc}0& 1\\ -\left(\alpha +\beta p\left(t\right)\right)& 0\end{array}\right]}^{\text{T}}\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]=\left[\begin{array}{cc}\left(\alpha +\beta p\left(t\right)\right)& 0\\ 0& -1\end{array}\right]$

$JA\left(t\right)=\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]\left[\begin{array}{cc}0& 1\\ -\left(\alpha +\beta p\left(t\right)\right)& 0\end{array}\right]=\left[\begin{array}{cc}-\left(\alpha +\beta p\left(t\right)\right)& 0\\ 0& 1\end{array}\right]$

$A^{\text{T}}\left(t\right)J+JA\left(t\right)=0$

Corollary 2. The Monodromy Matrix of the Hill’s equation

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

is a symplectic matrix.

As a consequence of this Corollary, Hill’s equation can not be asymptotically stable, it may be stable (bounded) or unstable.

Definition 7. The eigenvalues of the Monodromy matrix $M$ , equivalently the roots of its characteristic polynomial $p_M(\; \lambda \; )$

$p_M\left(\lambda \right)={\lambda }^2-tr\left(M\right)\lambda +1$ (11)

are called multipliers of the Hill’s Equations (1) or (2), denoted by $\lambda$

Definition 8. Associated to every multiplier $\lambda$ , there exist (an infinite) numbers called characteristic exponents $\mu$ related to a multiplier by $\lambda =e^{\left(\mu +\frac{2\text{π}jk}{T}\right)T}$ , $k\in \mathbb{Z}$ .

The roots of $p_M\left(\lambda \right)$ are the multipliers of (1):

${\lambda }_{1,2}=\frac{tr\left(M\right)\pm \sqrt{t{r}^2\left(M\right)-4}}{2}$ (12)

Ÿ If $t{r}^2\left(M\right)<4$ the multipliers are complex conjugates with modulus =1. This case corresponds to a stable system.

Ÿ If $t{r}^2\left(M\right)>4$ , the multipliers are real and reciprocal and as a consequence one of the eigenvalues will be greater than one therefore this case corresponds to an unstable system.

Ÿ $t{r}^2\left(M\right)=4$ the multipliers are real and repeated $=+1$ or $=-1$ . In this case Hill Equation is stable if only if $M$ is a diagonal matrix, otherwise the Hill Equation (1) is unstable.

The boundaries between stability-instability correspond to this last case, i.e. when $\left|tr\left(M\right)\right|=2$ . It is clear that $M$ depends on the parameters $\alpha$ , $\beta$ . It is customary to define [1] $\varphi \left(\alpha \mathrm{,}\beta \right)\triangleq tr\left(M\right)$ , Hochstadt [24] was the first to recognize the important properties of $\varphi \left(\alpha \mathrm{,}\beta \right)$ .

Theorem 6. Hochstadt [24] The function $\varphi \left(\alpha \mathrm{,}\beta \right)$ for any $\beta$ constant. The functions $\varphi \left(\alpha ,\beta \right)\pm 2=0$ have an infinite number of roots. For any ${\beta }_0$ , and for ${\alpha }_0$ sufficiently negative, $\varphi \left({\alpha }_0\mathrm{,}{\beta }_0\right)$ is positive, therefore increasing $\alpha$ appears the first root for the equation $\varphi \left(\alpha ,\beta \right)-2=0$ , which corresponds to a double multiplier at $+1$ , and from there appear two roots (not necessarily different) at $-1$ , then two roots $+1$ , up to infinity.

Due to the Hochstadt Theorem 6, there are two infinite sequences:

$\begin{array}{l}{\lambda }_0\mathrm{,}{\lambda }_1\mathrm{,}{\lambda }_2\mathrm{,}{\lambda }_3\mathrm{,}{\lambda }_4\mathrm{,}{\lambda }_5\mathrm{,}\cdots \\ {\overline{\lambda }}_1\mathrm{,}{\overline{\lambda }}_2\mathrm{,}{\overline{\lambda }}_3\mathrm{,}{\overline{\lambda }}_4\mathrm{,}{\overline{\lambda }}_5\mathrm{,}\cdots \end{array}$ (13)

Moreover they interlace as:

${\lambda }_0\mathrm{,}{\overline{\lambda }}_1\mathrm{,}{\overline{\lambda }}_2\mathrm{,}{\lambda }_1\mathrm{,}{\lambda }_2\mathrm{,}{\overline{\lambda }}_3\mathrm{,}{\overline{\lambda }}_4\mathrm{,}{\lambda }_3\mathrm{,}{\lambda }_4\mathrm{,}{\overline{\lambda }}_5\mathrm{,}{\overline{\lambda }}_6\mathrm{,}\cdots$ (14)

This fact is illustrated in Figure 2

For the values of $\alpha$ in which $\left|\varphi \left(\alpha ,1.5\right)\right|>2$ are shown in grey color and correspond to the unstable zones.

Remark 2. The unstable regions in Figure 2 are also known as Arnold Tongues [25] , and they are labeled from left to right with $\mathrm{0,1,2,3,}\cdots$ .

Note that the boundary of the 0-th Arnold tongue has a T-periodic solution since the discriminant $\varphi \left(\alpha ,\beta \right)-2=0$ has an $\alpha$ root on its boundary surface, the first Arnold tongue has 2T-periodic solutions, the second has T-periodic solutions and so on. In such a way that there is an alternation between T-periodic solutions and 2T-periodic solutions, in consecutive boundaries of the Arnold tongues.

Theorem 7. The Arnold tongues associated to the Hill’s Equation (1) begin at $\alpha ={\left(\frac{k\text{π}}{T}\right)}^2$ , $k=0,1,2,\cdots$ .

Proof. Let’s consider the case when $\beta =0$ in (2)

$\dot{y}\left(t\right)=A\left(t\right)y(\; t\; )$

then $A\left(t\right)=\left[\begin{array}{cc}0& 1\\ -\left(\alpha +\beta p\left(t\right)\right)& 0\end{array}\right]$ will be $A=\left[\begin{array}{cc}0& 1\\ -\alpha & 0\end{array}\right]$ , with the general solution:

$y\left(t\right)=e^{At}y(\; 0\; )$

If we are looking for periodic solutions:

$y\left(t\right)=y\left(t+T\right)=e^{A\left(t+T\right)}y\left(0+T\right)=e^{At}{e}^{AT}t\left(0\right)\Rightarrow e^{AT}=I$

then:

$e^{AT}=\left[\begin{array}{cc}\cos \left(\sqrt{\alpha }T\right)& \frac{1}{\sqrt{\alpha }}\sin \left(\sqrt{\alpha }T\right)\\ -\frac{\alpha }{\sqrt{\alpha }}\sin \left(\sqrt{\alpha }T\right)& \cos \left(\sqrt{\alpha }T\right)\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$

from where:

$\cos \left(\sqrt{\alpha }T\right)=1$

$\sqrt{\alpha }T=k\text{π},k=0,1,2,\cdots$

$\therefore$ the Arnold tongues of any Hill’s Equation (1) begin at:

$\alpha ={\left(\frac{k\text{π}}{T}\right)}^2,k=0,1,2,3,\cdots$ (15)

When the periodic function $p\left(t\right)=\text{sign}\left(\cos \left(t\right)\right)$ of the Hill’s Equation (1) is known as the Meissner’s equation, with a period $T=2\text{π}$ the points where the

Arnold tongues are born given the previous theorem are $\alpha$ : 0, $\frac{1}{4}$ , 1, $\frac{9}{4}$ , 4, $\cdots$ as is shown in the Figure 3.

Remark 3. Notice that according with (15) and theorem 6 generically, there is an infinite number of Arnold tongues of the Hill’s Equation (1) in continuous time.

The classical problem of Hill is to find the values of $\alpha$ , $\beta$ for which (1) has periodic solutions, i.e., to determine the boundaries of the Arnold tongues of (1).

Definition 9. [4] A solution of a differential equation $\dot{y}\left(t\right)=A\left(t\right)y\left(t\right)$ as in (2), where $A\left(t\right)=\left[\begin{array}{cc}0& 1\\ -\left(\alpha +\beta p\left(t\right)\right)& 0\end{array}\right]$ , $A\left(t+T\right)=A\left(t\right)$ , is said to be nono-

scillatory, if it has at most one zero on $\left(-\infty \mathrm{,}+\infty \right)$ . If it has an infinite number of zeros without accumulation points of the zeros on $\left(-\infty \mathrm{,}+\infty \right)$ , is said to be oscillatory.

Lemma 8. [4] [5] Nonoscillation criteria for Hill’s equation For each fixed $\beta$ of (1), there is an ${\alpha }_0=\alpha \left(\beta \right)$ , for which if

Ÿ $\alpha \le {\alpha }_0$ all solutions of (1) are nonoscillatory.

Ÿ $\alpha >{\alpha }_0$ all solutions of (1) are oscillatory.

2.4. Discretization Methods

Consider the following system

$\dot{x}\left(t\right)=A\left(t\right)x\left(t\right)$ (16)

When we apply a discretization method we have to change $t$ by $kh$ , where $h$ is the sampling interval and $k\in \mathbb{Z}$ is the discrete time. By notation we will omit the sampling interval i.e. $x_k$ instead of $x\left(kh\right)$ .

Definition 10. The discretization of (16) by the method of forward Euler is defined by [26] [27]

$x_{k+1}-x_k=h{A}_k{x}_k$ (17)

where $A_k=A(\; k\; h\; )$

Definition 11. The discretization of (16) by the method of backward Euler is defined by [26] [27]

$x_k-x_{k-1}=h{A}_k{x}_k$ (18)

Definition 12. Tustin method [26]

The discretization of (16) by the Tustin is defined by

$x_{k+1}-x_k=\frac{1}{2}h{A}_{k+1}{x}_{k+1}+\frac{1}{2}h{A}_k{x}_k$

$x_{k+1}={\left[I-\frac{1}{2}h{A}_{k+1}\right]}^{-1}\left[I+\frac{1}{2}h{A}_k\right]x_k$ (19)

Definition 13. Pole-Zero matching equivalence³

³Even if the method mention zero matching we will only care about the pole matching because in our linear Hamiltonian system we don’t have the matrix $B$ , $C$ that are related with the zeros.

Let $a$ be a pole of the system

$\dot{x}\left(t\right)=Ax(\; t\; )$

the technique of pole zero Matching Equivalence [26] consists on setting $e^{ah}$ as a pole of the discrete system.

The Zero-order hold consists in sampling a continuous signal and hold it during the sampling interval so it is defined by:

Definition 14. Zero-order hold [26]

The discrete system of (16) is defined by⁴

$x_{k+1}=\Phi \left(k+1,k\right)x_k$ (20)

where $\Phi \left(t\mathrm{,}{t}_0\right)$ is the state transition matrix of (16).

2.5. Solution of Linear Difference Equations

Since we are dealing with discrete systems we mention a theorem that gives the structure of the solution of a difference equation.

Consider the following difference linear equation

$u\left(k+n\right)+p_{n-1}u\left(k+n-1\right)+\cdots +p_{0}u\left(k\right)=0$ (21)

where $p_0,\cdots ,p_{n-1}$ are constants and $p_0\ne 0$ .

Theorem 9. [28]

Ÿ Suppose that (21) has characteristic roots ${\lambda }_1,\cdots ,{\lambda }_{n_1}$ with multiplicities ${\alpha }_1,\cdots ,{\alpha }_{n_1}$ respectively. Then (21) has the n independent solutions ${\lambda }_1^k$ , $\cdots$ , $k^{{\alpha }_1-1}{\lambda }_1^k$ , ${\lambda }_2^k$ , $\cdots$ , $k^{{\alpha }_2-1}{\lambda }_2^k$ , $\cdots$ , ${\lambda }_{n_1}^k$ , $\cdots$ , $k^{{\alpha }_{n_1}-1}{\lambda }_{n_1}^k$ .

Ÿ If the characteristic roots include a complex pair $\lambda =a\pm b$ , then real-valued solutions of (21) can be found by using polar form $\lambda =r{e}^{\pm i\theta }=r\left(\cos \theta \pm i\sin \theta \right)$ then ${\lambda }^k=r^k{e}^{\pm i\theta k}=r^k\left(\cos \theta \pm i\sin \theta \right)$ .

3. Discretization of Linear Hamiltonian Systems

⁴Remember that even when we write $\Phi \left(k+\mathrm{1,}k\right)$ actually is $\Phi \left(h\left(k+1\right)\mathrm{,}hk\right)$ .

In this section we will mention some methods that preserves the Hamiltonian structure, i.e. that after discretizing a continuous time Hamiltonian system give us as a result a discrete time Hamiltonian system, according with definitions 5 and 6. Opposite to what we could think not all the discretization methods preserve the Hamiltonian structure.

Consider the following Hamiltonian system

$\dot{x}\left(t\right)=\lambda JHx\left(t\right)=\lambda Ex\left(t\right)$ (22)

where $x=\left[\begin{array}{c}y\\ z\end{array}\right]$ , $J=\left[\begin{array}{cc}0& I\\ -I& 0\end{array}\right]$ , $\lambda \in {ℝ}_{+}$ , $H^{\text{T}}=H=\left[\begin{array}{cc}A& B^{\text{T}}\\ B& D\end{array}\right]$ , $E=\left[\begin{array}{cc}B& D\\ -A& -B^{\text{T}}\end{array}\right]$ .

3.1. Rasvan’s Procedure [29]

Using the discretization method of Euler with step $h$ , but using forward Euler in the first equation and backward Euler in the second equation [29] , [30]

$\begin{array}{l}\frac{y\left(\left(k+1\right)h\right)-y\left(kh\right)}{h}=\lambda By\left(kh\right)+\lambda Dz\left(kh\right)\\ \frac{z\left(kh\right)-z\left(\left(k-1\right)h\right)}{h}=-\lambda Ay\left(kh\right)-\lambda B^{\text{T}}z\left(kh\right)\end{array}$ (23)

will result into the system

$x_{k+1}=C\left(\lambda \right)x_k$

where $x=\left(\begin{array}{c}y\\ z\end{array}\right)$ , $C\left(\lambda \right)={\left(\begin{array}{cc}I& -\lambda D\\ 0& I+\lambda B^T\end{array}\right)}^{-1}\left(\begin{array}{cc}I+\lambda B& 0\\ -\lambda A& I\end{array}\right)$ .

Lemma 10. The matrix $C\left(\lambda \right)$ is symplectic [29] [30] .

3.2. Comparison of Discretization Methods

Consider the following Hamiltonian system

$\dot{x}\left(t\right)=JHx\left(t\right)=Ax\left(t\right)$ (24)

where $H=H^{\text{T}}\in {ℝ}^{2n\times 2n}$ and $J$ as in (22).

When a discretization method is applied to a Hamiltonian system (24) it will take the following structure

$x_{k+1}=R{x}_k$ (25)

In the Table 1 we compare some discretization method, looking for those who preserve the Hamiltonian structure.

The proof can be found in [31] or by verifying that $R^{\text{T}}JR=J$ for the given discretization method.

3.3. Example

Consider the following differential equation

$\dot{x}\left(t\right)=Ax\left(t\right)$ (26)

where $A$ is a Hamiltonian matrix.

Example 1. Let the matrix $A$ of the system (26) have the following spectrum $\sigma \left(A\right)=\left\{\pm 10i,\pm 50i,\pm 110i\right\}$ note here that different to the previous example we have just “stable” eigenvalues with zero real part.

However in Figure 4 the forward Euler method gives us an unstable discrete system and the backward Euler method gives us an asymptotically stable discrete system. Furthermore as we expected the Rasvan’s procedure, the Tustin method and

the Pole matching⁵ give us as a result a “stable” discrete system but with some differences between Rasvan, Tustin and Pole matching method related with the step $h$ [27] . The Figure 5 shows the response to some initial conditions of the Hamiltonian system of example 1, and also the different responses to the same initial conditions for the different methods analyzed previously.

⁵When we mention Pole matching we mean Pole Zero matching and we are also referring to the Zero order hold as it gives us the same discrete system.

In the discrete time approximation the forward Euler method gives us an unstable system and the backward Euler method gives us an asymptotically stable system, only the Rasvan’s procedure, the Tustin method and pole matching give us a stable system as the continuous time approach moreover this information is consistent with the discrete eigenvalue stability information.

4. Discretized Hill’s Equation

In this section we will discretize the Hill’s equation in the integral equation form

$y\left(T\right)={\displaystyle {\int }_0^T}A\left(\tau \right)y\left(\tau \right)\text{d}\tau +y\left(0\right)$ (27)

and then analyze the result as a discrete system without forgetting that we are discretizing a continuous time system therefore the resultants Arnold’s tongues will correspond to the continuous time Hill’s equation. We can use any discretization method as long as it preserves the Hamiltonian structure of the Hill equation.

By chosing the step $h=\frac{T}{K}$ , such that we can set a discrete period $K$ to discretize the continuous time Hill’s Equation (2) and obtain a discrete system of the following form

$y_{\left(k+1\right)h}=A_{kh}{y}_{kh}$ (28)

where $A_{\left(k+K\right)h}=A_{kh}$ , $T$ is the continuous time period of the Hill’s Equation (2). For notation we will drop the $h$ and write only $y_{k+1}=A_k{y}_k$ , $A_{k+K}=A_k$ , $\forall k\in \mathbb{Z}$ .

Theorem 11 (Lifting technique). Consider the following discrete time periodic system

$y_{k+1}=A_k{y}_k,A_{k+K}=A_k$ (29)

the system (29) is equivalent to a linear invariant time discrete system

$x_{k+1}=\overline{A}{x}_k$

Proof. Let’s define the following discrete vector state $x_k=\left[\begin{array}{cccc}{y}_{Kn}& y_{Kn+1}& \cdots & y_{Kn+K-1}\end{array}\right]$ , then we can transform (29) in this⁶

$x_{\left(n+1\right)}=\overline{A}{x}_n$ (30)

where $\overline{A}$ is a constant diagonal blocks matrix,

$\overline{A}=\text{diag}\left(A_{K-1}\cdots A_1{A}_0,\cdots ,A_{K-2}\cdots A_0{A}_{K-1}\right)$ .

Corollary 3. Let $M=\Phi \left(K,0\right)=A_{K-1}{A}_{K-2}\cdots A_1{A}_0$ be the Monodromy matrix of the system (29), then the system (29) is:

Ÿ stable if $\sigma \left(M\right)\subset {\overline{D}}_1\triangleq \left\{z\in ℂ:\left|z\right|=1\right\}$ & $\mu \in \sigma \left(M\right)$ & $\left|\mu \right|=1$ , $\mu$ is a simple root of the minimal polynomial $M$ .

Ÿ Unstable if $\exists \mu \in \sigma \left(M\right):\left|\mu \right|>1$ or if $\sigma \left(M\right)\subset {\overline{D}}_1$ , $\exists \mu :\left|\mu \right|=1$ & is not a simple root of the minimal polynomial.

Proof.

$\det \left(\lambda I-\overline{A}\right)={\displaystyle \prod _{i=1}^K}\det \left(\lambda I-{\overline{A}}_{ii}\right)$

⁶This procedure is known as lifting technique [32] .

on the other hand

$\begin{array}{c}\det \left(\lambda I-A_{K-1}\cdots A_1{A}_0\right)=\det \left(\lambda I-A_0{A}_{K-1}\cdots A_1\right)=\cdots \\ =\det \left(\lambda I-A_{K-2}\cdots A_0{A}_{K-1}\right)\end{array}$

which implies that

$\sigma \left(A_{K-1}\cdots A_1{A}_0\right)=\sigma \left(A_0{A}_{K-1}\cdots A_1\right)=\cdots =\sigma \left(A_{K-2}\cdots A_0{A}_{K-1}\right)$

therefore the stability of (29) is given by $M=A_{K-1}\cdots A_1{A}_0$ .

Example

Example 2. Consider the following system

$\dot{y}\left(t\right)=A\left(t\right)y\left(t\right)$ (31)

where $A\left(t\right)=\left[\begin{array}{cc}0& 1\\ -\left(\alpha +\beta p\left(t\right)\right)& 0\end{array}\right]$ , $\alpha$ , $\beta$ are constant, with the periodic

function $p\left(t\right)=\text{sign}\left(\cos \left(t\right)\right)$ we obtain the Figure 6, where the white zones shows the stability zones and the grey zones denote instability zones.

By fixing $\left(\alpha \mathrm{,}\beta \right)$ and using the theorem 11 to discretize the system (31) and using the corollary 3 to obtain the Monodromy matrix and to analyze its spectrum, plotting a grey dot when the pair $\left(\alpha \mathrm{,}\beta \right)$ is unstable, and by repeating this process while doing a sweep on the pair $\left(\alpha \mathrm{,}\beta \right)$ we obtain the Arnold tongues (Figure 6).

Remark 4. If we have in mind that we need to calculate the Monodromy matrix and its eigenvalues for a grid of points $\alpha$ and $\beta$ in order to plot the Arnold’s tongues. So let’s say that computing the Monodromy matrix for one point $\left(\alpha \mathrm{,}\beta \right)$ takes 0.005 seconds, for a resolution of points on alpha and beta of 1000 point this would take approximately 13.8 hours to compute and plot the Arnold tongues, however since this computational process can be parallelized we use GPUs to reduce the computational time from hours to seconds.

Notice that the objective of discretizing the Hill’s equation was to recover partially some properties of the continuous time system.

5. Hill’s Equation: Discrete Time Case

Unlike the discretized case, the discrete Hill’s equation has some properties completely different, that we will describe for the first time.

In this section we will study the Discrete Hill’s equation presented among others by Chulaevsky in 1989 [33] ⁷ and we will give it a graphical interpretation of parametric stability for the first time i.e. discrete Arnold tongues [25] . We will give a formula of where the discrete Arnold’s tongues begin and proof the nonoscilation criteria for discrete Hill’s equation.

Consider the following difference equation

$q_{k+2}+\left(\alpha +\beta p\left(k\right)\right)q_{k+1}+q_k=0$ (32)

named here as the Discrete Hill’s equation [33] , where $k\mathrm{,}K\in \mathbb{Z}$ , $p\left(k+K\right)=p\left(k\right)$ , is a sequence of period $K$ and of zero average, i.e. ${\displaystyle \sum _{k=0}^{K-1}}p\left(k\right)=0$ . It may be obtained by discrete variational methods [34] .

As in the case of continuous time, we can define the state vector $y_k=\left[\begin{array}{cc}{q}_{k+1}& q_k\end{array}\right]$ to rewrite the Equation (32) as

$y_{k+1}=A_k{y}_k$ (33)

where $A_k=\left[\begin{array}{cc}-\left(\alpha +\beta p\left(k\right)\right)& -1\\ 1& 0\end{array}\right]$ , $A_{k+K}=A_k$ .

Remark 5. Notice that $A_k$ is a symplectic matrix, i.e.

$A_k^{\text{T}}J{A}_k=J$

where $J=\left[\begin{array}{cc}0& I\\ -I& 0\end{array}\right]$ $\therefore$ (33) is a discrete hamiltonian system [22] .

Remark 6. Notice that since $A_{k+K}=A_k$ on (33) is periodic we can use the theorem 11 and the Corollary 3 to analyze the stability conditions of (33).

5.1. Examples

⁷It may be proven that this Equation (32) results to apply a discrete Variational Method to a discrete Lagrangian through a discrete Euler-Lagrange equation, this will be reported in another paper.

Example 3. Consider the discrete Hill equation $q_{k+2}+\left(\alpha +\beta p\left(k\right)\right)q_{k+1}+q_k=0$ as in (32) for some sequences, going from $K=1$ up to $K=6$ , the following Figure 7 and Figure 8 shows the parametric stability or Arnold Tongues, of the Discrete Hill’s equation.

Remark 7. For $K$ even there is a symmetry with respect to both axes; for $K$ odd the symmetry with respect to $\beta =0$ is preserved, but appear a symmetry with respect to the origin.

Where the white zones shows the stability zones and the grey zones denote instability zones or Arnold tongues labeled from left to right with $0,1,\cdots ,K-1$ and $K$ .

Remark 8. Notice that in continuous time generically there is an infinite number of tongues, this however is not true for the discrete case, in the discrete

case we have generically, only $K+1$ tongues.

Theorem 12. The Arnold’s tongues of the discrete Hill’s equation $y_{k+1}=A_k{y}_k$ as (33), where $A_k=\left[\begin{array}{cc}-\left(\alpha +\beta p\left(k\right)\right)& -1\\ 1& 0\end{array}\right]$ , $A_{k+K}=A_k$ , begin at

$\left|l_0{\left(-\alpha \right)}^K+l_1{\left(-\alpha \right)}^{K-2}+l_2{\left(-\alpha \right)}^{K-4}+\cdots +l_{\frac{K-i}{2}}{\left(-\alpha \right)}^i\right|=2$ (34)

where

$i=\{\begin{array}{ll}1\hfill & \text{if}K\text{is}\text{even}\hfill \\ 0\hfill & \text{if}K\text{is}\text{odd}\hfill \end{array}$

$l_0=1$

$l_1=-(\; K\; )$

$l_2=-\frac{\left(K\right)\left(K-1\right)}{2!}-\frac{\left(K-2\right)}{1!}{l}_1$

$l_3=-\frac{\left(K\right)\left(K-1\right)\left(K-2\right)}{3!}-\frac{\left(K-2\right)\left(K-3\right)}{2!}{l}_1-\frac{\left(K-4\right)}{1!}{l}_2$

$⋮$

$\begin{array}{c}{l}_{\frac{K-i}{2}}=-\frac{\left(K\right)\left(K-1\right)\cdots \left(K-\frac{K-i}{2}-1\right)}{\frac{K-i}{2}!}\\ -\frac{\left(K-2\right)\cdots \left(K-\frac{K-i}{2}-2\right)}{\left(\frac{K-i}{2}-1\right)!}{l}_1-\cdots -\frac{2+i}{1!}{l}_{\left(\frac{K-i}{2}-1\right)}\end{array}$

Proof. We are looking for K-periodic solutions of $y_{k+1}=A_k{y}_k$ , i.e.

$y_k=y_{k+K}=\Phi \left(k+K,0\right)y_0=\Phi \left(k+K,K\right)\Phi \left(K,0\right)y_0=\Phi \left(k,0\right)M{y}_0$

since the tongues begin at $\beta =0$ , let’s analyze $M$ for $\beta =0$ .

$M={A_{K-1}{A}_{K-2}\cdots A_0|}_{\beta =0}={\tilde{A}}^K$

where

$\tilde{A}=\left[\begin{array}{cc}-\alpha & -1\\ 1& 0\end{array}\right]$ ,

$\sigma \left(\tilde{A}\right)=\left\{{\lambda }_1=-\frac{1}{2}\alpha -\frac{1}{2}i\sqrt{4-{\alpha }^2},{\lambda }_2=-\frac{1}{2}\alpha +\frac{1}{2}i\sqrt{4-{\alpha }^2}\right\}$ ,

Since $M$ is a symplectic matrix of dimensions $2\times 2$ , its characteristic polynomial is as follows

$P_M\left(\lambda \right)={\lambda }^2-Tr\left({\tilde{A}}^K\right)\lambda +1$

so the Arnold tongues will begin at the stability boundary $\left|Tr\left({\tilde{A}}^K\right)\right|=2$

Since $\sigma \left({\tilde{A}}^K\right)=\left\{{\lambda }_1^K,{\lambda }_2^K\right\}$

$Tr\left({\tilde{A}}^K\right)={\lambda }_1^K+{\lambda }_2^K$ (35)

notice that $-\alpha ={\lambda }_1+{\lambda }_2$ , then by substituting in (35) we obtain the polynomial in (34).

Example 4. For a period $K=5$

$Tr\left({\tilde{A}}^5\right)=-{\alpha }^5-l_1{\alpha }^3-l_2\alpha$ (36)

with $l_1=-5$ , $l_2=-\frac{5\times 4}{2}-3\left(-5\right)=5$ , then

$-{\alpha }^5+5{\alpha }^3-5\alpha =2$

has the periodic roots on $\alpha =\left\{-2,-0.618,-0.618,1.618,1.618\right\}$ and

$-{\alpha }^5+5{\alpha }^3-5\alpha =-2$

has the anti-periodic roots on $\alpha =\left\{-1.618,-1.618,0.618,0.618,2\right\}$ .

Remark 9. It may be observed that the boundaries of the Arnold Tongues are K-periodic or 2K-periodic; alternating and starting with K-periodic (Figure 9).

Remark 10. Worth notice that all the Arnold tongues with the same period $K$ begin at the same points, since the Equation (34) is independent of the sequence $p\left(k\right)$ (Figure 10).

5.2. Nonoscillation Criteria for Discrete Hill’s Equation

We have found that the Nonoscillation Criteria for continuous time Hill’s equation is also fulfilled for the discrete Hill Equation (32).

Definition 15. A solution of a difference equation $y_{k+1}=A_k{y}_k$ as in (33), where $A_k=\left[\begin{array}{cc}-\left(\alpha +\beta p\left(k\right)\right)& -1\\ 1& 0\end{array}\right]$ , $A_{k+K}=A_k$ , is said to be nonoscillatory if

there is at most one change of sign in the solution of (33). If there are an infinite number of zeros of the solutions, they are called Oscillatory.

Theorem 13. For each fixed $\beta$ of (33), there is an ${\alpha }_0=\alpha \left(\beta \right)$ , for which if

• $\alpha <{\alpha }_0$ all solutions of (32) are nonoscillatory.

• $\alpha \ge {\alpha }_0$ all solutions of (32) are oscillatory.

Even more all the nonoscillatory solutions of (32) fall into the 0-th Arnold tongue.

Proof. Consider first the case with $\beta =0$

$q_{k+2}+\alpha q_{k+1}+q_k=0$ (37)

its characteristic polynomial is

$\left({\lambda }^2+\alpha \lambda +1\right)q_k=0$ (38)

with the following general solutions, using the Theorem 9:

• for $\alpha <-2$ the general solution of (37) is Nonoscillatory

$q_k=C{\left(\frac{\alpha }{2}-\frac{\sqrt{{\alpha }^2-4}}{2}\right)}^k+D{\left(\frac{\alpha }{2}+\frac{\sqrt{{\alpha }^2-4}}{2}\right)}^k$

• for $\alpha =-2$ the general solution of (37) is Nonoscillatory

$q_k=C{1}^k+Dt{1}^k$

• for $\left|\alpha \right|<2$ the general solution of (37) is oscillatory

$q_k=C{r}^k{e}^{i{\theta }_k}+D{r}^k{e}^{-i{\theta }_k}$

where $r$ is the magnitude of $\lambda$ , $\therefore r=1$ .

• for $\alpha =2$ the general solution of (37) is oscillatory

$q_k=C{\left(-1\right)}^k+Dk{\left(-1\right)}^k$

• for $\alpha >2$ the general solution of (37) is oscillatory

$q_k=C{\left(-\frac{\alpha }{2}-\frac{\sqrt{{\alpha }^2-4}}{2}\right)}^k+D{\left(-\frac{\alpha }{2}+\frac{\sqrt{{\alpha }^2-4}}{2}\right)}^k$

where $C$ , $D$ are constants, here we see that ${\alpha }_0=-2$ , for $\beta =0$ .

Let us go back to the case with $\beta \ne 0$ .

Due to corollary 3 the Equation (32) is equivalent to an invariant linear system with the same spectrum of the Monodromy matrix $M$ therefore the solutions of (32) are as follow

Let $\lambda \in \sigma \left(M\right)$ and $\alpha +\beta p\left(k\right)$ be such that:

• for $\lambda >1$ by hamiltonian properties the general solution of (32) will be Nonoscillatory

$q_k=C{\lambda }^k+D{\left({\lambda }^{-1}\right)}^k$

• for $\lambda =1$ by hamiltonian properties the general solution of (32) is Nonoscillatory

$q_k=C{1}^k+Dk{1}^k$

• for $\lambda =a+ib$ the general solution of (37) is oscillatory

$q_k=C{r}^k{e}^{i{\theta }_k}+D{r}^k{e}^{-i{\theta }_k}$

• for $\lambda =-1$ by hamiltonian properties the general solution of (32) is oscillatory

$q_k=C{\left(-1\right)}^k+Dt{\left(-1\right)}^k$

• for $\lambda <-1$ by hamiltonian properties and since $\lambda$ is negative, the general solution of (32) is oscillatory

$q_k=C{\left(\lambda \right)}^k+D{\left({\lambda }^{-1}\right)}^k$

where $C$ , $D$ are constants, notice here that the 0-th Arnold tongue lies in the case where $\lambda \ge 1$ i.e. all the nonoscillatory solutions of (32) fall into the 0-th Arnold tongue.

We have proved the next result:

Theorem 14. Given Equation (32), a discrete Hill equation, with a K-periodic coefficient $p\left(k\right)$ ; $k=0,1,2,\cdots ,K-1$ , then (32) is nonoscillatory if and only if $\left(\alpha \mathrm{,}\beta \right)$ belong to the 0-th Arnold tongue.

Example 5. To illustrate this result, let us use the Arnold tongues obtained in the Figure 8(d) of the example 3, then the Figure 11 shows the solutions $q_k$

over the time, for a different pair of $\left(\alpha \mathrm{,}\beta \right)$ with the initial condition $y_0={\left[\begin{array}{cc}1& 0\end{array}\right]}^{\text{T}}$ .

In the Figure 11 we can see that only the solution of the 0-th Arnold tongue has a nonoscillatory solution.

Let’s summarize the properties discussed in this work. So in the Table 2 we compare the discrete Hill’s equation and the continuous time Hill’s equation side by side.

6. Conclusions

In this paper we found that the Monodromy matrix is a symplectic matrix for both continuous Hill’s equation and discrete Hill’s equation.

It is important to use a discretization method that preserves the Hamiltonian structure otherwise the stability of the discretized system won’t correspond to the original system.

We found that unlike the case of continuous time Hill’s equation where there is an infinite number of tongues, in the discrete time case there is only a finite

number of them, this is mainly because of the nature of their state transition matrix, while one can have infinite roots of $\alpha$ the other can only have a finite number of roots of $\alpha$ . We found that the solutions of the Discrete time Hill’s equation are nonoscillatory if and only if the corresponding $\left(\alpha \mathrm{,}\beta \right)$ parameters belongs to the 0-th Arnold tongue.

We are leaving for future work the study of higher order discrete Hill’s equation since here we were focused on the one degree of freedom case.

Cite this paper

Servín, J.G.R. and Collado, M.J. (2017) Arnold Tongues for Discrete Hill’s Equation. Applied Mathematics, 8, 1859-1882. https://doi.org/10.4236/am.2017.812133

References