Non-smooth stability analysis of the parametrically excited impact oscillator

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 archetype 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 coefﬁcient is presented based on Hill’s inﬁnite determinant and this approxima-tion can be made arbitrarily accurate. A comparison of numerical and theoretical results is presented for the unilaterally constrained Mathieu equation.


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 where W is the angle of the pendulum with the downward vertical and aðtÞ ¼ € 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 W ¼ 0 and W ¼ p is described by Hill's equation [6,17,30] where gðtÞ ¼ gðt þ pÞ is a real piece-wise continuous function.
Harmonic excitation leads to the so-called Mathieu equation being a special case of Hill's equation with gðtÞ ¼ aþ2b 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 constraint at W ¼ 0 and the impact law 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 W % WÀ 1 6 W 3 , that is to say, the 'nonlinear Mathieu equation' 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 b and near a chosen resonance frequency ffiffiffi a p ¼ 1; 2,3 . ... 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 where e is the restitution coefficient.The unilateral constraint limits the dynamics to xðtÞ Z 0 and imposes a Newtonian impact law More physically correct would be to formulate the unilaterally constrained system (4) in the framework of nonsmooth dynamics [7,10,15,16] and to let a contact force l and an impulsive contact force L (both per unit mass) appear in the equation of motion and impact equation as Lagrange multipliers, i.e.
together with the set-valued contact law (Signorini's law) and the generalized Newtonian impact law where 0 r a ?bZ 0 denotes the inequality complementarity condition a Z 0, bZ 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 l is for this system always zero since the external force gðtÞxðtÞ and inertial force € x vanish during persistent contact for which xðtÞ ¼ _ xðtÞ ¼ € 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 _ x þ ðt i Þþe _ 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 W to 0 r WðtÞrp, 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 € x ðtÞþa _ x ðtÞþ g ðtÞxðtÞ ¼ 0 can easily be transformed to the standard form (4) by using the transformation xðtÞ ¼ e ð1=2Þat xðtÞ such that gðtÞ ¼ g ðtÞÀ 1 4 a 2 [17].The unilaterally constrained Hill's equation ( 4) is therefore an archetype 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.

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) € yðtÞþgðtÞyðtÞ ¼ 0, gðtÞ ¼ gðt þ pÞ has two continuously differentiable solutions y 1 ðtÞ and y 2 ðtÞ with the initial conditions 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Þ ¼ RðtÞ cos cðtÞ, y 2 ðtÞ ¼ RðtÞ sin cðtÞ, with RðtÞ40 for all t and the differential equations 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 with the state vector yðtÞ ¼ ðyðtÞ _ yðtÞÞ T and the time-dependent system matrix AðtÞ.
More generally, the fundamental solution matrix Uðt 1 ,t 0 Þ is defined as the mapping holds from which one can deduce the inverse of the fundamental solution matrix to be Uðt 1 ,t 0 Þ ¼ Uðt 0 ,t 1 Þ À1 .The fundamental solution matrix U T ¼ Uðp,0Þ is referred to as the monodromy matrix.The trace of the monodromy matrix is known as the discriminant of Hill's equation [17].The monodromy matrix U T has the characteristic equation because detðU T Þ ¼ 1, and the eigenvalues which are called characteristic multipliers (or Floquet multipliers).The characteristic multipliers are reciprocal in the sense that l 1 ¼ 1=l 2 .The discriminant D ¼ l 1 þ l 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 9D9 o2, then the characteristic multipliers l 1 ¼ l 2 are complex conjugated and located on the unit circle because 9l 1;2 9 ¼ 1.
If it holds that 9D9 42, then the characteristic multipliers l 1 and l 2 are real and distinct and we define the order by The normalization of the eigenvectors v 1 and v 2 is chosen such that v 1;2 lie in the first or fourth quadrant and Jv 1;2 J ¼ 1.If y 2 ðpÞ ¼ 0 then either l 1 Ày 1 ðpÞ ¼ 0 or l 2 Ày 1 ðpÞ ¼ 0 and the expression ( 16) for v 1 or v 2 degenerates.In this case there still exist two linearly independent real eigenvectors of which either v 1 or v 2 is equal to ð0 1Þ T .
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 halfopen interval ð0,p.The floor function bÁc and fractional part fÁg will be used to count the zeros.The definition and properties of the floor function and fractional part are given in Appendix A.
Proof.According to Proposition 1, one may write y 2 ðtÞ ¼ RðtÞ sin cðtÞ with RðtÞ40 for all t.The zeros of y 2 ðtÞ are therefore given by the condition sin cðtÞ ¼ 0, i.e. cðtÞ ¼ p,2p,3p, . . .and using (90) it therefore holds that Consider a solution y(t) of (2) with initial condition yð0Þ ¼ r 0 cos y 0 Z0 and _ yð0Þ ¼ r 0 sin y 0 4 0, where r 0 40 and y 0 A ðÀp=2,p=2.With mðtÞ A N 0 we will denote the rightcontinuous monotonically increasing step function which counts the number of reflections in the interval ð0,t such that ðÀ1Þ mðtÞ yðtÞ ¼ 9yðtÞ9.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ðpÞ and use the short-hand notation m ¼ mðpÞ.By definition, it holds that m¼n for y 0 ¼ p=2, because n is defined as the number of zeros of y 2 ðtÞ on the interval ð0,p and yðtÞ ¼ r 0 y 2 ðtÞ for y 0 ¼ p=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 y 0 Z 0 and _ yð0Þ ¼ r 0 sin y 0 , where r 0 4 0 and Àp=2 o y 0 r p=2.Let m denote the number of zeros of y(t) on the interval ð0,p.It holds that 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Þ _ yð0Þ which we write in polar coordinates using Proposition 1 as Remark.Some care needs to be taken for the case when y 2 ðpÞ ¼ 0, i.e. cðpÞ ¼ np, n A N 0 .In this case it holds that cotðcðpÞÞ ¼ þ 1 and the inequality ( 23) is therefore fulfilled which implies that m ¼n.Moreover, if cðpÞ ¼ np, then it can immediately be verified from ( 22) that m¼n.

The unilaterally constrained Hill's equation
In this section the unilaterally constrained Hill's equation is analyzed which consists of Hill's differential equation which holds for almost all t and xðtÞ Z 0, and the Newtonian impact law we denote the pre-and post-impact velocity and with eA½0; 1 the restitution coefficient.The velocity _ xðtÞ of the unilaterally constrained Hill's equation is considered to be right-continuous by convention, i.e. _ xðtÞ ¼ _ x þ ðtÞ, and the initial condition ðxð0Þ, _ 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 y 0 , _ xð0Þ ¼ r 0 sin y 0 , then it necessarily holds that r 0 Z 0 and y 0 A ðÀp=2,p=2.The solution x(t) of the unilaterally constrained Hill's equation is therefore confined to the domain The impact law (25) becomes active when The impact law can therefore formally be written for position and velocity as a homogeneous map 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 _ xð0Þ ¼ _ yð0Þ such that ðxð0Þ _ xð0ÞÞ T A D. It holds that xðtÞ ¼ yðtÞ during the interval t A ½0,t 1 , where t 1 is the time-instant for which impact occurs in the unilaterally constrained Hill's equation.The solution x(t) after t ¼ t 1 is reflected and scaled with e, i.e.

xðtÞ ¼ ÀeyðtÞ ð 28Þ
for t A ½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Þ ¼ ðÀeÞ mðtÞ yðtÞ ð 29Þ and, using m ¼ mðpÞ, we obtain the relationship 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].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Þ, _ xð0ÞÞ such that m is infinite.
However, the angle cðpÞ ¼ R p 0 RðtÞ À2 dt is finite because RðtÞ is continuous and strictly positive for all t, see Proposition 1.The number n ¼ bcðpÞ=pc is therefore finite and the proof follows 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 xðtÞ ¼ ðxðtÞ _ xðtÞÞ T , the system is written in first-order form as with the system matrix AðtÞ and impact map G given by The homogeneity of the impact conditions ( 27) allows us to write the impact map G as ÀeI and the matrix G therefore commutes with any arbitrary matrix.The solution of the unilaterally constrained Hill's equation can be obtained by concatenation of nonimpulsive 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 _ xðtÞ ¼ AðtÞxðtÞ.The fun- and it necessarily holds that n T x The impulsive motion at each collision time-instant t i is described by the impact map G. Concatenation of non-impulsive and impulsive motion gives The impact map G ¼ ÀeI commutes with all the fundamental solution matrices UðÁ,ÁÞ in (34) and Eq.(34) therefore simplifies to by a half-hyperplane S ¼ fx A D9y ¼ y c g, where y c is defined by (21), The switching of the number m on the half-hyperplane S leads together with (36) to a piece-wise linear, or, more precisely, a cone-wise linear Poincare ´map.
Proposition 4 (Poincare ´map).Let x j denote the solution xðpjÞ, j A N 0 , of the unilaterally constrained Hill's equation with The discrete time-evolution of x j is given by the cone-wise linear Poincare ´map Proof.From (35) it is clear that x j þ 1 ¼ ðÀeÞ m Uðp,0Þx j where m depends on x j .Write x j in polar coordinates as x j ¼ ðr j cos y j r j sin y j Þ T .If x j A D 1 , then it holds that y j 4 y c which implies that m¼n.Similarly, if The cone-wise linearity of the Poincare ´map introduced in Proposition 4 suggests to analyze the map by using polar coordinates x j ¼ ðr j cos y j r j sin y j Þ T , i.e.
Evaluation of the Poincare ´map for e40 yields where as y j þ 1 is not defined if r j þ 1 ¼ 0 for e ¼ 0. The map y j /y j þ 1 is therefore autonomous as it does not depend on r j which expresses the cone-wise linearity of the Poincare ´map.The map y j /y j þ 1 can be simplified even further by using a non-linear transformation.
Proposition 5. Let q j ¼ y 1 ðpÞþy 2 ðpÞtan y j and y 2 ðpÞ a0.It holds that where Proof.Substitution of (39) in q j ¼ y 1 ðpÞþy 2 ðpÞ tan y j gives 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 n ¼ Q ðq n Þ, and the fixed points therefore satisfy in which we recognize the characteristic equation ( 14) of the monodromy matrix U T .Hence, the fixed points q n agree with the real characteristic multipliers l 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]: p and we define 9q n 1 9 4 9q n 2 9 and equivalently 9l 1 9 4 9l 2 9. Furthermore, because l 1 l 2 ¼ 1, it holds that 9q n 1 9 41 4 9q n 2 9.The stability of the fixed points q n 1;2 is determined by the derivative of the map Q 0 ðqÞ ¼ 1=q 2 from which we see that q n 1 is asymptotically stable (9Q 0 ðq n 1 Þ9 o 1) whereas with initial condition q 0 has a closed form solution given by as we can simply verify or derive with the z-transformation [9].If q 0 aq n 2 ¼ l À1 1 ¼ l 2 , then the solution q j is attracted to q n 1 which follows from lim j-1 The map has a single fixed point q n ¼ 1 2 D ¼ 71.The map with initial condition q 0 has the closed form solution The fixed point q n is therefore unstable but globally attractive.Furthermore, if q j q n o 0 then it holds that q j þ 1 q n 4 0.
9D9 o 2: The map has no fixed points and the solution is quasiperiodic, wandering between R À 0 and R þ (see [9]).The number of iterations that the solution q j remains in R À 0 or R þ is bounded.
The above results on the map q j þ 1 ¼ Q ðq j Þ for 9D9 4 2 can be understood by noting that the values of q n 1;2 define for y 2 ðpÞ a 0 the angles y n 1 and y n 2 of the eigenvectors v 1 and v 2 (see ( 16)) of the monodromy matrix U T : tan y n If x j is chosen in the direction of an eigenvector, say x j ¼ r j v 1 , then the next iteration point x j þ 1 will be again in the direction of that eigenvector, i.e.
In other words, if y j ¼ y n 1 then also or, equivalently, that q n 1 is a fixed point of the map q j þ 1 ¼ Q ðq j Þ.The fixed points q n 1;2 ¼ l 1;2 are either both positive if D42 or both negative if Do2 and the sign of y 2 ðpÞ is given by ðÀ1Þ n .Hence, if 9D9 42 and y 2 ðpÞ a0, then it holds that v 1;2 A D 1 for ðÀ1Þ n D42 and v 1;2 A D 2 for ðÀ1Þ n Do2.If 9D9 42 and y 2 ðpÞ ¼ 0 then it holds that y 1 ðpÞ _ y 2 ðpÞ ¼ 1 and 9y 1 ðpÞþ _ y 2 ðpÞ9 42 from which we deduce that _ y 2 ðpÞ a 0 and signð _ y 2 ðpÞÞ ¼ signðDÞ.If y 2 ðpÞ ¼ 0 then it holds that n 4 0 and n is even for _ y 2 ðpÞ 40 and n is odd for _ y 2 ðpÞ o0.Consequently, the condition 9D9 4 2 and y 2 ðpÞ ¼ 0 implies that ðÀ1Þ n D42.Moreover, it holds that D 1 ¼ D for y 2 ðpÞ ¼ 0. We conclude that, if 9D9 42, then the location of the eigenvectors is determined by the condition for all values of y 2 ðpÞ.
The attractivity of the fixed point q n 1 implies that the solution x j will be drawn towards the eigenvector v 1 when j-1, because v 1 belongs to the characteristic multiplier which is largest in magnitude.The long-term behaviour of the dynamics is therefore determined by 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 D. In the following propositions, the spectral decomposition of the monodromy matrix U T is used, where V ¼ ðv 1 v 2 Þ is the matrix of eigenvectors and K ¼ diagðl i Þ is the diagonal matrix of eigenvalues of U T .Such a spectral decomposition can be made if the characteristic multipliers l 1;2 are distinct, implying 9D9 a2.Let c j ¼ ðc 1j c 2j Þ T be the coordinates on the basis of v 1 and v 2 such that with P ¼ V ÀT V À1 being dependent on the matrix of eigenvectors V.The positive definiteness of P ¼ P T 4 0 implies that V is a positive definite function.It holds that x j þ 1 ¼ ðÀeÞ m U T x j and, using (44), the increment of V therefore gives If e n 9l 1 9 o 1, then it holds that e 2m l 2 1;2 o 1 for m¼n and m ¼ n þ 1 which implies that V j þ 1 o V j for x j a 0. Hence, V is a quadratic time-independent Lyapunov function.The origin is therefore globally uniformly asymptotically stable if e n 9l 1 9o 1.
We now prove the last part of the proposition which states that if v 1;2 A D 1 and e n 9l 1 9 4 1, then the origin is unstable.A solution x j which starts in the direction of the eigenvector v 1 A D 1 , i.e. x 0 ¼ c 10 v 1 , will remain along v 1 and therefore in D 1 for all j Z 0. It therefore holds that m¼ n during all iterations of the Poincare ḿap.The solution x j ¼ c 1j v 1 is therefore given by c 1j ¼ ðÀeÞ n l 1 c 1,jÀ1 ¼ ððÀeÞ n l 1 Þ j c 10 which grows unbounded for e n 9l 1 9 41.The value of c 10 can be chosen arbitrarily small and the origin is therefore unstable.& Proposition 7 (Stability for the real case v 1;2 A D 2 ).If v 1;2 A D 2 then the equilibrium x n ¼ 0 of the unilaterally constrained Hill's equation is globally uniformly asymptotically stable if e n þ 1 9l 1 9 o 1 and unstable if e n þ 1 9l 1 9 4 1.
Proof.If v 1;2 A D 2 then the eigenvalues are real and distinct such that l 2 1 4 1 and l 2 2 ¼ l À2 1 o 1.Moreover, if v 1;2 A D 2 then it holds that y 2 ðpÞ a 0 and therefore cos y n 1;2 40.Consider the discrete Lyapunov candidate function where y n 1 and y n 2 are the angles of v 1 and v 2 defined by (41).The positive definiteness of V is assured through P ¼ P T 4 0. It holds that x j þ 1 ¼ ðÀeÞ m U T x j and the increment of V yields The decrease in the Lyapunov function V depends on the location of x j : and it holds that y c o y j r p=2.
Hence, it holds that y n 1 oy n 2 o y j , i.e. x j lies inside the cone spanned by Àv 1 and v 2 , which implies that c 1j o0 and c 2j 4 0. Furthermore, from the positiveness of together with c 1j we conclude that 9c 1j =c 2j 9 r 9cos y n 2 =cos y n 1 9.We now define Using l 1 ¼ l À1 2 4 1 we obtain: The expressions between the brackets are strictly negative for all m Z 0 which implies that V j þ 1 ÀV j o 0.
Hence, the quadratic positive definite function V is strictly decreasing for x j a 0. The origin is therefore globally uniformly asymptotically stable if e n þ 1 9l 1 9 o 1. Instability for e n þ 1 9l 1 9 41 can be proven by considering an initial condition in the direction of v 1 (see the proof of Proposition 6).& Proposition 8 (Stability for the complex conjugated case).If 9D9 o2 and eo1 then the equilibrium x n ¼ 0 of the unilaterally constrained Hill's equation is globally uniformly asymptotically stable.
Proof.If 9D9o 2 then the eigenvalues l 1;2 are complex conjugated (and distinct) with magnitude 9l 1;2 9 ¼ 1 and v 1 ¼ v 2 .The coordinates c 1j and c 2j are in this case complex conjugated, i.e. c 1j ¼ c 2j and 9c 1j 9 ¼ 9c 2j 9, because x j is real.Consider the discrete Lyapunov candidate function being a hermitian matrix.Clearly, it holds that P ¼ P n 4 0 and V is positive definite.Using The increment is non-positive, V j þ 1 ÀV j r 0, because m Z 0 which proves that the equilibrium is uniformly stable.It holds that m¼n for x j A D 1 and m ¼ n þ1 for . However, the dynamics of the map q j þ 1 ¼ Q ðq j Þ proves that the cones D 1 and D 2 do not have a positively invariant subset (or sub-cone) other than the origin if 9D9 o2, i.e. no limit set exists which lies exclusively in D 1 \f0g or exclusively in D 2 \f0g.This implies that the solution x j will wander between D 1 and D 2 and the Lyapunov function strictly decreases whenever x j A D 2 \f0g.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 D 1 and the largest positively invariant set in D 1 is the origin.Hence, it holds that lim j-1 V j ¼ 0, which proves that the equilibrium is globally uniformly attractive.& Proposition 9.If D ¼ À2, n¼0 and eo1 then the equilibrium x n ¼ 0 of the unilaterally constrained Hill's equation is globally uniformly asymptotically stable.
Proof.The condition n ¼0 implies that y 2 ðpÞ 40.If D ¼ À2 then the eigenvalues of U T are l 1;2 ¼ À1.The matrix A 1 ¼ ðÀeÞ n U T is not a stable matrix for D ¼ À2 and n¼ 0. The matrix A 2 ¼ ÀeA 1 is a stable matrix for eo1.
unique fixed point q n ¼ À1 which is globally attractive.Furthermore, if q 0 4 0 then q j oÀ1 for all j 4 0. The cone D 2 has therefore a globally attractive sub-cone.The long-term dynamics is therefore governed by the stable matrix A 2 .& Proposition 10.If D ¼ 2 and n ¼0 then the equilibrium x n ¼ 0 of the unilaterally constrained Hill's equation is not attractive and unstable for all eZ0.
Proof.The condition n ¼0 implies that y 2 ðpÞ 4 0. If D ¼ 2 then the eigenvalues of U T are l 1;2 ¼ 1 and the map q j þ 1 ¼ Q ðq j Þ has a unique fixed point q n ¼ 1 which is globally attractive.Furthermore, if q 0 o0 then q j 41 for all j 40.The cone D 1 has therefore a globally attractive sub-cone.The long-term dynamics is therefore governed by the matrix A 1 ¼ U T which is not a stable matrix.The initial condition x 0 ¼ ðr 0 cos y 0 r 0 sin y 0 Þ T with tan y 0 ¼ ð1Ày 1 ðpÞÞ= y 2 ðpÞ will lead to a solution x j ¼ x 0 for all r 0 Z0.Hence the equilibrium is not attractive.Moreover, A 1 is non-diagonalizable and A j 1 will diverge for j-1 which implies unboundedness of solutions.& Proposition 11.If 9D9 r2, n 4 0 and eo1 then the equilibrium x n ¼ 0 of the unilaterally constrained Hill's equation is globally uniformly asymptotically stable.
Proof.The spectral radius of U T equals unity for 9D9 r2.Hence, if in addition n Z 1, then the spectral radii of A 1 and A 2 are strictly smaller than unity, i.e. the matrices A 1 and A 2 are stable matrices.The stability of A 1 implies that there exist symmetric positive definite matrices P and Q 1 such that Using P we define the matrix Q 2 such that and express Q 2 using A 2 ¼ ÀeA 1 and Q 1 as With P 4 0, Q 1 4 0 and eo1 we infer that Q 2 4 0. Hence, there exists a common quadratic Lyapunov function Vðx j Þ ¼ x T j Px j which proves that the origin is globally uniformly asymptotically stable.& The previous propositions are summarized in the following theorem, being the main result of the paper.
If 9D9 42 and e4e c , then the equilibrium x n ¼ 0 is unstable.
Proof.Propositions 6 and 10 prove that, if n ¼0 and DZ2, then the equilibrium is not attractive for all eZ0.The equilibrium is therefore in this case not globally uniformly asymptotically stable which is expressed by e c ¼ 0. The proof of e c ¼ 9l 1 9 À1=n for n 4 0 and ðÀ1Þ n D42 follows from (42) and Proposition 6.The proof of e c ¼ 9l 1 9 À1=ðn þ 1Þ for n Z 0 and ðÀ1Þ n DoÀ2 follows from (42) and Proposition 7. The proof of e c ¼ 1 for n 4 0 and À2 r Dr2 follows from Proposition 11 and for n¼0 and À2 rD o2 from Propositions 8 and 9. & For 9D9 42 the results of Theorem 2 can be put in a more tangible form by introducing the number m 1 as The number m 1 is the value of m, i.e. the number of zeros of y(t) on the interval ð0,p, of a solution curve yðtÞ which starts in the direction of the first eigenvector, i.e. yð0Þ ¼ v 1 .

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 D 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 where g Àk ¼ g k .In this section we will assume that ffiffiffiffiffi g 0 p a 0; 2, 4; 8, . .., i.e. the even parametric resonances are avoided.Similarly, the first eigensolution f 1 ðtÞ ¼ e st p 1 ðtÞ, see (17), can be written as a complex Fourier series where the characteristic exponent s is related to the discriminant through D ¼ 2 coshðpsÞ.Substitution of the Fourier representations ( 55) and ( 56) in Hill's equation ( 2) yields the condition Reordering terms gives the equality The requirement that (57) has to be fulfilled for all t leads to an infinite set of homogeneous equations for the Fourier coefficients c k .The linear homogeneous system of equations has only a non-trivial solution if the infinite determinant vanishes, where 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ðsÞ can be expressed as Using the identity sin 2 1 2 ips , the determinant condition is written as which agrees with [17], Theorem 2.9.The discriminant D can therefore be expressed as Furthermore, it holds that 9l p , see (15).If 9D9 4 2, then the critical restitution coefficient is given by e c ¼ 9l 1 9 À1=m 1 , see Corollary 1.The critical restitution coefficient can therefore be calculated from under the assumption that ffiffiffiffiffi g 0 p a0; 2,4; 8, . ... The value of Dð0Þ can be approximated by the determinant of a central k Â k block as will be shown in Section 5.3.

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 D 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.

Stability boundaries
The Mathieu equation ( 3) 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 _ y 1 ðtÞ is odd and _ y 2 ðtÞ is even.Using the transition property (12)  (see [17]).Hence, it holds that D ¼ traceðU T Þ ¼ 2y 1 ðpÞ.
The stability of the unconstrained Mathieu equation (3) depends on the value of D, being a function of the parameters a and b.The stability boundaries in the parameter plane ða,bÞ are given by Dða,bÞ ¼ 7 2, i.e. 9y 1 ðpÞ9 ¼ 1.The unity of the determinant, detðU T Þ ¼ y 1 ðpÞ 2 Ày 2 ðpÞ _ y 1 ðpÞ ¼ 1, implies that either y 2 ðpÞ ¼ 0 and/or _ y 2 ðpÞ ¼ 0 at a stability boundary.In other words, one can distinguish between stability boundaries for which y 2 ðpÞ ¼ 0 and stability boundaries for which _ y 2 ðpÞ ¼ 0. The value of the discriminant Dða,bÞ and the number nða,bÞ, i.e. the number of zeros of the function y 2 ðtÞ on the interval ð0,p, have been computed using direct numerical integration on a grid of 1000 Â 1000 points for the intervals a ¼ À8 . . .The stability of the equilibrium of the unilaterally constrained Mathieu equation is dependent on the number nða,bÞ and the discriminant Dða,bÞ, which both depend on the system parameters a and b, and the restitution coefficient e.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 e 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 e 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 40, 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 e c ¼ 0) and D42.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 D, which characterize the unconstrained Mathieu equation, one can use common approximation methods to investigate the stability properties of the unilaterally constrained Mathieu equation with e ¼ 0.

Averaging
A standard averaging technique [28] can be used to give an approximation for Dða,bÞ if b=a 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 dynamics in comoving coordinates [27] which yields a linear set of differential equations allowing for a closed form solution of the averaged equations.
Let o ¼ 1; 2,3, . . .be a resonant frequency of the Mathieu equation ( 3) and let b ¼ Eo 2 where E is a small parameter measuring the relative intensity of the parametric excitation.
Due to symmetry we know that Dða,ÀbÞ ¼ Dða,bÞ and it therefore suffices to consider EZ0.Furthermore, we consider the value of a to be close to resonance and set where d is a detuning parameter.Following [27] Substitution of ( 64) and ( 65) in (3) yields the Mathieu equation in comoving coordinates: The linear planar system (68) has the general solution 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 Near the first resonance, it therefore holds that DoÀ2 if 9d9 o1, or, correspondingly, if 91Àa9 o 9b9.From (72) we calculate the largest characteristic multiplier 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 Clearly, if E ¼ 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 ðpÞ gives  which is slightly larger than zero for 9d9 o1.We therefore infer that n¼0 holds around the resonance frequency o ¼ 1 for small values of E.
The stability criterion for n¼0 and DoÀ2 reads as (see Theorem 2) Using the approximation (74) for l 1 , the critical coefficient of restitution is approximated near the first resonance by e c ¼ e ÀðEp=2Þ ffiffiffiffiffiffiffiffi ffi Inversely, for a given value of a and e one can calculate the critical from which we see that a small value of the restitution coefficient e 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 e c , which has been chosen to be 0.4, 0.6 and 0.8, the value of b has been computed using (79) and is shown as dashed lines in Fig. 7.The approximation agrees fairly well with the numerical results for e c ¼ 0:8.Significant differences can be seen for e c ¼ 0:4 and e c ¼ 0:6 because E ¼ b can no longer considered to be small in the upper half of Fig. 7.

Approximation using Hill's determinant
A much better approximation of the discriminant Dða,bÞ and the critical restitution coefficient can be obtained by using Hill's infinite determinant as discussed in Section 4.
The function gðtÞ ¼ aþ2b cos 2t of the Mathieu equation can be represented by gðtÞ ¼ aþbðe 2it þ e À2it Þ and the Fourier coefficients are therefore g 0 ¼ a and g 1 ¼ g À1 ¼ b whereas all other Fourier coefficients are zero.The determinant Dð0Þ, see (59), therefore reads as The value of Dð0Þ can be approximated by the determinant of the central 5 Â 5 block Hence, the critical restitution coefficient is determined by the equation or equivalently by (62).Inversely, for a given value of a and e one can calculate an approximation of the critical value of b, based on the terms in (82) up to order b 2 , as This approximation is valid (for small b) 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 a and b, which cannot (easily) be done by using the method of Hill's determinant.

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 D 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 e.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 Poincare ´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 ¼ pþOðEÞ.In Section 5.2, the averaging method in comoving coordinates is used to obtain an approximation of the two-dimensional Poincare ´map.This approximate Poincare ´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 E 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 nonzero position, which does not agree with the equilibrium of the unconstrained system, i.e. xðtÞ Z x c 40, then accumulation points are possible.The unilaterally constrained Hill's equation with x c 40 does not have an equilibrium and the contact force l 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 Poincare ´map, is lost.
The present paper gives more insight in the stability properties of Hill's equation with unilateral constraint, being an archetype 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-offreedom autoparametric systems with unilateral constraints.
The fractional part can be expressed in trigonometric functions as where cotðpkÞ ¼ þ1 for k A Z.
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 kp with k A Z and the number of zeros on the interval ð0,a therefore amounts to The function f 2 ðxÞ ¼ sinðx þ bÞ has the number of zeros

Proposition 2 . 1 p Z p 0 dt y 1 ðtÞ 2 þ y 2 ðtÞ 2
Let n denote the number of zeros of y 2 ðtÞ on the interval ð0,p.It holds that n ¼

2 40 ð23Þ which we write as cotðcðpÞÞ 4
The zeros of y(t) are therefore determined by the condition cosðcðtÞÀy 0 Þ ¼ 0 or equivalently sinðcðtÞÀy 0 þp=2Þ ¼ 0. Using (91), the number of zeros of y(t) therefore amounts to gives together with(19) the inequality n r m rn þ 1, i.e. m ¼n or m ¼ n þ1.From property (92) follows that m¼n if and only if cotðcðpÞÞ þ cot Ày 0 þ p Àtanðy 0 Þ, i.e. y 0 4 y c ¼ Àarctanðcot ðcðpÞÞÞ where cotðcðpÞÞ ¼ y 1 ðpÞ=y 2 ðpÞ.If the inequality does not hold, then m a n and m must therefore be equal to n þ 1. &

Fig. 2 .
Fig.2.Solution x(t) of the unilaterally constrained Hill's equation and y(t) of Hill's equation.
35Þin which the transition property(12) of the fundamental solution matrix and the notation U T ¼ Uðp,0Þ for the monodromy matrix has been used.Let yðtÞ denote the solution of Hill's equation with yð0Þ ¼ xð0Þ.From (35) it becomes apparent that one can relate the solution xðpÞ of the unilaterally constrained Hill's equation to yðpÞ and m xðpÞ ¼ ðÀeÞ m U T yð0Þ ¼ ðÀeÞ m yðpÞ, ð36Þ which is the same result as (30) but now in first-order form.The number m depends on the angle y 0 of the initial condition xð0Þ and m¼n if y 0 4 y c and m ¼ n þ 1 else, see Proposition 3. The domain D, defined in (26), is therefore split into two cones D 1 and D 2 D 1 ¼ fx A D9y 4 y c g,

Fig. 3 .
Fig. 3. Domain D being in the cones D 1 and D 2 by S.

Proposition 6
discusses the real case v 1;2 A D 1 and Proposition 7 the real case v 1;2 A D 2 using a slightly different Lyapunov function.Proposition 8 deals with the complex conjugated case 9D9 o 2 whereas the cases D ¼ 7 2 are discussed in Propositions 9-11.

Theorem 2 . 9 9
Let D be the discriminant and l 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,p.The equilibrium x n ¼ 0 of the unilaterally constrained Hill's equation (4) is globally uniformly asymptotically stable if 0 re o e c , where the critical restitution coefficient is given by e c ¼ 0 if n ¼ 0 and DZ2 9l 1 À1=n if n 40 and ðÀ1Þ n D42, 9l 1 À1=ðn þ 1Þ if n Z 0 and ðÀ1Þ n DoÀ2, 1 if n 40 and À2 r Dr2 or if n ¼ 0 and À2 r Do2: cos 2tÞyðtÞ ¼ 0 is a Hill's equation with symmetry of the function gðtÞ ¼ aþ2b cos 2t, i.e. gðtÞ ¼ gðÀtÞ, gðtÞ ¼ gðt þ pÞ: 32 and b ¼ 0 . . .12. The stability boundaries Dða,bÞ ¼ 72 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 ða,bÞ if y 2 ðpÞ 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 ðpÞ ¼ 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.

Fig. 4 .
Fig. 4. Ince-Strutt diagram and discriminant D of the unconstrained Mathieu equation (3).In the grey stability domains holds 9D9 o 2. The instability domains are white.

y 1 ðtÞ ¼ cosh mt cos tÀ ffiffiffiffiffiffiffiffiffiffiffi ffi 1Àd 2 p 2 ffiffiffiffiffiffiffiffiffiffiffi ffi 1Àd 2 q 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b 2
1þ d sinh mt sin t ð71Þ for the first fundamental solution and the approximation D ¼ 2y 1 ðpÞ ¼ À2 cosh Ep ð72Þ for the discriminant.Using E ¼ b and Ed ¼ 1Àa for o ¼ 1, we can express the determinant as a function of a and b Dða,bÞ ¼ À2 cosh p Àð1ÀaÞ 2 q :

Fig. 6 .
Fig. 6.Ince-Strutt diagram of the unilaterally constrained Mathieu equation with critical value ec of the restitution coefficient.

Fig. 5 .
Fig. 5. Diagram with the value of n of the Mathieu equation.

2 ðaÀ16Þ 2 b 4 :aðaÀ4Þ 2 ðaÀ16Þ 2 b 4 !
ð81Þ This approximation can be improved by calculating larger central k Â k blocks in (80).However, the coefficients of b 2 and b 4 in (81) slightly change if larger central blocks are considered.Using (61) the discriminant D can be approximated by Dða,bÞ ¼ 2À4 sin 2 :
Appendix A. Floor function and fractional part Definition 1 (Floor function and fractional part).With bxc we denote the floor function defined by bxc ¼ fmax k A Z9krxg ð 85Þ and with fxg the fractional part defined by fxg ¼ xÀbxc: 19 Z 9l 2 9.There exist two linearly independent real eigenvectors v 1 and v 2 .If y 2 ðpÞ a0, then the eigenvectors are given by Proposition 6 (Stability for the real case v 1;2 A D 1 ).If v 1;2 A D 1 then the equilibrium x n ¼ 0 of the unilaterally constrained Hill's equation is globally uniformly asymptotically stable if e n 9l 1 9 o1 and unstable if e n 9l 1 9 41.Proof.If v 1;2 A D 1 then the characteristic multipliers are real and distinct and the decomposition (43) exists.Consider the discrete Lyapunov candidate function we deduce that UðÀp,0Þ ¼ Uð0,pÞ ¼ Uðp,0Þ À1 .Evaluation of UðÀp,0Þ ¼ Uðp,0Þ À1 the identity y 1 ðpÞ ¼ _ cos 2tÞðz 1 cos ot þz 2 sin otÞ sin ot, _ z 2 ¼ EoðdÀ2 cos 2tÞðz 1 cos ot þz 2 sin otÞ cos ot:ð66ÞThe system (66) can be averaged over one period of oscillation, keeping z 1 and z 2 constant, i.e.
The floor function has the follow-