# Existence of periodic orbits in grazing bifurcations of impacting mechanical oscillators

Arne Nordmark

► To cite this version:

Arne Nordmark. Existence of periodic orbits in grazing bifurcations of impacting mechanical oscillators. Nonlinearity, 2001, 10.1088/0951-7715/14/6/306. hal-01297283

HAL Id: hal-01297283

https://hal.science/hal-01297283

Submitted on 4 Apr 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>
# Existence of periodic orbits in grazing bifurcations of impacting mechanical oscillators

Arne B Nordmark

Department of Mechanics, Royal Institute of Technology, SE-100 44 Stockholm, Sweden

## Abstract

Grazing bifurcations are local bifurcations that can occur in dynamical models of impacting mechanical systems. The motion resulting from a grazing bifurcation can be complex. In this paper we discuss the creation of periodic orbits associated with grazing bifurcations, and we give sufficient conditions for the existence of a such a family of orbits. We also give a numerical example for an impacting system with one degree of freedom.

## 1. Introduction

When modelling impacting mechanical systems, one common approach is to use an impact law where velocities change discontinuously at the moment of impact. One effect of using such a model is that a trajectory that undergoes a low-velocity impact becomes sensitive to changes in initial conditions. The sensitivity is inversely proportional to the impact velocity. We will also find a new type of local bifurcation, the *grazing bifurcation* of a periodic orbit, see Nordmark [1], Foale and Bishop [2], Budd and Dux [3]. These bifurcations happen when a periodic orbit, with zero or more impacts during the period, is displaced by a parameter change in such a way that it encounters a new impact, which then will take place with zero impact velocity (a *grazing impact*). Grazing bifurcations are associated with rich and complex dynamical behaviour, such as the sudden loss of stability or existence of the orbit, the creation of a large number of periodic orbits, and the possibility of localized attracting chaotic motion. They are also typical bifurcations under the change of a single bifurcation parameter in these impact models, and for these models they are thus just as important to understand as other typical bifurcations of smooth dynamical systems.

Of the rich set of phenomena involved in grazing bifurcations, here we will study one particular aspect: the creation of a set of periodic orbits branching off from the grazing bifurcation point. Several of these orbits have been observed in various particular systems.
Foale and Bishop [4] studied linear impact oscillators with one degree of freedom. Ivanov [5] has studied a particular two-parameter grazing bifurcation for an equivalent system. Chin *et al* [6,7] have made a thorough study of maximal periodic orbits of a two-dimensional (2D) mapping, which is similar to the Poincaré mappings found for single-degree-of-freedom impact oscillators.

In this paper, we present a general scheme for finding the periodic orbits that branch off
in a grazing bifurcation. The analysis can be applied to systems with one or several degrees of
freedom, provided that an impact law is used. We give a set of sufficient conditions involving
only linear equations and linear inequalities for the existence of a periodic orbit of a given
impact pattern.

The paper is divided as follows. In section 2 we state the form of the Poincaré mapping near a grazing periodic orbit. From this we calculate periodic orbits and their stability. Section 3 is a parameter study covering two-dimensional Poincaré mappings. Selected results for a linear driven impact oscillator are shown in section 4.

## 2. Grazing bifurcation of a periodic orbit

We will study the dynamics of a class of mappings that have a particular form containing a square root term. These mappings are interesting since they occur as local Poincaré mappings for mechanical systems of the type considered by Fredriksson and Nordmark [8]. There, the mechanical system is assumed to have the following properties:

* there is a condition for the occurrence of impact that is a smooth surface of codimension one in state space;

* the equations of motion are also smooth slightly beyond the impact surface;

* there is a smooth impact law mapping that takes the system to a new state when reaching the impact surface, and this impact law becomes the identity mapping as the impact velocity approaches zero.

We will assume that our mechanical system has these characteristics.

### 2.1. Form of the local Poincaré mapping

Suppose now that our system has a periodic orbit with a single grazing impact. We may well allow for an additional finite number of non-grazing impacts. To study the stability and bifurcations of this orbit we introduce a transversal local Poincaré section through a point on the trajectory that is not on the impact surface, so the periodic orbit becomes a fixed point of the Poincaré mapping. A section coordinate $x \in \mathbb{R}^N$ and a bifurcation parameter $\mu \in \mathbb{R}^M$ are introduced in such a way that the fixed point is at $x = 0$ when $\mu = 0$ and that $\mu = 0$ is the parameter value for which this orbit is grazing.

From the assumed characteristics of the system it follows that non-grazing impact does not lead to any discontinuity or non-smoothness in the Poincaré mapping, but a grazing impact introduces a non-differentiability (but not a discontinuity). Also, trajectories starting in the Poincaré section close to the grazing trajectory may encounter a low-velocity impact, a grazing impact or no impact near the point where there is a grazing impact of the original trajectory. On the other hand, they will always encounter non-grazing impacts close to the points where the original trajectory had such impacts. This leads us to the idea that the possible low-velocity impact should be treated separately from the rest of the motion. Following Fredriksson and Nordmark [8] (where the results stated in the following are derived), we can study a modified system where impacts with impact velocities below a threshold (safely lower than any of the
non-grazing impacts of the original trajectory) are ignored and do not imply a state change according to the impact law. In this case the trajectories will penetrate a little bit beyond the impact surface, but above we have assumed that the equations of motion are still valid and smooth here. This modified system has a smooth local Poincaré mapping and we will view the effects of the real, low-velocity impact as a modification of this system.

Making things more precise, the Poincaré mapping near $x = 0, \mu = 0$ of the modified system (where low-velocity impacts are ignored) is smooth and described by the function

$$f(x, \mu): \mathbb{R}^N \times \mathbb{R}^M \to \mathbb{R}^N$$

with $f(0, 0) = 0$. For later reference we define

$$A = D_x f(0, 0). \qquad (1)$$

Next, we incorporate the effect of the possible low-velocity impact as a correction to $f$. This is done by writing the Poincaré mapping of the real system as $g \circ f$, with

$$g(x, \mu): \mathbb{R}^N \times \mathbb{R}^M \to \mathbb{R}^N.$$

Points near $x = 0$ may or may not undergo a low-velocity impact. We introduce a smooth function $h$ that has the value of zero on a set coinciding with the set of points leading to grazing impacts, has negative values for points leading to low-velocity impacts and has positive values for points that do not impact (with low velocity). The function $h$ should have the properties

$$h(x, \mu): \mathbb{R}^N \times \mathbb{R}^M \to \mathbb{R}$$

with $h(0, 0) = 0$. We write

$$C = D_x h(0, 0) \qquad (2)$$

and we should chose $h$ such that $C \neq 0$.

A non-negative value of $h$ means no low-velocity impact, so then $g$ should be the identity mapping. The function $g$ takes the particular form

$$g(x, \mu) = \begin{cases} x & \text{if } h(x, \mu) \ge 0 \\ b(x, y, \mu)y + x & \text{if } h(x, \mu) \le 0 \end{cases} \qquad (3)$$

where

$$y = \sqrt{-h(x, \mu)} \qquad (4)$$

and $b$ is a smooth function

$$b(x, y, \mu): \mathbb{R}^N \times \mathbb{R} \times \mathbb{R}^M \to \mathbb{R}^N.$$

We write

$$B = b(0, 0, 0) \qquad (5)$$

and assume $B \neq 0$.

We see that $g$ is a continuous mapping, but $D_x g$ is unbounded as $h(x, \mu)$ approaches zero from the negative side. Still, the state space volume must not increase if the system is conservative or dissipative, and, in particular, it must remain bounded, so $\det(D_x g)$ should be bounded, which implies $CB = 0$.

When $\mu = 0$, there is a fixed point at $x = 0$ of the Poincaré mapping $g \circ f$. We now want to study the dynamics of this mapping near $x = 0, \mu = 0$, and more specifically the possibility of other periodic orbits.

For an example of these mappings for a particular system, see section 3.
## 2.2. Changing the Poincaré section

We will briefly discuss the effect of the particular choice of Poincaré section. Let $p(x, \mu)$ be a mapping that maps points from the original section to a new one according to the flow (in the flow direction that avoids the part of the trajectory with possible low-velocity impacts). This will be a smooth mapping, and we will assume that it is invertible with respect to the first argument with inverse $p^{-1}$. Denoting quantities in the new section with a hat, we find

$$ \hat{x} = p(x, \mu) \qquad (6) $$

$$ \hat{f}(\hat{x}, \mu) = p(f(p^{-1}(\hat{x}, \mu), \mu), \mu) \qquad (7) $$

$$ \hat{g}(\hat{x}, \mu) = p(g(p^{-1}(\hat{x}, \mu), \mu), \mu) \qquad (8) $$

$$ \hat{h}(\hat{x}, \mu) = h(p^{-1}(\hat{x}, \mu), \mu) \qquad (9) $$

$$ \hat{b}(\hat{x}, y, \mu) = p(b(p^{-1}(\hat{x}, \mu), y, \mu), \mu). \qquad (10) $$

Furthermore, defining

$$ P = D_x p(0, 0) \qquad (11) $$

we find

$$ \hat{A} = PAP^{-1} \qquad (12) $$

$$ \hat{B} = PB \qquad (13) $$

$$ \hat{C} = CP^{-1} \qquad (14) $$

which means that quantities such as the eigenvalues of $A$ and the values of $CA^k B$ that are used in the following are invariant with respect to changes of the Poincaré section.

By the same reasoning, these results also apply to smooth changes of the coordinate system on the same Poincaré section.

## 2.3. Equations for periodic orbits

A periodic point $\bar{x}$ of period $n$ should satisfy

$$ \bar{x} = (g \circ f)^n(\bar{x}, \mu) \qquad (15) $$

but this equation is awkward to work with because of the two different expressions of $g$, and because of the derivative singularity at $h = 0$. Instead we will try to find periodic orbits with a fixed pattern of visits in $h > 0$ and $h < 0$, respectively. Using a terminology where only low-velocity impacts are important, we will say that a point with $h < 0$ will lead to an impact, but a point with $h > 0$ will not.

We will first classify periodic points as being impacting or not, depending on whether some point in the orbit will lead to an impact or not. A non-impacting periodic point satisfies the equation

$$ \bar{x} - f^n(\bar{x}, \mu) = 0 \qquad (16) $$

with the restrictions

$$ h(f^k(\bar{x}, \mu), \mu) > 0 \quad \text{for } 1 \le k \le n. \qquad (17) $$

Linear stability is given by

$$ (D_x(g \circ f)^n)(\bar{x}, \mu) = (D_x f^n)(\bar{x}, \mu). \qquad (18) $$
A way to characterize an impacting periodic point is by its impact sequence $(n_1, n_2, \dots, n_m)$, which denotes the situation of $n_1 - 1$ non-impacting iterations followed by one impacting, followed by $n_2 - 1$ non-impacting iterations followed by one impacting, and so on. The period is $n = \sum n_i$. We will look for periodic points with a specified impact sequence. To do this we rewrite (15) as a system of equations by introducing $x_i$ for the x value immediately following impacting iteration number $i - 1$, (thus $\bar{x} = x_1$), and $y_i$ for the y value associated with impact iteration number $i$, and using all these $x_i$ and $y_i$ values as unknowns. Furthermore, equation (4) is squared to remove the singularity. This gives us the set

$$
\begin{align}
x_2 - b(f^{n_1}(x_1, \mu), y_1, \mu)y_1 - f^{n_1}(x_1, \mu) &= 0 \notag \\
h(f^{n_1}(x_1, \mu), \mu) + y_1^2 &= 0 \notag \\
\vdots \notag \\
x_1 - b(f^{n_m}(x_m, \mu), y_m, \mu)y_m - f^{n_m}(x_m, \mu) &= 0 \notag \\
h(f^{n_m}(x_m, \mu), \mu) + y_m^2 &= 0
\end{align}
$$

of equations to be solved for $x_i$ and $y_i$, subject to the restrictions

$$h(f^k(x_i, \mu), \mu) > 0 \quad \text{for} \quad 1 \le k \le n_i - 1 \quad (20)$$

and

$$y_i > 0. \quad (21)$$

We need (20) since we have assumed that these iterations are non-impacting, and (21) since we have squared (4).

To check the stability of a periodic orbit we also need the linearization of the Poincaré mapping. In terms of the $x_i$, $y_i$ and $n_i$, this matrix is

$$(D_x(g \circ f)^n)(\bar{x}) = \prod_{i=1}^{m} \left\{ \begin{aligned} & \left[I - \frac{b(f^{n_i}(x_i, \mu), y_i, \mu)}{y_i}\right. \\ & \left. \qquad + (D_y b)(f^{n_i}(x_i, \mu), y_i, \mu)\right] (D_x h)(f^{n_i}(x_i, \mu), \mu) / 2 \right] (D_x f^{n_i})(x_i, \mu) \\ & \left. + (D_x b)(f^{n_i}(x_i, \mu), y_i, \mu)y_i \right\}. \end{aligned} \right. \quad (22)$$

## 2.4. Single-parameter bifurcations in the generic case

In the following we will study grazing bifurcations under the change of a single parameter. Thus we take the parameter space dimension $M = 1$. We will assume that there is an isolated fixed point $\tilde{x}(\mu)$ of $f$ for small $\mu$ with $\tilde{x}(0) = 0$. This is guaranteed by the implicit function theorem if we have

$$\det(I - A) \neq 0 \quad (23)$$

and the point satisfies

$$\tilde{x}(\mu) = f(\tilde{x}(\mu), \mu). \quad (24)$$

We also assume that the fixed point crosses the surface $h = 0$ transversally when $\mu$ is changed, such that

$$h(\tilde{x}(\mu), \mu) = e\mu + \mathcal{O}(\mu^2) \quad (25)$$
where

$$e = D_{\mu}(h(\tilde{x}(\mu), \mu))(0) = C(I - A)^{-1}(D_{\mu}f)(0, 0) + (D_{\mu}h)(0, 0) \quad (26)$$

and $e \neq 0$. When $e\mu$ is small and positive, $\tilde{x}$ is also a non-impacting fixed point for the full Poincaré mapping $g \circ f$, since $g$ is the identity mapping for positive values of $h$. We will now consider the existence of impacting periodic points near $x = 0$, $\mu = 0$ with $\mu \neq 0$.

Defining

$$z = \begin{bmatrix} x_1 \\ y_1 \\ \vdots \\ x_m \\ y_m \end{bmatrix} \quad (27)$$

we can summarize the set of equations (19) as $F(z, \mu) = 0$ with the known particular solution $F(0, 0) = 0$, and for small $\mu$ this defines $z(\mu)$ with $z(0) = 0$ uniquely by the implicit function theorem unless $(D_z F)(0, 0)$ is singular. This $z(\mu)$ will of course also have to be checked for compliance with (20) and (21).

We use the notation $D_{(n_1,...,n_m)}$ for $(D_z F)(0, 0)$. This matrix has a banded structure. For example, with $m = 3$ we have

$$D_{(n_1,n_2,n_3)} = \begin{bmatrix} -A^{n_1} & -B & I & 0 & 0 & 0 \\ CA^{n_1} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -A^{n_2} & -B & I & 0 \\ 0 & 0 & CA^{n_2} & 0 & 0 & 0 \\ I & 0 & 0 & 0 & -A^{n_3} & -B \\ 0 & 0 & 0 & 0 & CA^{n_3} & 0 \end{bmatrix}. \quad (28)$$

We can make the form of (19)–(21) a bit clearer by introducing rescaled versions $X$ and $Y$ of the offsets from the fixed point: $x = \tilde{x}(\mu) + e\mu X$, $y = e\mu Y$. Then we find using (24)

$$ \begin{aligned} f^k(x, \mu) &= f^k(\tilde{x}(\mu) + e\mu X, \mu) \\ &= f^k(\tilde{x}(\mu), \mu) + e\mu A^k X + \mathcal{O}(\mu^2) = \tilde{x}(\mu) + e\mu A^k X + \mathcal{O}(\mu^2) \end{aligned} \quad (29) $$

and using (25)

$$ \begin{aligned} h(f^k(x, \mu), \mu) &= h(\tilde{x}(\mu) + e\mu A^k X + \mathcal{O}(\mu^2), \mu) \\ &= h(\tilde{x}(\mu), \mu) + e\mu C A^k X + \mathcal{O}(\mu^2) = e\mu + e\mu C A^k X + \mathcal{O}(\mu^2). \end{aligned} \quad (30) $$

The first and second lines of (19) become

$$e\mu(X_2 - BY_1 - A^{n_1}X_1 + \mathcal{O}(\mu)) = 0$$

$$e\mu(1 + CA^{n_1}X_1 + \mathcal{O}(\mu)) = 0$$

and making $Z$ the vector

$$Z = \begin{bmatrix} X_1 \\ Y_1 \\ \vdots \\ X_m \\ Y_m \end{bmatrix} \quad (31)$$
the full system after dividing by $e\mu$ is

$$D_{(n_1, \dots, n_m)} Z(\mu) = \begin{bmatrix} 0 \\ -1 \\ \vdots \\ 0 \\ -1 \end{bmatrix} + \mathcal{O}(\mu) \quad (32)$$

with conditions

$$\operatorname{sign}(e\mu)(CA^k X_i(\mu) + 1) + \mathcal{O}(\mu) > 0 \quad \text{for} \quad 1 \le k \le n_i - 1 \quad (33)$$

$$\operatorname{sign}(e\mu)Y_i(\mu) > 0. \quad (34)$$

These equations have well defined one-sided limits as $\mu \to 0^+$ or $\mu \to 0^-$. The result can be summarized as follows.

**Theorem 1.** Suppose $\det(I - A) \ne 0$, $e \ne 0$, $\det(D_{(n_1, \dots, n_m)}) \ne 0$, and let $Z(0)$ be the solution of (32). If (33) and (34) are fulfilled for $\mu = 0^+$ or $\mu = 0^-$, then there exists a unique periodic orbit with impact pattern $(n_1, \dots, n_m)$ for all small $\mu$ of that sign, that approaches $x = 0$ as $\mu \to 0$. If any of the quantities in (33) and (34) are negative for $\mu = 0^+$ or $\mu = 0^-$, then no periodic orbit with impact pattern $(n_1, \dots, n_m)$ can be continued to $x = 0$ for small $\mu$ of that sign.

Of course, equations (33) and (34) can be satisfied for at most one sign of small $\mu$, so all branches are one-sided. Uniqueness tells us that there exists at most one branch per impact pattern.

Returning to the question of stability, equation (22) for these solutions has a leading term that is $\mathcal{O}(\mu^{-m})$:

$$(D_x(g \circ f)^n)(\bar{x}) = \prod_{i=1}^{m} \left[ -BCA^{n_i}/(2y_i) \right] + \mathcal{O}(\mu^{1-m}) \quad (35)$$

(with matrix multiplication from the left) since $y_i = \mathcal{O}(\mu)$ as $\mu \to 0$. In particular, if all $CA^{n_i}B \ne 0$, then there is an eigenvector near $B$ with eigenvalue near

$$\prod_{i=1}^{m} \left[ -CA^{n_i}B/(2y_i) \right] \quad (36)$$

which is also $\mathcal{O}(\mu^{-m})$ as $\mu \to 0$. Thus in general, these impacting periodic solutions have eigenvalues that approach infinity as $\mu \to 0$, and consequently they are not expected to be stable, but stability cannot be ruled out for non-generic cases.

## 2.5. Comments on non-generic cases

If theorem 1 applies, then the results are unchanged for small enough perturbations of the system (at least for a fixed impact pattern). We have put emphasis on a system where symmetries or similar phenomena do not put any important constraints on the form of $f$, $b$ or $h$. Then we can always perturb any bifurcation situation into one where theorem 1 applies, using arbitrarily small perturbations. This could, in principle, be accomplished using more than one bifurcation parameter.

There are a few different ways in which the bifurcation situation is not fully covered by theorem 1. If the non-impacting fixed point is not hyperbolic, we must also investigate ordinary
smooth bifurcations of the point. If $e = 0$ it is not guaranteed that the fixed point of $f$ crosses the line $h = 0$, and it is a sign of the degeneration of the choice of bifurcation parameter.

A more interesting situation is when $\det(D_{(n_1, \dots, n_m)}) = 0$. The implicit function theorem does not apply, and there may be several branches of (19). The branches may be found using standard techniques of smooth bifurcation theory, and some of them are probably pruned by (20) and (21). The branches may also be stable. An example of this is examined in section 4.

Another interesting situation is when some the quantities of (33) and (34) are zero and the rest positive. Higher-order terms are needed to determine existence, and a zero in (33) means that the bifurcating orbit is itself close to grazing.

Lastly, we have studied the situation for a fixed impact pattern. If there are an infinite number of orbits of different impact patterns for some sign of $\mu$, the intersection of the $\mu$ intervals may well be zero. In this case, we cannot be certain that more than a finite number of the orbits exists for given small non-zero $\mu$.

### 3. Analysis when $N = 2$

A class of models showing grazing bifurcations that has been the subject of much recent study is periodically forced oscillators with one degree of freedom. Here we present a general model.

The state of such a system is described by a position $q$, a velocity $u$ and the time $t$, where $t$ is only significant modulo the driving period $T$. Impacting is specified by position constraints $Q_1(t) \le q \le Q_2(t)$ and by impact laws, giving the velocity after impact $u_+$ from the velocity before impact $u_-$ at the two boundaries: $u_+ = U_i(u_-, t)$, $i = 1, 2$. The dynamics when not impacting is specified by an acceleration function $a$:

$$ \dot{q} = u \qquad (37) $$

$$ \dot{u} = a(q, u, t). \qquad (38) $$

All of $Q_i$, $U_i$ and $a$ are periodic in $t$ with period $T$. To keep the system within bounds, we require that $u_+ > \dot{Q}_1(t)$ if $u_- < \dot{Q}_1(t)$ (and $u_+ < \dot{Q}_2(t)$ if $u_- > \dot{Q}_2(t)$). To ensure that zero relative incoming velocity is mapped to zero relative outgoing velocity we require $u_+ = u_-$ if $u_- = \dot{Q}_i(t)$.

#### 3.1. Poincaré mappings near a grazing periodic orbit

We now assume that there is a periodic orbit with a single grazing impact. Furthermore, it is assumed that it occurs at the lower bound and at $t = t_G$. At the grazing impact we have $q = Q_1(t_G)$, $u = \dot{Q}_1(t_G)$, and we define the relative acceleration at grazing $a_G = a(Q_1(t_G), \dot{Q}_1(t_G), t_G) - \dot{Q}_1(t_G)$, with $a_G > 0$.

We now introduce a Poincaré section as a plane of constant forcing phase (which is also constant time modulo the period $T$). The simplest form of the function $g$ would be obtained as we let the section phase approach $t_G$ from above, but we have assumed that the section must not lie at an impact so we cannot use $t = t_G$ as a Poincaré section. Instead we will use a section at an arbitrary time a bit later than $t_G$ (but well before any other impact) and then transform back to $t_G$ using the flow of (37) and (38) without any impacts. The resulting mapping is almost the Poincaré mapping at $t_G$, except that all low-velocity impacts have taken place, regardless of whether they occur just before or just after $t_G$.
<figure><img src="image_2.png" /><figcaption>Figure 1. The definition of the mapping g.</figcaption></figure>

In this section we use coordinates

$$x = \begin{bmatrix} q - Q_1(t_G) \\ u - \dot{Q}_1(t_G) \end{bmatrix}. \qquad (39)$$

The dynamics without the possible low-velocity impact near $t = t_G$ is a mapping $f(x)$ which is smooth, since all other impacts of the periodic orbit are assumed to occur at non-zero velocity (Nordmark [1]). Still disregarding the possibility of a low-velocity impact, any trajectory of (37) and (38) passing through the Poincaré section at the point $x$ close to 0 defines a local minimum value $h(x)$ of the function $q(t) - Q_1(t)$, since $a_G > 0$. When $h(x) < 0$, the trajectory has two intersections with $q = Q_1(t)$, so let $t_C(x)$ denote the earlier of these time points.

Now we can define the mapping $g$. For $h(x) \ge 0$ we let $g(x) = x$. For $h(x) < 0$ we start at $x$ with $t = t_G$, integrate (37) and (38) forwards or backwards in time until $t = t_C$, use the impact law $U_1$, and integrate back to $t = t_G$ again. This resulting point in the Poincaré section defines $g(x)$ if $h(x) < 0$.

The effect of the mapping $g$ is shown in figure 1, where the $x$ coordinate is not shown. The incoming and outgoing trajectories are broken, and we see how following the flow to $t_C$, applying $U_1$ and following the flow from $t_C$ forwards, is equivalent to following the flow to $t_G$, applying $g$, and following the flow from $t_G$ forwards.

Let $\tau = -(D_u U_1)(\dot{Q}_1(t_G), t_G) \ge 0$. We then find

$$h(x) = x_1 + \mathcal{O}(x)^2 \qquad (40)$$

$$g(x) = x + \begin{bmatrix} 0 \\ 1 \end{bmatrix} \sqrt{2a_G(1 + \tau)y + \mathcal{O}(x, y)^2} \quad \text{if } h(x) < 0 \qquad (41)$$

where $x_1$ is the first component of $x$ and $y = \sqrt{-h(x)}$. Using the notation of (2) and (5), we
have

$$C = \begin{bmatrix} 1 & 0 \end{bmatrix} \qquad (42)$$

$$B = \begin{bmatrix} 0 \\ 1 \end{bmatrix} \sqrt{2a_G(1 + \tau)}. \qquad (43)$$

We also set

$$A = (D_x f)(0) = \begin{bmatrix} \alpha & \gamma \\ \beta & \delta \end{bmatrix}. \qquad (44)$$

In a real system, the somewhat artificial Poincaré section used here cannot be directly observed. To use the theory to interpret data from a simulation with a real Poincaré section that is away from the grazing phase, the formulae (7)–(14) should be used.

## 3.2. Coordinate transformation

A parameter study of the existence of bifurcating periodic orbits is simplified by using transformed coordinates. To begin with, if $\gamma = 0$, all the $A^{n_i}$ are lower triangular, and $A^{n_i} B$ is parallel to $B$, and we can easily construct null vectors for $D_{(n_1, ..., n_m)}$. Since this matrix was assumed to be non-singular, we must have $\gamma \neq 0$. Then $A, B$ and $C$ can be transformed into (Nordmark [1])

$$A = \begin{bmatrix} a_1 & 1 \\ -a_2 & 0 \end{bmatrix} \qquad (45)$$

$$B = \begin{bmatrix} 0 \\ 1 \end{bmatrix} \sqrt{2a_G(1 + \tau)\gamma} \qquad (46)$$

$$C = \begin{bmatrix} 1 & 0 \end{bmatrix} \qquad (47)$$

where $a_1 = \alpha + \delta$ is the trace of $A$, and $a_2 = \alpha\delta - \beta\gamma$ is the determinant. $a_2 > 0$ since it is the determinant of the Poincaré mapping of a flow. For a parameter study, we may also assume that $a_2 \le 1$, since otherwise we can let time run backwards, which makes $a_1 \leftarrow a_1/a_2$ and $a_2 \leftarrow 1/a_2$. We also note that $\gamma$ only influences the $Y_i$ values, and that the signs of the $Y_i$ are what matters.

All in all, the relevant parameters in the case $N = 2$ are then $a_1, a_2, \text{sign}(\gamma)$ and $\text{sign}(e\mu)$. It should be noted that the sign of $\gamma$ is the sign of $CAB$.

It should be remarked that this coordinate transformation is made to simplify the parameter study, and is not a very useful tool when studying a particular system. For example, if $\gamma$ is very small, the transformation is almost singular and large deformations are introduced. This case is more properly studied as a two-parameter bifurcation.

## 3.3. Single-impact bifurcating orbits

First, we study orbits with the simple impact pattern (n). These orbits are termed maximal orbits in Chin et al [6], and since there is only one low-velocity impact in the period, one can hope that these orbits are those most likely to become stable as we move away from the
bifurcation point. Equations (32)–(34) become (in the limit $\mu \to 0$, with $X = X_1$ and $Y = Y_1$)

$$ \begin{bmatrix} I - A^n & -B \\ CA^n & 0 \end{bmatrix} \begin{bmatrix} X \\ Y \end{bmatrix} = \begin{bmatrix} 0 \\ -1 \end{bmatrix} \quad (48) $$

$$ \operatorname{sign}(e\mu)(CA^k X + 1) > 0 \quad \text{for} \quad 1 \le k \le n-1 \quad (49) $$

$$ \operatorname{sign}(e\mu)Y > 0. \quad (50) $$

We can eliminate $Y$ and simplify to obtain

$$ \begin{bmatrix} [1 & 0](I - A^n) \\ [1 & 0] \end{bmatrix} X = \begin{bmatrix} 0 \\ -1 \end{bmatrix} \quad (51) $$

$$ \operatorname{sign}(e\mu)b_{n,k} > 0 \quad \text{for} \quad 1 \le k \le n-1 \quad (52) $$

$$ \operatorname{sign}(e\mu) \operatorname{sign}(\gamma)c_n > 0. \quad (53) $$

We have introduced the notation

$$ b_{n,k} = \begin{bmatrix} 1 & 0 \end{bmatrix} A^k X + 1 \quad (54) $$

and

$$ c_n = \begin{bmatrix} 0 & 1 \end{bmatrix} (I - A^n) X \quad (55) $$

where $X$ is the solution of (51). To denote the different sign combinations of $e\mu$ and $\gamma$ we will use the notation (sign of $e\mu$, sign of $\gamma$). Thus $(+, -)$ means $e\mu > 0$ and $\gamma < 0$.

If we introduce the two eigenvalues of $A$: $\lambda_1$ and $\lambda_2$, we can write $a_1 = \lambda_1 + \lambda_2$ and $a_2 = \lambda_1\lambda_2$. If $(a_1/2)^2 \ne a_2$, then $\lambda_1$ and $\lambda_2$ are distinct, and we find

$$ A^n = \frac{1}{\lambda_1 - \lambda_2} \begin{bmatrix} \lambda_1^{n+1} - \lambda_2^{n+1} & \lambda_1^n - \lambda_2^n \\ -\lambda_1\lambda_2(\lambda_1^n - \lambda_2^n) & -\lambda_1\lambda_2(\lambda_1^{n-1} - \lambda_2^{n-1}) \end{bmatrix}. \quad (56) $$

The curve $(a_1/2)^2 = a_2$, where eigenvalues are equal, is, in fact, not special in this discussion of periodic orbits. All the relevant formulae have well defined limits when $\lambda_1 \to \lambda_2$, and no new phenomena are connected with this limit. Thus we will skip the case of equal eigenvalues in the following.

We can solve (51) to obtain

$$ X = \begin{bmatrix} -1 \\ X_2 \end{bmatrix} \quad (57) $$

where

$$ X_2 = \frac{(\lambda_1^{n+1} - \lambda_2^{n+1}) - (\lambda_1 - \lambda_2)}{\lambda_1^n - \lambda_2^n} \quad (58) $$

and compute

$$ b_{n,k} = \frac{-(\lambda_1\lambda_2)^k(\lambda_1^{n-k} - \lambda_2^{n-k}) - (\lambda_1^k - \lambda_2^k) + (\lambda_1^n - \lambda_2^n)}{\lambda_1^n - \lambda_2^n} \quad (59) $$

and

$$ c_n = \frac{-(\lambda_1 - \lambda_2)(\lambda_1^n - 1)(\lambda_2^n - 1)}{\lambda_1^n - \lambda_2^n}. \quad (60) $$

The asymptotic value of the eigenvalue of largest modulus is given by (36)

$$ \begin{aligned} \Lambda &= -CA^n B/2y + \mathcal{O}(1) \\ &= \frac{a_G(1+\tau)^2\gamma^2}{e\mu} \left( \frac{\lambda_1^n - \lambda_2^n}{\lambda_1 - \lambda_2} \right)^2 \frac{1}{(\lambda_1^n - 1)(\lambda_2^n - 1)} + \mathcal{O}(1). \end{aligned} \quad (61) $$
### 3.4. Existence boundaries of single-impact orbits

The existence of a bifurcating periodic orbit is determined by the values of $a_1, a_2$, sign($e\mu$), sign($\gamma$) and $n$. Thus the parameter set where single-impact orbits exist is a subset of $\mathbb{R} \times [0, 1] \times \{+, -\}^2 \times \mathbb{N}^+$. Visualizing this subset in meaningful ways presents a bit of a problem. For a particular bifurcation, we have $a_1, a_2$ and $\gamma$ given and want a table for the different combinations of sign($e\mu$) and $n$. Another way is to give a diagram in the $a_1-a_2$ plane for given $n$. We will use either of these methods as seems best fitted to the situation at hand.

We see that if $n = 1$, equation (52) does not come into play. Thus only sign($e\mu\gamma$) is important, so the $(+, +)$ and $(-, -)$ cases coincide, as do the $(+, -)$ and $(-, +)$ cases. For given $n > 1$, all four combinations of sign($e\mu$) and sign($\gamma$) lead to different inequalities, and an orbit exists for at most one of these sign combinations. In the $a_1-a_2$ plane, there will be a set of non-overlapping regions corresponding to different sign combinations. The boundaries have to do with the violation of one or more of the inequalities (52) and (53). There are three types of boundaries.

* $\lambda_1^n - \lambda_2^n = 0$ with $\lambda_1 - \lambda_2 \neq 0$: this is when the system matrix is singular, and will be called a singular boundary. It can only happen when $(a_1/2)^2 < a_2$, so eigenvalues must be complex.

* $b_{n,k} = 0$ for some $k$. This will be called an iterate boundary of order $k$.

* $c_n = 0$. This will be called a sign boundary.

It will be convenient to label regions of the $a_1-a_2$ plane (figure 2):

1. $a_2 + 1 < a_1$, where $0 < \lambda_2 < 1 < \lambda_1$;

2. $2\sqrt{a_2} < a_1 < a_2 + 1$, where $0 < \lambda_2 < \lambda_1 < 1$;

3. $-2\sqrt{a_2} < a_1 < 2\sqrt{a_2}$, where $\lambda_{1,2}$ are complex;

4. $-(a_2 + 1) < a_1 < -2\sqrt{a_2}$, where $-1 < \lambda_2 < \lambda_1 < 0$;

5. $a_1 < -(a_2 + 1)$, where $\lambda_2 < -1 < \lambda_1 < 0$.

There are three regions that are simple, since they do not contain any internal boundaries (see appendix A).

In region 1, both $b_{n,k}$ and $c_n$ are positive for all $k$ and $n$. Thus orbits of all periods exist in the $(+, +)$ case and additionally period 1 exists in the $(-, -)$ case. The longer orbits first move towards the origin along the stable eigenvector and then exit along the unstable one, and orbits are bounded as $n \to \infty$.

In region 2, both $b_{n,k}$ and $c_n$ are negative for all $k$ and $n$. Period 1 exists for $(+, -)$, and all periods exist for $(-, +)$. The longer orbits start far from the origin and come in along the eigenvector corresponding to $\lambda_1$. Since $X_2 \to -\infty$ exponentially as $n \to \infty$, the smallness assumptions are not fulfilled for the longer periods when $\mu$ has a non-zero value and we should not expect to find more than a finite number of them.

In region 5, we have $b_{n,k} > 0$ and $c_n < 0$. Period 1 exists for $(-, +)$ and all periods exist for $(+, -)$. As in region 1, longer periods tend to follow the stable and unstable eigenvectors and orbits are bounded.

The complete orbit when $n = 15$ is shown in figure 3 for the parameter values: $\lambda_1 = 1.5$, $\lambda_2 = 0.5$ (region 1), $\lambda_1 = 0.8$, $\lambda_2 = 0.5$ (region 2) and $\lambda_1 = -0.5$, $\lambda_2 = -1.5$ (region 5). The points plotted are $A^k X$ for $0 \le k \le 14$. The forbidden area for the iterates is shaded.
<figure><img src="image_3.png" /><figcaption>Figure 2. Regions in the $a_1-a_2$ plane.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. Single-impact orbits of period 15 for parameters in region 1 (top left), region 2 (top right) and region 5 (bottom).</figcaption></figure>

### 3.5. Boundary structure in regions 3 and 4

Regions 3 and 4 are much more complicated, each containing an infinite number of boundaries. Period 1 exists in the $(+, -)$ and $(-, +)$ cases. Period 2 exists in the $(+, +)$ case for $a_1 < 0$ and in the $(-, +)$ case for $a_1 > 0$. The line $a_1 = 0$ is a singular boundary.
When $n \ge 3$ iterate boundaries start to come into play. In region 3, where eigenvalues are complex, we will write $\lambda_{1,2} = re^{\pm i\theta}$ with $0 < \theta < \pi$ (giving $a_1 = 2r \cos(\theta)$, $a_2 = r^2$). Then

$$ \operatorname{sign}(b_{n,k}) = \operatorname{sign} \left( \frac{-r^n \sin((n-k)\theta) + r^{n-k} \sin(n\theta) - \sin(k\theta)}{\sin(n\theta)} \right) \quad (62) $$

$$ \operatorname{sign}(c_n) = -\operatorname{sign}(\sin(n\theta)). \quad (63) $$

The analysis is easiest when $r = 1$ ($a_2 = 1$). We can show:

**Theorem 2.** Let $n$ be an integer with $n \ge 3$, $r = 1$ and denote $s = \theta/(2\pi)$. Then period $n$ exists in the $(-, +)$ case for $0 < s < 1/(2n)$. Furthermore, let $m$ be an integer relative prime to $n$ with $0 < m/n < \frac{1}{2}$. Let $m_-, n_-, m_+$ and $n_+$ be the unique integers given by $mn_- - m_-n = -1$, $mn_+ - m_+n = 1$ and $0 < n_-, n_+ < n$. Then period $n$ exists in the $(+, -)$ case for $m/n < s < m_-/n_-$. For the $(+, +)$ case, period $n$ exists for $m_+/n_+ < s < m/n$, unless $m = 1$ when the interval $1/(2n) < s < 1/n$ instead. These intervals are the only places where period $n$ exists.

For proof, see appendix B. We note that for given $n \ge 3$, the number of intervals is $\phi(n) + 1$, where $\phi(n)$ is the value of the Euler $\phi$ function.

Take $n = 8$ as an example. We first have the $(-, +)$ interval $(0, \frac{1}{16})$. Integers $m$ relative prime to 8 and smaller than 4 are 1 and 3. With $m = 1$, we find $m_+/n_+ = 0/1$ and $m_-/n_- = 1/7$. This gives the $(+, -)$ interval $(\frac{1}{8}, \frac{1}{7})$ and the $(+, +)$ interval $(\frac{1}{16}, \frac{1}{8})$. Using $m = 3$ gives us the $(+, -)$ interval $(\frac{3}{8}, \frac{2}{5})$ and the $(+, +)$ interval $(\frac{1}{3}, \frac{3}{8})$.

We can find all of these intervals through the Farey tree construction starting with the sequence of the two rational numbers $0/1$ and $1/2$. That is, between a pair of adjacent rationals $m_1/n_1$ and $m_2/n_2$ in the sequence, we insert the Farey mediant $(m_1+m_2)/(n_1+n_2)$. This will eventually give us all rationals between $0$ and $1/2$, and in a form where the numerator and denominator are relative prime. The relevance of this construction to the existence intervals, is that when we insert a Farey mediant $(m_1+m_2)/(n_1+n_2)$ between $m_1/n_1$ and $m_2/n_2$, we have found for the period $n = n_1 + n_2$ a $(+, -)$ interval $(m_1+m_2)/(n_1+n_2) < s < m_2/n_2$ and a $(+, +)$ interval $m_1/n_1 < s < (m_1+m_2)/(n_1+n_2)$ (unless $m_1/n_1 = 0/1$, when instead we have a $(-, +)$ interval $0 < s < 1/(2(1+n_2))$ and a $(+, +)$ interval $1/(2(1+n_2)) < s < 1/(1+n_2)$).

The intervals on the line $r = 1$ can be extended into areas of regions 3 and 4, with $r < 1$. Is also seems that all existence areas include one of the intervals. We make the following hypothesis, based on numerical experiments:

**Conjecture 1.** When $n \ge 3$ the areas in regions 3 and 4 where period $n$ exists are associated with the intervals in the following way:

* an interval $(0, \frac{1}{2n})$ extends to a $(-, +)$ area bounded by $0 < s < 1/(2n)$;

* an interval $(\frac{1}{2n}, \frac{1}{n})$ extends to a $(+, +)$ area bounded to the left by $b_{n,n-1} = 0$ and to the right by $s = 1/(2n)$;

* an interval $(\frac{(n-1)/2}{n}, \frac{1}{2})$ (where $n$ is necessarily odd) extends to a $(+, -)$ area bounded to the left by $-(a_2+1) < a_1$ and to the right by $b_{n,n-2} = 0$;

* an interval $(\frac{m_+}{n_+}, \frac{m_-}{n_-})$ with $\frac{m_+}{n_+} \ne \frac{0}{1}$ extends to a $(+, +)$ area bounded by $b_{n,n-n_+} = 0$;

* an interval $(\frac{m_-}{n_-}, \frac{m_+}{n_+})$ with $\frac{m_-}{n_-} \ne \frac{1}{2}$ extends to a $(+, -)$ area bounded by $b_{n,n-n_-} = 0$.

The different active iterate boundaries $b_{n,k} = 0$ for given $n$ do not seem to cross in the $a_1-a_2$ plane, and all iterate boundaries with $3 \le n \le 15$ are shown in figure 4.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Iterate boundaries $b_{n,k} = 0$ for $n \le 15$. Boundaries with $n \le 5$ are labelled $n, k$.</figcaption></figure>

### 3.6. Eigenvalue of largest modulus

For all the orbits considered here, one eigenvalue approaches infinity as $\mu \to 0$, as follows from (61). From this formula we can also see the sign of $\Lambda$:

* in region 1, $\operatorname{sign}(\Lambda) = -\operatorname{sign}(e\mu)$;

* in regions 2–4, $\operatorname{sign}(\Lambda) = \operatorname{sign}(e\mu)$;

* in region 5, $\operatorname{sign}(\Lambda) = (-1)^{n+1} \operatorname{sign}(e\mu)$.

### 3.7. Summary by period

The results obtained can be summarized as follows.

Period 1 exists for $(+, +)$ and $(-, -)$ when $a_2 + 1 < a_1$ and for $(+, -)$ and $(-, +)$ when $a_1 < a_2 + 1$. There is a sign boundary at $a_1 = a_2 + 1$.

Period 2 exists for $(+, +)$ when $a_2 + 1 < a_1$, for $(-, +)$ when $0 < a_1 < a_2 + 1$, for $(+, +)$ when $-(a_2 + 1) < a_1 < 0$, and for $(+, -)$ when $a_1 < -(a_2 + 1)$. The line $a_1 = a_2 + 1$ is a sign boundary and an iterate boundary of order 1. At $a_1 = 0$ there is a singular boundary and $a_1 = -(a_2 + 1)$ is a sign boundary.

When $n \ge 3$, period $n$ exists for $(+, +)$ when $a_2 + 1 < a_1$, and for $(-, +)$ when $2\sqrt{a_2} \cos(\pi/n) < a_1 < a_2 + 1$. For each integer $m$ relative prime to $n$ and with $0 < m < n/2$ there is a $(+, +)$ area and a $(+, -)$ area in regions 3 and 4. If $n$ is odd the leftmost $(+, -)$ area continues into region 5. For even $n$ region 5 is a separate $(+, -)$ area. At $a_1 = a_2 + 1$, all $b_{n,k}$ and $c_n$ are zero, and $a_1 = 2\sqrt{a_2} \cos(\pi/n)$ is a singular boundary. In regions 3 and 4 there is an iterate boundary $b_{n,k} = 0$ for each $k$ relative prime to $n$. If $n$ is even, $c_n = 0$ and $b_{n,k} = 0$ for even $k$ at $a_1 = -(a_2 + 1)$.

Figure 5 shows the different areas for $n = 8$, and it can be seen how each of the intervals at $r = 1$ given above is connected with an area.
<figure><img src="image_6.png" /><figcaption>Figure 5. Existence areas for period 8.</figcaption></figure>

### 3.8. Significance of the boundaries

There are some cases where a boundary in parameter space is associated with a branch switching from existing for one sign of $\mu$, to the other sign. For $n = 1$, this happens when $\gamma = 0$. For $n > 1$, it happens in region 3 when $\theta = \pi/n$. An example is given in the next section.

An iterate boundary shows the possibility of having the bifurcating orbit itself grazing.

At the boundary $a_1 = a_2 + 1$, the existence of the non-impacting fixed point is not guaranteed. The fixed point is not hyperbolic when $a_1 = a_2 + 1$, $a_1 = -(a_2 + 1)$ or $a_2 = 1$ in region 3. In these cases there might be other non-impacting solutions involved in the bifurcation.

### 3.9. Multiple-impact orbits

When $N = 2$, solving for impact patterns with more than one impact becomes particularly simple. Using the fact that $CB = 0$, it is easy to show that solutions must satisfy (51) and (52) with $X = X_i$ and $n = n_i$ for all $i$, and the sign condition (53) changes into

$$ \operatorname{sign}(e\mu) \operatorname{sign}(\gamma) \begin{bmatrix} 0 & 1 \end{bmatrix} (X_{i+1} - A^{n_i} X_i) > 0 \quad \text{for } 1 \le i \le m. \quad (64) $$

This is then a condition telling whether period $n_{i+1}$ can follow period $n_i$. Since the $X_i$ are the same as for the single-impact situation, we see that multiple-impact orbits tend to follow $\mu^2$ close to single-impact orbits between impacts, but it can switch orbit at impacts. In region 1 for example, any pair of periods can follow one another in the $(+, +)$ case, so we can immediately conclude that a unique periodic orbit exists for any given impact pattern.
## 4. Some numerical examples

A particularly simple case to analyse is the grazing impact of a non-impacting orbit of a linear periodically driven impact oscillator with a coefficient of restitution. This system has become the most widely studied example of grazing bifurcations. Since we can express the motion between impacts in closed form, it is sometimes possible to find closed-form impacting periodic solutions by pasting non-impacting solutions together. In particular, single-impact orbits can be found in this way. In the following, we will not use this technique, but rather derive results using the local theory presented in earlier sections.

As the local theory does not predict any stable bifurcating periodic orbits (except perhaps for degenerate cases, see below), the question arises as to where motion will go after the bifurcation. In the generic case, the grazing periodic orbit at $\mu = 0$ is not stable unless we are in region 2 with $\operatorname{sign}(\gamma) > 0$, see Nordmark [1]. Then the motion after bifurcation grows continuously out of the non-impacting orbit, and the scenario is determined by the value of the largest eigenvalue. In the other cases, a jump to an attracting motion that is not necessarily close must take place. This analysis is then not local, and we do not attempt to give any analysis here. These cases come up in the numerical example given below, though.

The linear oscillator can be written in a number of equivalent forms. Here we will use (in the notation of section 3)

$$a(q, u, t) = -\frac{2d}{w}u - \frac{1}{w^2}q \qquad (65)$$

$$Q_1(t) = (1 + \mu) \cos(t) - 1 \qquad (66)$$

$$U_1(u, t) = -\tau u - (1 + \tau)(1 + \mu) \sin(t) \qquad (67)$$

with $d, w > 0$. There is a non-impact solution $q(t) = 0, u(t) = 0$ when $-1 < \mu < 0$. At $\mu = 0$, it undergoes a grazing bifurcation with $t_G = 0, a_G = 1$ and $e = -1$. The non-impact Poincaré mapping $f$ becomes

$$f(x) = \exp \left( 2\pi \begin{bmatrix} 0 & 1 \\ -\frac{2d}{w} & -\frac{1}{w^2} \end{bmatrix} \right) \left( x - \begin{bmatrix} -\mu \\ 0 \end{bmatrix} \right) + \begin{bmatrix} -\mu \\ 0 \end{bmatrix}. \qquad (68)$$

For $d > 1$ we are in region 2 with

$$a_1 = 2 \exp(-2\pi d/w) \cosh(2\pi \sqrt{d^2 - 1/w}) \qquad (69)$$

$$a_2 = \exp(-4\pi d/w) \qquad (70)$$

$$\operatorname{sign}(\gamma) > 0. \qquad (71)$$

By varying $d$ and $w$, all $a_1$ and $a_2$ in this region are accessible, but the sign of $\gamma$ is fixed. For $0 < d < 1$ we are in region 3 with

$$a_1 = 2 \exp(-2\pi d/w) \cos(2\pi \sqrt{1 - d^2}/w) \qquad (72)$$

$$a_2 = \exp(-4\pi d/w) \qquad (73)$$

$$\operatorname{sign}(\gamma) = \operatorname{sign} \left( \sin \left( 2\pi \sqrt{1 - d^2}/w \right) \right). \qquad (74)$$

All combinations of $a_1, a_2$ and $\operatorname{sign}(\gamma)$ are available here.

Apart from the grazing bifurcation of the non-impacting solution at $\mu = 0$, the system displays lots of other grazing bifurcations involving impacting periodic orbits. For these
bifurcations, any combination of $a_1, a_2 < 1$ and $\text{sign}(\gamma)$ should be possible. We will probably have to determine the orbit and bifurcation point numerically, and thus obtain numerical values of the derivatives needed, but in other respects the analysis method is as outlined in this work.

### 4.1. Bifurcation at a resonance

When $d = 0.6$ and $w = 4.8$, we have $\text{sign}(\gamma) > 0$ and complex eigenvalues with $r = \exp(-\pi/4)$ and $\theta/(2\pi) = \frac{1}{6}$. We find that single-impact orbits of periods 1 and 2 should exist for small positive $\mu$ (negative $e\mu$), and period 4 exists for $\mu < 0$. Period 3 is right at the singular boundary. This type of degenerate bifurcation was treated by Ivanov [5] using exact equations for the single-impact orbits. To study it by local techniques we need an expansion to order two. We find

$$h = x_1 - x_2^2/2 + \mathcal{O}(x, y, \mu)^3 \quad (75)$$

$$b = (1 + \tau) \left[ \sqrt{2} - \frac{4}{w} \frac{d}{y} - 3\sqrt{2} \frac{d}{w} x_2 + \frac{\sqrt{2}}{2} \frac{w^2 - 1}{w^2} \mu \right] + \mathcal{O}(x, y, \mu)^2. \quad (76)$$

Seeking a solution of (19) of the form $x(\epsilon)$, $y(\epsilon)$ and $\mu(\epsilon)$, where $\epsilon$ is a parameter along the solution curve, we find

$$x_1 = \mathcal{O}(\epsilon)^2 \quad (77)$$

$$x_2 = \sqrt{2}(1 + \tau)\epsilon + \mathcal{O}(\epsilon)^2 \quad (78)$$

$$y = (1 + r^3)\epsilon + \mathcal{O}(\epsilon)^2 \quad (79)$$

$$\mu = (1 - \tau r^3)^2 \epsilon^2 + \mathcal{O}(\epsilon)^3 \quad (80)$$

$$\epsilon > 0. \quad (81)$$

Eliminating $\epsilon$, we see that period 3 exists for $\mu > 0$, but branches off like a square root of $\mu$, instead of linearly in $\mu$ as the other periodic solutions do. The eigenvalues can be determined using (22), and are both found to be $\tau r^3 + o(1)$. The branch is thus stable as $\mu$ approaches 0, in contrast to the other branches that are not associated with existence borders. These results were also derived by Ivanov.

### 4.2. Some bifurcations near the resonance

When $w$ is slightly smaller than 4.8, the period 3 branch first goes towards negative $\mu$ but soon turns in a saddle-node bifurcation and becomes stable. For $w$ a bit larger than 4.8, it exists for positive $\mu$ but soon becomes stable in a flip bifurcation. We will look at some numerical results for $w = 4.9$, $\tau = 1$.

The large eigenvalues are asymptotically $-32.1/\mu$ for period 1, $-4.37/\mu$ for period 2, $-0.00504/\mu$ for period 3 and $0.203/(-\mu)$ for period 4. We can follow the different orbits out from the bifurcation point, to see how they affect the dynamics of the system. We will follow the orbits until they bifurcate, and say something about what happens after these (secondary) bifurcation points. The non-impact and single-impact branches are plotted in figure 6. Stable branches are solid and unstable broken. The coordinate $x$ used is $q - (Q_1 + 1)$ sampled when $t = \pi \mod 2\pi$. Instead of plotting all $n$ $x$ values for an orbit of period $n$, only one $x$ value (the sample taken about half a driving period after the low-velocity impact) per
<figure><img src="image_7.png" /><figcaption>Figure 6. Non-impact and single-impact branches for $d = 0.6$, $w = 4.9$, $τ = 1$.</figcaption></figure>

period is plotted. Following each branch from the bifurcation and outwards, we find the following.

Periods 1 and 2 will stay unstable with an eigenvalue < -1 for all positive $\mu$.

Period 3 will go through a flip bifurcation at $\mu_{pd} = 0.005479$ and become stable. The flip is subcritical, however, so an unstable double-impact orbit of period 6 exists for a small interval above $\mu_{pd}$. This double-impact orbit then encounters a grazing bifurcation at $\mu_{gr3} = 0.006100$, with $\text{sign}(\gamma) > 0$ and $\lambda_1 = 2.15$ in region 1. Thus unstable orbits of all periods that are multiples of 6 exist below $\mu_{gr3}$.

Period 4 will go through a saddle-node bifurcation at $\mu_{sn} = -0.084338$ and turn back towards more positive $\mu$ on a stable branch. At $\mu_{gr4} = -0.024447$, it will encounter a grazing bifurcation, with $\lambda_1 = 0.42$, $\lambda_2 = 0.0050$, $\gamma > 0$ and $e < 0$. This is in region 2, and we predict the creation of orbits with periods that are any multiple of 4. This situation can be more thoroughly analysed (see, for example, Chin et al [6] or Nordmark [9]), and it has been shown that, as $\mu$ increases from $\mu_{gr}$, all orbits of sufficiently long period first becomes stable in a subcritical flip bifurcation and then disappears in a grazing bifurcation in region 4 with $\gamma < 0$. Since $\frac{1}{4} < \lambda_1 < \frac{2}{3}$, the windows where the orbits are stable do not overlap, and in between the windows there are intervals of chaos. In figure 7 a numeric bifurcation diagram of the motion is shown. We also show the stable period 4 branch for $\mu < \mu_{gr4}$, and the unstable period 4 branch. Only one of the four components (the same as in figure 6) is shown, and it should be noted that the $\mu$ scale is much enlarged over that figure. Although an infinite number of periodic windows accumulate upon $\mu = \mu_{gr4}$, only the two of lowest periods are clearly visible in the figure, with a thin chaotic band in between. The periods are $4 \times 7$ and $4 \times 8$, respectively. Since each window has a length that is $\lambda_1^2$ times the length of the previous one, they narrow very quickly with increasing period. The asymptotic results obtained near the bifurcation point fail to predict what happens after the largest periodic window. Here instead we have a large chaotic band that extends up to $\mu_{cr} = -0.022577$. At this point, the chaotic attractor collides with the stable manifold of the unstable period 4 orbit that comes
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Attracting motion between $\mu_{gr4}$ and $\mu_{cr}$, together with single-impact orbit of period 4.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Stable non-impacting period 1 (square), unstable single-impact period 4 (stars), and stable single-impact period 4 (rings) when $\mu = -0.05$. Points to the left of the curve will encounter a low-velocity impact during the next iteration.</figcaption></figure>

from the original grazing bifurcation. The collision is also visible in figure 7. The stable manifold of the unstable period 4 also marks the boundary of the basin of attraction for the non-impacting solution. Thus chaotic motion is no longer attracting and iterates converge to the non-impacting solution.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Stable non-impacting period 1 (square), unstable single-impact period 4 (stars), and noisy period 4 attractor when $\mu = -0.024$. Points to the left of the curve will encounter a low-velocity impact during the next iteration.</figcaption></figure>

When $\mu$ is increased from negative values, the following attracting solutions will then be seen.

* For $\mu < \mu_{sn}$, only the non-impacting period 1 orbit exists.

* For $\mu_{sn} < \mu < \mu_{gr4}$, a stable single-impact period 4 coexists with the non-impacting orbit. Figure 8 shows the Poincaré section when $\mu = -0.05$. Note how the stable period 4 orbit is close to grazing.

* For $\mu_{gr4} < \mu < \mu_{cr}$, either a stable $n+1$ impact period $4n$ with $n \ge 7$, or a chaotic (noisy period 4) motion coexists with the non-impacting orbit. Figure 9 shows the Poincaré section when $\mu = -0.024$. We see how the attractor now contains points with grazing trajectories.

* For $\mu_{cr} < \mu < 0$ only the non-impacting orbit is attracting.

* For $0 < \mu < \mu_{pd}$ we have not identified any candidate in the preceding section. Figure 10 shows the Poincaré section when $\mu = 0.003$. Numerically, it looks as if the chaotic motion existing for $\mu < \mu_{cr}$ becomes attracting again as the non-impacting orbit loses stability at $\mu = 0$. This time it is connected instead of being noisy and of period 4. The unstable single-impact orbits of periods 1–3 appear to be part of this attractor. This $\mu$ interval also contains periodic windows.

* For $\mu_{pd} < \mu$, we have a stable period 3 orbit.

### 4.3. Some comments on the numerical methods

Single-impact orbits for the linear impact oscillator are known to be found by solving a quadratic equation for the sine of the phase at impact, together with a check that the incoming velocity has the right sign. This will also give the eigenvalues, the point of the saddle-node bifurcation for
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Unstable single-impact period 3 (rings), and attractor when $\mu = 0.003$. Points to the left of the curve will encounter a low-velocity impact during the next iteration.</figcaption></figure>

period 4, and the flip bifurcation for period 3. Multiple-impact orbits were solved by Newton's method, which also gives the stability. Points of grazing periodic orbits were solved by using the grazing phase and $\mu$ as unknowns, integrating backwards and forwards from the assumed grazing point to a Poincaré section, and using the equality of the forward and backward points in the sections as two equations for the two unknowns. The resulting system of two equations was solved using Newton's method, and this procedure also gives the matrix $A$ and the sign of $\gamma$ and $e$. The point $\mu_{cr}$ was found using an interactive search.

**5. Review of results and discussion**

We have studied the problem of finding periodic orbits that are created in grazing bifurcations occurring in models of impacting systems. The main result is that for a given impact pattern, we can formulate a system of equations that is smooth, and we can find sufficient conditions for existence in the form of a finite set of linear equations combined with a finite set of linear inequalities. We also find a stability estimate for these orbits, and generally at least one eigenvalue becomes unbounded as the bifurcation point is approached.

For systems with a two-dimensional Poincaré mapping, such as a periodically driven oscillator with one degree of freedom, we have also made a parameter study for the existence of single-impact orbits. Especially when eigenvalues are complex and the damping is low, a large number of orbits is possible, and existence is connected with number theory.

As a local theory, we are not able to say much about what will happen to the periodic orbits as we become further from the bifurcation point. We obtain an estimate of the dependence of the largest eigenvalue on the bifurcation parameter, but the points where (secondary) bifurcations take place are usually not close. There is also always the possibility that the orbits will
encounter another grazing bifurcation before becoming stable. There is, however, one case (corresponding to (-, +) in region 2 for a 2D system) where points with long periods become stable in a flip bifurcation for small values of the bifurcation parameter. This phenomenon, as well as other fates of the bifurcating orbits, are shown in a numerical example. By following the branches, we will also encounter several other grazing bifurcations, and this emphasizes the importance of grazing bifurcations in the study of these systems.

Another important aspect of the unstable periodic orbits in the two-dimensional system is that they are saddle points. We have not stressed this aspect here, but it is clearly of importance. For example, before becoming unstable, the non-impacting fixed point has a basin of attraction, and numerically it seems that the border of the basin consists of the stable manifold of the closest of the bifurcating orbits. The basin often has a typical shape that is polygon with thin attached fibres (for example, in Foale and Bishop [2]), and the polygon is formed by successive preimages of the line $h = 0$ under $f$.

The fact that grazing bifurcations usually produce a large number of unstable periodic orbits clearly has importance for the global dynamics of the system. A conjecture is that chaotic motion is also usually produced in these bifurcations, if not necessarily attracting. The fact that a grazing bifurcation need not involve any stable periodic orbits can make them go unnoticed in numerical simulations. In a situation close to that of figure 6, the orbit of period 4 might very well go through a grazing bifurcation before having a chance to become stable. In that situation, an invisible grazing bifurcation will be the first bifurcation to occur as the driving amplitude is increased, and complex motion might well be created at this point.

## Appendix A. Signs of $b_{n,k}$ and $c_n$ in selected regions

### A.1. Sign of $b_{n,k}$ in regions 1, 2 and 5

The expression (59) may be written as

$$b_{n,k} = \frac{(\lambda_1^n - 1)(1 - \lambda_2^k) - (\lambda_1^k - 1)(1 - \lambda_2^n)}{\lambda_1^n - \lambda_2^n} \quad (A1)$$

or

$$b_{n,k} = \frac{(\lambda_1 - 1)(1 - \lambda_2)(\lambda_1 - \lambda_2) \sum_{i=k}^{n-1} \sum_{j=0}^{k-1} \sum_{m=j}^{i-1} \lambda_1^m \lambda_2^{i+j-m-1}}{(\lambda_1 - \lambda_2) \sum_{i=0}^{n-1} \lambda_1^i \lambda_2^{n-i-1}} \quad (A2)$$

In regions 1 and 2, where both $\lambda_1$ and $\lambda_2$ are positive, the expression (A2) shows that the sign of $b_{n,k}$ is the sign of $(\lambda_1 - 1)(1 - \lambda_2)$, so it immediately follows that $b_{n,k} > 0$ in region 1 and $b_{n,k} < 0$ in region 2.

In region 5, we separate into cases. If both $n$ and $k$ are even integers, the value of $b_{n,k}$ is unchanged if we change the signs of $\lambda_1$ and $\lambda_2$, so the result from region 1 gives $b_{n,k} > 0$. If $n$ is odd and $k$ even, we obtain from (A1)

$$b_{n,k} = \frac{(1 + |\lambda_1|^n)(|\lambda_2|^k - 1) + (1 - |\lambda_1|^k)(|\lambda_2|^n + 1)}{|\lambda_2|^n - |\lambda_1|^n} > 0. \quad (A3)$$

If $n$ is even and $k$ odd,

$$b_{n,k} = \frac{(1 - |\lambda_1|^n)(|\lambda_2|^k + 1) + (1 + |\lambda_1|^k)(|\lambda_2|^n - 1)}{|\lambda_2|^n - |\lambda_1|^n} > 0. \quad (A4)$$
If both $n$ and $k$ are odd, the expression (59) becomes

$$b_{n,k} = \frac{(|\lambda_2|^n - |\lambda_2|^k) + (|\lambda_1|^k - |\lambda_1|^n) + (|\lambda_1||\lambda_2|)^k (|\lambda_2|^{n-k} - |\lambda_1|^{n-k})}{|\lambda_2|^n - |\lambda_1|^n} > 0 \quad (A5)$$

and we have shown that $b_{n,k} > 0$ in region 5 for all values of $n$ and $k$.

## A.2. Sign of $c_n$

The expression (60) immediately gives $c_n > 0$ in region 1 and $c_n < 0$ in region 2. In region 4, we write

$$c_n = (-1)^n \frac{(|\lambda_2| - |\lambda_1|)(1 - \lambda_1^n)(1 - \lambda_2^n)}{|\lambda_2|^n - |\lambda_1|^n} \quad (A6)$$

so $c_n > 0$ for even $n$ and $c_n < 0$ for odd $n$. In region 5, we write

$$c_n = - \frac{(|\lambda_2| - |\lambda_1|)(1 - \lambda_1^n)(|\lambda_2|^n - (-1)^n)}{|\lambda_2|^n - |\lambda_1|^n} \quad (A7)$$

so $c_n < 0$ for all $n$.

In region 3, we write $\lambda_{1,2} = re^{\pm i\theta}$ with $0 < \theta < \pi$, and

$$c_n = - \frac{\sin(\theta)}{r^{n-1} \sin(n\theta)} (1 - 2r^n \cos(n\theta) + r^{2n}) \quad (A8)$$

so $\text{sign}(c_n) = -\text{sign}(\sin(n\theta))$.

# Appendix B. Proof of theorem 2

When $r = 1$ we can write

$$b_{n,k}(s) = -2 \frac{\sin(k\pi s) \sin((n-k)\pi s)}{\cos(n\pi s)} \quad (B1)$$

$$c_n(s) = -2 \sin(2\pi s) \tan(n\pi s) \quad (B2)$$

where the variable $s = \theta/(2\pi) \in (0, \frac{1}{2})$. We should assume $n \ge 3$. For given $n$, we are interested in open subintervals $\sigma_i \subset (0, \frac{1}{2})$ in the $s$ variable, where all $b_{n,k}$ are of the same sign when $1 \le k \le n-1$ and where $c_n$ also has a definite sign.

Firstly, $c_n$ changes sign at all multiples of $\frac{1}{2n}$, so let $\tau_j = (\frac{j-1}{2n}, \frac{j}{2n})$ for $1 \le j \le n$. Then the sign of $c_n$ in $\tau_j$ is $(-)^j$. Any $\sigma_i$ must be a subinterval of some $\tau_j$.

The possible zeros of $b_{n,k}$ are located at $s = \frac{m}{k}$ or $s = \frac{m}{n-k}$ for non-negative integer $m$. One notes that the denominator is always smaller than $n$. If $s$ is a double zero, and thus $s = \frac{m_1}{k} = \frac{m_2}{n-k}$, then $s = \frac{m_1+m_2}{n}$ (since if $\frac{m_1}{n_1} \le \frac{m_2}{n_2}$, then $\frac{m_1}{n_1} \le \frac{m_1+m_2}{n_1+m_2} \le \frac{m_2}{n_2}$). Thus a zero of $b_{n,k}$ is double if and only if it is located at an even multiple of $\frac{1}{2n}$. We also find that all poles are simple and located at odd multiples of $\frac{1}{2n}$. Thus, within the intervals $\tau_j$, $b_{n,k}$ is smooth with only simple zeros. In fact, for given $n$ and $k$, at most one such zero exists within a $\tau_j$, since between $\frac{m_1}{k} \ne \frac{m_2}{n-k}$ lies $\frac{m_1+m_2}{n}$, which is a multiple of $\frac{1}{2n}$. Likewise, if a $b_{n,k}$ has a (double) zero at an even multiple of $\frac{1}{2n}$, there are no zeros in the two adjoining $\tau_j$.

Again, we want to locate subintervals $\sigma_i$ of the $\tau_j$, where $b_{n,k}$ is of the same sign for all $k$. The strategy is to check $b_{n,k}$ values close to multiples of $\frac{1}{2n}$, and then show that any $\sigma_i$ must have at least one endpoint that is a multiple of $\frac{1}{2n}$.

We start by checking the values when $s$ is an even multiple of $\frac{1}{2n}$.
• If $s = \frac{m}{n}$, then $b_{n,k} = 2 \sin^2(\pi km/n)$.

• If $s = \frac{m}{n} + \epsilon$ and $\frac{km}{n}$ is an integer, $b_{n,k} = -2\pi^2 k(n-k)\epsilon^2 + \mathcal{O}(\epsilon^4)$.

We will say that two positive integers $m$ and $n$ are relative prime if and only if there are integers $m_1$ and $n_1$ such that $mn_1 - m_1n = 1$. If $m$ and $n$ are relative prime, then $\frac{km}{n}$ is never an integer for the allowed $k$ values, and $b_{n,k}$ is always positive at and near $s = \frac{m}{n}$. If $m$ and $n$ are not relative prime, there is an allowed value of $k$ such that $\frac{km}{n}$ is an integer, and $b_{n,k}$ is zero at $s = \frac{m}{n}$ and negative on either side of this point. Since $b_{n,1}$ is always positive, all the $b_{n,k}$ do not have the same sign near $\frac{m}{n}$ if $m$ and $n$ are not relative prime.

Now suppose $m$ and $n$ are relative prime. Since zeros of $b_{n,k}$ are rational numbers with denominator smaller than $n$, we take advantage of the fact that the largest rational $\frac{m_+}{n_+}$ smaller than $\frac{m}{n}$ and with $0 < n_+ < n$ is uniquely given by $mn_+ - m_+n = 1$. If $n_+ > 2$, then $\frac{2m-1}{2n} < \frac{m_+}{n_+}$ and all the $b_{n,k}$ are positive for $s \in (\frac{m_+}{n_+}, \frac{m}{n}) \subset \tau_{2m}$. $n_+ = 2$ is not possible. If $n_+ = 1$, which only happens for $m_+ = 0$ and $m = 1$, in the interval $(0, \frac{1}{n})$ all the $b_{n,k}$ have a pole at $s = \frac{1}{2n}$, so $b_{n,k}$ is negative in $(\frac{0}{1}, \frac{1}{2n}) = \tau_1$ and positive in $(\frac{1}{2n}, \frac{1}{n}) = \tau_2$.

The closest zero above $\frac{m}{n}$ is found similarly since the smallest rational $\frac{m_-}{n_-}$ larger than $\frac{m}{n}$ and with $0 < n_- < n$ is uniquely given by $mn_- - m_-n = -1$. Here, $n_-$ is always at least 2, so all the $b_{n,k}$ are positive in the interval $(\frac{m}{n}, \frac{m_-}{n_-}) \subset \tau_{2m+1}$. When $n_- = 2$, which only happens for odd $n$ with $m = (n-1)/2$ and $m_- = 1$, the interval is $(\frac{(n-1)/2}{n}, \frac{1}{2}) = \tau_n$.

Next we look at limit values near odd multiples of $\frac{1}{2n}$.

• If $s = \frac{2m+1}{2n}$ and $\frac{k(2m+1)}{n}$ is an odd integer, then $b_{n,k} = \frac{2(n-k)}{n}$.

• If $s = \frac{2m+1}{2n}$ and $\frac{k(2m+1)}{n}$ is an even integer, then $b_{n,k} = \frac{2k}{n}$.

• If $s = \frac{2m+1}{2n} + \epsilon$ and $\frac{k(2m+1)}{n}$ is not an integer, then $b_{n,k} = \frac{\sin(\pi k (2m+1)/n)}{\pi n \epsilon} + \mathcal{O}(1)$.

We find only three cases where $b_{n,k}$ has the same sign for all $k$ when $s$ is near an odd multiple of $\frac{1}{2n}$: below $\frac{1}{2n}$ it is negative, above $\frac{1}{2n}$ it is positive and below $\frac{1}{2}$ (odd $n$) it is positive. All of these were already considered above.

Now we have to show that a $\sigma_i$ must have an endpoint on a multiple of $\frac{1}{2n}$. Study a particular $\tau_j$. Suppose that all $b_{n,k}$ are positive at some interior point. There can be no zero of any $b_{n,k}$ between this point and the nearest even multiple of $\frac{1}{2n}$, since that would make $b_{n,k}$ negative there. Thus a $\sigma_i$ with positive values must extend to an even multiple of $\frac{1}{2n}$. Suppose instead that all $b_{n,k}$ are negative at some interior point. All $b_{n,k}$ must have a zero between (possibly at) the nearest even multiple of $\frac{1}{2n}$ and this point. Then there can be no zeros towards the other endpoint, and a $\sigma_i$ with negative values must extend to an odd multiple of $\frac{1}{2n}$. Since the signs near all endpoints are established above, all the $\sigma_i$ are already listed.

To summarize:

• when $s \in (0, \frac{1}{2n})$, $b_{n,k}$ is negative and $c_n$ is negative: this gives a $(-, +)$ case;

• when $s \in (\frac{1}{2n}, \frac{1}{n})$, $b_{n,k}$ is positive and $c_n$ is positive: this gives a $(+, +)$ case;

• when $m$ and $n$ are relative prime, $0 < m < n/2$, and $s \in (\frac{m}{n}, \frac{m_-}{n_-})$, $b_{n,k}$ is positive and $c_n$ is positive: this gives a $(+, -)$ case;

• when $m$ and $n$ are relative prime, $1 < m < n/2$, and $s \in (\frac{m_+}{n_+}, \frac{m}{n})$, $b_{n,k}$ is positive and $c_n$ is negative: this gives a $(+, +)$ case.
References

[1] Nordmark A B 1991 Non-periodic motion caused by grazing incidence in an impact oscillator *J. Sound Vib.* **145** 279–97

[2] Foale S and Bishop S R 1992 Dynamical complexities of forced impacting systems *Phil. Trans. R. Soc. A* **338** 547–56

[3] Budd C and Dux F 1994 Intermittency in impact oscillators close to resonance *Nonlinearity* **7** 1191–224

[4] Foale S and Bishop S R 1994 Bifurcations in impact oscillators *Nonlinear Dynam.* **6** 285–99

[5] Ivanov A P 1993 Stabilization of an impact oscillator near grazing incidence owing to resonance *J. Sound Vib.* **162** 562–5

[6] Chin W, Ott E, Nusse H E and Grebogi C 1994 Grazing bifurcations in impact oscillators *Phys. Rev. E* **50** 4427–44

[7] Chin W, Ott E, Nusse H E and Grebogi C 1995 Universal behaviour of impact oscillators near grazing incidence
*Phys. Lett. A* **201** 197–204

[8] Fredriksson M H and Nordmark A B 1997 Bifurcations caused by grazing incidence in many degree of freedom impact oscillators *Proc. R. Soc. A* **453** 1261–75

[9] Nordmark A B 1996 Universal limit mapping in grazing bifurcations *Phys. Rev. E* **55** 266–70