Research Article
Lyapunov Exponents Estimation of the Impact Oscillation Using Observed Time Series
Department of Electrical Engineering, Far East College, 49 Chung-Hua Road, Hsin-Shih Town, Tainan 74404, Taiwan
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
(1) |
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 Nb=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 non-periodic 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
(2) |
This is based on m-dimension system with the mapping I:Rm → Rm,
(3) |
For a small perturbation of this mapping tn ⇒ tn+δtn, the Taylor expansion of Eq. 3 and the linearization model is
(4) |
The stability at point tn is decided by the eigenvalues of the Jacobian matrix DI(n)(tn). 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
(5) |
where DI(n)(tn) is the Jacobian matrix for impact map
(6) |
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 Nb = 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 m-dimensional system, an equivalent dynamical system can be constructed using delay embedding from time series (Darbyshire and Broomhead, 1996; Babrook et al., 1997)
(7) |
and the embedding representation of Eq. 7 using the delay coordinate embedding techniques is
(8) |
or in matrix form
(9) |
Thus the estimated Jacobian matrix DIn(tn) of Eq. 6 is accomplished by estimation the difference vector δtn and δtn+1 in Eq. 4. For example, DI(t1) is calculated by examining how a local area evolves from one vector t = [t1t2t3 tm] to another t = [t2t3t4 tm+1] in Eq. 9, which is achieved by finding several nearby vectors Nb and forming a neighborhood matrix B1 (Nbxm) and its evolved matrix B2 (Nbxm) . Thus the Eq. 4 is expressed by
(10) |
and the pseudo-inverse solution of the estimated Jacobian matrix is given by
(11) |
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
(12) |
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
(13) |
where R is a square right upper triangular matrix and Q has orthogonal columns. Consider an arbitrary orthogonal matrix Q0 and apply the orthogonalisation procedure to the impact map gives
(14) |
The columns Q1 are now an orthogonal basis for the impact map at n = 1 and R0 gives the relationship between this map and the image of the Q0 under the Jacobian matrix DI(1). This process can be repeated recursively and gives
(15) |
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
(16) |
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 B2. |
3. | Evolve one step and locate nearest embedded vector and form a neighbourhood matrix B2. |
4. | Calculate the estimated Jacobian matrix using the pseudo-inverse method in Eq. 11. |
5. | Perform QR decomposition. |
6. | Calculate the average exponents using Eq. 16. |
7. | Repeat the processes (2-6) until the exponents converge to fixed values. |
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 non-positive 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.
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.