Evaluating noise sensitivity on the time series determination of lyapunov exponents applied to the nonlinear pendulum

This contribution presents an investigation on noise sensitivity of some of the most disseminated techniques employed to estimate Lyapunov exponents from time series. Since noise contamination is unavoidable in cases of data acquisition, it is important to recognize techniques that could be employed for a correct identification of chaos. State space reconstruction and the determination of Lyapunov exponents are carried out to investigate the response of a nonlinear pendulum. Signals are generated by numerical integration of the mathematical model, selecting a single variable of the system as a time series. In order to simulate experimental data sets, a random noise is introduced in the signal. Basically, the analyses of periodic and chaotic motions are carried out. Results obtained from mathematical model are compared with the one obtained from time series analysis, evaluating noise sensitivity. This procedure allows the identification of the best techniques to be employed in the analysis of experimental data.


Introduction
The analysis of chaotic behavior is becoming common in many different fields of science as engineering [1][2][3], medicine [4], ecology [5], biology [6] and economy [7,8]. The study of chaos employs proper mathematical and geometrical aspects and, therefore, new analytical, computational and experimental methods are developed to analyze the response of nonlinear dynamical systems. Alligood et al. [9] say that "of course, the idea of a real experiment being governed by a set of equations is a fiction. A set of differential equations, or a map, may model the process closely enough to achieve useful goals". An approach to deal with the response of dynamical system is based on the analysis of data derived from an experiment [10].
The first problem on the experimental analysis is that data acquisition furnishes a time series of the observable measurements and it is necessary to convert observations into state vectors. Therefore, state space reconstruction needs to be employed and, basically, there are two different methods for this aim: derivative coordinates and delay coordinates [11][12][13]. The method of delay coordinates has proven to be a powerful tool to analyze chaotic behavior of dynamical system. Ruelle [14], Packard et al. [11] and Takens [12] have introduced the basic ideas of this method and one of the drawbacks in its application is the determination of delay parameters.
The other problem in the experimental data is the noise contamination, which is unavoidable in cases of data acquisition. Many studies are devoted to evalu-ate noise suppression and its effects in the analysis of chaos, however, there are a small number of reports devoted to the effects of the system noise on chaos [15].
Nonlinear analysis also involves the determination of quantities, known as dynamical invariants, which are important to identify chaotic behavior. Lyapunov exponent is an example that evaluates the sensitive dependence to initial conditions estimating the exponential divergence of nearby orbits. These exponents have been used as the most useful dynamical diagnostic tool for chaotic system analysis. The signs of the Lyapunov exponents provide a qualitative picture of the system's dynamics and any system containing at least one positive exponent presents chaotic behavior. Lyapunov exponents can also be used for the calculation of other invariant quantities as the attractor dimension, which may be determined by the Kaplan-Yorke conjecture [16]. The determination of Lyapunov exponents of dynamical system with an explicitly mathematical model, which can be linearized, is well established from the algorithm proposed by Wolf et al. [17]. On the other hand, the determination of these exponents from time series is quite more complex. Basically, there are two different classes of algorithms: Trajectories, real space or direct method [17,19,19,20]; and perturbation, tangent space or Jacobian matrix method [21][22][23][24][25][26][27][28].
This article is concerned with the analysis of nonlinear dynamics from time series, and the main purpose is to evaluate noise sensitivity of some of the most disseminated procedures employed either to state space reconstruction or the determination of Lyapunov exponents. Signals are generated by numerical integration of the nonlinear pendulum mathematical model, selecting a single variable of the system as a time series. In order to simulate experimental data sets a random noise is introduced in the signal. State space reconstruction and the determination of Lyapunov exponents are carried out regarding periodic and chaotic signals. The number of data points is chosen as the minimum required for a correct estimation of the desirable measure. Results obtained from mathematical model are compared with the ones obtained from time series analysis, evaluating noise sensitivity. This procedure allows the identification of the best techniques to be applied in the analysis of experimental data.

Nonlinear pendulum
Consider the motion of a nonlinear pendulum where θ defines its position, α is the linear viscous damper pa-rameter and ω n is associated with the natural frequency of the system. Furthermore, a harmonic forcing with amplitude ρ and frequency Ω is considered. With these assumptions, the dynamical system is governed by the well-known equation of motion This equation may be rewritten as a system,u = f (u), u ∈ R 3 , where u 1 = x = θ, u 2 = y =θ and u 3 = Ωt. Numerical simulations are carried out employing the fourth-order Runge-Kutta method with time steps less than ∆t = 2π/100Ω. For all simulations, parameters α = 0.2 and ω n = Ω = 1.0 are considered. In order to simulate experimental data sets, a signal s = x+η is defined where η = AR(−1, +1) represents noise, with A being the amplitude, and R(−1, +1) is related to random number within the interval (−1, +1). If η = 0, there is no noise and an ideal experimental data is simulated. In this article, two other noise levels are contemplated: A = 0.314 and A = 0.628, representing, respectively, 5% and 10% of the maximum signal amplitude. The number of data points, N , is chosen as the minimum required for a correct estimation of the desirable measure.
Basically, periodic and chaotic signals are analyzed. When the forcing amplitude is ρ = 2.56, the pendulum presents a period-2 motion. Figure 1 shows the steady state orbit on phase space for this motion projected on a cylindrical space and on the plane. When ρ = 2.50, the motion becomes chaotic. Figure 2 presents the strange attractor of the motion projected on cylindrical and plane spaces.

State space reconstruction
The basic idea of the state space reconstruction is that a signal contains information about unobserved state variables, which can be used to predict the present state. Therefore, a scalar time series, s(t), may be used to construct a vector time series that is equivalent to the original dynamics from a topological point of view. The state space reconstruction needs to form a coordinate system to capture the structure of orbits in state space. The method of delay coordinates could be done using lagged variables, s(t + τ ), where τ is the time delay. Then, considering an experimental signal, s(t), where t = t 0 + (n − 1)∆t with n = 1, 2, 3, . . . , N, it is possible to use a collection of time delays to create a vector in a D e -dimensional space, u(t), which represents the reconstructed dynamics of the system.
The method of delays was first proposed by Ruelle [14] and Packard et al. [11] and then by Takens [12] and Sauer et al. [29]. Its use has become popular for dynamical reconstruction, however, the choice of the delay parameters, τ -time delay, and D e -embedding dimension, has not been fully developed. This article employs the average mutual information method to determine time delay [30] and the method of false nearest neighbors to estimate embedding dimension [31].
In order to start the analysis of the nonlinear pendulum state space reconstruction employing the method of delay coordinates, a period-2 signal with N = 20,000, is considered. Figure 3 presents results of the mutual information and the false nearest neighbors analysis, for different noise levels. Concerning the time delay determination, there is a difficulty to determine the first minimum of the information curve when A = 0, ideal signal. Nevertheless, time delay may be estimated Following the determination of delay parameters, the method of delay coordinates can be applied in order to reconstruct the state space. The numerical state space of the motion is presented in Fig. 4 together with the reconstructed spaces for different noise levels. The comparison among numerical and reconstructed state spaces allows one to observe just a small coordinate change from one to another.
The forthcoming analysis regards a chaotic signal with N = 20,000. Figure 5 considers results of the mutual information and the false nearest neighbor analysis, for different noise levels. The same procedure applied in the determination of time delay of periodic signal can be employed here, resulting on the following values: τ = 2.262 s when A = 0; τ = 1.885 s when A = 0.314; τ = 1.822 s when A = 0.628. Once again, the analysis of the embedding dimension estimates D e = 3 and the noise does not have significant influence on the results. Figure 6 presents the Poincaré map obtained either by numerical simulation or by reconstruction using the method of delay coordinates for three noise levels: A = 0, A = 0.314 and A = 0.628. A strange attractor is clearly identified, presenting a fractal-like structure. The comparison among numerical and reconstructed state spaces allows one to observe just a small coordinate change from one to another.

Lyapunov exponents
Lyapunov exponents have been used as the most useful dynamical diagnostic tool for chaotic system analysis. These exponents evaluate the sensitive dependence to initial conditions estimating the exponential divergence of nearby orbits. The dynamics of the system transform a D-sphere of states in a D-ellipsoid and, when there is a chaotic motion, a complex evolution exists. Lyapunov exponents are related to the expanding and contracting nature of different directions in phase space and the signs of these exponents provide a qualitative picture of the system's dynamics. The existence of positive Lyapunov exponents defines directions of local instabilities in the dynamics of the system and any system containing at least one positive exponent presents chaotic behavior.
The determination of Lyapunov exponents of dynamical system with an explicitly mathematical model, which can be linearized, is well established from the algorithm proposed by Wolf et al. [17]. On the other hand, the determination of these exponents from time series is quite more complex. Basically, there are two different classes of algorithms: Trajectories, real space or direct method; and perturbation, tangent space or Jacobian matrix method.
Trajectories method has been originally developed by Wolf et al. [17] and the basic idea associated with it is analyzing the evolution of two nearby orbits in the tangent space. This method calculates only the largest exponent, which is sufficient to identify chaotic behavior but, on the other hand, cause some difficulties when the determination of other quantities is needed. Other similar algorithms have also been developed exploiting the same idea [18,19]. These algorithms consider that the divergence rate trajectories fluctuates along the trajectory, with the fluctuation given by the spectrum of effective Lyapunov exponents. The average of effective Lyapunov exponent along the trajectory is the true Lyapunov exponent and, therefore, it may be calculated from the slope of a curve associated with a greater instability direction. Tangent space methods seem to be most promising for the calculation of the Lyapunov spectrum from time series. The product of the Jacobians along the trajectory can determine the spectrum of Lyapunov exponents from the evaluation of eigenvalues of this matrix. To make use of this, the map x(t + 1) = f (x(t)) must at least approximately be known and there are several approaches to extract this map from a time series. Sano and Sawada [21] and Eckmann et al. [22] have developed similar algorithms where the Jacobian matrix is evaluated with a least square error algorithm. Brown et al. [27] and Briggs [24] consider a high order polynomial approximation to define the Jacobian matrix. Kruel et al. [25] improve the algorithm due to Sano and Sawada with a different form to evaluate one of the covariant matrix that generates the Jacobian matrix.
In the present contribution, Lyapunov exponents are determined employing the algorithms proposed by Wolf et al. [17], Kantz [18], Rosenstein et al. [19] and Sano and Sawada [21]. Algorithms developed by Hegger et al. [32] are employed for all simulations except the ones due to Wolf et al. [17]. Results are compared with reference values calculated from the algorithm for differential equations proposed by Wolf et al. [17].

Algorithm due to Wolf et al.
Trajectories method proposed by Wolf et al. [17] to determine Lyapunov exponents from time series considers the reconstructed attractor and examines orbital divergence on length scales, working in tangent space. The method monitors the long-term evolution of a single pair of nearby orbits and is able to estimate non-negative Lyapunov exponents. In principle, this method permit to compute all Lyapunov spectrum but in reality it is limited to the maximum one [18]. Wolf et al. adopts the following definition of Lyapunov exponents: where M is the total number of replacement steps. The distance L(t) is defined as the Euclidean norm between two points. Besides reconstruction parameters, the algorithm has three parameters to be determined, which have a great influence on the exponent calculation: the signal length, N ; the smaller distance between the trajectories, d min , which is associated with noise; and a constant propagation time, t EV . In order to start the analysis employing the algorithm due to Wolf et al., a period-2 signal with N = 10,000 is considered. For this situation, the reference value calculated from the differential equation [17] is λ max = 0. Figure 7(a) shows the maximum Lyapunov exponent for three different noise levels assuming d min = 0.0001 and t EV = 1. The time series algorithm for the ideal signal furnishes λ max = 0. For noisy signals, however, the algorithm is not able to identify periodic response meaning that noise and chaos cannot be distinguished ( Fig. 7(a)). Nevertheless, it should be pointed out that the alteration of parameters t EV and d min , improve the results and the exponent values converge to zero, as it is desirable. To elucidate the influence of this alteration, Fig. 7(b) shows results when t EV = 15 for A = 0.314 and t EV = 30 for A = 0.628. Even though a convenient variation of these parameters can reduce the noise effect, notice that this procedure is difficult to be employed on experimental data when a reference value is not known.
At this point, a chaotic signal with N = 30,000 is considered using d min = 0.0001 and t EV = 1. Lyapunov exponents predicted by the algorithms for differential equation and time series are presented in Fig. 8. The first algorithm furnishes λ max = +0.16, which is considered as a reference value. On the other hand, the algorithm for time series applied to the ideal signal furnishes λ max = +0.39. The maximum Lyapunov exponent assumes greater values when noise level grows, exhibiting the difficulty to its estimation. The alteration of the parameter t EV allows one to reduce the discrepancy among the previous results. By considering t EV = 5 when A = 0.314 and t EV = 10 when A = 0.628, values of maximum exponents converge to the reference value. Again, this procedure is not satisfactory to identify chaotic motion on experimental data when a reference value is not known.  Rosenstein et al.

Algorithm due to Kantz and due to
The algorithm proposed by Kantz [18] considers the same idea of the one proposed by Wolf et al. [17], establishing that the divergence rate trajectories fluctuates along the trajectory, with the fluctuation given by the spectrum of effective Lyapunov exponents. The average of effective Lyapunov exponent along the trajectory is the true Lyapunov exponent and the maximum exponent is given by where |u(0) − u ε (0)| = ε and u(t) − u ε (t) = εv u (t), with v u (t) representing the eigenvectors associated with the maximum Lyapunov exponent, λ max . From an arbitrary point u t = (s t−m+1 , . . . , s t ), all delay vectors of the series falling into the εneighborhood of u t will be regarded as the beginning of the neighboring trajectories, which are simply given by the points of the time series consecutive in time. If the distance between neighboring trajectories is measured in their true phase space, it is possible to observe the fluctuations of the divergence rate described by the distribution of effective Lyapunov exponents. In order to evaluate the distance between the reference trajectory, u t , and a neighbor, u i , after a relative time δ referring to the time index of the point where the distance begin to be greater than ε, δ(0), one defines a distance, dist(u t , u i , δ) = |u t+δ − u i+δ |, which represents the modulus of the δth scalar component of the two trajectories. These distances are projections of the difference vectors in the true phase space onto a one-dimensional subspace spanned by the observable. Therefore, they are modulated with | cos φ|, where φ is the angle between the eigenvector corresponding to λ max and the local direction of the subspace on which the observable lives. Taking into account the time average over the full length of the time series, and also considering the logarithm of the average distance, the Lyapunov exponents are computed by, Initially, the difference vectors in the true phase space are pointing in a general direction, however, for intermediate range of δ, S(δ) increases linearly with the slope λ which corresponds an estimation of the Lyapunov exponent.
Rosenstein et al. [19] have proposed a similar algorithm where the distance between the trajectories is defined as the Euclidean norm in the reconstructed phase space and, also, they have used only one neighbor trajectory. Both algorithms present good results for discrete systems (maps), nevertheless, results are not so good for continuous systems (differential equations). Notice, however, that a continuous system can be converted to a discrete system regarding a time series that is the Poincaré map of the signal. Both algorithms needs the following parameters, besides reconstruction parameters: the signal length, N ; and ε related to the trajectory neighbor. Furthermore, Kantz's algorithm considers a range where parameter ε is varied.
At this point, results predicted by algorithms proposed by Kantz and by Rosenstein et al. are focused. The analysis is done regarding time series that represent Poincaré maps of signals where ∆T = 6.28 s, resulting in N = 2, 000. Since embedding theorems [29] establishes that embedding dimension needs to be greater than 2D e + 1 in order to obtain an appropriate reconstruction, one consider different values of this parameter in order to obtain better results. After several tests, D e = 6 is adopted and ε varies from 20 to 25 in Kantz's algorithm. Firstly, period-2 signals with different noise levels are contemplated. The reference value associated with these signals is λ max = 0, calculated from the differential equation. Figure 9 shows the S(δ) curves predicted by both algorithms presenting the same horizontal curves. The slope of these curves defines the Lyapunov exponent, meaning that λ max = 0, the same value predicted by the algorithm due to Wolf et al. for ideal signal. Notice, however, that noise does not influence results, which is an important improvement to the algorithm due to Wolf et al. In addition, it should be emphasized that either the number of data points, N , or the embedding dimension, D e , may alter the slope of the curve S(δ). Therefore, the determination of these parameters may introduce difficulties to evaluate Lyapunov exponents when a reference value is not known.
A chaotic signal is now in focus. Figure 10 presents the curves S(δ) for an ideal chaotic signal (A = 0), with N = 2, 000, D e = 6 and ε varying from 20 to 25 in Kantz's algorithm. The S(δ) curve presents a linear range which tends to reach a stabilized value. The slope of the curve in this linear range estimates the maximum Lyapunov exponent and may be computed employing a linear regression. The algorithm due to Kantz predicts λ max = +0.184 ± 0.004, while the algorithm due to Rosenstein   +0. 16.
The S(δ) curves for noisy signals are presented in Fig. 11. Employing linear regression, Kantz's algorithm predicts λ max = +0.202 ± 0.002 for both noise levels, while Rosenstein's predicts λ max = +0.185 ± 0.001. These results are again close to the reference value λ max = +0. 16 and one can conclude that noise does not have a significant influence on the results.

Algorithm due to Sano and Sawada
The algorithm proposed by Sano and Sawada [21] is a tangent space method that is based on the multiplicative ergodic theorem by Oseledec [33]. One of the   advantages of tangent space methods is the estimation of the Lyapunov spectrum, in contrast of the maximum exponent calculated by trajectories methods. Before proceeding to a description of the algorithm due to Sano and Sawada, one considers a trajectory s(t), in a Ddimensional space, which is the solution of an ordinary differential equation, The evolution of a tangent vector ξ, in a tangent space at s(t), is represented by linearizing the differential equation,ξ = Bξ, where B is the Jacobian matrix of f . The solution of this equation is expressed as ξ(t) = B t ξ(0), where B t is the linear operator which maps tangent vector ξ(0) to ξ(t). Sano and Sawada have established a procedure for estimating a linearized flow map B t from the observed data. Hence, Lyapunov exponents can be computed as whereB is an approximation of the flow map B t and {e j i } is a set of basis vectors of tangent space. Be-sides signal and reconstruction parameters, the algorithm needs the following parameters: R 0 , radius of the sphere and N NEIGH , number of neighboring points. In order to start the analysis of the algorithm due to Sano and Sawada, an ideal period-2 signal (A = 0) with N = 70, 100 is considered, assuming N NEIGH = 30 and R 0 = N/1000. Notice that this algorithm needs large data points to estimate the Lyapunov spectrum. Figure 12 shows the spectrum for this signal where all the exponents tend to zero, meaning that the motion is periodic. Noisy signals are now considered (A = 0.314 and A = 0.628). Figure 13 presents Lyapunov spectra calculated with N NEIGH = 30 and R 0 = N/1000. Notice that as the noise level increases, values of Lyapunov exponents decrease. Since results converge to negative values, the algorithm is capable to identify a periodic motion. Nevertheless, this characteristic causes problems on the determination of other dynamical invariants.
The forthcoming analysis considers chaotic signals and, at first, an ideal signal (A = 0) with N = 70, 100, is contemplated. Figure 14 shows the Lyapunov spectrum calculated with N NEIGH = 30 and R 0 = N/1000. The maximum exponent converges to a positive value, λ 1 = +0.016, which is ten times less than the reference value (λ max = +0.16). Again, this is not a problem to identify chaotic motion but for the determination of the dynamical invariants. Figure 15 shows Lyapunov spectra for noisy chaotic signals (A = 0.314 and A = 0.628) with N = 70, 100, and assuming N NEIGH = 30 and R 0 = N/1000. Once again, as the noise level increases, values of Lyapunov exponents decrease and, since results converge to negative values, the algorithm is not capable of distinguish noise from chaos, presenting noise sensitivity.

Conclusions
This contribution evaluates noise sensitivity of some of the most disseminated techniques employed on the state space reconstruction and the determination of Lyapunov exponents from time series. Signals are generated by numerical integration of the nonlinear pendulum mathematical model, selecting a single variable of the system as a time series. In order to simulate experimental data sets, a random noise is introduced in the time series. Results obtained from mathematical model are compared from the one obtained from time series analysis, evaluating noise sensitivity. State space reconstruction is done employing the method of delay coordinates. The determination of delay parameters, time delay and embedding dimension, are made employing, respectively, the method of average mutual information and the false nearest neighbors. Both methods present good results and are not noise sensitive. Lyapunov exponents are calculated employing four different algorithms: Wolf et al., Kantz, Rosenstein et al. and Sano and Sawada. Algorithms due to Wolf et al. and due to Sano and Sawada are noise sensitive and are not capable to distinguish noise from chaos. On the other hand, algorithms due to Kantz and due to Rosenstein et al. are not noise sensitive. Nevertheless, the use of these algorithms present difficulties on the determination of parameters, especially when there is no reference value. The authors agree that this contribution is useful to identify the best techniques that may be applied on experimental analysis, however, the inves-tigation of other physical systems and different kinds of noise are necessary to confirm these conclusions. Recently, the authors apply these conclusions in time series obtained from an experimental pendulum [34] showing that Kantz's algorithm allows one to establish a difference between periodic and chaotic motion. The algorithm due to Rosenstein et al., on the other hand, does not present good results for periodic motion.