Spectrum of Lyapunov exponents in non-smooth systems evaluated using orthogonal perturbation vectors

This paper covers application of the novel method of Lyapunov exponents (LEs) spectrum estimation in non smooth mechanical systems. In the presented method, LEs are obtained from a Poincaré map. By analysing the map instead of the full trajectory, problems with transition of perturbations through discontinuities can be avoided. However, the explicit formula of the map is usually not known. Therefore, the Jacobi matrix of the map is estimated using small perturbations of the initial point. In such a manner, direct calculation of the Jacobi matrix can be avoided. The article provides a detailed description of the method accompanied by clear schemes. The algorithm of Jacobi matrix estimation is elaborated and an example is given. Efficiency of the method is confirmed by a numerical experiment. The mechanical oscillator with impact has been simulated. Bifurcation diagrams and Lyapunov exponents graphs have been generated. It has been shown that the method provides values of the whole Lyapunov exponents spectrum with high accuracy.


Introduction
Lyapunov exponents (LEs) are one of the most useful criteria for determining the stability of solutions of dynamical systems in local [1] and also global sense [2]. They are extremely useful instruments for identification of the motion character in these systems. Determination of these numbers is one of the fundamental tasks in nonlinear dynamics. Definition of LEs was established by V. I. Oseledec [3], in a form suited for the theory of dynamical systems. These numbers were named after A. M. Lyapunov, due to the fact that their values are qualitative and quantitative illustration of his criterion of dynamical systems stability [1]. LEs can be treated as the exponential measure of divergence or convergence of orbits in the phase space that start close to each other. Mathematically, LEs are eigenvalues of the state transition matrix that determine evolution of an infinitesimal perturbation in the dynamical system. Hence, they can be treated as a measure of the sensitivity to initial conditions in the phase space.
For preliminary assessment of the system dynamics, it is enough to know the largest Lyapunov exponent. Positive value of the largest LE indicates high sensitivity to initial conditions of the system and is a sign of an irregular response, i.e. chaotic or even hyperchaotic. The largest LE equal to zero means regular, periodic or quasi-periodic system dynamics. If all the LEs are negative, then the limit set is a stable fixed point. However, for more detailed and accurate analysis of the system, the knowledge of the complete spectrum of LEs is essential. For any dynamical system, such spectrum is a set of real numbers sorted in the non-increasing order. Number of LEs is equal to the dimension of the phase space. Knowing all the LE values in the spectrum enables to estimate the fractal dimension of the system attractor (according to the Kaplan-Yorke formula [4]), determine the dimension of the system quasiperiodicity (defined by the number of LEs of zero value) or the scale of hyperchaos (defined by number of positive LEs). On the other hand, Lyapunov exponents also provide a criterion for global stability of the system [2]. From properties of the LEs one can conclude that the sum of all of them is equal to the divergence of the phase flow [5]. If the sum of all components in the spectrum is negative, then the dynamical system under consideration is dissipative and its attractor is globally stable. This sum is equal to zero for any conservative system. In the case of positive sum of all the LEs, an exponential divergence of the phase flow is observed, so the system is globally unstable.
The analytical determination of Lyapunov exponents is possible only in the considerations of simple linear dynamical systems. Owing to this fact, numerical methods are used to calculate LEs in most cases. The first numerical characteristics of chaotic behaviour of the dynamical system, which presented a divergence of neighbouring trajectories, was conducted by Henon and Heiles [6]. An efficient algorithm for calculating the complete spectrum of the Lyapunov exponents based on the Oseledec theorem was independently formulated by both Benettin et al. [7] and Shimada, Nagashima Other technique, which can be also applied in the case of systems with discontinuities, is an algorithm of estimation of the LEs from the scalar time signal. This method is based on a procedure of reconstruction of the system attractor from the time dependence of one of the coordinates that was introduced by Takens [34]. The first numerical algorithm for estimation of the largest LE, derived from it, was formulated by Wolf [35]. In the following years, this algorithm has been improved by other researchers in order to estimate the entire spectrum of the LEs and to recognize the so-called spurious Lyapunov exponents [36].
In this paper we demonstrate a method for calculating the spectrum of LEs which is based on the idea of estimation of the Jacobian matrix using small perturbations of the initial orthogonal vectors. This technique is especially useful for non-smooth oscillators, because such estimation procedure allows us to avoid direct calculation of the Jacobi matrix. We have successfully applied this procedure in order to test the

The method
The method under consideration is based on the classical algorithm by Parker & Chua [12]. However, calculation of linearized equations is replaced with the estimation of Jacobi matrices using small orthogonal perturbation vectors. Scheme of the Jacobi matrix estimation method is presented in the Fig. 1.   Fig. 1. Scheme of the Jacobi matrix estimation method  [fi(x1,... ,xN-1

where DU(xn) is a Jacobi matrix composed of (N-1) column vectors DUi(xn) =
Consequently, each column vector can be estimated as: Finally, in accordance with the Eq.(4), estimated Jacobi matrix of the Henon map is obtained in the following form: The example above implies that the precision of the Jacobi matrix estimation depends on the magnitude of the small parameter  Note that, in general, this procedure requires high accuracy of the trajectory simulation, especially in the vicinity of the discontinuity. Then, such numerically generated Jacobi matrix can be implemented in a classical algorithm for calculating Lyapunov exponents spectrum, including the Gramm-Schmidt orthonormalization [12].

Numerical example
In this paper, the method presented in the previous section has been applied to a simple mechanical system with impacts. Scheme of the system is presented in the Fig. 2. The system is a typical linear mass-spring-damper system. However, motion of the mass is constrained by the wall placed in the position Xw. The equation of motion of the system is as follows: where X is the coordinate of the oscillator, m is the mass of the oscillator, k is the spring stiffness, c is the damping coefficient, F is the amplitude of forcing, Ω is the angular frequency of forcing and Xw is the position of the wall. Obviously, this equation is not correct when a contact with the wall occurs. If tc is the moment of time in which collision with the wall takes place, then transition of the oscillator's velocity is defined by the following equation: (12) where e is the coefficient of restitution (CoR).

Fig. 2. Scheme of the analysed system with impacts
By simple rearrangements, Eq. (11-12) can be transformed to the form: where  is the damping ratio, ω 2 =k/m is the squared natural frequency

Observations and conclusions
It can be noticed that the bifurcation diagram presented in the Fig. 3 is consistent with the Lyapunov exponents graph. Moreover, it can be shown that the divergence of the system (13-14) is equal to  2 =-0.1. As one can notice in the graph, the sum of both Lyapunov exponents from the spectrum is also equal to -0.1 independently of the value of the control parameter η, which confirms validity of calculations.