Fast and simple Lyapunov Exponents estimation in discontinuous systems

Typically, to estimate the whole spectrum of n Lyapunov Exponents (LEs), it is necessary to integrate n perturbations and to orthogonalize them. Recently it has been shown that complexity of calculations can be reduced for smooth systems: integration of (n-1) perturbations is sufficient. In this paper authors demonstrate how this simplified approach can be adopted to non-smooth or discontinuous systems. Apart from the reduced complexity, the assets of the presented approach are simplicity and ease of implementation. The paper starts with a short review of properties of LEs and methods of their estimation for smooth and non-smooth systems. Then, the algorithm of reduced complexity for smooth systems is shortly introduced. Its adaptation to non-smooth systems is described in details. Application of the method is presented for an impact oscillator. Implementation of the novel algorithm is comprehensively explained. Results of simulations are presented and validated. It is expected that the presented method can simplify investigations of non-smooth dynamical systems and support research in this field.


Introduction
Lyapunov Exponents (LEs) are the measures of the sensitivity of a dynamical system to a perturbation of its initial conditions: they indicate average, exponential rates of expansion or contraction of an infinitesimal disturbance of the initial state [1]. As such, they are widely used in analysis of systems' dynamics. The negative sum of all the LEs in the spectrum is a necessary and sufficient condition for existence of an attractor. Positive value of the largest Lyapunov Exponent (LLE) indicates that the attractor is chaotic, whereas two positive LEs are the sign of hyper-chaos. If the whole spectrum of LEs is negative, then the limit set is a stable equilibrium point. When one LE equals 0 and the rest of them are negative, a stable limit cycle is present. Finally, existence of K LEs equal 0, with the remaining ones negative, confirms that the attractor is a K -torus [1]. Last but not least, values of the LEs are connected with the fractal dimension of the attractor [2]. Note that the Lyapunov Exponents being the topic of this paper are often called global [3], as they predict behavior of the system as the time goes to infinity. Their values remain the same for any Nonlinear and Complex Physics 2169 2 The method 2.1 Evolution of an initial perturbation of a dynamical system Consider a dynamical system described by an ODE [1]: where x ∈ R n is a state vector,ẋ = dx dt , n is the order of the system, t is the time, f (x, t):R n ×R → R n is a vector field, x 0 is a vector of initial conditions and t 0 is an initial time. The solution x (t) :R → R n , t ≥ t 0 , x (t 0 ) = x 0 is called a trajectory of the system (1) starting from the initial conditions x 0 and the initial time t 0 . If the system under consideration is discontinuous, it may be necessary to augment the equation (1) with appropriate transition rules, which define behavior of the trajectory as it crosses the discontinuity.
Assume that the initial conditions x 0 of the system (1) are perturbed by an infinitesimal, initial perturbation vector z 1 0 (the lower index "0" means that the perturbation is an initial value, i.e. z 1 0 = z 1 (t 0 ), and the upper index "1" is an ordinal number, as it is going to be explained later). Then, the system (1) generates a disturbed trajectory y(t) such thatẏ = f (y, t) , y (t 0 ) = x 0 + z 1 0 . The perturbation vector z 1 (t) = y (t) − x(t) is the difference between the perturbed solution and the undisturbed one. Evolution of z 1 is described by the following equation: where J ij = ∂fi ∂xj is the Jacobian matrix of the vector field f with respect to the state vector x. The linearization f x + z 1 , t − f (x, t) = J (x, t) z 1 is valid as long as the perturbation vector z 1 is infinitesimal and the vector field f is differentiable at the point x with respect to all the state variables x i [33].

Estimation of the largest Lyapunov Exponent (LLE)
The precise definition of Lyapunov Exponents (LEs) is based on the solution of the variational equation of the system (1). Its detailed explanation is beyond the scope of this paper. Readers interested in this topic can refer, for example, to the book [1]. Intuitively, the largest Lyapunov Exponent (LLE) can be regarded as an exponential rate of increase or decrease of the length of almost any initial, infinitesimal perturbation z 1 0 , according to the following approximate formula: where λ 1 is the LLE and e is the Euler's number. The formula (3) becomes exact as t → ∞. Solving the equation (3) with respect to λ 1 yields: The formula (4) provides a method for estimation of the LLE. One can solve the system (1) and the perturbation equation (2) for a long enough time t and substitute results to the formula (4). However, it has to be taken into account that z 1 (t) may

2170
The European Physical Journal Special Topics approach 0 or ∞ for large values of the time t. In such case, to overcome numerical problems, it is necessary to periodically normalize the perturbation vector z 1 . Consequently, the formula (4) has to be slightly modified before its application in a numerical procedure. For more details on implementation of this LEs estimation method, one can refer to the book [1]. An alternative approach to calculation of the LLE [11,12] results from the following transformation of the formula (4): where the operator " · " denotes the scalar (dot) product of vectors. The equation (5) can be proven by differentiation of ln |z 1 (t) | with respect to the time t. Please note that the last expression in (5) is just the mean value of z 1 (τ ) · dz dτ / z 1 (τ ) 2 over the time interval t 0 ≤ τ ≤ t. Therefore, to estimate the LLE, one can simply evaluate the following formula in each time step: (6) and the average of all the values λ * 1 approaches the LLE: where λ 1 is the LLE, ∆t ∈ R + is a small time step and k ∈ N + is the number of iterations of the procedure, so that t − t 0 = k * ∆t. The appropriate value of k is not known a priori. The calculation can be stopped when the obtained value of λ 1 stabilizes. This method has the following assets: -it requires only the simplest mathematical operations, such as addition, multiplication and division, the logarithmic function is no longer necessary [11]; -application of the formula (6) results in faster calculations of the LEs [12]; -the formula (6) is independent of the length of the perturbation vector z 1 , and therefore its value does not change when z 1 is normalized.
The latter can be proven by substituting dz 1 (t) dt = J (x, t) z 1 (t) from the equation (2) to the formula (6).

Estimation of the LEs spectrum
Any dynamical system of the order n is characterized by n LEs: λ 1 ≥ λ 2 ≥ · · · ≥ λ n . To calculate the whole spectrum of n LEs, it has to be taken into account that each LE is connected with a direction in the phase space [1]. It means that almost any initial perturbation z 1 0 has a component along the direction of each LE: λ 1 , λ 2 , . . . , λ n . The formula (3) is valid due to the fact that, as the perturbation z 1 evolves, its component along the direction of λ 1 becomes relatively larger than other components. Consequently, the perturbation z 1 (t) aligns with the direction of λ 1 as the time goes by. This fact can be used to calculate other LEs. Consider another random, initial perturbation vector z 2 0 . Its evolution z 2 (t) can be predicted by means of the equation (2). In order to calculate the second LE, i.e. λ 2 , the vector z 2 must have no component in the direction of λ 1 . Under such condition, it is going to align with the direction of λ 2 and the rate of change of its length is going to be connected with the value of λ 2 . To achieve this, the perturbation z 2 must be continually orthogonalized with respect to the first perturbation z 1 , which aligns with the direction of λ 1 . Such orthogonalization can be performed using the following formula: where z 2 is the second perturbation after orthogonalization with respect to the first one, i.e. z 1 . Analogously, to estimate the i th LE λ i , it is necessary to observe a perturbation z i (t) which has no components along the directions of λ 1 , λ 2 , . . . , λ i−1 .
Then, it aligns with the direction of λ i and changes its length according to the value of λ i . For this purpose, the perturbation z i must be orthogonalized with respect to z 1 , z 2 , . . . , z i−1 . Consequently, z i becomes orthogonal to the directions of λ 1 , λ 2 , . . . , λ i−1 . This process is performed by means of the following orthogonalization formula [1]: where z i is the i th perturbation after orthogonalization with respect to z 1 , z 2 , . . . , z i−1 . In a popular approach to estimation of the LEs spectrum [1], the orthogonalization procedure is performed together with normalization. Such process is referred to as "Gram-Schmidt Orthonormalization". However, as it was mentioned before, the formula (6) for LEs estimation is independent of the length of the perturbation, so the normalization process does not have to be synchronized with estimation of the LEs. After orthogonalizing the perturbations z 2 , . . . , z n , the whole spectrum of LEs can be computed by means of the formula (10), analogous to the equation (6).
where z 1 = z 1 , as the first perturbation does not require orthogonalization. Similarly as in the case of λ 1 , the final value of the i th LE λ i is obtained from an average of all the values λ * i (11): At this point, it is worth to mention an important feature of the presented method: in order to estimate the whole spectrum of n LEs λ 1 , λ 2 , . . . , λ n , one needs to integrate only n−1 perturbations [12]. This is a result of the following reasoning. The first perturbation z 1 is a vector in n-dimensional state space. The second one, z 2 , is kept orthogonal to z 1 , so it has no component along z 1 . Therefore, the subspace in which the perturbation z 2 exists has the dimension n − 1. In general, the i th perturbation z i is contained in a subspace of the dimension n − i + 1. Consequently, the subspace of the last perturbation z n is one-dimensional, which corresponds to a straight line. As a result, any random vector of the nth initial perturbation z n 0 ,

2172
The European Physical Journal Special Topics after orthogonalization with respect to z 1 , z 2 , . . . , z n−1 , falls into a 1-dimensional subspace, which is connected with the direction of the last LE, i.e. λ n . Now, taking into account the fact that the formula (10) is independent from the length of the perturbation, it becomes clear that the last LE can be calculated using any random vector orthogonalized with respect to z 1 , z 2 , . . . , z n−1 .

Adaptation of the method to non-smooth and discontinuous systems
Estimation of the LLE using the formula (6), or the whole LEs spectrum by means of the equation (10), is possible only if the time derivative of a perturba- when the Jacobian matrix J (x, t) can be calculated at any point of the trajectory. However, at the points where the vector field f is non-smooth or discontinuous, the Jacobian matrix J (x, t) does not exist. Therefore, the procedure of LEs estimation described above must be modified in order to apply it to non-smooth systems.
Assume that a solution of a dynamical system passes through a discontinuity of the vector field f at the times t i , i ∈ 1, . . . , m, such that t 0 < t 1 < t 2 < · · · < t m < t, where t 0 is the initial time and t is the final time. Then, the trajectory is piecewise smooth and the Jacobian matrix J exists in the following time intervals: the initial one [t 0 , t 1 − ], the intermediate ones [t i + , t i+1 − ], i = 1, . . . , m and the final one [t m + , t] for a small number ∈ R + . Between these intervals, contact of the trajectory with the discontinuity occurs in the segments (t i − , t i + ) , i = 1, . . . , m. Although the time spans of these segments can be arbitrarily small, an instantaneous change of the state and the perturbation vectors can occur in them. This division of the trajectory is illustrated in Figure 1.
Please note that, although the scheme and the description involve the first perturbation z 1 , the same applies to the others: z 2 , . . . , z n .
Using the notation proposed in the previous paragraph, the formula (5) can be expanded in the following manner: In the last expression, the terms describing changes of the perturbation's length in the time intervals where the trajectory is smooth, i.e.
alternate with the segments in which the discontinuity appears, i.e. |z 1 (ti+ )| |z 1 (ti− )| . The former ones can be transformed into integrals, according to the formula (5). However, the latter cannot, unless the theory of distributions is involved [33]. Therefore, taking t m+1 = t, the expression (12) takes the following form (13): The first sum in (13), which contains integrals, can be interpreted as the mean value of the expression under the integral, i.e. z 1 (τ ) · dz 1 dτ (τ )/ z 1 (τ ) 2 , at all the points at which this expression can be evaluated. In other words, to calculate the first sum, one can simply use the formula (7) and skip the points in which a discontinuity occurs. The influence of the discontinuities is described by the second sum. In order to estimate it, it is necessary to apply the method described in the paper [28].
In the work [28] a general algorithm, which enables to estimate the transition rule from z 1 (t i − ) to z 1 (t i + ), is established. The method takes into account the first-order approximations of the vector field beyond the discontinuity, as well as the first-order approximations of the functions which define transition of the state vector through the discontinuity. As a result of application of the method, the transition rule in the form of a matrix T (x, t) ∈ R n×n is obtained, so that Although the detailed explanation of the algorithm is beyond the scope of this paper, a practical, illustrative example is provided in the next chapter. For more details on this method, please refer to the original paper [28].
For now, assume that the transition rule has already been established. Its application in the equation (13) yields the expression (14).
Again, the first sum, containing the integrals, can be approximated as the mean value of the expression z 1 (τ ) · dz 1 dτ (τ ) / z 1 (τ ) 2 over the time t, as long as the instances of time in which a discontinuity is encountered t 1 , . . . , t m are skipped. Let t − t 0 = k * ∆t, where ∆t is a small time step and k is the number of iterations of the procedure. Then, the expression (14) can be transformed into (15): The formula (15) constitutes the novel, simplified algorithm of estimation of the LLE in a system with discontinuity. The method works as follows: -integrate the system (1) and the perturbation (2) from the initial time t 0 to the final time t; -after each time step ∆t, in which a discontinuity is not encountered, estimate the value λ * 1 (6); -if the discontinuity is encountered, estimate the value λ * * 1 (16):

2174
The European Physical Journal Special Topics -the final value of the LLE is calculated as a sum of all the λ * 1 and λ * * 1 divided by the number of time steps (17): where the index j is the ordinal number of the LE (for the LLE j = 1).
To calculate the remaining LEs, values λ * j must be estimated by means of the formula (10) instead of (6), as the latter is valid for the LLE only. Moreover, the equation defining λ * * 1 (16) requires its counterpart (18) for the other LEs in the spectrum: where the index i corresponds to the number of contacts with the discontinuity and the index j corresponds to the ordinal number of the LE. Note that the transition rule for the perturbations z 2 , . . . , z n is the same as for the first one, i.e. z j (t i + ) = T (x, t) z j (t i − ), but the perturbation after the transition z j (t i + ) requires orthogonalization by means of the formula (9) before calculating λ * * j (18). Summing up, the final values of all the LEs are estimated using the expression (17) with use of the definitions of λ * j , λ * * j from the equations (10), (18) respectively and the orthogonalization formula (9).

Analysis of the convergence rate and the execution speed
Selection of an appropriate algorithm for an arbitrary purpose requires analysis of its convergence rate and execution speed, so that the result is not only correct, but also obtained in a reasonable time. Therefore, in this subsection, the new method of LEs estimation and the classical one, described by the formula (4) and widely explained in the book [1], are compared from the point of view of calculations efficiency.
Although the classical algorithm [1] was designed for smooth systems, it can be adopted to non-smooth or discontinuous ones in a similar manner as described in the Section 2.4 for the novel method, i.e. the matrix T (x, t) can be applied to predict transition of perturbations through a discontinuity of the system, which enables to use the formula (4) for estimation of LEs. The difference is that the classical algorithm requires integration of n perturbations to calculate the whole spectrum of n LEs, whereas in the new one it is sufficient to integrate n − 1 perturbations. Nevertheless, in both cases the influence of such modification on execution speed of LEs estimation procedures is negligible. This is due to the fact that the only additional computational cost, which appears in discontinuous systems, is multiplication of n perturbation vectors by the n × n matrix T (x, t) at each time when the trajectory encounters a discontinuity of the vector field. The computational cost of such multiplication is not larger than cost of performing one integration step of the equation (2) for each perturbation. Moreover, such transformation of perturbations is exact [28], so it does not impact the convergence rate negatively. Therefore, in the analysis of the convergence rate and the execution speed, it is enough to investigate algorithms for smooth systems, as their adaptation to discontinuous equations hardly increases the computational cost and does not influence the convergence rate.
In order to compare convergence rates of the new algorithm and the classical one [1], it is important to analyze the formulas (5-7), which show a connection between the classical method [1] and the new one. The formula (5) shows that the result obtained from the classical algorithm (4) is equal to the average value of the expression (6), i.e. the time integral of (6) divided by the total time of calculations. The crucial difference between the novel method (7) and the classical one (4) is that the former estimates the integral (5) numerically, as the sum (7), whereas in the latter the integral (5) is computed directly from the antiderivative (4). Therefore, as long as estimation of the integral (5) by means of the sum (7) is accurate enough, both methods should return the same result at each iteration. Consequently, their convergence rates are expected to be equal.
Comparison of execution speed of the novel method and the classical one [1], which does not take into account adaptations to discontinuous systems, has been presented in the publication [12]. It has been shown that the novel method is faster than the original one in most trials. In particular, the novel method proved to be always faster when the order of the system is n ≤ 4. As adaptations to discontinuous systems hardly influence the execution speed of both methods, similar dependencies hold when LEs of non-smooth equations are investigated.

Numerical example 3.1 Description of the system
In order to demonstrate application of the novel method of LEs spectrum estimation for non-smooth and discontinuous systems, an impact oscillator is analyzed. A scheme of the system is presented in Figure 2 below. The oscillator is described by the following set of equations: where X is the coordinate of the oscillator, m is its mass, k is the stiffness of the spring, c is the damping coefficient, F and Ω are the amplitude and the angular frequency of the external forcing respectively, X w is the coordinate of the wall, with which the oscillator can collide elastically, t is the time and t i are the times at which collisions take place. The first relationship defines evolution of the system when X < X w and the second one describes the transition of the state through the discontinuity at X = X w . Introducing dimensionless parameters and transformation of the system to the form (1) yield: T is the state vector and x is its derivative with respect to the dimensionless time τ = t k/m , f is the vector field, x 1 = kX/F ,    (20). position x 1 remains unchanged, whereas the sign of the velocity x 2 is reversed). It is assumed that x 1 ≤ x w for all τ . Any perturbation z 1 of the system (20) evolves in accordance with the equation: where z 1 is the perturbation and J is the Jacobian matrix. The equation (21) is valid when x 1 < x w .

Determining the transition matrix of a perturbation
In order to determine the transition rule, which shows how infinitesimal perturbations change when a trajectory crosses the discontinuity at x 1 = x w in the system (20), the method described in the paper [28] is applied. It is founded on the first-order approximations of both: the vector field near the discontinuity and the formula defining transition of the state. A scheme presenting application of the algorithm to the system (20) is depicted in Figure 3. In Figure 3 x(τ ) is a trajectory, lim T is the state just before collision, τ i is the time at which the current collision occurs, e 1 , e 2 are the unit vectors in the directions of x 1 , x 2 respectively and δ ∈ R + is an arbitrarily small number. Consider the unperturbed state x − and the two states perturbed in subsequent directions of the phase space: The minus sign in the former has been selected to prevent from violating the assumption x 1 ≤ x w . Now, using first order approximations, one can estimate the time at which all the three trajectories, starting from the unperturbed state x − and the perturbed states y − , w − , will cross the discontinuity, and what will be the final states after that time. As Figure 3 indicates, the final state starting from x − is called x + . The perturbed state y − transforms into y + = x + − z 1 and w − changes into w + = x + + z 2 after crossing the discontinuity. The task is to find such matrix T (x − , τ i ) which can transform an arbitrary perturbation just before the crossing the discontinuity z − into the perturbation just after transition through the discontinuity z Obviously, as the unperturbed state T is the state just before collision, it crosses the discontinuity immediately, as well as the second perturbed state w − = x − + δe 2 . However, taking into account the equations of the system (20), it can be noticed that the first perturbed state needs the time equal ∆t = δ/x − 2 to cross the discontinuity. Now, using the first order approximations, the states x + , y + and w + , obtained by integration of (20) over the time ∆t = δ/x − 2 from the initial conditions x − , y − , w − respectively, are estimated in the following manner: Recall that the diagonal matrix D = diag(1, −1) represents transformation of the state vector when a collision occurs. The expressions (22), (24) describe the situations in which an impact occurs immediately, therefore the initial conditions are firstly multiplied by D and then integrated over the time ∆t. On the other hand, the formula (23) concerns the case in which the impact takes place after the time ∆t. Consequently, multiplication by D takes place after integration.
The perturbations after the collision z 1 , z 2 are calculated as follows: Assuming that the transformation of perturbations through the discontinuity is linear, the matrix T (x − , τ i ) can be found. It is expected to have the following property: if z − is a perturbation just before the collision, then the perturbation immediately after the impact is z + = T (x − , τ i ) z − . This should hold for the perturbations analyzed in this subchapter as well, i.e. z 1 = T (x − , τ i ) [−δ, 0] T and T . The former expression shows that the left column of T is z 1 / − δ, whereas the latter informs that the right column of T is z 2 /δ, which results directly from properties of the product of matrices. Therefore, the transition matrix of a perturbation T (x − , τ i ) can be calculated as follows: The limit δ → 0 + is introduced due to the fact the perturbations to be transformed are infinitesimal. The matrix T x − , τ i presented in the expression (27) can be applied to calculate the quantities λ * * j according to the formulas (16) and (18), and thus to estimate the values of the LEs of the system (20) by means of the equation (17).

Technical details of the numerical simulation and its results
A simulation of the system (20) has been created using C++ programming language by means of Code::Blocks IDE. For ODE integration, the odeint library contained in the Boost package has been applied. The RK4 method with the fixed time step equal 5 × 10 −4 has been used. The following values of parameters have been adopted: β = 0.05, x w = 2. The non-dimensional forcing frequency η has been selected as a control parameter in the range η ∈ [0.7, 0.76]. For precise localization of impacts, the linear interpolation method has been implemented [1]. The system (20) has been integrated starting from zero initial conditions. Along with the system (20), a perturbation of its trajectory z 1 has been simulated using the equation (21). Initial values of both components of z 1 have been selected randomly in the interval [0, 1].
As it was explained in the introduction, the proposed method of LEs estimation requires only n − 1 perturbations to calculate the whole spectrum of LEs. Therefore, for the system (20) of the order 2, only one perturbation z 1 needs to be integrated according to the formula (21). In order to calculate the second LE, another vector z 2 has been selected randomly. However, it has not been integrated using the formula (21). It was enough to orthogonalize it with respect to z 1 , which made it adopt the direction of the second LE in the phase space, as it was explained in the introduction.
Before starting calculation of the LEs, the system (20) and the perturbation (21) were simulated over the time τ = 500 × 2π/η to skip transient phenomena. If the norm of the perturbation z 1 .dropped below 10 −8 or exceeded 10 8 , it was normalized. After skipping the transient periods, computation of LEs started. In each integration step, the following actions were made. Firstly, it was checked whether a collision occurred, i.e. whether the expression x 1 − x w changed its sign during the integration step. If not, the program proceeded to calculation of λ * 1 and λ * 2 . The norm of z 1 was checked. The perturbation was normalized when its length dropped below 10 −8 or exceeded 10 8 . Then, the value λ * 1 was computed according to the formula (6) and it was saved. After that, the vector z 2 was orthogonalized with respect to z 1 by means of the equation (8). If its norm dropped below 10 −8 or exceeded 10 8 , it was normalized. Then, the vector z 2 , obtained by orthogonalization of z 2 with respect to z 1 , was used for calculation of λ * 2 from the formula (10). The computed value λ * 2 was kept. On the other hand, if the collision occurred, the precise time of its occurrence τ i , the state just before the collision x − and the perturbation z 1 were found using the linear interpolation method [1]. Using these values, the perturbation