INTRODUCTION
Impact dynamics is considered to be one of the most important problems in mechanical vibrating systems. Such impacting oscillation may occur in the motion with amplitude constraining stops (Fig. 1). This problem has been researched in detail using mathematical model (Shaw and Holmes, 1983; Budd et al., 1995; Luo and Xie, 2001) and the simulation results showed in Fig. 2 (Lee, 2005). Different types of impacting response can be predicted from bifurcation diagram for varying the control parameter. For the purpose of condition monitoring of the system, Lee (2006) calculated the spectrum of Lyapunov exponents for varying the driving frequency. The results showed an agreement with the bifurcation diagram of Fig. 2. In this study, the estimation of Lyapunov exponents using the impacting time series is considered, where the delay embedding techniques proposed by Darbyshire and Broomhead (1996) and Babrook et al. (1997) are applied. In the following study, a complete procedure of computing Lyapunov exponents based on time series is developed and the results are compared to the previous theoretical results.
IMPACT OSCILLATION
An impact oscillator is considered in Fig. 1, where the system is under a harmonic excitation with an amplitude constraint. The equation of motion, when not impacting, is
where
are the acceleration, the velocity, the displacement, the mass, the damping,
the stiffness, the forcing amplitude, the forcing frequency and the amplitude
constraint, respectively.

Fig. 1: 
Impact oscillator 

Fig. 2: 
Bifurcation diagram of the impact oscillator for the control
parameters m = 1, c= 0.01, k = 1, p = 1, r = 0.8, g = 0 and ω = 0 to
5 
The impact occurs whenever x = g and the velocity is modeled as v(t^{+})
= rv(t^{}) where r is restitution coefficient, (t^{+}) is
the time after impact and (t^{}) is the time before impact.

Fig. 3: 
The impact oscillator for the control parameters m = 1, c
= 0.01, k = 1, p = 1, r = 0.8, g = 0 and ω = 2 (a) Displacement (b)
Phase trajectory (c) Impact map (d) Estimated Lyapunov exponents using the
observed time series of impact velocity, where the nearby vectors N_{b}=5
and the embedding dimension m = 2 
The result shows in Fig. 2 is a bifurcation diagram with
driving frequency range from 0 to 5. The regions of the stable and unstable
impact oscillations are clearly demonstrated. For the case of the driving frequency
at ω = 2, the displacement of the impact oscillation is shown in Fig.
3a. This is an example of a periodic impact oscillation and the phase trajectory
demonstrates a close orbit in Fig. 3b. For the same control
parameters but at driving frequency ω=2.8, the result shows in Fig.
4a and b demonstrate a nonperiodic impact oscillation.
LYAPUNOV EXPONENTS ESTIMATION
The Lyapunov exponents are measures of the average rate of the convergence or divergence of the system and defined as
This is based on mdimension system with the mapping I:R^{m} → R^{m},
For a small perturbation of this mapping t_{n} ⇒ t_{n}+δt_{n}, the Taylor expansion of Eq. 3 and the linearization model is
The stability at point t_{n} is decided by the eigenvalues of the Jacobian matrix DI^{(n)}(t_{n}). Thus m exponents of the average rate λ^{m} can be calculated if the Jacobian matrix is available.
In the following investigation, the estimation of the Jacobian matrix is considered using the simulated impacting time series. Based on the definition of the Lyapunov exponents in Eq. 2, the exponents for impact oscillation are
where DI^{(n)}(t_{n}) is the Jacobian matrix for impact map

Fig. 4: 
The impact oscillator for the control parameters m = 1, c
= 0.01, k = 1, p = 1, r = 0.8, g = 0 and ω = 2.8 (a) Displacement (b)
Phase trajectory (c) Impact map (d) Estimated Lyapunov exponents using the
observed time series of impact velocity, where the nearby vectors N_{b
}= 5 and the embedding dimension m = 2 
Now consider a practical situation that the Jacobian matrix of the impact system is unknown and only a time series of impacting signals is known. Thus, the aim is to estimate Jacobian matrix using observed time series.
For mdimensional system, an equivalent dynamical system can be constructed using delay embedding from time series (Darbyshire and Broomhead, 1996; Babrook et al., 1997)
and the embedding representation of Eq. 7 using the delay coordinate embedding techniques is
or in matrix form
Thus the estimated Jacobian matrix DI^{n}(t_{n}) of Eq. 6 is accomplished by estimation the difference vector δt_{n} and δt_{n+1} in Eq. 4. For example, DI(t_{1}) is calculated by examining how a local area evolves from one vector t = [t_{1}t_{2}t_{3}…t_{m}] to another t = [t_{2}t_{3}t_{4}…t_{m+1}] in Eq. 9, which is achieved by finding several nearby vectors N_{b} and forming a neighborhood matrix B_{1} (N_{b}xm) and its evolved matrix B_{2} (N_{b}xm) . Thus the Eq. 4 is expressed by
and the pseudoinverse solution of the estimated Jacobian matrix is given by
where ( )^{T} is transposition matrix and ( )^{1} is inverse matrix.
In computation of the eigenvalues of the matrix J in Eq. 5, 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.
12. To overcome this problem, the QR factorisation techniques in proving
the multiplicative ergodic theory (Eckmann and Rulle, 1985; Johnson et al.,
1987) are applied. 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 gives
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
This method removes the need to compute products of large numbers of matrices.
Thus the Lyapunov exponents have the form of a simple ergodic average and can
be expressed by
Summary of the processes of estimating Lyapunov exponents from time series:
1. 
Embed the data using Eq. 9. 
2. 
Locate nearest embedded vector and form a neighbourhood
matrix B_{2}. 
3. 
Evolve one step and locate nearest embedded vector and form a
neighbourhood matrix B_{2}. 
4. 
Calculate the estimated Jacobian matrix using the pseudoinverse method
in Eq. 11. 
5. 
Perform QR decomposition. 
6. 
Calculate the average exponents using Eq. 16. 
7. 
Repeat the processes (26) until the exponents converge to fixed values. 
RESULTS
The result shows in Fig. 3c is the points set in impact map at driving frequency ω = 2 and the estimated Lyapunov exponents are shown in Fig. 3d. The Lyapunov exponents converge to a nonpositive value after 5000 iterations. This indicates that the impact oscillation is stable and is in consistent with Fig. 2. For the case of driving frequency at ω=2.8, the results show in Fig. 4c. The “strange” point sets in impact map demonstrate the impact oscillator is chaos. One positive Lyapunov exponents shown in Fig. 4d confirms that the impact oscillation is unstable. Both results in Fig. 3d and 4d are also in consistent with Fig. 2.
CONCLUSIONS
The algorithms of computing Lyapunov exponents were investigated using observed impacting series. The calculating processes were (1) estimated Jacobian matrix by embedding techniques and (2) averaged the exponents by applying QR decomposition. Finally, the results of the estimated Lyapunov exponents showed an agreement with the bifurcation diagram.