Research Article  Open Access
On the Use of Interval Extensions to Estimate the Largest Lyapunov Exponent from Chaotic Data
Abstract
A method to estimate the (positive) largest Lyapunov exponent (LLE) from data using interval extensions is proposed. The method differs from the ones available in the literature in its simplicity since it is only based on three rather simple steps. Firstly, a polynomial NARMAX is used to identify a model from the data under investigation. Secondly, interval extensions, which can be easily extracted from the identified model, are used to calculate the lower bound error. Finally, a simple linear fit to the logarithm of lower bound error is obtained and then the LLE is retrieved from it as the third step. To illustrate the proposed method, the LLE is calculated for the following wellknown benchmarks: sine map, Rössler Equations, and MackeyGlass Equations from identified models given in the literature and also from two identified NARMAX models: a chaotic jerk circuit and the tent map. In the latter, a Gaussian noise has been added to show the robustness of the proposed method.
1. Introduction
One of the most applied techniques to confirm whether a system is chaotic or not is based on the estimation of the Lyapunov exponents [1]. The calculation of Lyapunov exponents is usually based on the average divergence or convergence of close trajectories along certain directions in state space. In chaotic systems, two trajectories separate exponentially with time despite having very similar initial conditions. Since the work of [2] several numerical methods to estimate Lyapunov exponents were proposed [3–8], just to mention a few. A comparison between estimation methods can be found in [9], for instance. The authors in [5] developed the first practical algorithm to calculate the Lyapunov exponents by estimating the growth of the corresponding set of vectors as the system evolves. This method allows the estimation of the complete spectrum of Lyapunov exponents. Amongst these exponents, the (positive) largest Lyapunov exponent (LLE) is the exponent considered to be the main reason for the separation rate. Therefore the estimation of such an exponent is used to build up the chaotic nature of the data under scrutiny. The authors in [1, 7] independently used statistical properties of the local divergence rates of nearby trajectories of the systems under investigation to estimate the LLE. These properties come in handy when the goal is to get high exactness of numerical estimates of the LLE, as shown in the work [10]. Recent and different applications using the LLE can be found in [11–14].
In [15] a straightforward method to compute the LLE using interval extensions and the original equations of motion was proposed. The exponent is estimated from the slope of the line derived from the lower bound error when considering two interval extensions of the original system [16, 17]. The authors of [15] demonstrated that the method is quick and simple to be used and could be considered as an alternative to other algorithms available in the literature. However, it requires the original equations of motion of the system and it can not be applied directly to time series. In order to overcome this limitation, a system identification methodology based on the wellknown polynomial NARMAX (Nonlinear AutoRegressive Moving Average model with eXogenous inputs) [18] approach is introduced. Such an approach (system representation) has been shown to adequately represent a large number of chaotic systems [18, 19]. A similar method using the same system representation was proposed in [20] but, in their procedure, the Jacobian matrix must be found and its values estimated. For polynomial NARMAX it was argued that the Jacobian matrix is easily obtained but that may not be the case when other types of representation, such as Neural Networks [20], are considered. Using the interval extensions, it is shown that it is possible to put together the advantages pointed out in [20] without the need to calculate the Jacobian matrix to estimate the LLE. Moreover, as NARMAX is based on a polynomial structure, interval extensions are easily generated by simple arithmetic manipulations, such as commutative and associative. In addition to that, as indicated in [15], the proposed method requires similar or less points to achieve convergence when compared to other methods available in the literature.
The rest of this paper is laid out as follows. Preliminary concepts are introduced in Section 2 and it gives the necessary background material to understand the main idea of the work. Then, the proposed method is presented in Section 3. Section 4 presents the results obtained by the proposed approach when applying it to three different systems that have been identified in the literature: sine map, Rössler Equations, and MackeyGlass Equations. Moreover, two identified polynomial NARMAX models have been considered: a chaotic jerk circuit [21] and the tent map. In the end, a Gaussian noise has been added to demonstrate the robustness of the proposed technique. Finally, the conclusions are given in Section 5.
2. Preliminary Concepts
In this section a brief review on the calculation of the LLE is given along with the recent result on the lower bound error and the interval extensions.
2.1. The LLE Calculation
To review the basic ideas laid out in [7], consider the following equation:where is the set of all other delay vectors in an neighbourhood of the vector (data from trajectories of the system under investigation) and is the number of elements in . The initial and final point of each window over the time series is given by and , respectively. The LLE can be estimated by searching for a linear scaling in plot versus . In order to find the linear scaling, it is necessary to tune the parameter , to find the neighbourhood, and to define the number of elements in . This will not be necessary for the method proposed in Section 3; however there is the need of finding the equations of motion of the system under investigation and this will be accomplished by using a system identification methodology as reviewed in the next section.
2.2. The Polynomial NARMAX
A NARMAX model can be written as follows [22]:where , , and are, respectively, the output, input, and noise terms at the discrete time . The parameters , , and are the maximum lag considered for output, input, and noise, respectively. Noise terms are frequently included in the parameter estimation process to avoid bias in the estimates. Here is assumed to be a polynomial with nonlinearity degree . Nonlinear systems are frequently modelled using particular cases of the polynomial NARMAX, such as NAR, NARX, and NARMA [16, 23–25]. In this paper, they all are treated under the name polynomial NARMAX. Further details of this mathematical representation can be found in [18].
The great advantage of using the polynomial NARMAX representation is that the problem of finding the relevant terms and the estimates of their coefficients can be cast into a linear regression problem and therefore the wellknown least squares method can be applied. In the eighties, Billings [18] devised an efficient method to take the advantage of the classical methods for solving the least squares problem to include the structure selection using the socalled ERR (Error Reduction Ratio) which attempts to measure the contribution of each term to the overall variance of the output signal. Such a method was called OLS (Orthogonal Least Squares method) and has been successfully used to obtain dynamically valid models from data generated by nonlinear systems behaving chaotically (see, e.g., [26, 27]). Additionally, it is straightforward to generate interval extensions from the polynomial NARMAX, as seen in Definition 1.
2.3. The Lower Bound Error
This subsection is an adaptation from [16, 17] on the lower bound error applied to continuous nonlinear systems.
Let , a metric space ; the relationwhere , is a recursive function or a map of a state space into itself and denotes the state at the discrete time . The sequence obtained by iterating (3) starting from an initial condition is called the orbit of [28, 29].
Definition 1 (see [16]). An interval extension of is an interval valued function of an interval variable , with the propertywhere an interval is a closed set of a real number , where
Example 2. Let us consider . Some examples of interval extensions are
Equations (6) are mathematically equivalent but they present different sequence of their basic arithmetic operations. On floating point standard neither commutative nor distributive properties hold [30, 31]. Therefore, (6) may exhibit different results when a computerbased approach is used to estimate their values. We have exploited these features in this work to estimate the LLE.
Definition 3. Let represent a pseudoorbit, which is defined by an initial condition, an interval extension of , some specific hardware and software, numerical precision standard, and discretization scheme. A pseudoorbit is an approximation of an orbit and can be represented by or , such thatwhere is the error and .
A pseudoorbit characterises an interval where the genuine orbit lies. Henceforth an interval related to each estimation of a pseudoorbit is given byFrom (7) and (8) it is clear that
Theorem 4 [16] establishes that at least one of two pseudoorbits must have an error greater than or equal to the lower bound error.
Theorem 4. Let two pseudoorbits and be derived from two interval extensions. Let be the lower bound error of a map ; then or .
3. The Proposed Method
The method proposed in this work can be summarized by the following steps:(1)Fit a model to the data series. Polynomial NARMAX models were the choice in this work.(2)Choose two different interval extensions of the system under investigation following Definition 1.(3)With exactly the same initial conditions simulate the two interval extensions.(4)Use the least squares method to fit a line to the slope of the logarithm curve of the absolute value of the lower bound error defined in Section 2.3. The slope of the line is the LLE.
To the interested reader, a detailed description and explanation of why this simple method works may be found in [15, 32].
4. Results
In this section, some examples are given to illustrate the usefulness of the proposed methods.
4.1. Example 1: Sine Map
A unidimensional sine map is defined as
The LLE from the original equation (10) was calculated using the method proposed in [15]. To this end, two interval extensions given by the following equations were considered:
The curve of the logarithm of the error bound is depicted in Figure 1(a) and slope of the linear regression, that is, the LLE, is shown to be 1.133. This value is in good agreement with the values found in the literature such as the one given in [33].
(a) Simulation of (11), with results for and and the same initial condition; that is, . LLE calculated is the inclination of the curve given by 1.133.
(b) Simulation of (13), with results for and and the same initial condition; that is, . LLE calculated is the inclination of the curve given by 1.114.
From the data generated by iterating the model in (10) the following polynomial NARMAX given in [34] was identified:Considering two equivalent interval extensions of the identified model (12) we haveFigure 1(b) shows the logarithm of the lower bound error for the two interval extensions given by (13). Table 1 summarizes the results. Note that the similarity of the LLE can be also used to validate the model.

4.2. Example 2: Rössler Equations
Rössler’s Equations [35] are given byThis system has been identified in [33] and the resulting polynomial NARMAX model is given by the following equation:
Equation (15) was considered as the first interval extension. The second interval extension is determined by rewriting the termin (15) asfor instance. Several other interval extensions could be determined but two of them would suffice for the purposes of finding the LLE.
As can be seen the interval extensions used to simulate the systems are mathematically equivalent, but they differ in the sequence of arithmetical operations. Figure 2 shows the logarithm of the lower bound error of the simulation using the two specified interval extensions. Table 2 summarizes the result of LLE for Rössler Equations. The literature suggests a value of 1.242 and the estimated value using the model procedure suggested in [33] is . The estimated value seems slightly different, but if the confidence interval to calculate the LLE using data is taken into consideration, one could consider instead and therefore the value of LLE calculated in [33] is in good agreement with the value in the literature. Similar reasoning can be applied to the value estimated using the method presented in this paper, that is, 1.530.
4.3. Example 3: MackeyGlass Equations
The MackeyGlass is an interesting system used in many papers as an example of a chaotic and infinite dimensional system, since it is a timedelay system [36, 37]. The system equation is given bywhere , , , and . The authors in [20] identified the following model to this system using a sequence of 234 data points sample at :
The interval extension used in this example is obtained just by replacing the term by which means that only the order of multiplication has been changed. Figure 3 presents the logarithm of the error for the MackeyGlass Equation. Using the method proposed here the LLE is found to be exactly the same as the value provided by [20], as reported in Table 3. It is worth emphasizing that there is no need to calculate the Jacobian matrix of the model in the proposed method.

4.4. Example 4: A Chaotic Jerk Circuit
A chaotic jerk circuit given in [21] was also used to further illustrate the proposed approach. The circuit is composed of resistors, capacitors, and operational amplifiers as shown in Figure 4. In this work, a polynomial NAR with nonlinearity degrees and was used to identify the variable of the chaotic jerk circuit, as shown in Figure 5. The data was produced by means of the simulation of the electronic circuit using LTspice. The obtained model has terms, which have been omitted in this manuscript for the sake of simplicity.
For this case, the interval extensions can be found by changing the order in which the model terms were computed. The two interval extensions were generated by the simple order commutation of elements in the polynomial, as illustrated in Example 2. The interval extensions used to simulate the systems are mathematically equivalent, but are not computationally equivalent. Figure 5 shows the logarithm of the lower bound error of the simulation using the two specified interval extensions and Table 4 summarizes the results of the LLE for the chaotic jerk circuit. For comparison purposes [21] shows that, by using experimental data, the value of LLE can be estimated as . On the other hand, the proposed methodology estimates the LLE value as which is in a good agreement with the value for the original system.

4.5. Example 5: The Tent Map
To verify the robustness of the proposed approach under the presence of noise, a tent map (see (22)) was used as a system to be identified, within an additive zero mean Gaussian noise with standard deviation .
A polynomial NAR was identified with and , yielding the model The interval extensions can be found by changing the order in which two terms were computed. In this example, the interval extensions were generated by the simple position change of the terms and Figure 6 presents the logarithm of the error for the identified model whereas Figure 7 shows the first return map of both identified and original tent map.
Despite the presence of noise affecting the system, the LLE found by the method proposed in this paper was in good agreement with the one analytically computed from the tent map, as shown in Table 5.

5. Conclusion
This paper presents a system identification approach to calculate the LLE based on the interval extensions. The method proposed here can be seen as an extension of the method developed in [15] in which the LLE is calculated by taking into consideration the properties of the exponential growing of the lower bound error [16]. The method is designed specifically for the calculation of the LLE from chaotic data, which was not possible to undertake with method proposed in [15]. A polynomial NARMAX model was used to identify the system, but other representations could also be used in the first step. Nevertheless, the polynomial NARMAX represents a very suitable mathematical representation to produce interval extensions, which can be obtained by simple manipulation upon the terms of the identified polynomial model. Compared to existing approaches, such as the one indicated in [20], the proposed method does not require the calculation of the Jacobian matrix and, more importantly, is very simple to use. We believe that the proposed method is also a good alternative for computing Lyapunov exponents from limited data. However, the computation of Lyapunovexponent spectrum, as performed in [38], is left for future investigations.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
References
 M. T. Rosenstein, J. J. Collins, and C. J. De Luca, “A practical method for calculating largest Lyapunov exponents from small data sets,” Physica D: Nonlinear Phenomena, vol. 65, no. 12, pp. 117–134, 1993. View at: Publisher Site  Google Scholar  MathSciNet
 V. I. Oseledec, “The multiplicative ergodic theorem: the lyapunov characteristic numbers of dynamical systems,” Transactions of the Moscow Mathematical Society, vol. 19, pp. 179–210, 1968. View at: Google Scholar  MathSciNet
 G. Benettin, L. Galgani, A. Giorgilli, and J.M. Strelcyn, “Lyapunov Characteristic Exponents for smooth dynamical systems and for hamiltonian systems; A method for computing all of them. Part 2: Numerical application,” Meccanica, vol. 15, no. 1, pp. 21–30, 1980. View at: Publisher Site  Google Scholar
 G. Benettin, L. Galgani, A. Giorgilli, and J.M. Strelcyn, “Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part 1: theory,” Meccanica, vol. 15, no. 1, pp. 9–20, 1980. View at: Publisher Site  Google Scholar
 A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano, “Determining Lyapunov exponents from a time series,” Physica D: Nonlinear Phenomena, vol. 16, no. 3, pp. 285–317, 1985. View at: Publisher Site  Google Scholar  MathSciNet
 M. Sano and Y. Sawada, “Measurement of the Lyapunov spectrum from a chaotic time series,” Physical Review Letters, vol. 55, no. 10, pp. 1082–1085, 1985. View at: Publisher Site  Google Scholar  MathSciNet
 H. Kantz, “A robust method to estimate the maximal Lyapunov exponent of a time series,” Physics Letters A, vol. 185, no. 1, pp. 77–87, 1994. View at: Publisher Site  Google Scholar
 G. Rangarajan, S. Habib, and R. D. Ryne, “Lyapunov exponents without rescaling and reorthogonalization,” Physical Review Letters, vol. 80, no. 17, pp. 3747–3750, 1998. View at: Publisher Site  Google Scholar
 K. Geist, U. Parlitz, and W. Lauterborn, “Comparison of different methods for computing Lyapunov exponents,” Progress of Theoretical and Experimental Physics, vol. 83, no. 5, pp. 875–893, 1990. View at: Publisher Site  Google Scholar  MathSciNet
 B. J. Kim and G. H. Choe, “High precision numerical estimation of the largest Lyapunov exponent,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 5, pp. 1378–1384, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 G. Wainrib and M. N. Galtier, “A local Echo State Property through the largest Lyapunov exponent,” Neural Networks, vol. 76, pp. 39–45, 2016. View at: Publisher Site  Google Scholar
 A. J. Steyer and E. S. Van Vleck, “A stepsize selection strategy for explicit RungeKutta methods based on Lyapunov exponent theory,” Journal of Computational and Applied Mathematics, vol. 292, pp. 703–719, 2016. View at: Publisher Site  Google Scholar  MathSciNet
 C. J. GavilánMoreno and G. EspinosaParedes, “Using Largest Lyapunov Exponent to Confirm the Intrinsic Stability of Boiling Water Reactors,” Nuclear Engineering and Technology, vol. 48, no. 2, pp. 434–447, 2016. View at: Publisher Site  Google Scholar
 W. Caesarendra, B. Kosasih, A. K. Tieu, and C. A. S. Moodie, “Application of the largest Lyapunov exponent algorithm for feature extraction in low speed slew bearing condition monitoring,” Mechanical Systems and Signal Processing, vol. 5051, pp. 116–138, 2015. View at: Publisher Site  Google Scholar
 E. M. A. M. Mendes and E. G. Nepomuceno, “A very simple method to calculate the (Positive) largest lyapunov exponent using interval extensions,” International Journal of Bifurcation and Chaos, vol. 26, no. 13, Article ID 1650226, 2016. View at: Publisher Site  Google Scholar
 E. G. Nepomuceno and S. A. M. Martins, “A lower bound error for freerun simulation of the polynomial NARMAX,” Systems Science & Control Engineering, vol. 4, no. 1, pp. 50–58, 2016. View at: Publisher Site  Google Scholar
 E. Nepomuceno, S. Martins, G. Amaral, and R. Riveret, “On the lower bound error for discrete maps using associative property,” Systems Science Control Engineering, vol. 5, no. 1, pp. 462–473, 2017. View at: Publisher Site  Google Scholar
 S. A. Billings, Nonlinear system identification: NARMAX methods in the time, frequency, and spatiotemporal domains, West Sussex: John Wiley & Sons, London, UK, 2013. View at: Publisher Site  MathSciNet
 L. A. Aguirre and C. Letellier, “Modeling nonlinear dynamics and chaos: A review,” Mathematical Problems in Engineering, vol. 2009, Article ID 238960, 2009. View at: Publisher Site  Google Scholar
 L. A. Aguirre and S. Billings, “Retrieving dynamical invariants from chaotic data using narmax models,” International Journal of Bifurcation and Chaos, vol. 05, no. 02, pp. 449–474, 1995. View at: Publisher Site  Google Scholar
 J. C. Sprott, “A new chaotic jerk circuit,” IEEE Transactions on Circuits and Systems II: Express Briefs, vol. 58, no. 4, pp. 240–243, 2011. View at: Publisher Site  Google Scholar
 S. Chen and S. A. Billings, “Representations of nonlinear systems: the {NARMAX} model,” International Journal of Control, vol. 49, no. 3, pp. 1013–1032, 1989. View at: Publisher Site  Google Scholar  MathSciNet
 L. Piroddi and W. Spinelli, “An identification algorithm for polynomial {NARX} models based on simulation error minimization,” International Journal of Control, vol. 76, no. 17, pp. 1767–1781, 2003. View at: Publisher Site  Google Scholar  MathSciNet
 J. R. Ayala Solares, H.L. Wei, and S. A. Billings, “A novel logisticNARX model as a classifier for dynamic binary classification,” Neural Computing and Applications, pp. 1–15, 2017. View at: Publisher Site  Google Scholar
 O. M. Mohamed Vall and R. M'hiri, “An approach to polynomial NARX/NARMAX systems identification in a closedloop with variable structure control,” International Journal of Automation and Computing, vol. 5, no. 3, pp. 313–318, 2008. View at: Publisher Site  Google Scholar
 E. M. Mendes and S. A. Billings, “On Identifying Global Nonlinear Discrete Models from Chaotic Data,” International Journal of Bifurcation and Chaos, vol. 07, no. 11, pp. 2593–2601, 1997. View at: Publisher Site  Google Scholar
 L. A. Aguirre, G. G. Rodrigues, and E. M. Mendes, “Nonlinear Identification and Cluster Analysis of Chaotic Attractors from a Real Implementation of Chua's Circuit,” International Journal of Bifurcation and Chaos, vol. 07, no. 06, pp. 1411–1423, 1997. View at: Publisher Site  Google Scholar
 W. Rudin, Principles of Mathematical Analysis, McGrawHill, 3rd edition, 1976. View at: MathSciNet
 R. Gilmore, M. Lefranc, and N. B. Tufillaro, The topology of chaos: Alice in stretch and squeezeland, Weinheim Chichester: WileyVCH John Wiley distributor, 2011. View at: Publisher Site
 M. L. Overton, Numerical computing with {IEEE} floating point arithmetic, Society for Industrial and Applied Mathematics, 2001. View at: MathSciNet
 IEEE Standards Committee and others, “7542008 IEEE standard for floatingpoint arithmetic,” IEEE Computer Society Std, 2008. View at: Publisher Site  Google Scholar
 E. G. Nepomuceno and E. M. A. M. Mendes, “On the analysis of pseudoorbits of continuous chaotic nonlinear systems simulated using discretization schemes in a digital computer,” Chaos, Solitons & Fractals, vol. 95, pp. 21–32, 2017. View at: Publisher Site  Google Scholar
 M. V. Correa, L. A. Aguirre, and E. M. A. M. Mendes, “Modeling chaotic dynamics with discrete nonlinear rational models,” International Journal of Bifurcation and Chaos, vol. 10, pp. 1019–1032, 2000. View at: Google Scholar
 E. G. Nepomuceno, R. H. Takahashi, G. F. Amaral, and L. A. Aguirre, “Nonlinear identification using prior knowledge of fixed points: a multiobjective approach,” International Journal of Bifurcation and Chaos, vol. 13, no. 5, pp. 1229–1246, 2003. View at: Publisher Site  Google Scholar  MathSciNet
 O. E. Rössler, “An equation for continuous chaos,” Physics Letters A, vol. 57, no. 5, pp. 397398, 1976. View at: Publisher Site  Google Scholar
 M. C. Mackey and L. Glass, “Oscillation and chaos in physiological control systems,” Science, vol. 197, no. 4300, pp. 287–289, 1977. View at: Publisher Site  Google Scholar
 J. Doyne Farmer, “Chaotic attractors of an infinitedimensional dynamical system,” Physica D: Nonlinear Phenomena, vol. 4, no. 3, pp. 366–393, 1982. View at: Publisher Site  Google Scholar
 X. Zeng, R. Eykholt, and R. A. Pielke, “Estimating the Lyapunovexponent spectrum from short time series of low precision,” Physical Review Letters, vol. 66, no. 25, pp. 3229–3232, 1991. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2018 Erivelton G. Nepomuceno et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.