INTRODUCTION
Impact dynamics is considered to be one of the most important problems which arise in vibrating systems. Such impacting oscillators may occur in the motion with amplitude constraining stops (Fig. 1). This problem has been researched in detail using bifurcation theory^{[17]}. Different types of impacting response due to different ranges of driving frequency or control parameters can be predicted from bifurcation diagrams. In this study, calculation of the Lyapunov exponents for impact oscillator is considered. Lyapunov exponents measure the average divergence or convergence rate of nearby trajectories on a particular return map in space. If the largest Lyapunov exponent is positive the trajectory is chaos, whereas nonpositive exponents indicate the trajectory is a stable motion. Thus the spectrum of Lyapunov exponents is one of the most useful diagnostics for systems.
Lyapunov exponents described the system’s behavior and provide the stability or instability of an equilibrium point of the nonlinear system of the differential equation. For mdimension system, with the mapping P: R^{m}, → R^{m},
a small perturbation of this mapping t_{n} → t_{n}+δ t_{n}, the Taylor expansion of Eq. 1 and the linearization model is
The stability at point t_{n} is decided by the eigenvalues of the Jacobian matrix DP^{(n)}(t_{n}). The Lyapunov exponents measure the average rate of the convergence or divergence of the system and are defined by many researchers^{[810]}.
Thus m exponents of the average rate λ_{m} can be calculated if the Jacobian matrix is available.
In general, if the underlying dynamic system is smooth, then the spectrum of Lyapunov exponents can be calculated using the Poincare return map^{[10]}. For example, a three dimensions flow in (t, x, v) space, the solutions can be projected onto a particular section t = ψ, where ψ is a constant in [0, T = 2π/ω]. Thus a twodimensional Poincare map is defined by a difference equation,
where, x, v are vectors and f, g are nonlinear transformations. The Jacobian matrix is expressed as:

Fig. 1: 
Impact oscillator 
In the current problem of the impact oscillator under consideration is nonsmooth due to the amplitude constraint. Thus the equation of the Poincare map is nonsmooth when impact occurs and the derivative in Eq. 5 is not defined^{[1114]}. To overcome this problem the impact map is considered to develop the algorithms.
MATHEMATICAL MODEL
The impact oscillator is considered in Fig. 1, where the system is under a harmonic excitation with an amplitude constraint. The impact oscillator is governed by the following equation
where, ,
x are the acceleration, the velocity, the displacement and m, c, k, p, ω,
g are the mass, the damping, the stiffness, the forcing amplitude, the forcing
frequency, the amplitude constraint, respectively. The impact occurs whenever
x = g and the velocity is modeled by ,
where the t^{+} represents the time after impact, the t¯ represents
the time before impact and r is restitution coefficient. For simplicity, we
set m = 1, c = 0, k = 1 and p = 1 then the system convert to
With initial conditions x(t_{0}) = x_{0} and ,
the solution of the displacement x and the velocity v between impacts are:
where, z = 1/(1ω^{2}).
JACOBIAN MATRIX FOR IMPACT MAP
To obtain the difference equation of Eq. 1, the variables x_{n}, v_{n} and t_{n} are considered as the impact n, then the variables x_{n+1}, v_{n+1} and t_{n+1} are obtained from Eq. 8 and 9 for the initial conditions
Thus the difference equations can introduce as below:
where, x_{n+1} = x_{n} = g. Thus the impact return map can
be constructed from a three dimensions flow in (t, x, v) space onto a two dimension
map (t,v) when x = g. Using Eq. 11 and 12,
the Jacobian matrix for impact map can be derived^{[1,2]},
The components in Eq. 13 are:
where,
are the accelerations at the two impacts.
LYAPUNOV EXPONENTS
Using the Jacobian matrix of Eq. 13 and the definition of Eq. 3, two Lyapunov exponents of impact oscillator can be obtained by:
In computation of the eigenvalues of the matrix J, that is given by:
Only a few iterations, the matrix J becomes very large for chaotic case and null for the periodic case because the product of the matrix DI in Eq. 20. To overcome this problem, the QRfactorisation technique applied^{[1517]}. Given a matrix A, there is a unique factorization
where, R is a square right upper triangular matrix and Q has orthogonal columns. Consider an arbitrary orthogonal matrix Q_{0} and apply the orthogonalisation procedure to the impact map
The columns Q_{1} are now an orthogonal basis for the impact map at n = 1 and R_{0} gives the relationship between this map and the image of the Q_{0} under the Jacobian matrix DI^{(1)}. This process can be repeated recursively and gives:
Thus the calculation of the Lyapunov exponents in Eq. 19 is converted to the summation of the matrix In(R_{j}) in Eq. 24.

Fig. 2: 
Bifurcation diagram of the impact oscillator under the parameters
m=1, c=0, k=1, p=1, r=0.8 and g=0 

Fig. 3: 
Estimated Lyapunov exponents of the impact oscillator under
the parameters m=1, c=0, k=1, p=1, r=0.8 and g=0 
The result shown in Fig. 2 is the plot of bifurcation diagram
for varying the driving frequency; the regions of the stable and unstable oscillator
are demonstrated by points set in impact map. The corresponding results of the
estimated Lyapunov exponents are shown in Fig. 3, where the
positive values indicate the stable motions and the nonpositive values indicate
unstable motions. For example, the results in Fig. 4a show
periodic points set in impact map. The corresponding results of the estimated
Lyapunov exponents are shown in Fig. 4b, where both exponents
converge to nonpositive values after 5000 iterations. This indicates that the
impact oscillator is stable at driving frequency ω = 2.

Fig. 4: 
Impact oscillator under the parameters m=1, c=0, k=1, p=1,
r=0.8, g=0 and ω=2 (a) impact map (b) Lyapunov exponents 

Fig. 5: 
Impact oscillator under the parameters m=1, c=0, k=1, p=1,
r=0.8, g=0 and ω=2.8 (a) impact map (b) Lyapunov exponents 
For the case of driving frequency at ω = 2.8, the results show strange
points set in Fig. 5a and one of the estimated Lyapunov exponents
converges to a positive value in Fig. 5b. This confirms that
the impact oscillator is unstable.
CONCLUSIONS
The algorithm of computing Lyapunov exponents is successfully applied to the impact oscillator. The impact oscillator is considered as nonsmooth dynamic system; therefore the Jacobian matrix for impact map is derived. The results show the plot of bifurcation diagram consists with the estimated Lyapunov exponents for varying the driving frequency.