## Abstract

Impact oscillators are non-smooth systems with such complex behaviours that their numerical treatment by traditional methods is not always successful. We design non-standard finite-difference schemes in which the intrinsic qualitative parameters of the system—the restitution coefficient, the oscillation frequency and the structure of the nonlinear terms—are suitably incorporated. The schemes obtained are unconditionally stable and replicate a number of important physical properties of the involved oscillator system such as the conservation of energy between two consecutive impact times. Numerical examples, including the Duffing oscillator that develops a chaotic behaviour for some positions of the obstacle, are presented. It is observed that the cpu times of computation are of the same order for both the standard and the non-standard schemes.

## 1. Introduction

Mechanical systems subject to impacts, and particularly impact oscillators, have become a major area of study because they arise in a variety of engineering and real-life problems such as rattling gears, vibrating absorbers, car suspension, impact print hammers, and so on. The main concerns are the excessive noise, wear and fatigue that these systems may induce. Thus, there is an essential need to study them properly in order to improve their performance.

We consider a one-dimensional mechanical system whose motion is described by the initial value-problem(1.1)The motion is subject to convex constraints, that is, *y*(*t*) remains in a convex set *K* of non-empty interior, which is assumed to be the interval *K*=[0,∞) in this introduction. Thus, *y*=0 represents the position of the obstacle. The loss of energy during each impact time *t* is modelled by a percussion law,(1.2)where *e*∈[0,1] is the restitution coefficient (cf. Brogliato 2000, p. 117). We also assume that the function *f* is uniquely decomposable as so that *ω* appears as the frequency of the oscillation. In order to include the effect of the normal reaction force that occurs during the collision between the system and the obstacle, it is convenient to write the problems (1.1) and (1.2) as a differential inclusion(1.3)where ∂*Ψ*_{[0,∞)}(*y*) is the subdifferential of the indicator function *Ψ*_{[0,∞)} of [0,∞) (cf. Goeleven *et al*. 2003).

In general, there are only a few cases for which analytical solutions to equations (1.1) and (1.2) or (1.3) can be found. Hence, the crucial need for numerical methods. Here are some of the methods. The first approach is to compute each time of impact and utilize classical schemes. However, this method works in the case when the impacts are isolated. Another method consists of approximating the solutions *y*_{λ} of a Yosida-like regularization of (1.3) at the cost of a severe restriction on the size of the time-step Δ*t* and of the parameter *λ*, viz. (e.g. Paoli 2001; Paoli & Schatzman 2002*a*). Our point of departure is an impact scheme studied by Paoli (1993, 2001) and Paoli & Schatzman (2002*a*,*b*), and referred to hereafter as PS scheme. In order to handle the constraint automatically without approximating the impact times, these authors propose a sequence of discrete solutions (*y*^{k}) given by(1.4)where *F* ^{k} is some approximation of *f*(*t*,*y*) at the time *t*_{k}=*k*Δ*t* or *t*_{k+1}.

The PS scheme is convergent. However, being based on the standard finite-difference method, it fails to replicate a number of essential physical properties of the impact problem between two consecutive impact times. In this paper, we fill this gap by using the non-standard finite-difference method. Since the pioneering works of Mickens in the mid-1980s, the non-standard approach has shown great potential in the design of reliable schemes that preserve significant properties of solutions of differential models in science and engineering (e.g. Mickens 1994, 1997, 2000; Al-Kahby *et al*. 2000; Anguelov & Lubuma 2001, 2003; Anguelov *et al*. 2003). Of particular interest are also the edited book of Mickens (2000), Mickens' special issue (Gumel 2003), as well as the papers by Gumel *et al*. (2003), Moghadas *et al*. (2003) and Anguelov *et al*. (2005), where real-life problems are handled by non-standard schemes. Unlike the results in these references, this work constitutes, to the authors' best knowledge, the first time that the non-standard approach is implemented for differential inclusions and for vibro-impact problems. The main feature of the non-standard method proposed in this paper is an explicit incorporation of the qualitative parameter *ω* of oscillation frequency as well as a suitable approximation, via *F* ^{k}, of nonlinear terms that occur in *f* (see definition 2.1 below). This yields the following new impact scheme:(1.5)

In the particular case of the harmonic oscillator (i.e. *f*(*t*, *y*)≡−*ω*^{2}*y*), the explicit solution of the impact problem in equations (1.1)–(1.3) with *y*_{0}=1 and *v*_{0}=0 is found to bewhere, for *n*=0,1,2,…, is the time when the impact number *n*+1 occurs. Assume that *t*^{*}=*t*_{k} is between two impact times, that is, *t*^{*}∈(*s*_{n},*s*_{n+1}) for some *n*. Then, with suitable initial value on this interval, the discrete solution obtained by the non-standard method (1.5) is such that *y*(*t*^{*})=*y*^{k}. On the contrary, the solution *y*^{k} given by the PS scheme (1.4) is not a good approximation of *y*(*t*^{*}) whenever the quantity *ω*Δ*t* is large. In fact, numerous illustrative numerical examples are provided throughout the relevant sections. In all the examples, we observe that the times of computation, for a given Δ*t*, are of the same order for both the non-standard and the standard schemes.

The paper is organized as follows. In §2, we give a short self-contained presentation of the non-standard finite-difference method. Section 3 is devoted to the design of energy-preserving non-standard finite-difference schemes for oscillator problems without impacts. The schemes obtained are also unconditionally stable in the sense of Lax–Richtmyer in the case of linear oscillators. A brief account on the theory of vibro-impact problems is presented in §4, and §5 presents the PS scheme and a new non-standard impact scheme. Section 6 is devoted to concluding remarks, which include possible extensions of this study.

## 2. The non-standard finite-difference method, quid?

The non-standard finite-difference approach is still in the early stages of its development. To make the paper relatively self-contained, in this section we give a short presentation of the non-standard finite-difference method. (For survey papers, we refer the reader to Mickens (1999, 2002) and to chapter 1 of Mickens (2000).) We consider the general case of an initial value problem for an autonomous first order differential equation,(2.1)where *T* can be infinite as in dynamical systems. We could, in fact, consider a system of *n* differential equations in *n* unknowns. However, we restrict ourselves to the scalar equation (2.1) to simplify the exposition. We assume that the problem (2.1) has a unique solution *y* (see standard references for sufficient conditions on *g*). In what follows, we assume implicitly that involved functions possess all needed smoothness properties.

For numerical purposes, we replace the interval (0,*T *) with the mesh {*t*_{k}=*k*Δ*t*}_{k≥0}, the parameter Δ*t*>0 being the time-step size. We denote, by *y*^{k}, an approximation to the solution *y* at the node *t*=*t*_{k}:(2.2)with superscripts referring to time-step, not to powers. The non-standard approach aims to design difference equations whose solution sequences (*y*^{k}) preserve essential physical properties of the exact solution *y*(*t*) without restriction on the size of Δ*t*. This philosophy clearly excludes the forward Euler scheme:(2.3)which is the oldest and simplest standard finite-difference method for (2.1). The ideal schemes would be the so-called exact schemes (Mickens 1994) or dynamically consistent schemes (Al-Kahby *et al*. 2000) defined as difference equations with solutions *y*^{k} satisfyinginstead of equation (2.2). Unfortunately, exact schemes of practical interest do not generally exist since few differential equations can be solved explicitly. It is normal practice in numerical analysis to define the quality of an algorithm by applying it to a test differential equation. For example, the absolute stability theory of the Euler method (2.3) and the exponentially fitted method for stiff systems are tailored from the decay equation (Lambert 1991). In line with this practice, we consider the model equation(2.4)where *m*≥1 is an integer and *λ*≠0 is a real number. (This reduces to the decay equation when *m*=1.) The solution of (2.4) at the time *t*=*t*_{k+1} isthis is an exact scheme of (2.4), which on setting *y*^{k} ≔*y*(*t*_{k}), may be rewritten as(2.5)

Mickens (1994, 2000) provided exact schemes for a large number of ordinary and partial differential equations. He singled out the important observation that the complex structure of the denominator of the discrete derivative in (2.5)_{1} and the non-local form of the approximation of the nonlinear term *λy*^{m} in (2.5)_{2} constitute a general property of these schemes. This observation motivates the following formal definition by Anguelov & Lubuma (2001) of Mickens' non-standard approach.

A difference equation to determine approximate solutions *y*^{k} to the solution *y*(*t*) of (2.1) is called a non-standard finite-difference method if at least one of the following conditions is met.

The classical denominator Δ

*t*of the discrete derivative*D*_{Δt}*y*^{k}is replaced by a non-negative function*ϕ*(Δ*t*) such that(2.6)Nonlinear terms that occur in

*g*(*y*) are discretized in a non-local way (i.e. by a suitable function of several points of the mesh).

The power of the non-standard finite-difference method is measured with the help of qualitative stability as depicted below (Anguelov & Lubuma 2001).

A difference equation in *y ^{k}* is called (qualitatively) stable with respect to some property

*P*of the differential equation (2.1) or of its solution

*y*whenever the discrete equation or its solution replicates the property

*P*for any value of Δ

*t*.

A point of order is necessary at this stage. Mickens (1994, p. 84) set five rules for the construction of discrete models that have the capability to replicate the properties of the exact solution. Our definition of a non-standard finite-difference scheme retains two only of Mickens' rules namely, parts (*a*) and (*b*) of definition 2.1. The other rules are expressed in terms of definition 2.2. For example, the schemes considered in this paper are stable with respect to the order of the involved differential equation, i.e. the difference equation and the differential equation have the same order.

At the end of this section, we show how to modify the forward Euler method (2.3) into a non-standard one, which is much better in the sense of definition 2.2. The simple examples (2.4) and (2.5) illustrate the need for the structure of the right-hand side of the differential equation to be intrinsically reflected in the discrete schemes if they are required to replicate the qualitative properties of the solution of the differential equation. In the general framework of equation (2.1), the properties of its solutions will be captured by a number(2.7)where traces all the fixed or critical points of the differential equation in (2.1). (We assume that the differential equation has a non-zero number of fixed points , which are all hyperbolic, that is, ; the situation of non-hyperbolic fixed points is handled by Anguelov & Lubuma (2003) and Anguelov *et al*. (2003).) In the spirit of (2.5)_{1} and definition 1.1, we renormalize the denominator of the discrete derivative in (2.3) through a non-negative function *ϕ* satisfying (2.6). This leads to the non-standard forward Euler method:(2.8)

Apart from being convergent and of order 1, as its classical analogue (2.3), the non-standard method (2.8) has the important qualitative property stated in the next result, which is proved by Mickens (1994) (for scalar equations) and Lubuma & Roux (2003) (for systems).

*Let ϕ in* *(2.6)* *be such that*(2.9)

*(A typical example is ϕ**(z)*=1−*e*^{−z}.*) Then the non-standard forward Euler method is elementary stable in the following sense*. *For every* Δ*t*>0*, the fixed points* *of the scheme* *(2.8)* *are exactly those of the differential equation in* *(2.1)*. *(This is the stability with respect to the number of fixed points in terms of* *definition 2.2.)* *Furthermore, each fixed point has the same Lyapunov stability*/*instability* *(cf*. *Stuart & Humphries 1998)* *for both the scheme* *(2.8)* *and the differential equation in* *(2.1)*. *(This is the stability with respect to the type of fixed points)*.

The order or the step one of the difference equation (2.8) coincides with the order of the differential equation (2.1). This requirement, which was mentioned earlier, is essential in the construction of non-standard schemes. Classical linear multi-step methods lack this property; it is shown in Anguelov & Lubuma (2001) that they are elementary unstable. It is also shown in this reference that the classical Runge–Kutta methods are elementary unstable. Furthermore, it was observed in Gumel *et al*. (2003) and Moghadas *et al*. (2003) that the classical Runge–Kutta method of order 4 fails to preserve the positivity property of solutions of epidemic models. We will therefore not consider classical schemes (for further specific motivation, see remark 5.1 and the comments in §5).

Each solution of the differential equation in (2.1) is either an increasing or a decreasing function on the interval [0,T). It follows from Anguelov & Lubuma (2003) that the non-standard forward Euler scheme (2.8) is stable with respect to the monotonicity property of solutions whenever for any Δ*t* and *y*. Of course, the classical forward Euler scheme fails to preserve all these properties.

## 3. Energy-preserving schemes

The physical problems investigated in this paper are in one way or another linked to conservative oscillators described by second-order differential equations of the form(3.1)

Equation (3.1) can be transformed into an equivalent first order autonomous system of two ordinary differential equations to which the analogue of the non-standard method in §2 could be applied. However, it is more interesting from the physical perspective to view this equation as a model of Hamiltonian systems in classical mechanics. Equation (3.1) is indeed equivalent to its first integral(3.2)with(3.3)where *H* represents the sum of kinetic energy and potential energy. Equation (3.2) is therefore a statement of conservation of energy. Our aim is to derive finite-difference methods, which are stable with respect to the principle of conversation of energy. From the analysis of §2 about the approximation of the derivative and with *ϕ* satisfying (2.6), we consider a discrete energy of the general form(3.4)where the discrete potential energy will be determined shortly. The requirement for discrete principle of conservation of energy means that(3.5)or, equivalently,(3.6)for *k*≥0. Using the well-known relation *a*^{2}−*b*^{2}=(*a*−*b*)(*a*+*b*), equation (3.6) takes the equivalent form(3.7)

At a fixed time *t*^{*}=*k*Δ*t*=*t*_{k}, it is clear that the first fraction in equation (3.7) approximates the second-order derivative in equation (3.1). On the other hand, for consistency of the scheme (3.7) with the differential equation (3.1), the second fraction in (3.7) must approximate the term *g*[*y*(*t*^{*})^{2}] in (3.1). In view of the mean value theorem, this will be true provided that we choose the potential energy as(3.8)

Instead of the second fraction in equations (3.7 ) and (3.8), an approximation of the formis in general utilized for the term *g*[*y*(*t*^{*})^{2}] in (3.1), where γ(.,.,.) is an appropriate function (cf. Mickens 1994; Anguelov & Lubuma 2001; Anguelov *et al*. 2005). If γ(.,.,.) is chosen such thatthen the associated non-standard scheme (3.7 ) is of order 2. However, we are not interested in such higher order schemes since the non-smoothness of the specific mechanical systems considered in this paper will result in a loss of the order of convergence (see remark 5.1 and the comments introducing §5).

We will often refer to two examples in this paper. The first one is the harmonic oscillator(3.9)With , application to equation (3.9) of the scheme (3.7 ), coupled with (3.8), yields the non-standard scheme(3.10)or, equivalently,(3.11)which is actually the exact scheme of the harmonic oscillator (Mickens 1994). The scheme (3.11) forms part of Gautschi trigonometric method (Gautschi 1961). In the case when the standard denominator (Δ*t*)^{2} is utilized in place of (4/ω^{2})sin^{2}(ωΔ*t*/2), the scheme (3.10) is also known as the Störmer–Verlet method (Hairer *et al*. 2002).

With reference to theorem 2.1, it is not clear in the present general context how the denominator of the discrete derivative in (3.4) should depend on the qualitative parameters related to the differential equation (3.2). Nevertheless, the fact that the discrete derivative in the exact scheme (3.10) of the harmonic oscillator includes explicitly the frequency *ω* of the oscillation motivates our specific choice below of the denominator. (This means that we implicitly assume that the function *g* in (3.1) satisfies *g*(0)=1 in its Maclaurin expansion.) We have the following result.

*The non-standard finite-difference scheme*(3.12)*for* *(3.1)* *is equivalent to the discrete principle of conservation of energy* *(3.6)* *where the discrete potential energy is given by* *(3.8)* *and* . *The scheme* *(3.12)* *reduces to the exact scheme* *(3.10)* *in the case when* *(3.1)* *is the harmonic oscillator* *(3.9)*.

The theorem results from the identity combined with the arguments which, from equations (3.5) and (3.6), led to equations (3.7) and (3.8). ▪

Apart from the discrete solution (position) *y*^{k} at the time *t*_{k}, we will sometimes need an approximation *v*^{k} of the velocity . The discrete velocity is computed from the following system, which is equivalent to the scheme (3.12):(3.13)

Instead of (3.13), it is possible, in the particular case of the harmonic oscillator, to utilize the explicit system(3.14)which is equivalent to the scheme (3.11) (see Mickens 1994).

The second main example for the illustration of our results is the Duffing oscillator(3.15)for which theorem 3.1 yields the non-standard scheme(3.16)

Note the non-local way in which the nonlinear term*g*(

*y*

^{2}) in (3.1) is discretized in the non-standard scheme (3.12) in general, and in (3.16) in particular. On the other hand, for the second order equation (3.1), the initial values

*y*(0)=

*y*

_{0}and are usually given. However, a value for

*y*

^{1}is needed in order to start the scheme (3.12). In analogy with a classical procedure, we utilize the approximation(3.17)

In (3.12), put *k*=0 and replace *y*^{−1} with its expression obtained from (3.17) and this gives the missing starting value whenever the structure of *G* makes (3.12) an explicit scheme. We prefer this procedure, which, in the case of equations (3.10) and (3.16), for example, leads to the values(3.18)and(3.19)respectively. When the scheme (3.12) is not explicit, some alternative procedures are necessary to determine *y*^{1}. A possible procedure is the less accurate approximation(3.20)which readily gives a starting value *y*^{1}.

We stress that the conservation of energy alone does not necessarily make a scheme good. This requirement should be viewed as an addition to the usual properties of convergence, stability, etc. In this regard, figures 1 and 2 produced for the values *ω*=10, Δ*t*=0.001, *y*_{0}=1 and *v*_{0}=0, show how bad the Störmer–Verlet method(3.21)for (3.9) is compared with the exact scheme (3.10). However, both schemes satisfy the principle of conservation of energy; the standard scheme (3.21) is known to be unstable in the sense of Lax–Richtmyer for *ω*Δ*t*≥2.

All the computations in this paper were carried out on a PIV-1.4 GHz personal computer, using Matlab for computing and graphing.

## 4. Vibro-impact problems

We focus now on a mechanical system of one degree of freedom whose motion is described by(4.1)(4.2)(4.3)where and are given initial conditions. The motion is subject to convex constraints in the sense that the position *y*(*t*) is constrained to remain inside a closed convex set, which is the interval *K*=[*z*_{0},∞) or *K*=(−∞,*z*_{0}] of is the position of the obstacle. Following Brogliato (2000, p. 117), the loss of energy during each impact time *t* is modelled by a percussion law of the form(4.4)where *e*∈[0,1] is the restitution coefficient. The constrained problem in (4.1)–(4.3) may be rewritten as a differential inclusion (cf. Goeleven *et al*. 2003):(4.5)

In (4.5), where *y*∈*K*,(4.6)is the subdifferential of the indicator functionof the convex set *K* with *N*_{K}(*y*) being the normal cone to *K* at *y*, that is,

If denotes the normal reaction force that occurs during the collision of the system with the obstacle, then we have the set-valued relationwhich means that =0 whenever (no contact between the system and the obstacle) whereas *y*(*t*)∈∂*K* in the case ≠0 of contact (cf. Goeleven *et al*. 2003, p. 120).

The following explicit existence result is proved by Paoli (1993) and Paoli & Schatzman (1993) in the general setting of multi-degree systems.

*Let K the above interval of* *be such that the initial position in equation* *(4.2)* *satisfies y*_{0}∈*K*. *Let the function* *in* *(4.1)* *be continuous in the three arguments and satisfy a Lipschitz condition in the last two arguments*. *Then*, *there exists a Lipschitz continuous function y* : [0,*T*]→*K with a distributional derivative* *that has a bounded variation*. *This function satisfies equation* *(4.2)* *and is a solution of* (4.1)–(4.4) *in the following sense of variational inequality*:

*p*<+∞, of a sequence extracted from a family of solutions of Yosida-like regularization of the problem (4.5) (cf. Paoli & Schatzman 2002

*a*,

*b*).

## 5. Non-standard impact scheme

There are three main approaches to the numerical solution of the impact problem (4.5). In the event driven scheme's approach, standard schemes (such as the Runge–Kutta method) are utilized to approximate the unconstrained motion, coupled with a procedure for approximating the impact times. The impact law permits the restarting of the algorithm. It can be shown that the order of convergence of the standard schemes is preserved or reduced accordingly as the impacts are isolated or there is an accumulation of impacts (i.e. infinite impacts in a finite time; cf. Janin & Lamarque 2001). This is the most convenient approach. However, in the case of accumulation of impacts, the numerical detection of all impacts is impossible or can be very expensive. The second approach is that of compliant models. In this approach, the impact law is ‘relaxed’ by a normal compliance condition. This leads to a very stiff ordinary differential equation. Standard schemes can then be utilized to approximate the regularization of (4.5) at the cost of a severe restriction on the step size Δ*t* and on the regularization parameter (cf. Paoli 2001; Paoli & Schatzman 2002*a*). Impact schemes constitute the third approach. Their philosophy consists in taking into account the unilateral constraint as well as the restitution coefficient *e*, without any determination or approximation of the impact times. This type of scheme is the point of departure of this paper. They were first proposed and studied by Paoli (1993, 2001) and Paoli & Schatzman (2002*a*,*b*), and are referred to hereafter as PS schemes. More precisely, Paoli & Schatzman proposed a convergent sequence of discrete solutions (*y *^{k}) constructed recursively as follows. The starting values, obtained according to the classical procedure described after example 3.2, are(5.1)where the impact law (4.4) permits to set(5.2)Thereafter, for *k*≥1, the algorithm is(5.3)where *F* ^{k} is some approximation of at time *t*_{k} or *t*_{k+1}. For example, *F* ^{k} can be defined by(5.4)where the function *F*, satisfying the natural consistency relationis continuous in all the arguments and Lipschitz continuous with respect to the last three arguments.

For the implementation of the algorithm, the following lemma is instrumental.

*Let K be a closed convex subset of* , *λ a positive number and* . *Then**where**is the projection of y onto K*.

Now, since the constraint or the inclusion in *K* is also satisfied by the discrete weighted average position , using lemma 5.1, the PS scheme (5.3) implies the ‘workable’ formulation(5.5)

We mention some important facts. The velocity of a vibro-impact system is a discontinuous function owing to the impact law (4.4). This discontinuity makes the (global) order of convergence of the PS scheme to be one for the position and zero for the velocity. However, no theoretical results are available in general about the order of convergence of the PS scheme. The dynamics of vibro-impact systems is often complex. In particular, these systems can be sensitive to initial data and/or have a chaotic behaviour. As a result, the order of convergence of the scheme is not that essential since any prescribed accuracy in the computation is lost in finite time.

As mentioned earlier, one of the key points in the Paoli–Schatzman procedure is that the impact qualitative parameter *e* is ‘suitably’ incorporated in the scheme, a fact that is consistent with the philosophy of the non-standard finite-difference approach presented in §§2 and 3. However, being based on a standard finite-difference method related in some way to (3.21), it is not surprising that its fails to replicate some essential physical properties of the system between two impact times. In an effort to improve it, we design a non-standard version, which notably takes into account the oscillation frequency and the energy property of the system.

For convenience, we analyse the impact problem in (4.1)–(4.5) with equation (4.1), being in the setting of the general conservative oscillator (3.1):(5.6)We then take *K*=[*z*_{0},∞) with *y*=*z*_{0} representing the position of the obstacle. Motivated by theorem 3.1 and the PS algorithm, we set(5.7)and propose the following non-standard impact scheme(5.8)which leads to the form(5.9)

The first starting value *y*^{0} of the non-standard impact scheme is as in equation (5.1). The second initial guess *y*^{1} is obtained in an indirect way, which is an adaptation of what was done for the non-standard scheme (3.12). More precisely, the first situation in (5.2) about *y*_{0} corresponds to absence of collision. Thus, as previously, we utilise the value of *y*^{−1} given by equation (3.17), put it in equation (3.12) or in (5.8) and (5.9), where *k*=0 and solve for *y*^{1}. On the contrary, the condition on *y*_{0} in (5.2)_{2} means that there is collision, i.e. *z*_{0}=*y*_{0}. In this case, the impact law (4.4) legitimates the replacement of equation (3.17) with(5.10)

From (5.10), we proceed as in the first case to obtain *y*^{1}. We recall that when the scheme (3.12) is not explicit; we must utilise an alternative way to determine *y*^{1}, for example, equation (3.20), in the case when there is no impact or its analogue, involving the left-hand side of (5.10), in the case of impact.

Owing to the condition (cf. equation (2.6))the results obtained by the non-standard impact scheme and the PS scheme display a similar convergence property for small values of Δ*t*. However, we emphasize that the non-standard scheme has no restriction on the step size, being unconditionally stable because of theorem 3.1. For the same reason, the new scheme preserves the principle of conservation of energy between two consecutive impact times. Here is a further important qualitative stability property, valid in the case when impact times are isolated. In the limit case *g*≡1, the non-standard scheme in (5.8) and (5.9) is ‘piecewise’ exact. This means that on each interval (*s*_{n},*s*_{n+1}) formed by two consecutive impact times, the non-standard scheme and equation (4.1) have the same general solution (Mickens 1994). Consequently, if denote all the nodes, which belong to [*s*_{n},*s*_{n+1}], then each discrete solution produced by the non-standard scheme is equal to the exact solution at the time (i.e. , *j*≥2) whenever the first two approximations satisfy the requirement and . (See example 5.1 below.) On the contrary, the discrete solutions given by the PS scheme (i.e. equations (5.8) and (5.9)), with (Δ*t*)^{2} as denominator of the discrete derivative, are not good approximations whenever the quantity *ω*Δ*t* is large. (See also remark 5.2 below.)

Consider the impact system (*g*≡1 in (4.1) and (5.6))modelling the position *y*(*t*) of a particle that acts as an harmonic oscillator with instantaneous impact against a fixed obstacle. (This oscillator was considered without impact in example 3.1.) It is clear that the explicit trajectory before the first impact time *s*_{0} isThus, the first impact occurs at the time *s*_{0}=*π*/2*ω*. Moreover, since , the impact law shows that at the second impact time *s*_{1}, we haveThe second impact occurs then at the time *s*_{1}=(*π*/2*ω*)+(*π*/*ω*). Proceeding by induction, we obtain the solutionwhere, for *n*=0,1,2,…, is the time when the impact number *n*+1 occurs.

The total energy of the system is . Thus, the system is conservative, with , for all times except at the impact times if, and only if, *e*=1. If 0≤*e*<1, then the system is conservative, with , on each interval of time (*s*_{n},*s*_{n+1}) for *n*≥0.

Application to this system of the non-standard impact scheme (5.8),(5.11)for *ω*=10, Δ*t*=0.001, *y*_{0}=1, *z*_{0}=0 and *v*_{0}=0 yields figure 3 for displacements and figure 4 for total energy corresponding to different values of *e*. The non-standard algorithm needs only 0.373 cpu time, whereas the standard algorithm needs 0.435 cpu time.

Consider again example 5.1. While the facts stated just before this example apply, we add another weakness of the standard scheme as well as a further important qualitative stability property of the non-standard impact scheme. Let *ω*≥1 and consider the standard scheme with Δ*t*=*π*/(4*ω*) so that the nodes *t*_{0}, *t*_{1} and *t*_{2} are in the interval [0,*s*_{0}] before the first impact. Then, we obtain *y*^{0}=1, and . Using (3.10), (3.11) and (3.18), similar computations with the non-standard scheme give *y*^{0}=1, and . This shows that the standard scheme gives a bad indication on the impact time.

Assume now that there is no constraint (impact) in example 5.1. The solution is then *y*(*t*)=cos *ωt* and it satisfies the boundedness property −1≤*y*(*t*)≤1 for all *t*≥0. The non-standard scheme is stable with respect to this property, being actually the exact scheme. However, the standard scheme does not replicate this property. The situation is even worse if we consider a mechanical system with two stops located at 1+*ϵ* and −(1+*ϵ*) for a very small *ϵ*>0: computation on a long period of time does not yield the relation −1≤*y*^{k}≤1 for the standard discrete solution and, thus, unlike the non-standard scheme, the standard one may generate spurious impacts.

Consider again the Duffing oscillator (3.15) (cf. example 3.2). We apply the non-standard scheme for the parameters in the table below:

Figure 5 shows the displacements for different values of *e*. This figure shows that the system oscillates between the fixed point *y*=0, even if the restitution coefficient is not equal to one: the system becomes periodic as displayed by the phase portraits in figure 6. Moreover, when 0≤*e*≤1, the total energy of the system converges to a constant. In figures 7 and 8, the total energy computed by the standard method shows some ‘bad’ oscillations while the total energy obtained by the non-standard method is constant. The computation requires 9.93 cpu time and 9.84 cpu time for the non-standard and standard algorithms, respectively.

We revisit example 5.2 on changing only the position of the obstacle as indicated in the following table:

The phase portraits in figures 9 and 10 show a possible chaotic behaviour of the Duffing oscillator. There is a clear distinction between the phase portraits computed by the standard and the non-standard schemes. The chaotic behaviour is also confirmed by figures 11 and 12, which show that the total energy computed by both methods (standard and non-standard) is not constant: the standard scheme is dissipative whereas the non-standard scheme is not. The non-standard scheme being more reliable, the authors believe that the chaotic behaviour is owing to the mechanical system itself and not to the numerical methods. This statement is consistent with a well known fact that vibro-impact systems can develop complex and chaotic behaviours (see, for example, Budd *et al*. 1995).

## 6. Concluding remarks

One of the main concerns in the study of impact oscillators is that these are non-smooth systems, which have complex and even chaotic behaviours. Typically, at the time of the impact, which is supposed to be instantaneous, the velocity of the oscillator changes instantaneously. As a result, this may induce uncomfortable phenomena such as noise, fatigue and so on. Theoretical analysis and simulations performed in laboratories on simple models have contributed to gain some understanding. For more complex systems, not only is it necessary to utilize numerical methods, but it is essential to construct schemes, which preserve, as much as possible, the significant physical properties of the systems.

A major contribution of this paper is the design of numerical schemes in which the intrinsic qualitative features of vibro-impact systems are incorporated in order for the schemes to replicate the dynamics of these non-smooth mechanical systems. The paper is the first application of Mickens' non-standard finite-difference method to vibro-impact problems (which mathematically translate into differential inclusions) and enriches, therefore, the list of real-life problems which have recently been handled by the non-standard approach (see Mickens 2000; Gumel 2003; Gumel *et al*. 2003; Moghadas *et al*. 2003; Anguelov *et al*. 2005 and references therein).

To be more precise, we have constructed a non-standard finite-difference method for impact oscillators with one degree of freedom. Motivated by the structure of the exact scheme of the harmonic oscillator with one frequency, the construction utilizes two of Mickens' rules. These are the renormalization of the denominator of the discrete derivative and the non-local approximation of the nonlinear terms involved in the equation (cf. Mickens 1994). Compared with Paoli & Schatzman impact schemes (Paoli 1993, 2001; Paoli & Schatzman 2002*a*,*b*), which are based on standard finite-difference discretizations, the schemes constructed in this work are unconditionally stable and they replicate a number of key physical properties (e.g. conservation of energy, etc.) of the oscillators between two consecutive impact times. The numerous numerical simulations considered in the paper show that, for a given mesh step size, the times of computation are of the same order for both the non-standard and the standard schemes.

Our interest for future research is to extend this study to multi-degree vibro-impact problems (e.g. Peterka 2000), especially those which involve partial differential equations. As a first task in this direction, we think of the vibrations of a beam between stops (Dumont 2002, 2003). Another topic of interest is the investigation of suitable non-standard finite-difference schemes when further impact laws by Brogliato (2000) are utilized.

## Acknowledgments

This work was carried out while the second author was a visiting professor at the Université de la Réunion. The hospitality of the Department of Mathematics and of the laboratory IREMIA are greatly appreciated. The authors would like to thank Prof. R. E. Mickens and the referees for their remarks and constructive criticisms.

## Footnotes

- Received June 24, 2004.
- Accepted October 6, 2004.

- © 2005 The Royal Society