# Non-smooth stability analysis of the parametrically excited impact oscillator

Remco Leine

► To cite this version:

Remco Leine. Non-smooth stability analysis of the parametrically excited impact oscillator. International Journal of Non-Linear Mechanics, 2012, 47 (9), pp.1020 - 1032. 10.1016/j.ijnonlinmec.2012.06.010. hal-01403578

HAL Id: hal-01403578

https://hal.science/hal-01403578

Submitted on 26 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Non-smooth stability analysis of the parametrically excited impact oscillator

R.I. Leine

Institute of Mechanical Systems, Department of Mechanical and Process Engineering, ETH Zurich, CH-8092 Zurich, Switzerland

The aim of this paper is to give a Lyapunov stability analysis of a parametrically excited impact oscillator, i.e. a vertically driven pendulum which can collide with a support. The impact oscillator with parametric excitation is described by Hill's equation with a unilateral constraint. The unilaterally constrained Hill's equation is an archetypel of a parametrically excited non-smooth dynamical system with state jumps. The exact stability criteria of the unilaterally constrained Hill's equation are rigorously derived using Lyapunov techniques and are expressed in the properties of the fundamental solutions of the unconstrained Hill's equation. Furthermore, an asymptotic approximation method for the critical restitution coefficient is presented based on Hill's infinite determinant and this approximation can be made arbitrarily accurate. A comparison of numerical and theoretical results is presented for the unilaterally constrained Mathieu equation.

Keywords: Parametric excitation, Non-smooth dynamics, Impulsive dynamics, Hill's equation

## 1. Introduction

The objective of the paper is to give more insight into the stability properties of non-smooth dynamical systems with parametric excitation. The impact oscillator with parametric excitation is studied which is described by Hill's equation with a unilateral constraint.

The theory of parametrically excited systems has applications in a wide range of disciplines, e.g. the quadrupole ion trap [19], the exact plane wave solutions in general relativity [3], parametric amplifiers [1], rotor dynamical instabilities [23,24,28], parametric resonance in power transmission belts [20] and celestial mechanics [25]. The importance of parametric excitation has led to a wealth of literature on the theoretical and experimental analysis of parametrically excited systems, see [8,28] and references therein. The attention in the literature is mainly focussed on the dynamics of the planar vertically driven pendulum and on its linearization which is described by Hill's equation. The vertically driven pendulum, of which the suspension point is driven periodically up and down (see Fig. 1a), is one of the simplest mechanical systems with parametric excitation. The system is also known as the parametrically excited pendulum, the vertically forced pendulum or the Kapitza pendulum. The dynamics of the vertically driven pendulum can (after some scaling) be expressed by

$$ \ddot{\theta}(t) + \left( \frac{g}{l} - \frac{a(t)}{l} \right) \sin \theta(t) = 0, \quad (1) $$

where $\theta$ is the angle of the pendulum with the downward vertical and $a(t) = \ddot{u}(t)$ is the vertical acceleration of the suspension point.

Extensive studies of the non-linear dynamics of the vertically driven pendulum can be found in [2,4,5,8,26]. The dynamics of the pendulum system in the vicinity of its equilibrium positions $\theta=0$ and $\theta=\pi$ is described by Hill's equation [6,17,30]

$$ \ddot{y}(t) + g(t)y(t) = 0, \quad (2) $$

where $g(t) = g(t+\pi)$ is a real piece-wise continuous function. Harmonic excitation leads to the so-called Mathieu equation

$$ \ddot{y}(t) + (\alpha + 2\beta \cos 2t)y(t) = 0, \quad (3) $$

being a special case of Hill's equation with $g(t) = \alpha + 2\beta \cos 2t$.

Systems with some degree of non-smoothness or switching behavior are often referred to as non-smooth dynamical systems [7,10,15]. Mechanical systems with impact and/or friction form an important subclass of non-smooth dynamical systems. The stability properties of non-smooth (mechanical) systems are currently receiving much attention, see [16] and references therein. However, the stability of equilibria of parametrically excited non-smooth systems has hardly been addressed.

Following [22], various contributions study the dynamics of the planar pendulum with impact but with a horizontally driven suspension point and the system is therefore not parametrically excited. Moreover, the existence of an equilibrium is lost under horizontal excitation.

Ivanov [11] studies the dynamics of Hill's equation with an impulsive parametric excitation in which the function $g(t)$ contains Dirac functions. Although impulsive action is added to Hill's equation, the linearity of the system is maintained.

Quinn [21] gives an in-depth study of two parametrically excited pendula under vertical harmonic excitation which can collide with each other. A symmetric response of the pendula agrees with the dynamics of a single vertically driven pendulum (1) with a unilateral
<figure><img src="image_2.png" /><figcaption>Fig. 1. Vertically driven pendulum. (a) Unconstrained and (b) with unilateral constraints.</figcaption></figure>

constraint at $\vartheta=0$ and the impact law $\dot{\vartheta}^+ = -\varepsilon\dot{\vartheta}^-$. The method of averaging in amplitude and phase coordinates (Lagrange standard form) is used to describe the dynamics in the vicinity of the equilibrium and near a resonance frequency under the assumption of a small amplitude of the excitation $a(t)$. The analysis takes the non-linearity of the pendulum into account using the approximation $\sin \vartheta \approx \vartheta - \frac{1}{6}\vartheta^3$, that is to say, the 'nonlinear Mathieu equation' $\ddot{y}(t) + (\alpha + 2\beta \cos 2t)(y(t) - \frac{1}{6}y(t)^3) = 0$ is considered. The averaged equations of motion in amplitude and phase coordinates are used to derive an approximate impact event map which maps the state of the system at a collision time instant to the state at the following collision time instant. The results in [21] can be used to derive an approximate stability criterion for the equilibrium point of the (linear or non-linear) Mathieu equation with unilateral constraint. The drawback of the approach taken in [21] is that the approximation is only valid for small values of the excitation parameter $\beta$ and near a chosen resonance frequency $\sqrt{\alpha} = 1, 2, 3, \dots$.

In this paper, a detailed Lyapunov stability analysis is presented of the equilibrium of Hill's equation in the presence of a unilateral constraint with restitution, hereafter called the (unilaterally) constrained Hill's equation. The aim of this paper is twofold. Firstly, the exact stability criteria of the unilaterally constrained Hill's equation are rigorously derived using Lyapunov techniques and are expressed in the properties of the fundamental solutions of the unconstrained Hill's equation (2). Secondly, an asymptotic approximation method for the critical restitution coefficient is presented based on Hill's infinite determinant and this approximation can be made arbitrarily accurate.

Hill's equation in the presence of a unilateral constraint with restitution is described by

$$\ddot{x}(t) + g(t)x(t) = 0,$$

$$x(t_i) = 0: \quad \dot{x}^{+}(t_i) = -\varepsilon\dot{x}^{-}(t_i), \qquad (4)$$

where $\varepsilon$ is the restitution coefficient. The unilateral constraint limits the dynamics to $x(t) \ge 0$ and imposes a Newtonian impact law $\dot{x}^+ = -\varepsilon\dot{x}^-$. More physically correct would be to formulate the unilaterally constrained system (4) in the framework of non-smooth dynamics [7,10,15,16] and to let a contact force $\lambda$ and an impulsive contact force $\Lambda$ (both per unit mass) appear in the equation of motion and impact equation as Lagrange multipliers, i.e.

$$\ddot{x}(t) + g(t)x(t) = \lambda \quad \text{for almost all } t,$$

$$x(t_i) = 0: \quad \dot{x}^{+}(t_i) - \dot{x}^{-}(t_i) = \Lambda, \qquad (5)$$

together with the set-valued contact law (Signorini's law)

$$0 \le \lambda \perp x \ge 0 \qquad (6)$$

and the generalized Newtonian impact law

$$x(t_i) = 0: \quad 0 \le \Lambda \perp \dot{x}^{+}(t_i) + \varepsilon\dot{x}^{-}(t_i) \ge 0, \qquad (7)$$

where $0 \le a \perp b \ge 0$ denotes the inequality complementarity condition $a \ge 0, b \ge 0, ab=0$ [7,10,16]. Such a (mechanical) system with impulsive effects due to unilateral constraints may be conveniently cast in terms of a measure differential inclusion, e.g. see [7,16]. However, the contact force $\lambda$ is for this system always zero since the external force $g(t)x(t)$ and inertial force $\ddot{x}$ vanish during persistent contact for which $x(t) = \dot{x}(t) = \ddot{x}(t) = 0$. Furthermore, the constraint is always actively participating in the impact process as there is only one constraint in the system which implies that the equality $\dot{x}^{+}(t_i) + \varepsilon\dot{x}^{-}(t_i) = 0$ holds at a collision time instant $t_i$. Moreover, it will be shown in this paper that accumulation points of impact events (Zeno behaviour) cannot occur in this particular system. Hence, the system (5) with set-valued force laws (6) and (7) simplifies to the unilaterally constrained Hill's equation in the form of (4). The simpler form (4) has been chosen in this paper to describe the dynamics of the unilaterally constrained Hill's equation instead of the more general framework of non-smooth dynamics (involving measure differential inclusions) in order to improve the readability for a heterogenous audience.

The unilaterally constrained Hill's equation (4) describes the dynamics in the vicinity of the equilibrium positions of a vertically driven pendulum with a vertical wall limiting the angle $\vartheta$ to $0 \le \vartheta(t) \le \pi$, see Fig. 1b. Similarly, (4) can be considered to be the linearization of a unilaterally constrained Euler column under dynamic axial loading which can only deflect in the unconstrained direction and of which only the first bending mode is considered. The case with linear damping $\ddot{x}(t) + \alpha\ddot{x}(t) + \bar{g}(t)\bar{x}(t) = 0$ can easily be transformed to the standard form (4) by using the transformation $x(t) = e^{(1/2)\alpha t}\bar{x}(t)$ such that $g(t) = \bar{g}(t) - \frac{1}{4}a^2$ [17]. The unilaterally constrained Hill's equation (4) is therefore an prototype of a parametrically excited non-smooth dynamical system with state jumps.

The paper is organized in the following way. First, basic properties of Hill's equation (2) are reviewed in Section 2 and some novel results on the number of zeros of fundamental solutions are derived. Subsequently, the stability properties of the unilaterally constrained Hill's equation are studied in Section 3 using Lyapunov stability techniques and are expressed in the properties of the fundamental solutions of the unconstrained Hill's equation. This theoretical result opens the way to use standard approximation methods for the stability analysis of the unilaterally constrained Hill's equation. An approximation technique for the critical restitution coefficient based on Hill's infinite determinant is presented in Section 4. Finally, numerical results are given for the unilaterally constrained Mathieu equation in Section 5 and the numerically obtained Ince-Strutt diagram is compared with the approximations using Hill's infinite determinant and the averaging method. The paper closes with conclusions and discussion in Section 6.

## 2. Properties of Hill's equation

In this section some properties of the unconstrained Hill's equation (2) are derived or reviewed, which will be useful when analyzing the unilaterally constrained Hill's equation in Section 3. Hill's equation (2)

$$\ddot{y}(t) + g(t)y(t) = 0, \quad g(t) = g(t + \pi)$$

has two continuously differentiable solutions $y_1(t)$ and $y_2(t)$ with the initial conditions

$$y_1(0) = 1, \quad \dot{y}_1(0) = 0, \quad y_2(0) = 0, \quad \dot{y}_2(0) = 1,$$
which are usually referred to as normalized solutions or fundamental solutions. The following proposition (see [17]) proves that the fundamental solutions can conveniently be described in polar coordinates.

**Proposition 1** (Magnus & Winkler [17]). The fundamental solutions of (2) can be expressed in polar coordinates by

$$y_1(t) = \varrho(t) \cos \psi(t), \quad y_2(t) = \varrho(t) \sin \psi(t), \qquad (8)$$

with $\varrho(t) > 0$ for all $t$ and the differential equations

$$\ddot{\varrho}(t) - \frac{1}{\varrho(t)^3} + g(t)\varrho(t) = 0, \quad \psi(t) = \int_0^t \frac{dt}{\varrho(t)^2}, \qquad (9)$$

with initial conditions $\varrho(0) = 1$, $\dot{\varrho}(0) = 0$, $\psi(0) = 0$ and $\dot{\psi}(0) = 1$.

**Proof.** Substitution of the fundamental solutions $y_1(t)$ and $y_2(t)$ expressed in polar coordinates (8) in Hill's equation (2) gives two differential equations

$$(\ddot{\varrho} - \varrho\dot{\psi}^2 + g(t)\varrho) \cos \psi - (2\dot{\varrho}\dot{\psi} + \varrho\ddot{\psi}) \sin \psi = 0,$$

$$(\ddot{\varrho} - \varrho\dot{\psi}^2 + g(t)\varrho) \sin \psi + (2\dot{\varrho}\dot{\psi} + \varrho\ddot{\psi}) \cos \psi = 0.$$

It can therefore be deduced that both the terms $\ddot{\varrho} - \varrho\dot{\psi}^2 + g(t)\varrho = 0$ and $2\dot{\varrho}\dot{\psi} + \varrho\ddot{\psi} = 0$ vanish. If the latter expression is multiplied with $\varrho(t) > 0$, then the equality $d(\varrho^2\dot{\psi})/dt = 0$ is obtained. The product $\varrho^2\dot{\psi}$ is therefore constant. Using the initial condition $\varrho(0) = 1$, $\psi(0) = 0$ and $\dot{\psi}(0) = 1$ gives $\psi(t) = \int_0^t \varrho(t)^{-2} dt$ and substitution of $\dot{\psi} = \varrho^{-2}$ in the former equation yields $\ddot{\varrho} - \varrho^{-3} + g(t)\varrho = 0$. $\square$

Hill's equation can be put in first-order form as

$$\dot{\mathbf{y}}(t) = \mathbf{A}(t)\mathbf{y}(t), \quad \mathbf{A}(t) = \begin{pmatrix} 0 & 1 \\ -g(t) & 0 \end{pmatrix}, \qquad (10)$$

with the state vector $\mathbf{y}(t) = (y(t) \dot{y}(t))^{\mathrm{T}}$ and the time-dependent system matrix $\mathbf{A}(t)$. The fundamental solutions constitute the fundamental solution matrix

$$\mathbf{\Phi}(t,0) = \begin{pmatrix} y_1(t) & y_2(t) \\ \dot{y}_1(t) & \dot{y}_2(t) \end{pmatrix},$$

which is therefore the solution of the matrix differential equation $\mathbf{\Phi}(t,0) = \mathbf{A}(t)\mathbf{\Phi}(t,0)$ with initial condition $\mathbf{\Phi}(0,0) = \mathbf{I}$. The system (10) has a unit Wronskian [12]

$$\det(\mathbf{\Phi}(t,0)) = e^{\int_0^t \operatorname{trace}(\mathbf{A}(t)) dt} = 1,$$

because $\operatorname{trace}(\mathbf{A}(t)) = 0$. The fundamental solution matrix $\mathbf{\Phi}(t,0)$ maps the initial condition $\mathbf{y}(0)$ to the state $\mathbf{y}(t)$

$$\mathbf{y}(t) = \mathbf{\Phi}(t,0)\mathbf{y}(0). \qquad (11)$$

More generally, the fundamental solution matrix $\mathbf{\Phi}(t_1, t_0)$ is defined as the mapping

$$\mathbf{y}(t_1) = \mathbf{\Phi}(t_1, t_0)\mathbf{y}(t_0),$$

which fulfills the matrix differential equation $\dot{\mathbf{\Phi}}(t, t_0) = \mathbf{A}(t)\mathbf{\Phi}(t, t_0)$ with $\mathbf{\Phi}(t_0, t_0) = \mathbf{I}$. Furthermore, the transition property

$$\mathbf{\Phi}(t_2, t_0) = \mathbf{\Phi}(t_2, t_1)\mathbf{\Phi}(t_1, t_0) \qquad (12)$$

holds from which one can deduce the inverse of the fundamental solution matrix to be $\mathbf{\Phi}(t_1, t_0) = \mathbf{\Phi}(t_0, t_1)^{-1}$. The fundamental solution matrix $\mathbf{\Phi}_T = \mathbf{\Phi}(\pi, 0)$ is referred to as the monodromy matrix. The trace of the monodromy matrix

$$\Delta = \operatorname{trace}(\mathbf{\Phi}_T) = y_1(\pi) + \dot{y}_2(\pi) \qquad (13)$$

is known as the discriminant of Hill's equation [17]. The monodromy matrix $\mathbf{\Phi}_T$ has the characteristic equation

$$\lambda^2 - \Delta\lambda + 1 = 0, \qquad (14)$$

because $\det(\mathbf{\Phi}_T) = 1$, and the eigenvalues

$$\lambda_{1,2} = \frac{1}{2} \Delta \pm \frac{1}{2} \sqrt{\Delta^2 - 4}, \qquad (15)$$

which are called characteristic multipliers (or Floquet multipliers). The characteristic multipliers are reciprocal in the sense that $\lambda_1 = 1/\lambda_2$. The discriminant $\Delta = \lambda_1 + \lambda_2$ plays an essential role in the analysis of Hill's equation and it is useful to distinguish between three cases:

* If it holds that $|\Delta| < 2$, then the characteristic multipliers $\lambda_1 = \bar{\lambda}_2$ are complex conjugated and located on the unit circle because $|\lambda_{1,2}| = 1$.

* If it holds that $|\Delta| > 2$, then the characteristic multipliers $\lambda_1$ and $\lambda_2$ are real and distinct and we define the order by $|\lambda_1| \ge |\lambda_2|$. There exist two linearly independent real eigenvectors $\mathbf{v}_1$ and $\mathbf{v}_2$. If $y_2(\pi) \ne 0$, then the eigenvectors are given by

$$\mathbf{v}_1 = \frac{\operatorname{sign}(y_2(\pi))}{\sqrt{y_2(\pi)^2 + (\lambda_1 - y_1(\pi))^2}} \quad \lambda_1 - y_1(\pi), \\ \mathbf{v}_2 = \frac{\operatorname{sign}(y_2(\pi))}{\sqrt{y_2(\pi)^2 + (\lambda_2 - y_1(\pi))^2}} \quad \lambda_2 - y_1(\pi). \qquad (16)$$

The normalization of the eigenvectors $\mathbf{v}_1$ and $\mathbf{v}_2$ is chosen such that $\mathbf{v}_{1,2}$ lie in the first or fourth quadrant and $\|\mathbf{v}_{1,2}\| = 1$. If $y_2(\pi) = 0$ then either $\lambda_1 - y_1(\pi) = 0$ or $\lambda_2 - y_1(\pi) = 0$ and the expression (16) for $\mathbf{v}_1$ or $\mathbf{v}_2$ degenerates. In this case there still exist two linearly independent real eigenvectors of which either $\mathbf{v}_1$ or $\mathbf{v}_2$ is equal to $(0, 1)^{\mathrm{T}}$.

* If it holds that $|\Delta| = 2$, then the characteristic multipliers are equal $\lambda_1 = \lambda_2 = \Delta/2 = \pm 1$.

Hill's equation has two (complex) eigensolutions

$$f_1(t) = e^{\sigma t} p_1(t), \quad f_2(t) = e^{-\sigma t} p_2(t), \qquad (17)$$

where $p_i(t) = p_i(t+\pi)$, $i=1,2$, are complex periodic functions and $\sigma$ is the characteristic exponent defined by $e^{\pm \pi \sigma} = \lambda_{1,2}$, i.e. $\Delta = \lambda_1 + \lambda_2 = 2 \cosh(\pi \sigma)$.

In Section 3 on the stability analysis of the unilaterally constrained Hill's equation, it will be of interest to know the number of zeros of the fundamental solution $y_2(t)$ on the half-open interval $(0, \pi]$. The floor function $\lfloor \cdot \rfloor$ and fractional part ${\cdot}$ will be used to count the zeros. The definition and properties of the floor function and fractional part are given in Appendix A.

**Proposition 2.** Let $n$ denote the number of zeros of $y_2(t)$ on the interval $(0, \pi]$. It holds that

$$n = \left\lfloor \frac{1}{\pi} \int_0^\pi \frac{dt}{y_1(t)^2 + y_2(t)^2} \right\rfloor. \qquad (18)$$

**Proof.** According to Proposition 1, one may write $y_2(t) = \varrho(t) \sin \psi(t)$ with $\varrho(t) > 0$ for all $t$. The zeros of $y_2(t)$ are therefore given by the condition $\sin \psi(t) = 0$, i.e. $\psi(t) = \pi, 2\pi, 3\pi, \dots$ and using (90) it therefore holds that

$$n = \left\lfloor \frac{\psi(\pi)}{\pi} \right\rfloor \qquad (19)$$

with $\psi(\pi) = \int_0^\pi \varrho(t)^{-2} dt$. $\square$

Consider a solution $\mathbf{y}(t)$ of (2) with initial condition $\mathbf{y}(0) = r_0 \cos \theta_0 \ge 0$ and $\dot{\mathbf{y}}(0) = r_0 \sin \theta_0 > 0$, where $r_0 > 0$ and $\theta_0 \in (-\pi/2, \pi/2]$. With $m(t) \in \mathbb{N}_0$ we will denote the right-continuous monotonically increasing step function which counts the number of reflections in the interval $(0, t]$ such that
($-1)^{m(t)}y(t) = |y(t)|$. It can easily be verified that $m(0) = 0$ and $m(t)$ equals the number of zeros of $y(t)$ on the interval $(0,t]$. In the following, we will be interested in the number $m(\pi)$ and use the short-hand notation $m = m(\pi)$. By definition, it holds that $m=n$ for $\theta_0 = \pi/2$, because $n$ is defined as the number of zeros of $y_2(t)$ on the interval $(0,\pi]$ and $y(t) = r_0y_2(t)$ for $\theta_0 = \pi/2$. The calculation of $m$ is given by the following proposition.

**Proposition 3.** Consider a solution $y(t)$ of (2) with initial condition $y(0) = r_0 \cos \theta_0 \ge 0$ and $\dot{y}(0) = r_0 \sin \theta_0$, where $r_0 > 0$ and $-\pi/2 < \theta_0 \le \pi/2$. Let $m$ denote the number of zeros of $y(t)$ on the interval $(0,\pi]$. It holds that

$$m = \begin{cases} n & \text{if } \theta_0 > \theta_c, \\ n+1 & \text{if } \theta_0 \le \theta_c, \end{cases} \qquad (20)$$

with

$$\theta_c = -\arctan\left(\frac{y_1(\pi)}{y_2(\pi)}\right). \qquad (21)$$

**Proof.** From (11) follows that the solution $y(t)$ is given by the linear combination $y(t) = y_1(t)y(0)+y_2(t)\dot{y}(0)$ which we write in polar coordinates using Proposition 1 as

$$y(t) = \varrho(t)r_0(\cos \theta_0 \cos \psi(t) + \sin \theta_0 \sin \psi(t)) = \varrho(t)r_0(\cos(\psi(t)-\theta_0))$$

with $\varrho(t) > 0$. The zeros of $y(t)$ are therefore determined by the condition $\cos(\psi(t)-\theta_0) = 0$ or equivalently $\sin(\psi(t)-\theta_0+\pi/2) = 0$. Using (91), the number of zeros of $y(t)$ therefore amounts to

$$m = \left[ \frac{\psi(\pi) - \theta_0 + \frac{\pi}{2}}{\pi} \right] - \left[ \frac{-\theta_0 + \frac{\pi}{2}}{\pi} \right], \quad -\frac{\pi}{2} < \theta_0 \le \frac{\pi}{2}. \qquad (22)$$

Property (88) gives together with (19) the inequality $n \le m \le n+1$, i.e. $m=n$ or $m=n+1$. From property (92) follows that $m=n$ if and only if

$$\cot(\psi(\pi)) + \cot\left(-\theta_0 + \frac{\pi}{2}\right) > 0 \qquad (23)$$

which we write as $\cot(\psi(\pi)) > -\tan(\theta_0)$, i.e.

$$\theta_0 > \theta_c = -\arctan(\cot(\psi(\pi)))$$

where $\cot(\psi(\pi)) = y_1(\pi)/y_2(\pi)$. If the inequality does not hold, then $m \ne n$ and $m$ must therefore be equal to $n+1$. $\square$

**Remark.** Some care needs to be taken for the case when $y_2(\pi) = 0$, i.e. $\psi(\pi) = n\pi$, $n \in \mathbb{N}_0$. In this case it holds that $\cot(\psi(\pi)) = +\infty$ and the inequality (23) is therefore fulfilled which implies that $m=n$. Moreover, if $\psi(\pi) = n\pi$, then it can immediately be verified from (22) that $m=n$.

### 3. The unilaterally constrained Hill's equation

In this section the unilaterally constrained Hill's equation is analyzed which consists of Hill's differential equation

$$\ddot{x}(t) + g(t)x(t) = 0, \qquad (24)$$

which holds for almost all $t$ and $x(t) \ge 0$, and the Newtonian impact law

$$\dot{x}^{+}(t_i) = -\varepsilon \dot{x}^{-}(t_i) \qquad (25)$$

at time-instant for which $x(t_i) = 0$. With the notation

$$\dot{x}^{-}(t_i) = \lim_{t_1 \uparrow t_i} \dot{x}(t), \quad \dot{x}^{+}(t_i) = \lim_{t_i \downarrow t_i} \dot{x}(t),$$

we denote the pre- and post-impact velocity and with $\varepsilon \in [0, 1]$ the restitution coefficient. The velocity $\dot{x}(t)$ of the unilaterally constrained Hill's equation is considered to be right-continuous by convention, i.e. $\dot{x}(t) = \dot{x}^{+}(t)$, and the initial condition $(x(0), \dot{x}(0))$ is

likewise considered to be a post-impact state. Hence, if the initial condition is written in polar coordinates

$$x(0) = r_0 \cos \theta_0, \quad \dot{x}(0) = r_0 \sin \theta_0,$$

then it necessarily holds that $r_0 \ge 0$ and $\theta_0 \in (-\pi/2, \pi/2]$. The solution $x(t)$ of the unilaterally constrained Hill's equation is therefore confined to the domain

$$\mathcal{D} = \{\mathbf{x} = (x, \dot{x})^{\mathrm{T}} \in \mathbb{R}^2 \mid x = r \cos \theta, \dot{x} = r \sin \theta, r \ge 0, -\frac{\pi}{2} < \theta \le \frac{\pi}{2}\}. \qquad (26)$$

The impact law (25) becomes active when $x(t_i) = 0$ and induces a jump $\dot{x}^{+}(t_i) = -\varepsilon \dot{x}^{-}(t_i)$ in the velocity whereas the position $x(t)$ remains continuous at the collision time-instant, i.e. $x^{+}(t_i) = x^{-}(t_i) = x(t_i) = 0$. The impact law can therefore formally be written for position and velocity as a homogeneous map

$$x^{+}(t_i) = -\varepsilon x^{-}(t_i), \quad \dot{x}^{+}(t_i) = -\varepsilon \dot{x}^{-}(t_i), \qquad (27)$$

because $x^{+}(t_i) = x^{-}(t_i) = 0$. The homogeneity of the impact conditions (27) is due to the fact that the unilateral constraint is located at $x=0$. A non-zero location of the unilateral constraint would give a completely different type of dynamics and is not studied in this paper but some remarks are given at the end of Section 6.

The homogeneity of the linear differential equation (24) and the homogeneity of the impact conditions (27) have important consequences for the solution $x(t)$ of the unilaterally constrained Hill's equation. Consider a solution curve $x(t)$ of the unilaterally constrained Hill's equation as well as the solution curve $y(t)$ of the unconstrained Hill's equation, see Fig. 2, both with the same initial condition $x(0) = y(0)$ and $\dot{x}(0) = \dot{y}(0)$ such that $(x(0), \dot{x}(0))^{\mathrm{T}} \in \mathcal{D}$. It holds that $x(t) = y(t)$ during the interval $t \in [0, t_1]$, where $t_1$ is the time-instant for which $x(t_1) = y(t_1) = 0$. At $t = t_1$ an impact occurs in the unilaterally constrained Hill's equation. The solution $x(t)$ after $t = t_1$ is reflected and scaled with $\varepsilon$, i.e.

$$x(t) = -\varepsilon y(t) \qquad (28)$$

for $t \in [t_1, t_2]$ due to the linearity and homogeneity of (24). The second impact occurs when $x(t)$ becomes again zero, which is the next zero of $y(t)$. The number of zeros of $y(t)$ in the interval $(0, t]$ is given by $m(t)$, see Proposition 3. The solution $x(t)$ is therefore directly related to $y(t)$ and $m(t)$ through

$$x(t) = (-\varepsilon)^{m(t)} y(t) \qquad (29)$$

and, using $m = m(\pi)$, we obtain the relationship

$$x(\pi) = (-\varepsilon)^m y(\pi). \qquad (30)$$

In non-smooth dynamical systems, the so-called 'Zeno-behaviour' can be present which is characterized by the occurrence of infinitely many non-smooth events in a finite time interval. For instance, in the bouncing ball system, the equilibrium may be reached through an infinite number of impacts in a finite time [14]. In non-smooth dynamics, this is referred to as an accumulation point of impact events. The occurrence of an accumulation point implies that uniqueness of the solution in backward time is lost [16].

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Solution $x(t)$ of the unilaterally constrained Hill's equation and $y(t)$ of Hill's equation.</figcaption></figure>
**Theorem 1.** Accumulation points of impact events do not occur in dynamics of the unilaterally constrained Hill's equation (4).

**Proof.** The presence of an accumulation point would imply that there exists an initial condition $(x(0), \dot{x}(0))$ such that $m$ is infinite. However, the angle $\psi(\pi) = \int_{0}^{\pi} q(t)^{-2} dt$ is finite because $q(t)$ is continuous and strictly positive for all $t$, see Proposition 1. The number $n = \lfloor \psi(\pi)/\pi \rfloor$ is therefore finite and the proof follows from $n \le m \le n+1$. $\square$

The absence of accumulation points of impact events implies that the solution of (4) is unique in forward and backward time. Moreover, if the equilibrium is attractive, then the attraction is asymptotic in the sense that attraction cannot occur in finite time [16].

The direct relationship between the solutions of the unilaterally constrained and unconstrained Hill's equation can be expressed in first-order form and be related to the fundamental solution matrix. Using $\mathbf{x}(t) = (x(t) \ \dot{x}(t))^{\mathrm{T}}$, the system is written in first-order form as

$$ \dot{\mathbf{x}}(t) = \mathbf{A}(t)\mathbf{x}(t) $$

$$ \mathbf{n}^{\mathrm{T}}\mathbf{x}^{-}(t_i) = 0: \quad \mathbf{x}^{+}(t_i) = \mathbf{G}\mathbf{x}^{-}(t_i) \qquad (31) $$

with the system matrix $\mathbf{A}(t)$ and impact map $\mathbf{G}$ given by

$$ \mathbf{A}(t) = \begin{pmatrix} 0 & 1 \\ -g(t) & 0 \end{pmatrix}, \quad \mathbf{G} = -\varepsilon\mathbf{I}, \quad \mathbf{n} = \begin{pmatrix} 1 \\ 0 \end{pmatrix}. \qquad (32) $$

The homogeneity of the impact conditions (27) allows us to write the impact map $\mathbf{G}$ as $-\varepsilon\mathbf{I}$ and the matrix $\mathbf{G}$ therefore commutes with any arbitrary matrix. The solution of the unilaterally constrained Hill's equation can be obtained by concatenation of non-impulsive motion given by Hill's differential equation and the impact equations. The non-impulsive motion between two consecutive collision time-instants $t_i$ and $t_{i+1}$ is described by the linear homogeneous differential equation $\dot{\mathbf{x}}(t) = \mathbf{A}(t)\mathbf{x}(t)$. The fundamental solution matrix $\mathbf{\Phi}(t_{i+1}, t_i)$ therefore maps the post-impact state $\mathbf{x}^{+}(t_i)$ to the pre-impact state $\mathbf{x}^{-}(t_{i+1})$,

$$ \mathbf{x}^{-}(t_{i+1}) = \mathbf{\Phi}(t_{i+1}, t_i)\mathbf{x}^{+}(t_i), \qquad (33) $$

and it necessarily holds that $\mathbf{n}^{\mathrm{T}}\mathbf{x}^{\pm}(t_i) = \mathbf{n}^{\mathrm{T}}\mathbf{x}^{\pm}(t_{i+1}) = 0$. The impulsive motion at each collision time-instant $t_i$ is described by the impact map $\mathbf{G}$. Concatenation of non-impulsive and impulsive motion gives

$$ \mathbf{x}(\pi) = \mathbf{\Phi}(\pi, t_m)\mathbf{G}\mathbf{\Phi}(t_m, t_{m-1})\mathbf{G}\dots\mathbf{\Phi}(t_2, t_1)\mathbf{G}\mathbf{\Phi}(t_1, 0)\mathbf{x}(0). \qquad (34) $$

The impact map $\mathbf{G} = -\varepsilon\mathbf{I}$ commutes with all the fundamental solution matrices $\mathbf{\Phi}(\cdot, \cdot)$ in (34) and Eq. (34) therefore simplifies to

$$ \begin{aligned} \mathbf{x}(\pi) &= \mathbf{G}^m \mathbf{\Phi}(\pi, t_m) \mathbf{\Phi}(t_m, t_{m-1}) \dots \mathbf{\Phi}(t_2, t_1) \mathbf{\Phi}(t_1, 0) \mathbf{x}(0) \\ &= \mathbf{G}^m \mathbf{\Phi}(\pi, 0) \mathbf{x}(0) \\ &= (-\varepsilon)^m \mathbf{\Phi}_T \mathbf{x}(0) \end{aligned} \qquad (35) $$

in which the transition property (12) of the fundamental solution matrix and the notation $\mathbf{\Phi}_T = \mathbf{\Phi}(\pi, 0)$ for the monodromy matrix has been used. Let $\mathbf{y}(t)$ denote the solution of Hill's equation with $\mathbf{y}(0) = \mathbf{x}(0)$. From (35) it becomes apparent that one can relate the solution $\mathbf{x}(\pi)$ of the unilaterally constrained Hill's equation to $\mathbf{y}(\pi)$ and $m$

$$ \mathbf{x}(\pi) = (-\varepsilon)^m \mathbf{\Phi}_T \mathbf{y}(0) = (-\varepsilon)^m \mathbf{y}(\pi), \qquad (36) $$

which is the same result as (30) but now in first-order form. The number $m$ depends on the angle $\theta_0$ of the initial condition $\mathbf{x}(0)$ and $m=n$ if $\theta_0 > \theta_c$ and $m=n+1$ else, see Proposition 3. The domain $\mathcal{D}$, defined in (26), is therefore split into two cones $\mathcal{D}_1$ and $\mathcal{D}_2$

$$ \mathcal{D}_1 = \{\mathbf{x} \in \mathcal{D} | \theta > \theta_c\}, $$

$$ \mathcal{D}_2 = \{\mathbf{x} \in \mathcal{D} | \theta \le \theta_c\}, \qquad (37) $$

by a half-hyperplane $\Sigma = \{\mathbf{x} \in \mathcal{D} | \theta = \theta_c\}$, where $\theta_c$ is defined by (21), such that $\mathcal{D} = \mathcal{D}_1 \cup \mathcal{D}_2$ and $\mathcal{D}_1 \cap \mathcal{D}_2 = \{0\}$, see Fig. 3. If $y_2(\pi) = 0$ then it holds that $\theta_c = -\pi/2$ and $\mathcal{D}_1 = \mathcal{D}$. The switching of the number $m$ on the half-hyperplane $\Sigma$ leads together with (36) to a piece-wise linear, or, more precisely, a cone-wise linear Poincaré map.

**Proposition 4 (Poincaré map).** Let $\mathbf{x}_j$ denote the solution $\mathbf{x}(\pi j)$, $j \in \mathbb{N}_0$, of the unilaterally constrained Hill's equation with $\mathbf{x}_0 = \mathbf{x}(0) \in \mathcal{D}$. The discrete time-evolution of $\mathbf{x}_j$ is given by the cone-wise linear Poincaré map

$$ \mathbf{x}_{j+1} = \begin{cases} \mathbf{A}_1 \mathbf{x}_j & \text{if } \mathbf{x}_j \in \mathcal{D}_1, \\ \mathbf{A}_2 \mathbf{x}_j & \text{if } \mathbf{x}_j \in \mathcal{D}_2, \end{cases} \qquad (38) $$

where $\mathbf{A}_1 = (-\varepsilon)^n \mathbf{\Phi}_T$ and $\mathbf{A}_2 = -\varepsilon \mathbf{A}_1 = (-\varepsilon)^{n+1} \mathbf{\Phi}_T$.

**Proof.** From (35) it is clear that $\mathbf{x}_{j+1} = (-\varepsilon)^m \mathbf{\Phi}(\pi, 0)\mathbf{x}_j$ where $m$ depends on $\mathbf{x}_j$. Write $\mathbf{x}_j$ in polar coordinates as $\mathbf{x}_j = (r_j \cos \theta_j \ r_j \sin \theta_j)^{\mathrm{T}}$. If $\mathbf{x}_j \in \mathcal{D}_1$, then it holds that $\theta_j > \theta_c$ which implies that $m=n$. Similarly, if $\mathbf{x}_j \in \mathcal{D}_2$, then it holds that $m=n+1$. $\square$

The cone-wise linearity of the Poincaré map introduced in Proposition 4 suggests to analyze the map by using polar coordinates $\mathbf{x}_j = (r_j \cos \theta_j \ r_j \sin \theta_j)^{\mathrm{T}}$, i.e.

$$ \tan \theta_j = \frac{\dot{x}(\pi j)}{x(\pi j)}, \quad r_j = \sqrt{x(\pi j)^2 + \dot{x}(\pi j)^2}. $$

Evaluation of the Poincaré map for $\varepsilon > 0$ yields

$$ \begin{aligned} \tan \theta_{j+1} &= \frac{\dot{x}(\pi(j+1))}{x(\pi(j+1))} = \frac{(-\varepsilon)^m (\dot{y}_1(\pi)\mathbf{x}(\pi j) + \dot{y}_2(\pi)\dot{x}(\pi j))}{(-\varepsilon)^m (y_1(\pi)\mathbf{x}(\pi j) + y_2(\pi)\dot{x}(\pi j))} \\ &= \frac{\dot{y}_1(\pi) + \dot{y}_2(\pi)\tan \theta_j}{y_1(\pi) + y_2(\pi)\tan \theta_j}, \end{aligned} \qquad (39) $$

where as $\theta_{j+1}$ is not defined if $r_{j+1} = 0$ for $\varepsilon = 0$. The map $\theta_{j+1} \mapsto \theta_{j+1}$ is therefore autonomous as it does not depend on $r_j$ which expresses the cone-wise linearity of the Poincaré map. The map $\theta_{j+1} \mapsto \theta_{j+1}$ can be simplified even further by using a non-linear transformation.

**Proposition 5.** Let $q_j = y_1(\pi) + y_2(\pi)\tan \theta_j$ and $y_2(\pi) \ne 0$. It holds that

$$ q_{j+1} = Q(q_j) = \Delta - \frac{1}{q_j}, \qquad (40) $$

where $\Delta = \mathrm{trace}(\mathbf{\Phi}_T) = y_1(\pi) + \dot{y}_2(\pi)$.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Domain $\mathcal{D}$ being in the cones $\mathcal{D}_1$ and $\mathcal{D}_2$ by $\Sigma$.</figcaption></figure>
**Proof.** Substitution of (39) in $q_j = y_1(\pi) + y_2(\pi) \tan \theta_j$ gives

$$q_{j+1} = y_1(\pi) + \frac{y_2(\pi)\dot{y}_1(\pi) + \dot{y}_2(\pi)(q_j - y_1(\pi))}{q_j} = \operatorname{trace}(\mathbf{\Phi}_T) - \frac{\operatorname{det}(\mathbf{\Phi}_T)}{q_j},$$

where

$$\operatorname{trace}(\mathbf{\Phi}_T) = y_1(\pi) + \dot{y}_2(\pi) \text{ and}$$

$$\det(\mathbf{\Phi}_T) = y_1(\pi)\dot{y}_2(\pi) - y_2(\pi)\dot{y}_1(\pi) = 1. \quad \square$$

The fixed points of the map $q_{j+1} = Q(q_j)$ are those values of $q$ which are mapped to themselves, i.e. $q^* = Q(q^*)$, and the fixed points therefore satisfy

$$q^{*2} - \Delta q^* + 1 = 0$$

in which we recognize the characteristic equation (14) of the monodromy matrix $\mathbf{\Phi}_T$. Hence, the fixed points $q^*$ agree with the real characteristic multipliers $\lambda_{1,2}$ of the monodromy matrix, see (15). The map $q_{j+1} = Q(q_j)$ is known as the Riccati difference equation (or first-order rational difference equation) and has been studied in detail in [9]. The dynamics can be considered for thee different cases [9]:

* $|\Delta| > 2$: The map $q_{j+1} = Q(q_j)$ has two distinct fixed points

$$q_{1,2}^* = \lambda_{1,2} = \frac{1}{2} \Delta \pm \frac{1}{2} \sqrt{\Delta^2 - 4}$$

and we define $|q_1^*| > |q_2^*|$ and equivalently $|\lambda_1| > |\lambda_2|$. Furthermore, because $\lambda_1\lambda_2 = 1$, it holds that $|q_1^*| > 1 > |q_2^*|$. The stability of the fixed points $q_{1,2}^*$ is determined by the derivative of the map $Q'(q) = 1/q^2$ from which we see that $q_1^*$ is asymptotically stable ($|Q'(q_1^*)| < 1$) whereas $q_2^*$ is unstable ($Q'(q_2^*) > 1$). The map $q_{j+1} = Q(q_j)$ with initial condition $q_0$ has a closed form solution given by

$$q_j = \frac{(\lambda_1^{j+1} - \lambda_2^{j+1})q_0 - (\lambda_1^j - \lambda_2^j)}{(\lambda_1^j - \lambda_2^j)q_0 - (\lambda_1^{j-1} - \lambda_2^{j-1})}$$

as we can simply verify or derive with the z-transformation [9]. If $q_0 \neq q_2^* = \lambda_1^{-1} = \lambda_2$, then the solution $q_j$ is attracted to $q_1^*$ which follows from

$$\lim_{j \to \infty} q_j = \lim_{j \to \infty} \frac{\lambda_1^{j+1} q_0 - \lambda_1^j q_0}{\lambda_1^j q_0 - \lambda_1^{j-1}} = \lambda_1 = q_1^*.$$

* $|\Delta| = 2$: The map has a single fixed point $q^* = \frac{1}{2}\Delta = \pm 1$. The map with initial condition $q_0$ has the closed form solution

$$q_j = \frac{q_0 \Delta(1+j) - 2j}{2q_0 j + \Delta(1-j)}$$

and limit

$$\lim_{j \to \infty} q_j = \frac{q_0 \Delta - 2}{2q_0 - \Delta} = q^*.$$

The fixed point $q^*$ is therefore unstable but globally attractive. Furthermore, if $q_j q^* < 0$ then it holds that $q_{j+1} q^* > 0$.

* $|\Delta| < 2$: The map has no fixed points and the solution is quasi-periodic, wandering between $\mathbb{R}_0^-$ and $\mathbb{R}^+$ (see [9]). The number of iterations that the solution $q_j$ remains in $\mathbb{R}_0^-$ or $\mathbb{R}^+$ is bounded.

The above results on the map $q_{j+1} = Q(q_j)$ for $|\Delta| > 2$ can be understood by noting that the values of $q_{1,2}^*$ define for $y_2(\pi) \neq 0$ the angles $\theta_1^*$ and $\theta_2^*$ of the eigenvectors $\mathbf{v}_1$ and $\mathbf{v}_2$ (see (16)) of the monodromy matrix $\mathbf{\Phi}_T$:

$$\tan \theta_1^* = \frac{q_1^* - y_1(\pi)}{y_2(\pi)} = \tan \theta_c + \frac{\lambda_1}{y_2(\pi)},$$

$$\tan \theta_2^* = \frac{q_2^* - y_1(\pi)}{y_2(\pi)} = \tan \theta_c + \frac{\lambda_2}{y_2(\pi)}. \quad (41)$$

If $\mathbf{x}_j$ is chosen in the direction of an eigenvector, say $\mathbf{x}_j = r_j \mathbf{v}_1$, then the next iteration point $\mathbf{x}_{j+1}$ will be again in the direction of that eigenvector, i.e. $\mathbf{x}_{j+1} = r_{j+1} \mathbf{v}_1$. In other words, if $\theta_j = \theta_1^*$ then also $\theta_{j+1} = \theta_1^*$ which implies that $\theta_1^*$ is a fixed point of the map $\theta_{j+1} \mapsto \theta_{j+1}$, or, equivalently, that $q_1^*$ is a fixed point of the map $q_{j+1} = Q(q_j)$.

The fixed points $q_{1,2}^* = \lambda_{1,2}$ are either both positive if $\Delta > 2$ or both negative if $\Delta < 2$ and the sign of $y_2(\pi)$ is given by $(-1)^n$. Hence, if $|\Delta| > 2$ and $y_2(\pi) \neq 0$, then it holds that $\mathbf{v}_{1,2} \in D_1$ for $(-1)^n \Delta > 2$ and $\mathbf{v}_{1,2} \in D_2$ for $(-1)^n \Delta < 2$. If $|\Delta| > 2$ and $y_2(\pi) = 0$ then it holds that $y_1(\pi)\dot{y}_2(\pi) = 1$ and $|y_1(\pi) + \dot{y}_2(\pi)| > 2$ from which we deduce that $\dot{y}_2(\pi) \neq 0$ and $\operatorname{sign}(\dot{y}_2(\pi)) = \operatorname{sign}(\Delta)$. If $y_2(\pi) = 0$ then it holds that $n > 0$ and $n$ is even for $\dot{y}_2(\pi) > 0$ and $n$ is odd for $\dot{y}_2(\pi) < 0$. Consequently, the condition $|\Delta| > 2$ and $y_2(\pi) = 0$ implies that $(-1)^n \Delta > 2$. Moreover, it holds that $D_1 = D$ for $y_2(\pi) = 0$. We conclude that, if $|\Delta| > 2$, then the location of the eigenvectors is determined by the condition

$$\mathbf{v}_{1,2} \in \begin{cases} \mathcal{D}_1 & (-1)^n \Delta > 2, \\ \mathcal{D}_2 & (-1)^n \Delta < -2. \end{cases} \quad (42)$$

for all values of $y_2(\pi)$.

The attractiveness of the fixed point $q_1^*$ implies that the solution $\mathbf{x}_j$ will be drawn towards the eigenvector $\mathbf{v}_1$ when $j \to \infty$, because $\mathbf{v}_1$ belongs to the characteristic multiplier which is largest in magnitude. The long-term behaviour of the dynamics is therefore determined by the matrix $\mathbf{A}_1 = (-\varepsilon)^n \mathbf{\Phi}_T$ if $\mathbf{v}_1 \in D_1$ or by the matrix $\mathbf{A}_2 = (-\varepsilon)^{n+1} \mathbf{\Phi}_T$ if $\mathbf{v}_1 \in D_2$. This insight leads directly to a stability condition. In the following, a more rigorous stability proof is given by using a discrete quadratic Lyapunov function. First a number of propositions are presented which discuss the various cases of the discriminant $\Delta$. Proposition 6 discusses the real case $\mathbf{v}_{1,2} \in D_1$ and Proposition 7 the real case $\mathbf{v}_{1,2} \in D_2$ using a slightly different Lyapunov function. Proposition 8 deals with the complex conjugated case $|\Delta| < 2$ whereas the cases $\Delta = \pm 2$ are discussed in Propositions 9–11.

In the following propositions, the spectral decomposition

$$\mathbf{\Phi}_T = \mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1} \quad (43)$$

of the monodromy matrix $\mathbf{\Phi}_T$ is used, where $\mathbf{V} = (\mathbf{v}_1 \ \mathbf{v}_2)$ is the matrix of eigenvectors and $\mathbf{\Lambda} = \operatorname{diag}(\lambda_i)$ is the diagonal matrix of eigenvalues of $\mathbf{\Phi}_T$. Such a spectral decomposition can be made if the characteristic multipliers $\lambda_{1,2}$ are distinct, implying $|\Delta| \neq 2$. Let $\mathbf{c}_j = (c_{1j} \ c_{2j})^\mathrm{T}$ be the coordinates on the basis of $\mathbf{v}_1$ and $\mathbf{v}_2$ such that $\mathbf{x}_j = c_{1j} \mathbf{v}_1 + c_{2j} \mathbf{v}_2$, i.e.

$$\mathbf{c}_j = \mathbf{V}^{-1} \mathbf{x}_j. \quad (44)$$

**Proposition 6** (Stability for the real case $\mathbf{v}_{1,2} \in D_1$). If $\mathbf{v}_{1,2} \in D_1$ then the equilibrium $\mathbf{x}^* = \mathbf{0}$ of the unilaterally constrained Hill's equation is globally uniformly asymptotically stable if $\varepsilon^n |\lambda_1| < 1$ and unstable if $\varepsilon^n |\lambda_1| > 1$.

**Proof.** If $\mathbf{v}_{1,2} \in D_1$ then the characteristic multipliers are real and distinct and the decomposition (43) exists. Consider the discrete Lyapunov candidate function

$$V(\mathbf{x}_j) = V_j = \mathbf{x}_j^\mathrm{T} \mathbf{P} \mathbf{x}_j \quad (45)$$

with $\mathbf{P} = \mathbf{V}^{-\mathrm{T}} \mathbf{V}^{-1}$ being dependent on the matrix of eigenvectors $\mathbf{V}$. The positive definiteness of $\mathbf{P} = \mathbf{V}^\mathrm{T} > 0$ implies that $\mathbf{V}$ is a positive definite function. It holds that $\mathbf{x}_{j+1} = (-\varepsilon)^m \mathbf{\Phi}_T \mathbf{x}_j$ and, using (44), the increment of $\mathbf{V}$ therefore gives

$$\begin{aligned} V_{j+1} - V_j &= \mathbf{x}_j^\mathrm{T} (\varepsilon^{2m} \mathbf{\Phi}_1^\mathrm{T} \mathbf{P} \mathbf{\Phi}_T - \mathbf{P}) \mathbf{x}_j \\ &= \mathbf{x}_j^\mathrm{T} \mathbf{V}^{-\mathrm{T}} (\varepsilon^{2m} \mathbf{\Lambda}^2 - \mathbf{I}) \mathbf{V}^{-1} \mathbf{x}_j \\ &= \mathbf{c}_j^\mathrm{T} (\varepsilon^{2m} \mathbf{\Lambda}^2 - \mathbf{I}) \mathbf{c}_j \\ &= c_{1j}^2 (\varepsilon^{2m} \lambda_1^2 - 1) + c_{2j}^2 (\varepsilon^{2m} \lambda_2^2 - 1). \end{aligned} \quad (46)$$

If $\varepsilon^n |\lambda_1| < 1$, then it holds that $\varepsilon^{2m} \lambda_{1,2}^2 < 1$ for $m=n$ and $m=n+1$ which implies that $V_{j+1} < V_j$ for $\mathbf{x}_j \neq \mathbf{0}$. Hence, $\mathbf{V}$ is a quadratic
time-independent Lyapunov function. The origin is therefore globally uniformly asymptotically stable if $\epsilon^n|\lambda_1| < 1$.

We now prove the last part of the proposition which states that if $\mathbf{v}_{1,2} \in \mathcal{D}_1$ and $\epsilon^n|\lambda_1| > 1$, then the origin is unstable. A solution $\mathbf{x}_j$ which starts in the direction of the eigenvector $\mathbf{v}_1 \in \mathcal{D}_1$, i.e. $\mathbf{x}_0 = c_{10}\mathbf{v}_1$, will remain along $\mathbf{v}_1$ and therefore in $\mathcal{D}_1$ for all $j \ge 0$. It therefore holds that $m=n$ during all iterations of the Poincaré map. The solution $\mathbf{x}_j = c_{1j}\mathbf{v}_1$ is therefore given by

$$c_{1j} = (-\epsilon)^n \lambda_1 c_{1,j-1} = ((-\epsilon)^n \lambda_1)^j c_{10}$$

which grows unbounded for $\epsilon^n|\lambda_1| > 1$. The value of $c_{10}$ can be chosen arbitrarily small and the origin is therefore unstable. $\square$

**Proposition 7 (Stability for the real case** $\mathbf{v}_{1,2} \in \mathcal{D}_2$). If $\mathbf{v}_{1,2} \in \mathcal{D}_2$ then the equilibrium $\mathbf{x}^* = \mathbf{0}$ of the unilaterally constrained Hill's equation is globally uniformly asymptotically stable if $\epsilon^{n+1}|\lambda_1| < 1$ and unstable if $\epsilon^{n+1}|\lambda_1| > 1$.

**Proof.** If $\mathbf{v}_{1,2} \in \mathcal{D}_2$ then the eigenvalues are real and distinct such that $\lambda_1^2 > 1$ and $\lambda_2^2 = \lambda_1^{-2} < 1$. Moreover, if $\mathbf{v}_{1,2} \in \mathcal{D}_2$ then it holds that $y_2(\pi) \ne 0$ and therefore $\cos \theta_{1,2}^* > 0$. Consider the discrete Lyapunov candidate function

$$V(\mathbf{x}_j) = V_j = \mathbf{x}_j^T \mathbf{P} \mathbf{x}_j \qquad (47)$$

with $\mathbf{P} = \mathbf{V}^{-T} \mathbf{B} \mathbf{V}^{-1}$ and

$$\mathbf{B} = \begin{pmatrix} \cos^2 \theta_1^* & 0 \\ 0 & \lambda_1^4 \cos^2 \theta_2^* \end{pmatrix}, \qquad (48)$$

where $\theta_1^*$ and $\theta_2^*$ are the angles of $\mathbf{v}_1$ and $\mathbf{v}_2$ defined by (41). The positive definiteness of V is assured through $\mathbf{P} = \mathbf{P}^T > 0$. It holds that $\mathbf{x}_{j+1} = (-\epsilon)^m \mathbf{\Phi}_T \mathbf{x}_j$ and the increment of V yields

$$\begin{aligned} V_{j+1} - V_j &= \mathbf{x}_j^T (\epsilon^{2m} \mathbf{\Phi}_T^T \mathbf{P} \mathbf{\Phi}_T - \mathbf{P}) \mathbf{x}_j \\ &= \mathbf{x}_j^T \mathbf{V}^{-T} (\epsilon^{2m} \mathbf{\Lambda} \mathbf{B} \mathbf{\Lambda}^{-1} \mathbf{B}) \mathbf{V}^{-1} \mathbf{x}_j \\ &= \mathbf{c}_j^T (\epsilon^{2m} \mathbf{\Lambda}^2 - \mathbf{I}) \mathbf{B} \mathbf{c}_j \\ &= c_{1j}^2 \cos^2 \theta_1^* (\epsilon^{2m} \lambda_1^{-2} - 1) + c_{2j}^2 \cos^2 \theta_2^* \lambda_1^4 (\epsilon^{2m} \lambda_2^2 - 1). \end{aligned} \qquad (49)$$

The decrease in the Lyapunov function V depends on the location of $\mathbf{x}_j$:

* If $\mathbf{x}_j \in \mathcal{D}_2 \setminus \{0\}$ then $m=n+1$ and $\epsilon^{2m}\lambda_1^2 < 1$ which implies $V_{j+1} - V_j < 0$.

* If $\mathbf{x}_j \in \mathcal{D}_1 \setminus \{0\}$ then $m=n \ge 0$ and it holds that $\theta_c < \theta_j \le \pi/2$. Hence, it holds that $\theta_1^* < \theta_2^* < \theta_j$, i.e. $\mathbf{x}_j$ lies inside the cone spanned by $-\mathbf{v}_1$ and $\mathbf{v}_2$, which implies that $c_{1j} < 0$ and $c_{2j} > 0$. Furthermore, from the positiveness of

$$x_{1j} = c_{1j} \cos \theta_1^* + c_{2j} \cos \theta_2^* > 0$$

together with $c_{1j}$ we conclude that $|c_{1j}/c_{2j}| \le |\cos \theta_2^*/\cos \theta_1^*|$. We now define $d_1 = c_{1j}^2 \cos^2 \theta_1^* + c_{2j}^2 \cos^2 \theta_2^*$ and $d_2 = -c_{1j}^2 \cos^2 \theta_1^* + c_{2j}^2 \cos^2 \theta_2^*$. It holds that $d_1 > 0$ and $d_2 \ge 0$ for $\mathbf{x}_j \in \mathcal{D}_1 \setminus \{0\}$.

Using $\lambda_1 = \lambda_2^{-1} > 1$ we obtain:

$$\begin{aligned} V_{j+1} - V_j &= \frac{1}{2}(d_1 - d_2)(\epsilon^{2m} \lambda_1^2 - 1) + \frac{1}{2}(d_1 + d_2)(\epsilon^{2m} \lambda_1^2 - \lambda_1^4) \\ &= \frac{1}{2}d_1(2\epsilon^{2m} \lambda_1^2 - 1 - \lambda_1^4) + \frac{1}{2}d_2(1 - \lambda_1^4) \\ &= \frac{1}{2}d_1(2(\epsilon^{2m} - 1)\lambda_1^2 - (\lambda_1^2 - 1)^2) + \frac{1}{2}d_2(1 - \lambda_1^4). \end{aligned} \qquad (50)$$

The expressions between the brackets are strictly negative for all $m \ge 0$ which implies that $V_{j+1} - V_j < 0$.

Hence, the quadratic positive definite function V is strictly decreasing for $\mathbf{x}_j \ne \mathbf{0}$. The origin is therefore globally uniformly asymptotically stable if $\epsilon^{n+1}|\lambda_1| < 1$. Instability for $\epsilon^{n+1}|\lambda_1| > 1$

can be proven by considering an initial condition in the direction of $\mathbf{v}_1$ (see the proof of Proposition 6). $\square$

**Proposition 8 (Stability for the complex conjugated case).** If $|\Delta| < 2$ and $\epsilon < 1$ then the equilibrium $\mathbf{x}^* = \mathbf{0}$ of the unilaterally constrained Hill's equation is globally uniformly asymptotically stable.

**Proof.** If $|\Delta| < 2$ then the eigenvalues $\lambda_{1,2}$ are complex conjugated (and distinct) with magnitude $|\lambda_{1,2}| = 1$ and $\mathbf{v}_1 = \bar{\mathbf{v}}_2$. The coordinates $c_{1j}$ and $c_{2j}$ are in this case complex conjugated, i.e. $c_{1j} = \bar{c}_{2j}$ and $|c_{1j}| = |c_{2j}|$, because $\mathbf{x}_j$ is real. Consider the discrete Lyapunov candidate function

$$V(\mathbf{x}_j) = V_j = \mathbf{x}_j^T \mathbf{P} \mathbf{x}_j \qquad (51)$$

with $\mathbf{P} = \bar{\mathbf{V}}^{-T} \mathbf{V}^{-1} = (\mathbf{V}^{-1})^* \mathbf{V}^{-1}$ being a hermitian matrix. Clearly, it holds that $\mathbf{P} = \mathbf{P}^* > 0$ and V is positive definite. Using $\mathbf{x}_{j+1} = (-\epsilon)^m \mathbf{\Phi}_T \mathbf{x}_j$ the increment of V yields

$$\begin{aligned} V_{j+1} - V_j &= \mathbf{x}_j^T (\epsilon^{2m} \mathbf{\Phi}_T^T \mathbf{P} \mathbf{\Phi}_T - \mathbf{P}) \mathbf{x}_j \\ &= \mathbf{c}_j^T (\epsilon^{2m} \mathbf{\Lambda}^* \mathbf{\Lambda}^{-1} \mathbf{I}) \mathbf{c}_j \\ &= |c_{1j}|^2 (\epsilon^{2m} |\lambda_1|^{-2} - 1) + |c_{2j}|^2 (\epsilon^{2m} |\lambda_2|^{-2} - 1) \\ &= (\epsilon^{2m} - 1) \|\mathbf{c}_j\|^2. \end{aligned} \qquad (52)$$

The increment is non-positive, $V_{j+1} - V_j \le 0$, because $m \ge 0$ which proves that the equilibrium is uniformly stable. It holds that $m=n$ for $\mathbf{x}_j \in \mathcal{D}_1$ and $m=n+1$ for $\mathbf{x}_j \in \mathcal{D}_2$. If $n=0$ then it follows that $V_{j+1} = V_j$ for $\mathbf{x}_j \in \mathcal{D}_1 \setminus \{0\}$ and $V_{j+1} < V_j$ for $\mathbf{x}_j \in \mathcal{D}_2 \setminus \{0\}$. However, the dynamics of the map $q_{j+1} = Q(q_j)$ proves that the cones $\mathcal{D}_1$ and $\mathcal{D}_2$ do not have a positively invariant subset (or sub-cone) other than the origin if $|\Delta| < 2$, i.e. no limit set exists which lies exclusively in $\mathcal{D}_1 \setminus \{0\}$ or exclusively in $\mathcal{D}_2 \setminus \{0\}$. This implies that the solution $\mathbf{x}_j$ will wander between $\mathcal{D}_1$ and $\mathcal{D}_2$ and the Lyapunov function strictly decreases whenever $\mathbf{x}_j \in \mathcal{D}_2 \setminus \{0\}$. We can therefore invoke LaSalle's invariance principle for discrete-time systems [13]. The vanishing of the increment $V_{j+1} - V_j = 0$ holds in $\mathcal{D}_1$ and the largest positively invariant set in $\mathcal{D}_1$ is the origin. Hence, it holds that $\lim_{j \to \infty} V_j = 0$, which proves that the equilibrium is globally uniformly attractive. $\square$

**Proposition 9.** If $\Delta = -2$, $n=0$ and $\epsilon < 1$ then the equilibrium $\mathbf{x}^* = \mathbf{0}$ of the unilaterally constrained Hill's equation is globally uniformly asymptotically stable.

**Proof.** The condition $n=0$ implies that $y_2(\pi) > 0$. If $\Delta = -2$ then the eigenvalues of $\mathbf{\Phi}_T$ are $\lambda_{1,2} = -1$. The matrix $\mathbf{A}_1 = (-\epsilon)^n \mathbf{\Phi}_T$ is not a stable matrix for $\Delta = -2$ and $n=0$. The matrix $\mathbf{A}_2 = -\epsilon \mathbf{A}_1$ is a stable matrix for $\epsilon < 1$. If $\Delta = -2$ then the map $q_{j+1} = Q(q_j)$ has a unique fixed point $q^* = -1$ which is globally attractive. Furthermore, if $q_0 > 0$ then $q_j < -1$ for all $j > 0$. The cone $\mathcal{D}_2$ has therefore a globally attractive sub-cone. The long-term dynamics is therefore governed by the stable matrix $\mathbf{A}_2$. $\square$

**Proposition 10.** If $\Delta = 2$ and $n=0$ then the equilibrium $\mathbf{x}^* = \mathbf{0}$ of the unilaterally constrained Hill's equation is not attractive and unstable for all $\epsilon \ge 0$.

**Proof.** The condition $n=0$ implies that $y_2(\pi) > 0$. If $\Delta = 2$ then the eigenvalues of $\mathbf{\Phi}_T$ are $\lambda_{1,2} = 1$ and the map $q_{j+1} = Q(q_j)$ has a unique fixed point $q^* = 1$ which is globally attractive. Furthermore, if $q_0 < 0$ then $q_j > 1$ for all $j > 0$. The cone $\mathcal{D}_1$ has therefore a globally attractive sub-cone. The long-term dynamics is therefore governed by the matrix $\mathbf{A}_1 = \mathbf{\Phi}_T$ which is not a stable matrix. The initial condition $\mathbf{x}_0 = (r_0 \cos \theta_0 \ r_0 \sin \theta_0)^T$ with tan $\theta_0 = (1-y_1(\pi))/y_2(\pi)$ will lead to a solution $\mathbf{x}_j = \mathbf{x}_0$ for all $r_0 \ge 0$. Hence the equilibrium is not attractive. Moreover, $\mathbf{A}_1$ is non-diagonalizable and $\mathbf{A}_1^j$ will diverge for $j \to \infty$ which implies unboundedness of solutions. $\square$
**Proposition 11.** If $|\Delta| \le 2$, $n > 0$ and $\varepsilon < 1$ then the equilibrium $\mathbf{x}^* = \mathbf{0}$ of the unilaterally constrained Hill's equation is globally uniformly asymptotically stable.

**Proof.** The spectral radius of $\Phi_T$ equals unity for $|\Delta| \le 2$. Hence, if in addition $n \ge 1$, then the spectral radii of $\mathbf{A}_1$ and $\mathbf{A}_2$ are strictly smaller than unity, i.e. the matrices $\mathbf{A}_1$ and $\mathbf{A}_2$ are stable matrices. The stability of $\mathbf{A}_1$ implies that there exist symmetric positive definite matrices $\mathbf{P}$ and $\mathbf{Q}_1$ such that

$$ \mathbf{A}_1^T \mathbf{P} \mathbf{A}_1 - \mathbf{P} = - \mathbf{Q}_1. $$

Using $\mathbf{P}$ we define the matrix $\mathbf{Q}_2$ such that

$$ \mathbf{A}_2^T \mathbf{P} \mathbf{A}_2 - \mathbf{P} = - \mathbf{Q}_2. $$

and express $\mathbf{Q}_2$ using $\mathbf{A}_2 = -\varepsilon \mathbf{A}_1$ and $\mathbf{Q}_1$ as

$$ \mathbf{Q}_2 = -\mathbf{A}_2^T \mathbf{P} \mathbf{A}_2 + \mathbf{P} = -\varepsilon^2 \mathbf{A}_1^T \mathbf{P} \mathbf{A}_1 + \mathbf{P} = (1-\varepsilon^2) \mathbf{P} + \mathbf{Q}_1. $$

With $\mathbf{P} > 0$, $\mathbf{Q}_1 > 0$ and $\varepsilon < 1$ we infer that $\mathbf{Q}_2 > 0$. Hence, there exists a common quadratic Lyapunov function $V(\mathbf{x}_j) = \mathbf{x}_j^T \mathbf{P} \mathbf{x}_j$ which proves that the origin is globally uniformly asymptotically stable. $\square$

The previous propositions are summarized in the following theorem, being the main result of the paper.

**Theorem 2.** Let $\Delta$ be the discriminant and $\lambda_1$ be the characteristic multiplier with largest magnitude of the unconstrained Hill's equation (2) and let $n$ denote the number of zeros of the fundamental solution $y_2(t)$ on the interval $(0, \pi]$. The equilibrium $\mathbf{x}^* = \mathbf{0}$ of the unilaterally constrained Hill's equation (4) is globally uniformly asymptotically stable if $0 \le \varepsilon < \varepsilon_c$, where the critical restitution coefficient is given by

$$ \varepsilon_c = \begin{cases} 0 & \text{if } n=0 \text{ and } \Delta \ge 2 \\ |\lambda_1|^{-1/n} & \text{if } n > 0 \text{ and } (-1)^n \Delta > 2, \\ |\lambda_1|^{-1/(n+1)} & \text{if } n \ge 0 \text{ and } (-1)^n \Delta < -2, \\ 1 & \text{if } n > 0 \text{ and } -2 \le \Delta \le 2 \text{ or if } n=0 \text{ and } -2 \le \Delta < 2. \end{cases} \quad (53) $$

If $|\Delta| > 2$ and $\varepsilon > \varepsilon_c$, then the equilibrium $\mathbf{x}^* = \mathbf{0}$ is unstable.

**Proof.** Propositions 6 and 10 prove that, if $n=0$ and $\Delta \ge 2$, then the equilibrium is not attractive for all $\varepsilon \ge 0$. The equilibrium is therefore in this case not globally uniformly asymptotically stable which is expressed by $\varepsilon_c = 0$. The proof of $\varepsilon_c = |\lambda_1|^{-1/n}$ for $n > 0$ and $(-1)^n \Delta > 2$ follows from (42) and Proposition 6. The proof of $\varepsilon_c = |\lambda_1|^{-1/(n+1)}$ for $n \ge 0$ and $(-1)^n \Delta < -2$ follows from (42) and Proposition 7. The proof of $\varepsilon_c = 1$ for $n > 0$ and $-2 \le \Delta \le 2$ follows from Proposition 11 and for $n=0$ and $-2 \le \Delta < 2$ from Propositions 8 and 9. $\square$

For $|\Delta| > 2$ the results of Theorem 2 can be put in a more tangible form by introducing the number $m_1$ as

$$ m_1 = \begin{cases} n & \text{if } (-1)^n \Delta > 2, \\ n+1 & \text{if } (-1)^n \Delta < -2. \end{cases} \quad (54) $$

The number $m_1$ is the value of $m$, i.e. the number of zeros of $y(t)$ on the interval $(0, \pi]$, of a solution curve $\mathbf{y}(t)$ which starts in the direction of the first eigenvector, i.e. $\mathbf{y}(0) = \mathbf{v}_1$.

**Corollary 1.** If $|\Delta| > 2$, then it holds that $\varepsilon_c = |\lambda_1|^{-1/m_1}$.

## 4. Approximation of the critical restitution coefficient using Hill's determinant

The results of the previous section, which are summarized in Theorem 2, show that the stability properties of the unilaterally

constrained Hill's equation are completely determined by the properties of the fundamental solutions of the unconstrained Hill's equation. More precisely, the critical restitution coefficient only depends on the value of the discriminant $\Delta$ and the number $n$. This insight suggests that standard approximation techniques for the unconstrained Hill's equation can be used to approximate the critical restitution coefficient of the unilaterally constrained Hill's equation. In this section the method of Hill's infinite determinant will be explored.

The function $g(t)$ in Hill's equation can be represented by a complex Fourier series as

$$ g(t) = \sum_{k=-\infty}^{\infty} g_k e^{2ikt}, \quad (55) $$

where $g_{-k} = \overline{g_k}$. In this section we will assume that $\sqrt{g_0} \ne 0, 2, 4, 8, \dots$, i.e. the even parametric resonances are avoided. Similarly, the first eigensolution $f_1(t) = e^{\sigma t} p_1(t)$, see (17), can be written as a complex Fourier series

$$ f_1(t) = e^{\sigma t} \sum_{k=-\infty}^{\infty} c_k e^{2ikt}, \quad (56) $$

where the characteristic exponent $\sigma$ is related to the discriminant through $\Delta = 2 \cosh(\pi\sigma)$. Substitution of the Fourier representations (55) and (56) in Hill's equation (2) yields the condition

$$ e^{\sigma t} \left[ \sum_{k=-\infty}^{\infty} (\sigma + 2ik)^2 + \sum_{s=-\infty}^{\infty} g_s e^{2ist} \right] c_k e^{2ikt} = 0. $$

Reordering terms gives the equality

$$ \sum_{k=-\infty}^{\infty} \left( (\sigma + 2ik)^2 c_k + \sum_{s=-\infty}^{\infty} g_s c_{k-s} \right) e^{2ikt} = 0. \quad (57) $$

The requirement that (57) has to be fulfilled for all $t$ leads to an infinite set of homogeneous equations

$$ (\sigma + 2ik)^2 c_k + \sum_{s=-\infty}^{\infty} g_s c_{k-s} = 0 \quad (58) $$

for the Fourier coefficients $c_k$. The linear homogeneous system of equations has only a non-trivial solution if the infinite determinant

$$ D(\sigma) = |H_{kl}|_{k,l = -\infty, \dots, \infty} \quad (59) $$

vanishes, where

$$ H_{kk} = \frac{g_0 + (\sigma + 2ik)^2}{g_0 - 4k^2}, \quad H_{kl} = \frac{g_{l-k}}{g_0 - 4k^2} \quad (l \ne k). \quad (60) $$

Each row in (58) has been divided by $g_0 - 4k^2$ to ensure convergence [18,29]. In [17,29] it is shown that the infinite determinant $D(\sigma)$ can be expressed as

$$ D(\sigma) = D(0) - \frac{\sin^2\left(\frac{1}{2}i\pi\sigma\right)}{\sin^2\left(\frac{1}{2}\pi\sqrt{g_0}\right)}. $$

Using the identity

$$ \sin^2\left(\frac{1}{2}i\pi\sigma\right) = \frac{1}{2}(1 - \cos(\pi\sigma i)) = \frac{1}{2}(1 - \cosh(\pi\sigma)) = \frac{1}{4}(2 - \Delta), $$

the determinant condition is written as

$$ D(\sigma) = D(0) - \frac{2 - \Delta}{4 \sin^2\left(\frac{1}{2}\pi \sqrt{g_0}\right)} = 0, $$
which agrees with [17], Theorem 2.9. The discriminant $\Delta$ can therefore be expressed as

$$ \Delta = 2 + 4D(0) \sin^2\left(\frac{1}{2}\pi\sqrt{g_0}\right). \qquad (61) $$

Furthermore, it holds that $|\lambda_1| = \frac{1}{2}|\Delta| + \frac{1}{2}\sqrt{\Delta^2 - 4}$, see (15).

If $|\Delta| > 2$, then the critical restitution coefficient is given by $\varepsilon_c = |\lambda_1|^{-1/m_1}$, see Corollary 1. The critical restitution coefficient can therefore be calculated from

$$ \varepsilon_c = \left( \left| 1 + 2D(0)\sin^2\left(\frac{1}{2}\pi\sqrt{g_0}\right) \right| + \sqrt{\left(1 + 2D(0)\sin^2\left(\frac{1}{2}\pi\sqrt{g_0}\right)\right)^2 - 1} \right)^{-1/m_1} \qquad (62) $$

under the assumption that $\sqrt{g_0} \neq 0, 2, 4, 8, \dots$. The value of $D(0)$ can be approximated by the determinant of a central $k \times k$ block as will be shown in Section 5.3.

# 5. The unilaterally constrained Mathieu equation

In order to illustrate the previous results, the stability properties of the unilaterally constrained Mathieu equation are studied in this section. First, the symmetry of the Mathieu equation is exploited in Section 5.1. A stability diagram for the unilaterally constrained Mathieu equation is obtained by using direct numerical integration to calculate the discriminant $\Delta$ and $n$. Subsequently, the method of averaging is employed in Section 5.2 to derive approximate expressions for the critical restitution coefficient in the vicinity of the first parametric resonance. Finally, the method of Hill's determinant is used in Section 5.2 to obtain an improved approximation for the critical restitution coefficient in this parameter region.

## 5.1. Stability boundaries

The Mathieu equation (3)

$$ \ddot{y}(t) + (\alpha + 2\beta \cos 2t)\dot{y}(t) = 0 $$

is a Hill's equation with symmetry of the function $g(t) = \alpha + 2\beta \cos 2t$, i.e.

$$ g(t) = g(-t), \quad g(t) = g(t + \pi). $$

This implies that, if $y_1(t)$ and $y_2(t)$ are solutions of (3), then $y_1(-t)$ and $y_2(-t)$ are also solutions of (3). The function $y_1(t)$ is therefore even and $y_2(t)$ is odd. For the same reason it holds that $\dot{y}_1(t)$ is odd and $\dot{y}_2(t)$ is even. Using the transition property (12) we deduce that $\Phi(-\pi, 0) = \Phi(0, \pi) = \Phi(\pi, 0)^{-1}$. Evaluation of

$$ \mathbf{\Phi}(-\pi, 0) = \begin{pmatrix} y_1(\pi) & -y_2(\pi) \\ -\dot{y}_1(\pi) & \dot{y}_2(\pi) \end{pmatrix} $$

gives together with $\mathbf{\Phi}(-\pi, 0) = \mathbf{\Phi}(\pi, 0)^{-1}$ the identity $y_1(\pi) = \dot{y}_2(\pi)$ (see [17]). Hence, it holds that $\Delta = \text{trace}(\mathbf{\Phi}_T) = 2y_1(\pi)$.

The stability of the unconstrained Mathieu equation (3) depends on the value of $\Delta$, being a function of the parameters $\alpha$ and $\beta$. The stability boundaries in the parameter plane $(\alpha, \beta)$ are given by $\Delta(\alpha, \beta) = \pm 2$, i.e. $|y_1(\pi)| = 1$. The unity of the determinant, $\det(\mathbf{\Phi}_T) = y_1(\pi)^2 - y_2(\pi)\dot{y}_1(\pi) = 1$, implies that either $y_2(\pi) = 0$ and/or $\dot{y}_2(\pi) = 0$ at a stability boundary. In other words, one can distinguish between stability boundaries for which $y_2(\pi) = 0$ and stability boundaries for which $\dot{y}_2(\pi) = 0$. The value of the discriminant $\Delta(\alpha, \beta)$ and the number $n(\alpha, \beta)$, i.e. the number of zeros of the function $y_2(t)$ on the interval $(0, \pi]$, have been computed using direct numerical integration on a grid of 1000 × 1000 points for the intervals $\alpha = -8... 32$ and $\beta = 0... 12$. The stability boundaries $\Delta(\alpha, \beta) = \pm 2$ of the unconstrained Mathieu equation are depicted in Fig. 4, which is often called the Ince-Strutt diagram.

The number $n$ changes its value in the parameter plane $(\alpha, \beta)$ if $y_2(\pi)$ changes sign. The boundary of the domains where $n$ is constant therefore agrees with those stability boundaries of the unconstrained Mathieu equation for which $y_2(\pi) = 0$, see Fig. 5. Fig. 4 indicates the number $m_1$ in the instability domains of the unconstrained Mathieu equation. Apparently, $m_1 = k$ in the $k$-th instability domain.

The stability of the equilibrium of the unilaterally constrained Mathieu equation is dependent on the number $n(\alpha, \beta)$ and the discriminant $\Delta(\alpha, \beta)$, which both depend on the system parameters $\alpha$ and $\beta$, and the restitution coefficient $\varepsilon$. The numerical results for the critical restitution coefficient are depicted in Fig. 6, being the Ince-Strutt diagram for the unilaterally constrained Mathieu equation. The level curves for $\varepsilon_c = 0, 0.2, 0.4, 0.6, 0.8$ and $1$ are shown in Fig. 6. The grey areas, being the stability domains of the unconstrained Hill's equation, have a critical restitution coefficient $\varepsilon_c = 1$. It can be seen that a decrease in the restitution coefficient enlarges the stability domain in those regions of the parameter space for which $n > 0$, especially when $n$ is large. The value of $n$ is zero in the so-called zeroth instability domain [30] (the lower left part of Fig. 6 labeled with $\varepsilon_c = 0$) and $\Delta > 2$. As follows from Theorem 2, the value of the restitution coefficient has no influence in the zeroth instability domain as the long-term behaviour is governed by non-impacting motion. The zeroth instability domain of the unconstrained Mathieu equation is therefore also unstable for the constrained Mathieu equation.

As the stability of the constrained Mathieu equation depends on the number $n$ and $\Delta$, which characterize the unconstrained Mathieu equation, one can use common approximation methods to investigate the stability properties of the unilaterally constrained Mathieu equation with $\varepsilon = 0$.

## 5.2. Averaging

A standard averaging technique [28] can be used to give an approximation for $\Delta(\alpha, \beta)$ if $\beta/\alpha$ is small. The averaging technique can be applied to the dynamics expressed in amplitude and phase coordinates (polar coordinates), see for instance [21], but this type of averaging results in an averaged system consisting of nonlinear differential equations. Here, the averaging is done on the

<figure><img src="image_5.png" /><figcaption>Fig. 4. Ince-Strutt diagram and discriminant $\Delta$ of the unconstrained Mathieu equation (3). In the grey stability domains holds $|\Delta| < 2$. The instability domains are white.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Diagram with the value of n of the Mathieu equation.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Ince-Strutt diagram of the unilaterally constrained Mathieu equation with critical value $\epsilon_c$ of the restitution coefficient.</figcaption></figure>

dynamics in comoving coordinates [27] which yields a linear set of differential equations allowing for a closed form solution of the averaged equations.

Let $\omega = 1, 2, 3, \dots$ be a resonant frequency of the Mathieu equation (3) and let $\beta = \epsilon\omega^2$ where $\epsilon$ is a small parameter measuring the relative intensity of the parametric excitation. Due to symmetry we know that $\Delta(\alpha, -\beta) = \Delta(\alpha, \beta)$ and it therefore suffices to consider $\epsilon \ge 0$. Furthermore, we consider the value of $\alpha$ to be close to resonance and set

$$ \alpha = \omega^2(1 - \epsilon\delta), \qquad (63) $$

where $\delta$ is a detuning parameter. Following [27], comoving coordinates $z_1(t)$ and $z_2(t)$ are introduced such that

$$ y(t) = z_1(t) \cos \omega t + z_2(t) \sin \omega t, \qquad (64) $$

Differentiation of $y(t)$ gives the relationship

$$ \dot{z}_1 \cos \omega t + \dot{z}_2 \sin \omega t = 0. \qquad (65) $$

Substitution of (64) and (65) in (3) yields the Mathieu equation in comoving coordinates:

$$ \begin{aligned} \dot{z}_1 &= -\epsilon\omega(\delta - 2 \cos 2t)(z_1 \cos \omega t + z_2 \sin \omega t) \sin \omega t, \\ \dot{z}_2 &= \epsilon\omega(\delta - 2 \cos 2t)(z_1 \cos \omega t + z_2 \sin \omega t) \cos \omega t. \end{aligned} \qquad (66) $$

The system (66) can be averaged over one period of oscillation, keeping $z_1$ and $z_2$ constant, i.e.

$$ \begin{aligned} \dot{z}_1 &= -\frac{\epsilon\omega}{\pi} \int_0^{\pi} (\delta - 2 \cos 2t) \cos \omega t \sin \omega t \, dt z_1 \\ &\quad -\frac{\epsilon\omega}{\pi} \int_0^{\pi} (\delta - 2 \cos 2t) \sin^2 \omega t \, dt z_2, \\ \dot{z}_2 &= \frac{\epsilon\omega}{\pi} \int_0^{\pi} (\delta - 2 \cos 2t) \cos^2 \omega t \, dt z_1 \\ &\quad + \frac{\epsilon\omega}{\pi} \int_0^{\pi} (\delta - 2 \cos 2t) \cos \omega t \sin \omega t \, dt z_2. \end{aligned} \qquad (67) $$

Consider the first resonance at $\omega = 1$. Evaluation of the averaged equation (67) gives

$$ \begin{aligned} \dot{z}_1 &= -\frac{\epsilon}{2}(1 + \delta)z_2, \\ \dot{z}_2 &= -\frac{\epsilon}{2}(1 - \delta)z_1. \end{aligned} \qquad (68) $$

The linear planar system (68) has the general solution

$$ z_1(t) = c_1 \cosh \mu t + c_2 \sinh \mu t, $$

$$ z_2(t) = -\frac{\sqrt{1-\delta^2}}{1+\delta}(c_1 \sinh \mu t + c_2 \cosh \mu t). \qquad (69) $$

with

$$ \mu = \frac{\epsilon}{2} \sqrt{1 - \delta^2}. \qquad (70) $$

In order to find an approximation for $y_1(t)$ we set $z_1(0) = 1$ and $z_2(0) = 0$ giving $c_1 = 1$, $c_2 = 0$. Hence, we obtain the approximation

$$ y_1(t) = \cosh \mu t \cos t - \frac{\sqrt{1-\delta^2}}{1+\delta} \sinh \mu t \sin t \qquad (71) $$

for the first fundamental solution and the approximation

$$ \Delta = 2y_1(\pi) = -2 \cosh \left( \frac{\epsilon\pi}{2} \sqrt{1-\delta^2} \right) \qquad (72) $$

for the discriminant. Using $\epsilon = \beta$ and $\epsilon\delta = 1 - \alpha$ for $\omega = 1$, we can express the determinant as a function of $\alpha$ and $\beta$

$$ \Delta(\alpha, \beta) = -2 \cosh \left( \frac{\pi}{2} \sqrt{\beta^2 - (1-\alpha)^2} \right). \qquad (73) $$

Near the first resonance, it therefore holds that $\Delta < -2$ if $|\delta| < 1$, or, correspondingly, if $|1 - \alpha| < |\beta|$. From (72) we calculate the largest characteristic multiplier

$$ \begin{aligned} \lambda_1 &= \frac{1}{2} \Delta - \frac{1}{2} \sqrt{\Delta^2 - 4} = -\cosh\left(\frac{\epsilon\pi}{2} \sqrt{1-\delta^2}\right) - \sinh\left(\frac{\epsilon\pi}{2} \sqrt{1-\delta^2}\right) \\ &= -e^{(\epsilon\pi/2)\sqrt{1-\delta^2}}. \end{aligned} \qquad (74) $$

Similarly, we obtain an approximation for $y_2(t)$ by setting $z_1(0) = 0$ and $z_2(0) = 1$, i.e. $c_1 = 0$ and $c_2 = -(1+\delta)/(\sqrt{1-\delta^2})$, which yields

$$ y_2(t) = - \frac{1+\delta}{\sqrt{1-\delta^2}} \sinh \mu t \cos t + \cosh \mu t \sin t. \qquad (75) $$

Clearly, if $\epsilon = 0$, then it holds that $y_2(t) = \sin t$ and the value of $n$ is on the verge of turning from zero to one. Evaluation of $y_2(\pi)$ gives

$$ y_2(\pi) = \frac{1+\delta}{\sqrt{1-\delta^2}} \sinh\left(\frac{\epsilon\pi}{2}\sqrt{1-\delta^2}\right), \qquad (76) $$
which is slightly larger than zero for $|\delta| < 1$. We therefore infer that $n=0$ holds around the resonance frequency $\omega = 1$ for small values of $\epsilon$.

The stability criterion for $n=0$ and $\Delta < -2$ reads as (see Theorem 2)

$$ \epsilon < \epsilon_c = |\lambda_1|^{-1}. \qquad (77) $$

Using the approximation (74) for $\lambda_1$, the critical coefficient of restitution is approximated near the first resonance by

$$ \epsilon_c = e^{-(\epsilon\pi/2)\sqrt{1-\delta^2}} = e^{-(\pi/2)\sqrt{\beta^2-(1-\alpha)^2}}. \qquad (78) $$

Inversely, for a given value of $\alpha$ and $\epsilon$ one can calculate the critical value of $\beta$ as

$$ \beta_c = \sqrt{(1-\alpha)^2 + \left(\frac{2}{\pi} \ln \epsilon\right)^2} \qquad (79) $$

from which we see that a small value of the restitution coefficient $\epsilon$ is beneficial for the stability of the unilaterally constrained Mathieu equation in the vicinity of the first resonance.

A comparison of the approximation (79) obtained with the averaging method and the (almost exact) numerical results is shown in Fig. 7. For a given value of $\epsilon_c$, which has been chosen to be 0.4, 0.6 and 0.8, the value of $\beta$ has been computed using (79) and is shown as dashed lines in Fig. 7. The approximation agrees fairly well with the numerical results for $\epsilon_c = 0.8$. Significant differences can be seen for $\epsilon_c = 0.4$ and $\epsilon_c = 0.6$ because $\epsilon = \beta$ can no longer considered to be small in the upper half of Fig. 7.

### 5.3. Approximation using Hill's determinant

A much better approximation of the discriminant $\Delta(\alpha, \beta)$ and the critical restitution coefficient can be obtained by using Hill's infinite determinant as discussed in Section 4.

The function $g(t) = \alpha + 2\beta \cos 2t$ of the Mathieu equation can be represented by $g(t) = \alpha + \beta(e^{2it} + e^{-2it})$ and the Fourier coefficients are therefore $g_0 = \alpha$ and $g_1 = g_{-1} = \beta$ whereas all other Fourier coefficients are zero. The determinant $D(0)$, see (59), therefore reads as

$$ D(0) = \begin{vmatrix} \dots & \dots & \dots & \dots & \dots \\ \dots & 1 & \frac{\beta}{\alpha-4^2} & 0 & 0 & 0 \\ \dots & \frac{\beta}{\alpha-2^2} & 1 & \frac{\beta}{\alpha-2^2} & 0 & 0 \\ \dots & 0 & \frac{\beta}{\alpha} & 1 & \frac{\beta}{\alpha} & 0 \\ \dots & 0 & 0 & \frac{\beta}{\alpha-2^2} & 1 & \frac{\beta}{\alpha-2^2} \\ \dots & 0 & 0 & 0 & \frac{\beta}{\alpha-4^2} & 1 \\ \dots & \dots & \dots & \dots & \dots & \dots \end{vmatrix}. \qquad (80) $$

The value of $D(0)$ can be approximated by the determinant of the central $5 \times 5$ block

$$ D(0) \approx 1 - \frac{4\alpha - 32}{\alpha(\alpha - 4)(\alpha - 16)} \beta^2 + \frac{3\alpha - 32}{\alpha(\alpha - 4)^2(\alpha - 16)^2} \beta^4. \qquad (81) $$

This approximation can be improved by calculating larger central $k \times k$ blocks in (80). However, the coefficients of $\beta^2$ and $\beta^4$ in (81) slightly change if larger central blocks are considered. Using (61) the discriminant $\Delta$ can be approximated by

$$ \Delta(\alpha, \beta) = 2 - 4 \sin^2 \left( \frac{1}{2} \pi \sqrt{\alpha} \right) \\ \cdot \left( 1 - \frac{4\alpha - 32}{\alpha(\alpha - 4)(\alpha - 16)} \beta^2 + \frac{3\alpha - 32}{\alpha(\alpha - 4)^2(\alpha - 16)^2} \beta^4 \right). \qquad (82) $$

<figure><img src="image_8.png" /><figcaption>Fig. 7. Approximate Ince-Strutt diagram around the first resonance by using averaging (dashed lines) and Hill's determinant (dash-dot lines).</figcaption></figure>

Hence, the critical restitution coefficient is determined by the equation

$$ \epsilon_c^{m_1} + \frac{1}{\epsilon_c^{m_1}} = -2 + 4 \sin^2 \left( \frac{1}{2} \pi \sqrt{\alpha} \right) \\ \cdot \left( 1 - \frac{4\alpha - 32}{\alpha(\alpha - 4)(\alpha - 16)} \beta^2 + \frac{3\alpha - 32}{\alpha(\alpha - 4)^2(\alpha - 16)^2} \beta^4 \right), \qquad (83) $$

or equivalently by (62). Inversely, for a given value of $\alpha$ and $\epsilon$ one can calculate an approximation of the critical value of $\beta$, based on the terms in (82) up to order $\beta^2$, as

$$ \beta_c^2 = \frac{\alpha(\alpha-4)(\alpha-16)}{4\alpha-32} \left( 1 - \frac{\epsilon^{m_1} + 2 + \frac{1}{\epsilon^{m_1}}}{4 \sin^2 \left( \frac{1}{2} \pi \sqrt{\alpha} \right)} \right). \qquad (84) $$

This approximation is valid (for small $\beta$) in each of the instability domains but the value of $m_1$ should be a priori known. For the Mathieu equation, however, it holds that $m_1 = k$ in the $k$-th instability domain.

A comparison of the approximation (84) in the first instability region ($m_1 = 1$) obtained with Hill's determinant (dash-dot lines) and the approximation (79) obtained with the averaging method and the (almost exact) numerical results is shown in Fig. 7. Clearly, the approximation (84) is much better than the approximation (79). However, using the averaging method one is able to estimate the value of $n$ (and therefore $m_1$) as a function of $\alpha$ and $\beta$, which cannot (easily) be done by using the method of Hill's determinant.

## 6. Conclusions and discussion

In this paper the stability conditions of the unilaterally constrained Hill's equation have been addressed in detail using Floquet theory and Lyapunov techniques. It has been shown that the stability of the equilibrium of the unilaterally constrained Hill's equation depends on the discriminant $\Delta$ and the number $n$ (i.e. the number of zeros of the second fundamental solution within one period) of the unconstrained Hill's equation and on the restitution coefficient $\epsilon$. The remarkable simplicity of the unilaterally constrained Hill's equation stems from the fact that, although the system can be considered to be strongly non-linear due to the presence of the unilateral constraint, its Poincaré map is cone-wise linear. The cone-wise linearity originates from the
homogeneity of the linear differential equation and the homogeneity of the impact map.

The practical merit of the paper is that a precise estimation of the critical restitution coefficient can be obtained by calculating the fundamental solutions of the unconstrained Hill's equation using direction numerical integration methods (ODE-solvers). In addition, two approximation methods are proposed which give closed form expressions for the critical restitution coefficient: the averaging method (Section 5.2) and the method of Hill's infinite determinant (Sections 4 and 5.3). A comparison of the approximation techniques applied to the unilaterally constrained Mathieu equation has been given in Section 5.

The averaging method in comoving coordinates, which is employed in Section 5.2, gives the same approximation of the critical restitution coefficient as obtained by the averaging method in [21], Section 2.2.2. In [21], a two-dimensional impact event map is constructed for the first instability domain using the averaging method in amplitude and phase coordinates and assuming that the time difference between consecutive impacts equals $t_{i+1}-t_i = \pi+\mathcal{O}(\epsilon)$. In Section 5.2, the averaging method in comoving coordinates is used to obtain an approximation of the two-dimensional Poincaré map. This approximate Poincaré map is cone-wise linear as opposed to the approximate impact event map of [21] which is fully non-linear. The use of higher-order averaging methods to improve the approximation for larger values of $\epsilon$ becomes cumbersome as it takes a much larger effort and, in addition, an improved approximation for $t_{i+1}-t_i$ needs to be obtained for averaging in amplitude and phase coordinates.

The approximation of the critical restitution coefficient using Hill's infinite determinant, see Sections 4 and 5.3, is very accurate and can easily be improved by considering larger central blocks for the determinant $D(0)$. However, the method using Hill's infinite determinant gives no direct way to determine the number $n$.

The analysis has shown that the impact time instants are defined by the zero-crossings of the solution of the unconstrained Hill's equation, which are always separated in time. Accumulation points of impacts can therefore not exist in the unilaterally constrained Hill's equation (4) as has been proven in Theorem 1. If, however, the location of the constraint is moved to a non-zero position, which does not agree with the equilibrium of the unconstrained system, i.e. $x(t) \ge x_c > 0$, then accumulation points are possible. The unilaterally constrained Hill's equation with $x_c > 0$ does not have an equilibrium and the contact force $\lambda$ in Eq. (5) does not vanish. The framework of non-smooth dynamics is therefore needed to describe and understand the dynamics of the unilaterally constrained Hill's equation with non-zero constraint position. Moreover, the non-linear dynamics becomes far more complicated because the homogeneity of the impact conditions (27), and also the cone-wise linearity of the Poincaré map, is lost.

The present paper gives more insight in the stability properties of Hill's equation with unilateral constraint, being an prototype of a parametrically excited non-smooth dynamical system with impulsive motion. Further research will focus on the application and extension of the obtained results to the stability analysis of multi-degree-of-freedom autoparametric systems with unilateral constraints.

## Appendix A. Floor function and fractional part

**Definition 1 (Floor function and fractional part).** With $[x]$ we denote the floor function defined by

$$[x] = \{\max k \in \mathbb{Z} | k \le x\} \qquad (85)$$

and with $\{x\}$ the fractional part defined by

$$\{x\} = x - [x]. \qquad (86)$$

The fractional part can be expressed in trigonometric functions as

$$\{x\} = \frac{1}{2} - \frac{1}{\pi} \arctan(\cot(\pi x)), \qquad (87)$$

where $\cot(\pi k) = +\infty$ for $k \in \mathbb{Z}$. The floor function has the following property:

$$[x] + [y] \le [x+y] \le [x] + [y] + 1, \qquad (88)$$

and the equality $[x] + [y] = [x+y]$ holds only if $\{x\} + \{y\} < 1$, or using (87), if $\cot(\pi x)+\cot(\pi y) > 0$. More precisely, the relations (88) can be formulated as

$$[x+y] = \begin{cases} [x] + [y] & \text{if } \cot(\pi x) + \cot(\pi y) > 0, \\ [x] + [y] + 1 & \text{if } \cot(\pi x) + \cot(\pi y) \le 0. \end{cases} \qquad (89)$$

The floor function can be used to count the number of zeros of sinusoidal functions. The function $f_1(x) = \sin(x)$ has the zeros $k\pi$ with $k \in \mathbb{Z}$ and the number of zeros on the interval $(0,a]$ therefore amounts to

$$m_1 = \lfloor \frac{a}{\pi} \rfloor. \qquad (90)$$

The function $f_2(x) = \sin(x+b)$ has the number of zeros

$$m_2 = \left\lfloor \frac{a+b}{\pi} \right\rfloor - \left\lfloor \frac{b}{\pi} \right\rfloor, \qquad (91)$$

which, using (89), can be expressed as

$$m_2 = \begin{cases} m_1 & \text{if } \cot(a) + \cot(b) > 0, \\ m_1 + 1 & \text{if } \cot(a) + \cot(b) \le 0. \end{cases} \qquad (92)$$

## References

[1] M.V. Bartuccelli, J.H.B. Deane, G. Gentile, S.A. Gourley, Global attraction to the origin in a parametrically driven nonlinear oscillator, Applied Mathematics and Computation 153 (1) (2004) 1–11.

[2] M.V. Bartuccelli, G. Gentile, K.V. Georgiou, On the stability of the upside-down pendulum with damping, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 458 (2018) (2002) 255–269.

[3] T. Birkandan, M. Hortaçsu, Examples of Heun and Mathieu functions as solutions of wave equations in curved spaces, Journal of Physics A: Mathematical and Theoretical 40 (2007) 1105–1116.

[4] H.W. Broer, I. Hoveijn, M. Van Noort, C. Simó, G. Végter, The parametrically forced pendulum: a case study in 1½ degree of freedom, Journal of Dynamics and Differential Equations 16 (4) (2004) 897–947.

[5] H.W. Broer, I. Hoveijn, M. Van Noort, G. Végter, The inverted pendulum: a singularity theory approach, Journal of Dynamics and Differential Equations 157 (1) (1999) 120–149.

[6] H.W. Broer, M. Levi, Geometrical aspects of stability theory for Hill's equations, Archive for Rational Mechanics and Analysis 131 (1995) 225–240.

[7] B. Brogliato, Absolute stability and the Lagrange–Dirichlet theorem with monotone multivalued mappings, Systems & Control Letters 51 (2004) 343–353.

[8] M. Cartmell (Ed.), Introduction to Linear, Parametric and Nonlinear Vibrations, Chapman and Hall, London, 1990.

[9] B.J. Daiuto, T.T. Hartley, S.P. Chicatelli, The Hyperbolic Map and Applications to the Linear Quadratic Regulator, Lecture Notes in Control and Information Sciences, vol. 110, Springer, Berlin, 1989.

[10] Ch. Glocker, Set-Valued Force Laws, Dynamics of Non-Smooth Systems, Lecture Notes in Applied Mechanics, vol. 1, Springer-Verlag, Berlin, 2001.

[11] A.P. Ivanov, The stability of mechanical systems subjected to impulsive actions, Journal of Applied Mathematics and Mechanics 65 (4) (2001) 617–629.

[12] D.W. Jordan, P. Smith, Nonlinear Ordinary Differential Equations, Oxford University Press, New York, 2007.

[13] J.P. LaSalle, The Stability and Control of Discrete Processes, Applied Mathematics Sciences, vol. 62, Springer, New York, 1961.

[14] R.I. Leine, T.F. Heimsch, Global uniform symptotic attractive stability of the non-autonomous bouncing ball system, Physica D (2011), http://dx.doi.org/10.1016/j.physd.2011.04.013.

[15] R.I. Leine, H. Nijmeijer, Dynamics and Bifurcations of Non-Smooth Mechanical Systems, Lecture Notes in Applied and Computational Mechanics, vol. 18, Springer Verlag, Berlin, 2004.
[16] R.I. Leine, N. van de Wouw, Stability and Convergence of Mechanical Systems, with Unilateral Constraints, Lecture Notes in Applied and Computational Mechanics, vol. 36, Springer Verlag, Berlin, 2008.

[17] W. Magnus, S. Winkler, Hill's Equation, Dover Publications, New York, 1979.

[18] A.H. Nayfeh, D.T. Mook, Nonlinear Oscillations, Wiley-Interscience, Chichester, 1979.

[19] W. Paul, Electromagnetic traps for charged and neutral particles, Reviews of Modern Physics 62 (1990) 531-540.

[20] F. Pellicano, A. Fregolent, A. Bertuzzi, F. Vestroni, Primary and parametric non-linear resonances of a power transmission belt: experimental and theoretical analysis, Journal of Sound and Vibration 244 (4) (2001) 669-684.

[21] D.D. Quinn, The dynamics of two parametrically excited pendula with impacts, International Journal of Bifurcations and Chaos 15 (6) (2005) 1975-1988.

[22] S.W. Shaw, R.H. Rand, The transition to chaos in a simple mechanical system, International Journal of Nonlinear Mechanics 24 (1) (1989) 41-56.

[23] A. Tondl, Quenching of Self-Excited Vibrations, Elsevier, Amsterdam, 1991.

[24] A. Tondl, T. Ruijgrok, F. Verhulst, R. Nabergoj, Autoparametric Resonance in Mechanical Systems, Cambridge University Press, New York, 2000.

[25] S.R. Valluri, R. Biggs, W. Harper, C. Wilson, The significance of the Mathieu-Hill differential equation for Newton's apsidal precession theorem, Canadian Journal of Physics 77 (5) (1999) 393-407.

[26] M. van Noort, The Parametrically Forced Pendulum, A Case Study in 1½ Degree of Freedom, Ph.D. thesis, Rijksuniversiteit Groningen, The Netherlands, 2001.

[27] F. Verhulst, Nonlinear Differential Equations and Dynamical Systems, Springer, Berlin, 1996.

[28] F. Verhulst, Parametric and auto-parametric resonance, Acta Applicandae Mathematicae 70 (2002) 231-264.

[29] E.T. Whittaker, G.N. Watson, A Course of Modern Analysis, Cambridge University Press, Cambridge, 1927.

[30] V.A. Yakubovich, V.M. Starzhinskii, Linear Differential Equations with Periodic Coefficients 1 and 2, Wiley, New York, 1975.