# Mappings of grazing-impact oscillators

Jaap Molenaar, John G de Weger, Willem van de Water

► To cite this version:

Jaap Molenaar, John G de Weger, Willem van de Water. Mappings of grazing-impact oscillators. Nonlinearity, 2001, 10.1088/0951-7715/14/2/307. hal-02520794

HAL Id: hal-02520794

https://hal.science/hal-02520794v1

Submitted on 26 Mar 2020

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Mappings of grazing-impact oscillators

Jaap Molenaar<fnref n="1" />, John G de Weger<fnref n="2" /> and Willem van de Water<fnref n="2" />

## Abstract

Impacting systems are found in a great variety of mechanical constructions and they are intrinsically nonlinear. In this paper it is shown how near-grazing systems, i.e. systems in which the impacts take place at low speed, can be described by discrete mappings. The derivation of this mapping for a harmonic oscillator with a stop is dealt with in detail. It is found that the resulting mapping for rigid obstacles is somewhat different from those presented earlier in the literature. The derivations are extended to systems with a compliant obstacle. We find that the map for impacts with a compliant obstacle is very similar to the one describing collisions with a rigid obstacle. A notable difference is a change of scale of the bifurcation parameter. We illustrate our findings in the limit of large damping, where the mechanism of period adding can be analysed exactly. The relevance of our results to experiments on practical impact systems is indicated.

## 1. Introduction

Impacting behaviour is found in a large number of mechanical systems. Examples are gear rattle, ships colliding against fenders, loosely fitting joints, suspension bridges and ball bearings. In atomic force microscopy, impacts occur on a mesoscopic scale and the interaction forces may have an intricate form (Berg and Briggs 1998, van de Water and Molenaar 2000). Although the dynamics in between impacts may be linear, the collision introduces an essential nonlinearity. We believe that impacts are the principal source of nonlinearity in discrete mechanical systems.

A special situation arises when the impacts are with zero velocity, so-called *grazing* impacts. For impacts that are close to grazing, it is possible to condense the mathematical description into a discrete mapping as was first done by Nordmark (1991, 1997); followed by Frederiksson and Nordmark (1997) and Frederiksson (1998). As the bifurcation properties of a mapping are analysed much more readily than those of a differential equation, the study

<footnote n="1">¹ Mathematics Department, Eindhoven University of Technology, PO Box 513, 5600 MB Eindhoven, The Netherlands</footnote>

<footnote n="2">² Physics Department, Eindhoven University of Technology, PO Box 513, 5600 MB Eindhoven, The Netherlands</footnote>
of these mappings led to the discovery of important organizing principles in the dynamics of impact oscillators.

The Nordmark map (Nordmark 1991) reads as

$$ \left. \begin{array}{l} x_{n+1} = \alpha x_n + y_n + \rho \\ y_{n+1} = -\gamma x_n \end{array} \right\} \quad \text{if } x_n \le 0 \qquad (1) $$

$$ \left. \begin{array}{l} x_{n+1} = -\sqrt{x_n} + y_n + \rho \\ y_{n+1} = -\gamma r^2 x_n \end{array} \right\} \quad \text{if } x_n > 0. \qquad (2) $$

Here, $(x_n, y_n)$ are the transformed coordinates of the two-dimensional phase space of the original impact oscillator taken at stroboscopic times $t = n \in \mathbb{Z}$. Equation (1) describes the non-impacting case with $x_n < 0$, whereas (2) is applied when the $x_n$-coordinate crosses the position of the wall at $x = 0$. The parameters $\alpha$ and $\gamma$ are related to the parameters of the ordinary differential equation, whereas the restitution coefficient $r$ gauges the energy loss at impact. The bifurcation parameter $\rho$ is proportional to the driving amplitude of the impact oscillator. At large negative values of $\rho$ no impacts occur, and equation (1) applies. When $\rho$ is increased in a quasi-stationary way from the non-impacting state, impacts first occur at $\rho = 0$. These impacts may be with a relatively large velocity and the transition that has taken place may be hysteretic, i.e. impacts may remain when $\rho$ is subsequently smoothly decreased and may only vanish at a negative value of $\rho$. The emergence of a square-root in (2) is characteristic for grazing collisions. Because the Jacobian of the map behaves as $x^{-1/2}$ at $x = 0$, we will refer to it as a square-root singularity. In the vicinity of the line $x = 0$, the map stretches phase space extremely.

Intimately tied to this square-root singularity is the phenomenon of period adding. Depending on the value of the parameters $\alpha$ and $\gamma$, this can come either as an infinite series of reverse period additions if $\rho$ is increased from 0, or as a series of (non-impacting) period-1 to period-$M$, $M = 1, 2, \dots$, transitions at given values of $\alpha$ and $\gamma$. In all cases, the $M$-periodic orbits reached are maximal periodic orbits that have a single impact per period. It is tempting to draw an analogy with period-doubling bifurcations, which are found in systems with a quadratic nonlinearity (Feigenbaum 1983). In both cases, the order of the nonlinearity determines a characteristic bifurcation scenario.

Figure 1 illustrates these period additions when the bifurcation parameter is reduced from a positive value downwards. As period additions are readily accessible in the experiment (de Weger et al 2000) we will concentrate on this characteristic feature of impact oscillators when we test the prediction of mappings.

The mapping (1) and (2) assumes that the wall is perfectly rigid so that impacts take place instantaneously. This idealization will not be fulfilled in real experiments where collisions occur with a more or less compliant wall which introduces a time delay at impact. The key question then is whether a finite compliance will 'smooth' the square-root singularity, thus essentially altering the typical bifurcation structure of grazing impact oscillators. The prime motivation of this work is a derivation for grazing impact mappings in the case of collisions with yielding obstacles. To this end we will first derive in sections 2 and 3 the mapping for a harmonic oscillator which undergoes grazing impacts with a hard wall. We will use a lucid and straightforward expansion that allows a precise association of the parameters of the resulting map with those of the differential equation. We will argue that (2) needs to be altered slightly but essentially. Having laid out our mathematical tools, the extension to the case of impacts with a non-rigid wall in section 6 appears straightforward.

Our procedure is to split the dynamics into a map $\mathcal{L}$ which describes the effect of the impact at stroboscopic times and a map which evolves the oscillator between collisions. Most
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Bifurcation diagram for $(\alpha, \gamma) = (0.25, 0.01)$ and $r = 0.9$ showing a series of period addings. We have used the mapping (37) and (38) with (48) which is derived for low-velocity impacts with a yielding wall. As bifurcation parameter we used, instead of $\rho$, the physical parameter $\sigma = F/F_g - 1$, where $F$ is the driving strength of the harmonic oscillator (3) for which the map was derived, and $F_g$ the driving strength (5) where impacts first occur. The relation between $\rho$ and $\sigma$ will be derived in this paper. Apart from a scale factor of the horizontal axis, the bifurcation diagram is almost identical to figure 3(b) in Chin *et al* (1994) which was computed using the Nordmark map (1) and (2).</figcaption></figure>

of the work is the derivation of $\mathcal{L}$. Such a procedure is not new (Frederikson and Nordmark 1997). In a recent paper Dankowicz and Nordmark (2000) discuss the general properties of a 'discontinuity map' (similar to $\mathcal{L}$, but now taken at the collision time itself) in the presence of a yielding wall.

In impacting systems there are two sources of discontinuity: the impact itself and the energy loss upon impact. The loss of energy is a crucial ingredient of impact oscillators. In an experimental study (de Weger *et al* 2000) of a mechanical oscillator we found collisions to be almost completely inelastic.

The study of Dankowicz and Nordmark (2000) is for perfectly elastic collisions. Depending on several small parameters, either a square-root ($x_n^{1/2}$) behaviour, or a dependence $x_n + c x_n^{3/2}$ of the impact map was found for generic systems that involve a yielding wall. Consequently, for some choice of the relative ordering of the parameters the square-root singularity of the Jacobian vanishes, thus essentially altering the dynamical properties of the system.

Practical systems suggest a natural ordering of the several small quantities of Dankowicz and Nordmark (2000), and a single parameter remains. Our conclusion will be that for most systems involving a yielding wall, an *effective* square-root behaviour survives. An exception is the case of perfectly elastic collisions where an *effective* square-root behaviour will only be felt when the impact velocity is large enough. Period-adding bifurcations are thus expected to also be found in experimental situations which involve compliant obstacles. In section 7 we will illustrate the predictive power of our mapping when we analyse period addings in the
limit of large friction.

## 2. Harmonic oscillator

A sketch of the physical system is shown in figure 2. We consider a periodically driven harmonic oscillator with mass $m$, friction coefficient $c$, spring constant $k$. The driving force has angular velocity $\omega$ and amplitude $f$. At an obstacle is placed a distance $\ell > 0$ from the rest position of the oscillator. In non-dimensional form the equation of motion reads as

$$ \ddot{u} + \nu\dot{u} + \Omega^2(u+1) = F \cos(2\pi t + \varphi) \quad (3) $$

with

$$ \nu = 2\pi c/m\omega $$

$$ \Omega^2 = 4\pi^2 k/m\omega^2 $$

$$ F = 4\pi^2 f/m\omega^2\ell. $$

The position has been scaled such that the rest position of the mass is at $u = -1$ and of the obstacle at $u = 0$. So, the distance $\ell$ is the unit of length. Time has been scaled such that the external forcing has frequency $2\pi$. For later convenience a phase difference $\varphi$ is introduced. With these conventions, the particular solution of (3) has the elegant form

$$ p(t) = (1 + \sigma) \cos(2\pi t) - 1 \quad (4) $$

with

$$ \sigma = (F - F_g)/F_g $$

$$ F_g = \sqrt{(\Omega^2 - 4\pi^2)^2 + 4\pi^2\nu^2} \quad (5) $$

$$ \varphi = \arccos[(\Omega^2 - 4\pi^2)/F_g]. $$

For $F = F_g$ the particular solution $p(t)$ grazes the wall at times $t = n \in \mathbb{Z}$ when it collides with zero velocity. We take the normalized excitation amplitude $\sigma$ as the bifurcation parameter. One of our tasks will be to derive the relation between $\sigma$ and the bifurcation parameter $\rho$ of the mapping. If we denote points in phase space by $u = (u, \dot{u})^T$, the solution $u_0(t)$ of the undriven oscillator (3) with $F = 0$ and with the obstacle ignored, is given by

$$ u_0(t) = \mathbf{P}(t - t_0)u_0(t_0) \quad (6) $$

<figure><img src="image_3.png" /><figcaption>**Figure 2.** The physical system studied in this paper. The mass $m$ of the oscillator can collide with a yielding wall. We assume that the massless wall is attached with frictionless springs to the fixed world. The resilience of the wall is determined by the stiffness of the second spring. We use normalized units with $m = 1$, the wall position at rest at $u = 0$, and the rest position of the oscillating mass at $u = -1$.</figcaption></figure>
with the evolution matrix $\mathbf{P}(t)$ of the free system

$$ \mathbf{P}(t) = \frac{1}{s_2 - s_1} \begin{pmatrix} s_2 e^{s_1 t} - s_1 e^{s_2 t} & e^{s_2 t} - e^{s_1 t} \\ s_1 s_2 (e^{s_1 t} - e^{s_2 t}) & s_2 e^{s_2 t} - s_1 e^{s_1 t} \end{pmatrix} \quad (7) $$

where $s_1 = \frac{1}{2}(-\nu + \sqrt{\nu^2 - 4\Omega^2})$ and $s_2 = \frac{1}{2}(-\nu - \sqrt{\nu^2 - 4\Omega^2})$. For later convenience we mention that $P_{11} - \nu P_{12} = P_{22}$. Note that $\mathbf{P}(0)$ equals the unit matrix. Given the initial conditions $u(t_0)$, the general solution of equation (3) is the sum of the solutions in (4) and (6) and can be written as

$$ u(t) = \mathbf{P}(t - t_0)[u(t_0) - p(t_0)] + p(t). \quad (8) $$

Expression (8) holds as long as no impact takes place, thus if $u < 0$. The time-one operator $\mathbf{P} \equiv \mathbf{P}(1)$ connects states at successive periods of the driving; its eigenvalues are $\exp(s_1)$ and $\exp(s_2)$. From this we directly conclude that $p(t)$ is asymptotically stable in the case of damping, i.e. if $\nu > 0$. So, the system tends to $p(t)$ in between impacts. In the overdamped case $\nu^2 > 4\Omega^2$ this convergence is fast and without oscillations, whereas in the underdamped case $\nu^2 < 4\Omega^2$ this approach is slow and $u(t)$ will oscillate around $p(t)$. The difference between under- and overdamped oscillations may lead to essentially different maps, as discussed in section 8. The trace $\alpha = \text{Tr}(\mathbf{P})$ and determinant $\gamma = \text{Det}(\mathbf{P})$ will play an important role in the following. They are given by the sum and the product of the eigenvalues:

$$ \begin{aligned} \alpha &= e^{s_1} + e^{s_2} = 2e^{-\frac{1}{2}\nu} \cosh\left(\frac{1}{2}\sqrt{\nu^2 - 4\Omega^2}\right) \\ \gamma &= e^{s_1} \times e^{s_2} = e^{-\nu}. \end{aligned} \quad (9) $$

### 3. Modelling the impact

To incorporate the effect of the obstacle, an impact model has to be introduced. We denote the velocity and time at collision by $(v_c, t_c)$ and the velocity and time at recoil by $(v_r, t_r)$. The duration of the impact, i.e. the time spent in the region $u > 0$, is denoted by $\Delta t$. An impact model then must specify the relations

$$ v_r = v_r(v_c, t_c) \quad \text{and} \quad t_r = t_c + \Delta t(v_c, t_c). \quad (10) $$

In the next sections we first derive expressions for the case of impacts with a perfectly rigid wall when the recoil is instantaneous ($\Delta t = 0$). After that the effect of a compliant wall is studied. In our approach the energy loss upon collision is done through a simple restitution model which allows for a useful comparison with the map derived by Nordmark (1991). Thus,

$$ v_r = -r v_c \qquad 0 \le r \le 1. \quad (11) $$

The case $r = 1$ corresponds to an elastic wall. If $r < 1$ energy is absorbed during the impact.

It is interesting to note that relation (11) still has a firm physical meaning if $r < 0$. Then, the mass does not bounce back, but penetrates the region $u > 0$ and instantaneously slows down upon passing the boundary $u = 0$. For $r = -1$ the obstacle is totally ignored, and we are effectively dealing with an unrestricted oscillator. The latter case can be conveniently used to check the formulae to be derived in the following.

Any realistic description of collisions should involve energy losses and (11) is the simplest thinkable model. We realize, however, that the choice (11) induces a discontinuity of the velocity, other than the one which is inherent to the impact. In our analysis these two sources of discontinuity remain entangled most of the time. As there are many ways to account for energy loss during collisions, the choice of energy loss models must be guided by situations of
practical interest. We believe that the model which is analysed here is relevant for macroscopic experiments (de Weger et al 2000). In another collision model studied by us in the context of atomic force microscopy (van de Water and Molenaar 2000), the energy loss is continuous. However, it appeared that the singularity structure of the dynamical system survived. The question of why this is so remains a subject of further study.

## 4. Instantaneous impacts

Since the dynamics between impacts (via (8)) and during impacts (via (10)) is explicitly known, the orbits can, in principle, be calculated. The dynamical equations (8) and (10) can be reduced to a discrete mapping for orbits which undergo near-grazing collisions. The reduction is possible for orbits which remain close to the particular solutions $p(t)$ in (4) with $\sigma \approx 0$. Possible impacts then occur in short time windows around $t = n \in \mathbb{Z}$. In the literature (Whiston 1987, 1992, Foale and Bishop 1992, Budd and Dux 1994, Budd et al 1995) one usually studies mappings which connect states of the system at successive impacts. We found it more convenient to focus on the Poincaré map connecting the states $u_n = (u_n, v_n)^T$ at stroboscopic times $t = n \in \mathbb{Z}$.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** The impact map $\mathcal{L}$. (a) If impacts occur before $t = t_n$, $\mathcal{L}$ maps the virtual position $u_n$ onto the true position $\hat{u}_n$. (b) If impacts occur after the stroboscopic time, $\mathcal{L}$ maps the true state $u_n$ onto the virtual state $\hat{u}_n$.</figcaption></figure>

A judicious choice of the state is a crucial step in the present derivation. For $u_n$ and $v_n$ we take the position and velocity the system would have at $t = n$ if the obstacle had been removed shortly before $t = n$. This is most conveniently understood from figure 3. If the impact takes place (shortly) after $t = n$ or if no impact occurs at all around $t = n$, the state involves the real values of position and velocity. If an impact occurs (shortly) before $t = n$, $u_n$ is a virtual state with $u_n > 0$. We account for a possible impact near $t = n$ through the local mapping $\mathcal{L}$:

$$ \hat{u}_n = \mathcal{L}(u_n). \qquad (12) $$

The action of $\mathcal{L}$ and the location of the points $u_n$ and $\hat{u}_n$ are illustrated in figure 3. The key idea is to calculate the orbit around $t = n$ with and without accounting for the presence of the obstacle. If an impact occurs before $t = n$ (figure 3(a)), $u_n$ represents the virtual state at $t = n$ obtained if no obstacle were present, while $\hat{u}_n$ represents the physical state obtained by incorporating the impact. If the impact takes place after $t = n$, the roles of $u_n$ and $\hat{u}_n$ are reversed, as shown in figure 3(b). If there is no impact or if the impact is grazing, $\mathcal{L}$ is the identity.

After application of $\mathcal{L}$, the system is propagated from state $\hat{u}_n$ at $t = n$ to state $u_{n+1}$ at $t = n + 1$ without accounting for the obstacle, thus by applying the evolution given by (8).
The full mapping is then given by

$$u_{n+1} = \mathbf{P}[\mathcal{L}(u_n) - \mathbf{p}(t_n)] + \mathbf{p}(t_{n+1}). \quad (13)$$

From (4) we have $\mathbf{p}(t_n) = \mathbf{p}(t_{n+1}) = (p(n), \dot{p}(n))^T = \sigma(1, 0)^T \equiv \sigma e_1$.

Using a parabolic approximation of the orbits around $t = n \in \mathbb{Z}$, one finds that an impact will take place if the quantity

$$\chi_n = v_n^2 - 2A u_n \quad (14)$$

is positive. Here, $A$ is the acceleration of the grazing orbit, i.e. $\ddot{p}(n)$ with $\sigma = 0$:

$$A = -4\pi^2. \quad (15)$$

If $\chi_n > 0$ and $u_n > 0$, an impact has occurred shortly before $t = n$, whereas if $\chi_n > 0$ and $u_n < 0$ the oscillator will hit it shortly after $t = n$. Both cases can be dealt with on an equal footing.

We now expand the dynamics of near-grazing orbits in terms of the small parameters $u_n$, $v_n$ and $\sigma$, assuming that they are all of the same order of magnitude. We are interested in a derivation up to and including the linear order, so terms involving (cross) products are omitted<fnref n="1" />. To find $\mathcal{L}$ explicitly, we express $\hat{u}_n$ in terms of $u_n$. In the derivation we use the impact characteristics $(v_c, t_c)$ as intermediates. To find approximations for $(v_c, t_c)$ that are accurate up to linear order in $(u_n, v_n, \sigma)$, the orbit $u(t)$ is expanded up to and including third order in time:

$$u(t) = u_n + v_n t + \frac{1}{2} A_n t^2 + \frac{1}{6} B_n t^3.$$

Here, we conveniently use $t = n$ as the origin of the time axis and introduce the notation $A_n = \ddot{u}(t = n)$ and $B_n = d^3u/dt^3(t = n)$. Setting $u(t_c) = 0$ we obtain the expression

$$t_c = \frac{1}{A_n} \left( \sqrt{\chi_n} - v_n + \frac{B_n}{3A_n} u_n \right). \quad (16)$$

Substitution of $t_c$ into the velocity $v \equiv \dot{u}$ yields

$$v_c = \sqrt{\chi_n} - \frac{2B_n}{3A_n} u_n. \quad (17)$$

From these expressions we see that in lowest order, $t_c$ and $v_c$ scale with $\sqrt{u_n}$. This is due to the parabolic character of the orbit near its maximum. According to (3) the acceleration of the oscillator is given by

$$\ddot{u} = -vv - \Omega^2(u+1) + (1+\sigma)F_g \cos(2\pi t + \phi) \quad (18)$$

which is a function of $(u, v, \sigma, t)$. Substituting $t = n$ we find that, in the linear approximation, $A_n$ can be written as

$$A_n = A - \Omega^2 u_n - vv_n + (\Omega^2 - 4\pi^2)\sigma \quad (19)$$

with the ‘grazing’ acceleration $A$ introduced in (15). In a similar way we obtain

$$\frac{d^3 u}{dt^3} = -v\ddot{u} - \Omega^2 v - 2\pi(1+\sigma)F_g \sin(2\pi t + \phi). \quad (20)$$

<footnote n="1">¹ Other assumptions about the relative ordering of $u_n$, $v_n$ and $\sigma$ are possible (de Weger 1999), but they do not affect our main conclusions.</footnote>
Substitution of (18) and $t = n$ yields

$$ B_n = \nu\Omega^2 u_n + (\nu^2 - \Omega^2)v_n - \nu\Omega^2\sigma. \quad (21) $$

We conclude that in (16) and (17) $A_n$ and $B_n$ can be replaced by their ‘grazing’ values $A$ and $B \equiv 0$, since terms of order higher than linear can be omitted:

$$ t_c(u_n, v_n) = \frac{1}{A}(\sqrt{\chi_n} - v_n) \qquad v_c(u_n, v_n) = \sqrt{\chi_n}. \quad (22) $$

Next, we expand $u(t)$ around the (still unknown) quantities $(\hat{u}_n, \hat{v}_n)$:

$$ u(t) = \hat{u}_n + \hat{v}_n t + \frac{1}{2} \hat{A}_n t^2 + \frac{1}{6} \hat{B}_n t^3. \quad (23) $$

Substituting $t = t_c$ we obtain the equation

$$ u(t_c) = 0 = \hat{u}_n + \hat{v}_n t_c + \frac{1}{2} \hat{A}_n t_c^2 + \frac{1}{6} \hat{B}_n t_c^3. \quad (24) $$

For the velocity we similarly find

$$ v(t_c) = -r v_c = \hat{v}_n + \hat{A}_n t_c + \frac{1}{2} \hat{B}_n t_c^2. \quad (25) $$

We note that $\hat{A}_n$ and $\hat{B}_n$ are given by expressions (19) and (21) after replacement of $(u_n, v_n)$ by $(\hat{u}_n, \hat{v}_n)$. If we substitute these expressions into (24) and (25) we meet with equations which are linear in $(\hat{u}_n, \hat{v}_n)$. After solving these, we obtain $(\hat{u}_n, \hat{v}_n)$ in terms of $(v_c, t_c)$. Substitution of (22) then yields $(\hat{u}_n, \hat{v}_n)$ in terms of $(u_n, v_n)$:

$$ \hat{u}_n = -(1+2r)u_n - \frac{1}{A}(r+1)v_n \left(v_n - \sqrt{\chi_n}\right) \quad (26) $$

and

$$ \hat{v}_n = v_n - (r+1) \left(\sqrt{\chi_n} - 2\nu u_n\right). \quad (27) $$

The latter two expressions describe the action of the map $\mathcal{L}$ introduced in (12). The map trivially reduces to the identity for the case $r = -1$ (no obstacle). Note that the second term on the right-hand side of (26) may be of linear order, although it contains products. This is the case if $v_n \gg u_n$, for which (26) reduces to $\hat{u}_n \approx -r u_n$. We emphasize that the expressions also yield the correct answer if the impact occurs after $t = n$. The second term on the right-hand side in (26) then ensures that $|\hat{u}_n| < |u_n|$, as is to be expected.

It is convenient to write the local map $\mathcal{L}$ in matrix notation:

$$ \mathcal{L}(u_n) = Nu_n + n(u_n) \quad (28) $$

with

$$ N(t) = \begin{pmatrix} -(1+2r) & 0 \\ 0 & 1 \end{pmatrix} \quad (29) $$

and

$$ n(u_n) = -(1+r) \left( \frac{1}{A} v_n \left( v_n - \sqrt{\chi_n} \right) \middle/ \sqrt{\chi_n} - 2\nu u_n \right). \quad (30) $$

The nonlinearity of $\mathcal{L}$ is contained in the vector $n$. Substitution of (28) into (13) yields for $\chi_n \ge 0$

$$ u_{n+1} = P(Nu_n + n) + \sigma(I - P)e_1 \quad (31) $$
while for $\chi_n \le 0$ we have the no-collision case

$$u_{n+1} = P u_n + \sigma (I - P) e_1. \qquad (32)$$

With equations (31) and (32) we have brought the dynamics of the harmonic impact oscillator to the form of a mapping. The mapping is approximate in that we have assumed near-grazing orbits. The mapping is not yet very transparent: it still lacks an impact criterion which depends on a single coordinate (such as the map (1) and (2)). In the next section we will reshape the mapping via a nonlinear transformation in order to obtain similarity with the map (1) and (2). It must be realized that this procedure might reduce the accuracy of the map (31) and (32) which was derived as a linear approximation to the near-grazing dynamics.

## 5. Natural coordinates

The map (31) and (32) can be put into a more concise form by realizing that the factor $\chi_n$, which we introduced in (14) and which discriminates between impacting and non-impacting orbits, plays a crucial role. The curve $\chi_n = 0$ in phase space represents the grazing orbit in the vicinity of $t = t_n$. This suggests choosing the family of curves given by $\chi_n = c$ with $c$ constant as level lines of the new coordinates:

$$x = l \chi. \qquad (33)$$

Here and in the following equation the constant factor $l > 0$ is introduced for later convenience. For the other new coordinate a convenient choice is the following linear combination of $\chi$ and $v$:

$$y = -l [P_{22} \chi + 2A P_{12} v + 2A \sigma (P_{22} - \gamma)]. \qquad (34)$$

This combination is chosen such that the second component of the two-dimensional map (31) and (32) gets the simple form given in (37) and (38). In figure 4 the level lines of the new coordinates $(x, y)$ are sketched. Inverting this transformation and retaining only terms linear in $(x, y, \sigma)$, we obtain the expressions

$$u = \frac{x}{l} \qquad v = -\frac{1}{2lA P_{12}} (y + P_{22} x + 2A l \sigma (P_{22} - \gamma)). \qquad (35)$$

The inverse transformation does not exist if $P_{12} = 0$, but this case corresponds to critical damping and can be trivially accommodated by adjusting the propagator $P(t)$ in equation (7).

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Level lines of the new x-coordinate, defined in (34). The level lines of the new y-coordinate (34) have a similar shape, but with respect to a shifted origin.</figcaption></figure>
Applying (33) and (34) at $t = n + 1$ we obtain $(x_{n+1}, y_{n+1})$ as functions of $(u_{n+1}, v_{n+1})$. The latter are expressed in terms of $(u_n, v_n)$ via the map (31) and (32). Application of the transformation (35) yields the map in transformed coordinates. As long as $r \neq -1$, which corresponds to the absence of the obstacle, an appropriate choice for the factor $l$ is

$$l = \frac{1}{4A^2(1+r)^2 P_{12}^2} \qquad (36)$$

After lengthy algebra, in which terms of order higher than linear in $(x, y, \sigma)$ are omitted, we find that the transformed map reads as

$$\left. \begin{aligned} x_{n+1} &= \alpha x_n + y_n + \rho \\ y_{n+1} &= -\gamma x_n \end{aligned} \right\} \quad \text{if} \quad x_n \le 0 \qquad (37)$$

$$\left. \begin{aligned} x_{n+1} &= -c_1 \sqrt{x_n} + c_2 x_n + y_n + \rho \\ y_{n+1} &= c_3 x_n \end{aligned} \right\} \quad \text{if} \quad x_n \ge 0. \qquad (38)$$

The coefficients $c_1, c_2$ and $c_3$, and the bifurcation parameter $\rho$, are given by

$$\begin{aligned} c_1 &= \operatorname{sign}(P_{12}) \\ c_2 &= \alpha - 2(1+r)P_{22} + (1+r)^2 P_{22}^2 \\ c_3 &= (1+2r)\gamma - (1+r)^2 P_{22}^3 \\ \rho &= \frac{1-\alpha+\gamma}{2|A|(1+r)^2 P_{12}^2} \sigma. \end{aligned} \qquad (39)$$

If $r = -1$, the obstacle is ignored and the dynamics is described by the map in (37) only, irrespective of the sign of $x_n$. The association of the parameters $\alpha$ and $\gamma$ in the collisionless map (37) with the parameters of the ordinary differential equation (3) can also be made by realizing that the eigenvalues of its Jacobian must be the same as those of the linear propagator $P$, as was also noted by Chin et al (1994). The form of the remaining parameters, however, can only be learned from a full nonlinear analysis as is done here.

In comparison with the Nordmark map (1) and (2), our map has several new features. First, we emphasize that the sign of the square-root term in equation (38) depends on the linear propagator $P_{12}$ and is not fixed, as in (2). We will explain in section 8 that a fixed negative sign prohibits period-1 maximal periodic orbits of the underdamped oscillator, in clear contradiction with both experiment and numerical simulation. In de Weger et al (2000) we show that the term proportional to $x_n$ in equation (38) is essential for a quantitative understanding of the destruction of maximal periodic orbits through an additional impact. This term is absent in the Nordmark map (1) and (2). Finally, we point out that the $c_3$ parameter differs from the corresponding prefactor in (2). In the next section we will demonstrate that with the groundwork now laid out, the derivation of the map for the case of non-instantaneous impacts can be done readily.

## 6. Non-instantaneous impacts

If the obstacle is not rigid but compliant, the impacting mass will penetrate the region $u > 0$ for some time interval, $\Delta t$ say. For $u < 0$ the motion of the oscillator very near collision is determined by the constant acceleration $A$ defined in (15). For $u > 0$ the motion is described by the dimensionless equation

$$\ddot{u} + v\dot{u} + \Omega^2(u+1) + \kappa^2 u = F \cos(2\pi t + \phi). \qquad (40)$$
The wall repels the mass harmonically with spring constant $κ^2$. The stiffness of the wall is gauged by the ratio $Γ = κ^2/Ω^2$ of the spring constants of wall and oscillator. An infinitely stiff wall has $Γ = ∞$. As was argued in section 3, we admit that energy may be instantaneously absorbed by the wall at the moment of collision. The collision velocity $v_c$ is thus reduced by a factor $r$ with $0 ≤ r ≤ 1^2$. Apart from this possibility of energy dissipation, energy may also be lost through the linear friction term $vú$ in (40). The solution of (40) with initial conditions $u(t_c) = 0$ and $v(t_c) = r v_c$ is known explicitly and the duration of the impact follows from the condition

$$u(t_c + \Delta t) = 0. \qquad (41)$$

Whilst $Δt$ as a function of $v_c$ at a given value of $κ$ may easily be computed numerically from (41), more insight is gained from approximate analytical expressions for $Δt$. These are most appropriately found by exploiting the physical characteristics of the system. Since the impacts are near grazing, the friction term $vu$ in (40) can be neglected in the first instance. As long as the wall is relatively stiff, we have $Γ ≫ 1$ and thus $κ^2 ≫ Ω^2$. Furthermore, we have $Δt ≪ 1$ and the stiffer the wall, the shorter the impact will be. This suggests expanding $Δt$ in powers of $κ^{-1}$, with the result

$$Δt = \frac{1}{κ} \arccos \frac{1 - β^2}{1 + β^2} + O(1/κ^2) \qquad (42)$$

with the dimensionless factor $β$ being given by

$$β = \frac{r v_c κ}{|\tilde{A}|}. \qquad (43)$$

The factor $\tilde{A}$ is defined as $\tilde{A} = (1 + σ)A + σΩ^2$, with the reduced excitation amplitude $σ$ defined in (5) and $A$ in (15). From (42) and (43) it is seen that the product $v_cκ$ determines the impact characteristics. In figure 5 $Δt$ is drawn as a function of $β$. Note that the time delay behaves linearly for low collision velocities, but converges to a constant limiting value for increasing values of $v_c$.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** The general behaviour of the residence time $Δt$ as a function of $β = rv_cκ/|\tilde{A}|$. For small $β$ the time delay grows linearly with $β$ and the impacts are soft. For large $β$ the time delay is constant and the impacts are hard. In both limits a map with a square-root singularity is derived.</figcaption></figure>

In the following we shall apply the result in (42) for the two limiting cases $β ≪ 1$ and $β ≫ 1$. We refer to these as low- and high-velocity impacts, respectively. Whether the collision velocity is low or high, not only depends on the collision velocity $v_c$, but also on the wall stiffness $κ$ and on the restitution coefficient $r$. As the magnitude of the collision velocity

<footnote n="2">² So, the minus sign present in equation (11) is not relevant in the case of a compliant wall. This implies that the role of the factor $r$ is slightly different now.</footnote>
depends on the dynamical state of the impact oscillator, it is not possible to give a *priori* criteria for wall hardness, for example as a criterion on the value of $\Gamma$. From an experimental point of view, an extremely interesting question is how the softness of the wall affects the bifurcation scenarios which are characteristic for impact oscillators.

## 6.1. Low-velocity impacts

For low-velocity impacts, we assume $\beta \ll 1$, i.e.

$$v_c \ll \frac{|\mathcal{A}|}{r \kappa}. \qquad (44)$$

To lowest order in $\beta$, the duration of the impact then is

$$\Delta t = \frac{2r}{|\mathcal{A}|} v_c = \frac{2\beta}{\kappa}. \qquad (45)$$

The stiffness of the wall has now dropped from (45), the orbit does penetrate the wall, but the time delay $\Delta t$ could have been obtained by setting $\kappa$ to zero in (40). The only remaining wall parameter is the restitution coefficient $r$. For purely elastic low-velocity collisions ($r = 1$), therefore, the presence of the wall can be ignored completely in our approximation. The dynamics is then effectively also free in the region $u > 0$ and is described by the regular map (37) everywhere. Our analysis does not go beyond terms which are first order in $x$, and in this case the presence of the wall will only be felt in the next higher order $x^{3/2}$. That this is indeed the case was shown by Dankowicz and Nordmark (2000).

To within the order of our approximation, the obstacle is only felt when it absorbs energy, thus if $r < 1$. To find the mapping for this case we derive the analogues of (26) and (27) starting from the impact law (10) with (45). To that end the derivation steps (23)–(25) have to be repeated. We eventually find that the local map $\mathcal{L}$ is still of the form (28)–(30), but with the replacement $r \to -r$. Thus in this case we have

$$\mathcal{L}(u_n) = \mathbf{N}u_n + n(u_n) \qquad (46)$$

with

$$\mathbf{N}(t) = \begin{pmatrix} -(1-2r) & 0 \\ 0 & 1 \end{pmatrix} \qquad n(u_n) = -(1-r) \begin{pmatrix} \frac{1}{\mathcal{A}} v_n (v_n - \sqrt{\chi_n}) \\ \sqrt{\chi_n} - 2\nu u_n \end{pmatrix}. \qquad (47)$$

Since equations (28)–(30) have the same structure as (46) and (47), the remainder of the evaluation does not change. In conclusion, if $r = 1$ the dynamics is given by the map in (37), irrespective of the sign of $x_n$. If $r \neq 1$, the map is given by equations (37) and (38) with the coefficients

$$\begin{align}
c_1 &= \operatorname{sign}(\mathbf{P}_{12}) \nonumber \\
c_2 &= \alpha - 2(1-r) \mathbf{P}_{22} + (1-r)^2 \mathbf{P}_{22}^2 \nonumber \\
c_3 &= (1-2r) \gamma - (1-r)^2 \mathbf{P}_{22}^3 \tag{48} \\
\rho &= \frac{1-\alpha+\gamma}{2|\mathcal{A}|(1-r)^2 \mathbf{P}_{12}^2} \sigma. \nonumber
\end{align}$$

It is a highly remarkable conclusion that for $r < 1$ the impact map is *the same* as the map for an infinitely stiff wall, but for a sign change of $r$. The consequence is that soft collisions with a soft wall at $r < 1$ should display *the same* bifurcation phenomena as collisions with a stiff wall,
but for a scale change $(1+r)^2/(1-r)^2$ of the bifurcation parameter $\rho$. The renormalization of the bifurcation parameter for soft collisions with energy loss is an important result of our analysis in which the discontinuity of the vector field remains intertwined with that of the velocity.

## 6.2. High-velocity impacts

High-velocity impacts correspond to the case $\beta \gg 1$, i.e.

$$ \frac{r \kappa v_c}{|A|} \gg 1. \qquad (49) $$

From equation (42) we conclude that for large enough $\beta$ the impact duration $\Delta t$ is given approximately by

$$ \Delta t \approx \frac{\pi}{\kappa} \qquad (50) $$

and is thus independent of the collision velocity $v_c$. This fact has interesting consequences. Combining equations (49) and (50) we find that for $\beta \gg 1$

$$ \Delta t \ll \frac{r \pi v_c}{|A|}. \qquad (51) $$

Repeating steps (23)–(25), but with the collision rule

$$ t_r = t_c + \frac{\pi}{\kappa} \qquad v_r = -r v_c \qquad (52) $$

we again arrive at (38) with (39), provided that terms which are much smaller than the terms already included in (38) are ignored. Effectively, we are then dealing with the case $\Delta t = 0$. This implies that, in leading order, high-velocity impacts at a compliant wall are described by the same map as instantaneous impacts at a rigid wall, i.e. by equations (37)–(39).

## 6.3. Summary and discussion

From the previous section we conclude the persistence of the square-root singularity in the mapping in almost all cases, except that of low-velocity, perfectly elastic collisions with a yielding wall. These elastic collisions highlight the discontinuity of the vector field, which for compliant walls was shown to involve an $x^{3/2}$ instead of a square-root behaviour of the map (Dankowicz and Nordmark 2000). In contrast, we find an $x^{1/2}$ behaviour for high-velocity impacts ($\beta \gg 1$). Using numerical simulations we will now show how this apparent paradox is resolved. The simulations are based on the ordinary differential equation<fnref n="3" />

$$ \ddot{u} + v\dot{u} + K(u) = F \cos(2\pi t + \phi) \qquad (53) $$

where

$$ K(u) = \begin{cases} \Omega^2(u+1) & u \le 0 \\ \Omega^2(u+1) + \kappa^2 u & u > 0 \end{cases} $$

<footnote n="3">³ An efficient numerical scheme for solving equation (53) in the presence of grazing impacts uses the analytical solutions for the motion between impacts. The (exact) positions $u(t)$ are computed in a small number of discrete points $t_1, \dots, t_k$ in each drive period. It is crucial not to miss boundary crossings of $u(t)$. These crossings are detected both directly by checking $u(t_1), \dots, u(t_k)$ and by computing the position $u(t)$ at the turning points of the velocity $\dot{u}$. The number $k$ of discrete points in the regions $u < 0$ and $u > 0$ is taken as being proportional to $\Omega$ and $\kappa$, respectively.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** (a) The (12) element of the Jacobian of a numerically computed map of an underdamped oscillator ($\nu, \Omega$) = (0.196, 7.368) that collides elastically ($r = 1$) with a yielding wall for stiffness ratios $\Gamma = 10, 10^2, 10^3$ and $10^4$, respectively. It was found from integrating the differential equation (53). (b) Bifurcation diagram computed from (53) for $\Gamma = 10$. (c) Same as (b), but now for $\Gamma = 1000$. The period-3 to period-1 transition at decreasing values of $\sigma$ is characteristic of a square-root singularity of the Jacobian.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** (a) Jacobian of a numerically computed map of an underdamped oscillator ($v, \Omega$) = (0.196, 7.368) that collides inelastically ($r$ = 0.9) with a yielding wall for stiffness ratios $\Gamma$ = 10, $10^2$, $10^3$ and $10^4$, respectively. It was found from integrating the differential equation (53). (b) Bifurcation diagram computed from (53) for $\Gamma$ = 10. (c) Same as (b), but now for $\Gamma$ = 1000. The horizontal scale of (b) has been expanded by a factor of $10^2$. The backward transition period-3 to period-1. The hysteretic period-3 to period-1 transition at decreasing values of $\sigma$ that is recognized in both cases is characteristic of a square-root singularity of the Jacobian.</figcaption></figure>
which combines (3) and (40). For any derived mapping, the central question is if it can
correctly predict bifurcation scenarios which are of practical relevance. One of those is the
phenomenon of period adding which is intimately tied to the square-root singularity of the
Jacobian. Inspired by the experiment of de Weger *et al* (2000) we have chosen the parameters
such that the dynamical state is a period-1 to period-3 transition which displays hysteresis.
In a quasi-stationary upward scan of *F* from the non-impacting state, impacts occur first at
*F* = *F*<sub>g</sub> (σ = 0). When σ is subsequently decreased, the non-impacting state is reached only
at a negative value σ = -σ<sub>H</sub>. At this point the orbit is closest to grazing and we take the
downward transition period-3 → period-1 as indicative of an effective square-root behaviour
of the map.

The singularity structure of the map is most clearly seen in its Jacobian. Because the
numerical map is not in natural coordinates, all elements of the Jacobian display the same
qualitative behaviour and we arbitrarily select the element $J_{12} = \partial u_{n+1}/\partial \dot{u}_n$. In figure 6(a) we
show the Jacobian as a function of $x = u_n$ of the stroboscopic map that was computed from
(53) for the case of perfectly elastic collisions ($r = 1$). For fairly compliant walls ($\Gamma = 10$), the
$x^{3/2}$ behaviour of the map can be clearly recognized. Although the Jacobian is never singular
at $x = 0$, it starts to develop an *effective* $x^{-1/2}$ behaviour for increasing $\Gamma$. This is in complete
agreement with our analysis and resolves the apparent paradox. Accordingly, for collisions
that are hard enough ($\beta \gg 1$), the bifurcation scenario is that of the square-root map. This
is illustrated in figure 6(c) that shows the bifurcation diagram for $\Gamma = 10^3$, where indeed the
characteristic period-3 → period-1 transition can be seen for decreasing $\sigma$, whereas a different
scenario is observed at $\Gamma = 10$ in figure 6(b). In figure 6(a) the parameter which gauges the
hardness of the collision $\beta = 21$ at $x = 0.02$ for $\Gamma = 10^3$, but it must be emphasized that
$\beta$ is not uniformly large as it varies as $x^{1/2}$. The single parameter $\beta$, therefore, controls the
cross-over between the two behaviours that were recognized by Dankowicz and Nordmark
(2000).

For collisions with energy loss ($r = 0.9$) the numerically computed map is shown in figure 7(a). In agreement with our analysis, the Jacobian now displays an $x^{-1/2}$ behaviour both for large and small $\beta$, but not for intermediate values. This is illustrated in the bifurcation diagrams of figure 7(b) for $\Gamma = 10$ and $10^3$. In both cases the period-3 to period-1 transitions which are typical for a square-root map can be recognized. As in the elastic case, we must realize that for $\Gamma = 10^3 \beta$ is not uniformly large.

We recall that both cases are described by the same map, but for a sign reversal of $r$
in the coefficients. This leads to a scale change of $\sigma$ with a factor of $(1+r)^2/(1-r)^2 \approx 10^2$. The $\sigma$-axis of the bifurcation diagram at $\Gamma = 10$ (figure 7(b)) has been scaled accordingly.

For inelastic collisions $r < 1$ the discontinuity of the vector field and the jump of
the velocity upon collision are interleaved, such that for $\beta \gg 1$ the map reflects that of
the vector field only and for $\beta \ll 1$ the singularity is due to the velocity change upon
collision. It is a remarkable finding that both cases are related through the scale change
$(r-1)^2/(r+1)^2$. Incidentally, as can be seen from the behaviour at $\beta = O(1)$, a
restitution-rule impact energy loss does not necessarily imply a square-root singularity of
the Jacobian.

**7. Period-adding bifurcations**

As mentioned in the introduction, one of the characteristic features of maps with a square-root
singularity is the phenomenon of period-adding bifurcations. The emergence of this type of
bifurcation can be understood in an extremely simple way from the map in the limit of large damping, $\nu \gg 1$. However, it will appear that this simplicity comes at a price, namely that the number of addings is not as large as those in figure 1 which was computed from the mapping (37), (38) and (48) $s_1 \approx 0$, $s_2 \approx -\nu$, and accordingly $\alpha \approx 1$, $\gamma \approx 0$, $P_{11} \approx 1$, $P_{12} \approx 1/\nu$ and $P_{22} \approx 0$. The two-dimensional map then tends to the one-dimensional map

$$x_{n+1} = x_n + \rho \quad \text{if} \quad x_n \le 0. \qquad (54)$$

$$x_{n+1} = -\sqrt{x_n} + x_n + \rho \quad \text{if} \quad x_n \ge 0 \qquad (55)$$

with $\rho$ given by

$$\rho \approx \frac{\nu\Omega^2}{2|A|(1-r)^2} \sigma. \qquad (56)$$

The influence of the restitution coefficient $r$ is now only felt through the scaling of the bifurcation parameter. The reduced map resembles those analysed by Nusse and Yorke (1992), Nusse et al (1994) and Lamba and Budd (1994).

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Sketch of a maximal period-$M$ orbit with $M = 6$. At $t = M - 1$ the orbit grazes the wall (position $u > 0$), so that this orbit will be destroyed as soon as the bifurcation parameter is increased further.</figcaption></figure>

The nature of maximal periodic orbits is illustrated in figure 8. Starting at some $x_1 < 0$, the orbit creeps in the direction of the origin with steps of size $\rho$ according to (54) until it collides and is thrown back onto the negative x-axis according to (55). Maximal M-periodic orbits are stable only in a given interval of the bifurcation parameter $\rho$. For too large $\rho$, a given M-periodic orbit is destroyed by an additional impact at time $M - 1$, whereas for too small $\rho$, the cycle loses stability as the orbit is not sufficiently thrown back at impact. We will now compute the stability interval $[\rho_M^L, \rho_M^U]$.

The cycle elements $x_1 \dots x_M$ of a maximal periodic orbit satisfy

$$x_1 < x_2 < \dots < x_{M-1} < 0 < x_M \qquad (57)$$

with $x_M$ mapped onto $x_1$. Since $x_{M+1}$ is found by applying (54) $M - 1$ times and (55) once,
the orbit satisfies

$$ 
\begin{aligned}
x_{M+1} &= -\sqrt{x_M} + x_M + \rho \\
&= -\sqrt{x_1 + (M-1)\rho} + x_1 + (M-1)\rho + \rho.
\end{aligned} 
\quad (58) $$

At the upper boundary $\rho_M^U$ of its existence interval, the M-periodic orbit grazes the wall at $t = M - 1$, as is illustrated in figure 8. Then we have

$$ x_1 = -(M-2)\rho_M^U \quad x_{M-1} = 0 \quad \text{and} \quad x_M = \rho_M^U. \quad (59) $$

From the condition $x_{M+1} = x_1$ it then follows that

$$ \rho_M^U = \frac{1}{M^2}. \quad (60) $$

It can be seen that maximal periodic orbits with increasing periods M live at decreasing values of the bifurcation parameter $\rho$.

The lower bound of $\rho_M^L$ of the $\rho$-interval of existence of the maximal period M orbit follows from the stability criterion

$$ \left| \frac{\partial x_{M+1}}{\partial x_1} \right| \le 1. \quad (61) $$

Using equation (58), the requirement becomes

$$ \left| -\frac{1}{2(x_1 + (M-1)\rho)^{1/2}} + 1 \right| < 1. \quad (62) $$

With the period condition $x_{M+1} = x_1$ we finally have

$$ \rho_M^L = \frac{1}{4M}. \quad (63) $$

In conclusion, for the strongly damped oscillator, maximal period-M orbits exist in parameter intervals

$$ \rho_M \in \left[ \frac{1}{4M}, \frac{1}{M^2} \right]. \quad (64) $$

Clearly, as M increases, the lower bound overtakes the upper bound and we conclude that maximal periodic orbits only exist up until a maximum period $M_{max} = 3$.

Our conclusion is beautifully illustrated in figure 9, where we compare iterations of the map with numerical simulations of the differential equation (53) for collisions with a yielding wall with stiffness ratio $\Gamma = 10$. In both cases a truncated period-adding sequence is seen with stability intervals which agree well with the prediction (64). A similar comparison was done for collisions with an infinitely hard wall. The appropriate map is the same as (54) and (55), but for a sign reversal of $r$, which only affects the scale of the bifurcation parameter $\rho$ (56). At the restitution coefficient used is $r = 0.9$ in figure 9, the scale change is $O(10^2)$. In the differential equation this change moves the windows of maximal periodic orbits to values of $\sigma$ which are so large that the map is no longer a valid approximation. Accordingly, no windows of periodic orbits were observed in the simulations of equation (53). The striking agreement between the simple mapping and the simulation results for the differential equation illustrates that allowance for a yielding wall is an essential ingredient of impact maps.

The map (37) and (38) reduces to a one-dimensional map if $\gamma = 0$, in all other cases the mapping is essentially two dimensional, which makes a bifurcation analysis much more
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Truncated series of inverse period addings for the strongly damped oscillator with ($\alpha, \gamma$) = (0.99, 0.01) and $r = 0.9$. (a) Iterations of the map (37) and (38) with (48) which is derived for low-velocity impacts with a yielding wall. The vertical lines indicate the stability intervals given by equation (64). (b) Integration of the differential equation (53) using a stiffness ratio $\Gamma = \kappa^2/\Omega^2 = 10$. The scale of the position $x$ is determined by the stroboscopic phase. In both cases the horizontal axis is the physical bifurcation parameter $\sigma = F/F_g - 1$.</figcaption></figure>
technical. The frictionless oscillator is studied in Lamba and Budd (1994). For the Nordmark map (1) and (2) an exhaustive analysis has been done by Chin et al (1994). We expect that many results of their work will carry over to our mapping. The characteristics of maximal periodic orbits in underdamped oscillators are compared with experimental results in de Weger et al (2000).

## 8. Conclusion

In this paper we have devised a mapping for orbits of a harmonic impact oscillator which are close to grazing. We have assumed that this occurs at small values of the control parameter $\sigma$ and small impact velocities. For period-adding bifurcations these quantities may not both be small at the same time. When in a quasi-stationary upward scan of $\sigma$ an orbit first starts impacting at $\sigma = 0$, it may do so at a rather high velocity. The impacts may remain at $\sigma < 0$ when $\sigma$ is subsequently lowered, only to vanish at some negative value $\sigma = -\sigma_H$. The transition therefore displays hysteresis. A quantitatively correct description of this phenomenon is a challenge for the derived mapping. An additional complication is that the collision velocity depends on the dynamical state of the impact oscillator, and it is difficult to decide a priori whether collisions are hard or soft. The predictive power of the maps (37)–(39) and (37), (38), (48) in situations of practical relevance is tested in de Weger et al (2000).

A notable difference between the Nordmark map (1) and (2) and the present one (37) and (38) is the presence of the sign factor in front of the square-root. In the overdamped case $\nu^2 - 4\Omega^2 > 0$ the sign factor has a constant positive sign, but it may change sign for underdamped oscillations since then

$$P_{12} = \gamma^{1/2} \frac{\sin \psi}{\psi} \quad \text{with} \quad \psi = \frac{1}{2} \sqrt{4\Omega^2 - \nu^2}.$$

We recall that the case $\psi = 2\pi$ corresponds to the situation where the oscillator is driven with a frequency that equals its own freely swinging frequency. It may readily be appreciated from figure 8 that maximal M-periodic orbits need to be thrown back sufficiently at impact, and therefore in general need a square-root term with a negative sign. The condition for the existence of maximal periodic orbits of period $M > 1$ is then $2\pi < \psi < 3\pi$ (modulo $2\pi$). It can also be appreciated that a map with a positive square-root term can only support period-1 impacting orbits. Further, it appears that the combination of a negative square-root sign with a negative $\alpha$ supports period-2 maximal periodic orbits. Since

$$\alpha = 2\gamma^{1/2} \cos \psi$$

period-2 orbits occur in $5\pi/2 < \psi < 3\pi$ (modulo $2\pi$). This simple observation agrees with both experiment and numerical simulation of the differential equation. We note that both the period-1 and the period-2 maximal periodic orbits were missed by Chin et al (1994) who considered only negative square-root terms and the case $\alpha > 0$.

A remarkable conclusion of this paper is that soft impacts with a yielding wall give rise to an impact map which differs from the hard-wall case, roughly through a rescaling of the bifurcation parameter $\rho$. For the compliant wall we have considered an admittedly simple collision model, but we believe that our main conclusion, namely the persistence of the square-root singularity in impact mappings, is robust.
Acknowledgments

We thank Arne Nordmark for valuable suggestions. We gratefully acknowledge financial
support by the ‘Nederlandse Organizatie voor Wetenschappelijk Onderzoek (NWO)’.

References

Berg J and Briggs G A D 1997 Nonlinear dynamics of intermittent-contact mode atomic force microscopy *Phys. Rev.* B **55** 14 899–908

Budd C J 1996 Non-smooth dynamical systems and the grazing bifurcations *Nonlinear Mathematics and its Applications* ed P J Aston (Cambridge: Cambridge University Press)

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

Budd C J, Dux F and Cliff A 1995 The effect of frequency and clearance variations on single-degree-of-freedom impact oscillators *J. Sound Vib.* **184** 475–502

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

Dankowicz H and Nordmark A B 2000 On the origin and bifurcations of stick-slip oscillations *Physica D* **136** 280–302

de Weger J 1999 The grazing bifurcation in impact oscillators *PhD Thesis* Eindhoven University of Technology

de Weger J, Binks D, Molenaar J and van de Water W 1996 Generic behaviour of grazing impact oscillators *Phys. Rev. Lett.* **76** 3951–4

de Weger J, van de Water W and Molenaar J 2000 Grazing impact oscillations *Phys. Rev.* E **62** 2030–41

Feigenbaum M J 1983 Universal behavior in non-linear systems *Physica D* **7** 16–39

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

Frederiksson M H 1998 Topics in impacting system dynamics *PhD Thesis* Royal Institute of Technology, Stockholm

Frederiksson M H and Nordmark A B 1997 Bifurcations caused by grazing incidence in many degrees of freedom impact oscillators *Proc. R. Soc. A* **453** 1261–76

Lamba H and Budd C J 1994 Scaling of Lyapunov exponents at nonsmooth bifurcations *Phys. Rev.* E **50** 84–90

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

—1997 Universal limit mapping in grazing bifurcations *Phys. Rev.* E **55** 266–70

Nusse H E, Ott E and Yorke J A 1994 Border-collision bifurcations: an explanation for observed bifurcation phenomena

*Phys. Rev.* E **49** 1073–6

Nusse H E and Yorke J A 1992 Border-collision bifurcations including period two to period three for piecewise smooth systems *Physica D* **57** 39–57

van de Water W and Molenaar J 2000 Dynamics of vibrating atomic force microscopy *Nanotechnology* **11** 192–9

Whiston G S 1987 Global dynamics of a vibro-impacting linear oscillator *J. Sound Vib.* **118** 395–429

Whiston G S 1992 Singularities in vibro-impact dynamics *J. Sound Vib.* **152** 427–60