Mathematical Problems in Engineering
Volume 2018 (2018), Article ID 6909151, 8 pages
https://doi.org/10.1155/2018/6909151

Research Article

On the Use of Interval Extensions to Estimate the Largest Lyapunov Exponent from Chaotic Data

Erivelton G. Nepomuceno,¹ Samir A. M. Martins,¹ Márcio J. Lacerda,¹ and Eduardo M. A. M. Mendes²

¹Control and Modelling Group (GCOM), Department of Electrical Engineering, Federal University of São João del-Rei, 36307-352 São João del-Rei, MG, Brazil
²Department of Electronic Engineering, School of Engineering, Federal University of Minas Gerais, 31270-901 Belo Horizonte, MG, Brazil

Correspondence should be addressed to Márcio J. Lacerda; lacerda@ufsj.edu.br

Received 28 November 2017; Revised 12 January 2018; Accepted 28 January 2018; Published 4 March 2018

Academic Editor: Gisele Mophou

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.

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 well-known benchmarks: sine map, Rössler Equations, and Mackey-Glass 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 well-known 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 Mackey-Glass 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 well-known 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 so-called 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 relation

where , 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 property

where 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 computer-based approach is used to estimate their values. We have exploited these features in this work to estimate the LLE.

Definition 3. Let represent a pseudo-orbit, which is defined by an initial condition, an interval extension of , some specific hardware and software, numerical precision standard, and discretization scheme. A pseudo-orbit is an approximation of an orbit and can be represented by or , such that

where is the error and .

A pseudo-orbit characterises an interval where the genuine orbit lies. Henceforth an interval related to each estimation of a pseudo-orbit is given by

From (7) and (8) it is clear that

Theorem 4 [16] establishes that at least one of two pseudo-orbits must have an error greater than or equal to the lower bound error.

Theorem 4. Let two pseudo-orbits 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.

Figure 1: LLE calculation of sine map: (a) original equations and (b) identified model. The -axis is the number of iterates and the -axis is the log of the lower bound error.

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 have

Figure 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.

Table 1: LLE for sine map investigated in [34]. The reference value is calculated using the original equations of the systems, as proposed in [15].

[15]
Our method

LLE
1.133
1.114

4.2. Example  2: Rössler Equations

Rössler’s Equations [35] are given by

This 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 term

in (15) as

for 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.

Table 2: LLE for Rössler’s Equations studied in [33].

[33] (data)
[33] (model)
Our method

LLE
1.242
1.566
1.530

Figure 2: LLE calculation for Rössler’s Equations. The -axis is the simulation time and the -axis is the log of the lower bound error. Using a polynomial NARMAX and interval extension, our method has estimated the LLE to be equal to 1.530, while the authors in [33] calculated it as 1.566.

4.3. Example  3: Mackey-Glass Equations

The Mackey-Glass is an interesting system used in many papers as an example of a chaotic and infinite dimensional system, since it is a time-delay system [36, 37]. The system equation is given by

where , , , 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 Mackey-Glass 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.

Table 3: LLE for Mackey-Glass Equation studied in [20].

Literature (data)
Literature (model)
Our method

LLE
0.007
0.006
0.006

Figure 3: LLE calculation for Mackey-Glass Equations. The -axis is the simulation time and the -axis is the log of the lower bound error. Using a polynomial NARMAX and interval extension, our method has estimated the LLE to be equal to 0.006, the same value found in the literature using the dynamical equations.

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.

Figure 4: The chaotic jerk circuit proposed by [21].

Figure 5: LLE calculation for chaotic jerk circuit. The -axis is the simulation time and the -axis is the log of the lower bound error. Using a polynomial NARMAX and interval extension, our method has estimated the LLE to be equal to 0.0787, which represents a good agreement with the value suggested by [21].

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.

Table 4: LLE for the variable of the chaotic jerk circuit investigated in [21]. The reference value was computed using circuit acquired data.

[21]
Our method

LLE
0.0735
0.0787

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.

Figure 6: LLE calculation for the tent map. The -axis is the number of iterates and the -axis is the log of the lower bound error. Using a polynomial NARMAX and interval extension, our method has estimated the LLE to be equal to 1.0461, while the correct value is 0.9993.

Figure 7: First return map for the tent map (, in black) and for the identified model (△, in red).

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.

Table 5: LLE for tent map.

Analytical
Our method

LLE
0.9993
1.0461

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 Lyapunov-exponent 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. 1-2, pp. 117–134, 1993.

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.

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.

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.

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.

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.

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.

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.

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.

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.

G. Wainrib and M. N. Galtier, “A local Echo State Property through the largest Lyapunov exponent,” Neural Networks, vol. 76, pp. 39–45, 2016.

A. J. Steyer and E. S. Van Vleck, “A step-size selection strategy for explicit Runge-Kutta methods based on Lyapunov exponent theory,” Journal of Computational and Applied Mathematics, vol. 292, pp. 703–719, 2016.

C. J. Gavilán-Moreno and G. Espinosa-Paredes, “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.

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. 50-51, pp. 116–138, 2015.

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.

E. G. Nepomuceno and S. A. M. Martins, “A lower bound error for free-run simulation of the polynomial NARMAX,” Systems Science & Control Engineering, vol. 4, no. 1, pp. 50–58, 2016.

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.

S. A. Billings, Nonlinear system identification: NARMAX methods in the time, frequency, and spatio-temporal domains, West Sussex: John Wiley & Sons, London, UK, 2013.

L. A. Aguirre and C. Letellier, “Modeling nonlinear dynamics and chaos: A review,” Mathematical Problems in Engineering, vol. 2009, Article ID 238960, 2009.

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.

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.

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.

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.

J. R. Ayala Solares, H.-L. Wei, and S. A. Billings, “A novel logistic-NARX model as a classifier for dynamic binary classification,” Neural Computing and Applications, pp. 1–15, 2017.

O. M. Mohamed Vall and R. M'hiri, “An approach to polynomial NARX/NARMAX systems identification in a closed-loop with variable structure control,” International Journal of Automation and Computing, vol. 5, no. 3, pp. 313–318, 2008.

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.

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.

W. Rudin, Principles of Mathematical Analysis, McGraw-Hill, 3rd edition, 1976.

R. Gilmore, M. Lefranc, and N. B. Tufillaro, The topology of chaos: Alice in stretch and squeezeland, Weinheim Chichester: Wiley-VCH John Wiley distributor, 2011.

M. L. Overton, Numerical computing with {IEEE} floating point arithmetic, Society for Industrial and Applied Mathematics, 2001.

IEEE Standards Committee and others, “754-2008 IEEE standard for floating-point arithmetic,” IEEE Computer Society Std, 2008.

E. G. Nepomuceno and E. M. A. M. Mendes, “On the analysis of pseudo-orbits of continuous chaotic nonlinear systems simulated using discretization schemes in a digital computer,” Chaos, Solitons & Fractals, vol. 95, pp. 21–32, 2017.

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.

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.

O. E. Rössler, “An equation for continuous chaos,” Physics Letters A, vol. 57, no. 5, pp. 397-398, 1976.

M. C. Mackey and L. Glass, “Oscillation and chaos in physiological control systems,” Science, vol. 197, no. 4300, pp. 287–289, 1977.

J. Doyne Farmer, “Chaotic attractors of an infinite-dimensional dynamical system,” Physica D: Nonlinear Phenomena, vol. 4, no. 3, pp. 366–393, 1982.

X. Zeng, R. Eykholt, and R. A. Pielke, “Estimating the Lyapunov-exponent spectrum from short time series of low precision,” Physical Review Letters, vol. 66, no. 25, pp. 3229–3232, 1991.