Next Article in Journal
Assessment of Polymer Concrete Sample Geometry Effect on Ultrasonic Wave Velocity and Spectral Characteristics
Next Article in Special Issue
Numerical Analysis of Stapes Prosthesis Constraining in the Case of Otosclerosis
Previous Article in Journal
Buckling Analysis of a Large Shelter with Composites
Previous Article in Special Issue
Dynamic Behavior of Aviation Polymer Composites at Various Weight Fractions of Physical Modifier
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Alternative Methods of the Largest Lyapunov Exponent Estimation with Applications to the Stability Analyses Based on the Dynamical Maps—Introduction to the Method

Division of Dynamics, Lodz University of Technology, Stefanowskiego 1/15, 90-001 Lodz, Poland
*
Author to whom correspondence should be addressed.
Materials 2021, 14(23), 7197; https://doi.org/10.3390/ma14237197
Submission received: 31 July 2021 / Revised: 26 October 2021 / Accepted: 12 November 2021 / Published: 25 November 2021

Abstract

:
Controlling stability of dynamical systems is one of the most important challenges in science and engineering. Hence, there appears to be continuous need to study and develop numerical algorithms of control methods. One of the most frequently applied invariants characterizing systems’ stability are Lyapunov exponents (LE). When information about the stability of a system is demanded, it can be determined based on the value of the largest Lyapunov exponent (LLE). Recently, we have shown that LLE can be estimated from the vector field properties by means of the most basic mathematical operations. The present article introduces new methods of LLE estimation for continuous systems and maps. We have shown that application of our approaches will introduce significant improvement of the efficiency. We have also proved that our approach is simpler and more efficient than commonly applied algorithms. Moreover, as our approach works in the case of dynamical maps, it also enables an easy application of this method in noncontinuous systems. We show comparisons of efficiencies of algorithms based our approach. In the last paragraph, we discuss a possibility of the estimation of LLE from maps and for noncontinuous systems and present results of our initial investigations.

1. Introduction

Lyapunov exponents are invariants characterizing numerous aspects of nonlinear systems’ dynamics from complexity, stability, loss of information about a system’s dynamical state, the type and structure of attractor—manifold to which the solution tends. The full spectrum of LE consists of a number of indicators equal to the analyzed system’s dimension. As Lyapunov exponents contain information about the limit of an exponential change of initial perturbation for infinite time range, procedures of LE estimation are very time intensive. Therefore, new methods which could increase the efficiency of LE estimation are still being developed. Even a comparatively minor improvement of a method means huge time savings. As far as investigations into the stability of dynamical systems are concerned, an application of the largest LE is warranted. Since analyzing the stability of dynamical systems is one of the most important challenges in science and engineering, we decided to attempt a development of the LLE estimation method. In the article, we try to demonstrate that our method is both simple and efficient. Additionally, we present the basics for its development, allowing further increase of the efficiency and potential for application for maps and systems with discontinuities.
While the first numerical study of the system behavior using LE dates back to 1964 and the work of Henon and Heiles [1], the LE estimation methods are still being improved [2,3,4,5,6,7,8,9], as LE’s are employed in many different areas of scientific and engineering research, including: mechanical systems [5,10,11,12,13], electrical systems [14,15,16] bioengineering [17], refs. [18,19] astronomy and astrophysics [20,21], materials [22], neuronal models investigations [23,24,25] optimal control [26], time series analyses [27], systems with different types of discontinuities [28,29,30,31], systems with parametric oscillations and fluctuating parameters [32], systems with time delay [12,33,34,35,36,37], systems with hidden attractors [38,39], chaotic fractional order derivative systems [40,41], hybrid-type systems [42,43] and synchronization phenomena analyses [44,45,46,47,48,49,50]. Since LE’s are applied in such a wide spectrum of scientific and engineering research, all the studies regarding properties of LE’s are highly justified.
Recently, we have studied different aspects of the nonlinear systems’ control with the use of different new nonlinear methods. We investigated the stability of continuous systems [51] and systems with discontinuities [9,28], control system’s optimization [52], synchronization phenomena of energy flow [48,53,54] and chaos-based control of energy flow [55,56,57]. We have also investigated efficiency of our novel method of Lyapunov spectrum estimation in [58] and showed that it allows for significant computation time savings.
As far as investigations into the stability of dynamical systems are concerned, application of the largest LE is warranted. Aproximately 60% of scientific research utilizes this simpler and faster indicator. In view of the above, we decided to extend studies of LLE’s properties and present the results of our new investigations.

2. The Method

Assume that a dynamical system is described by a set of differential equations in the form:
d x d t = f x ,   t
where x is a state vector, t is time and f is a vector field that (in general) depends on x and t . Consider a situation in which the state vector x is disturbed by an infinitesimal perturbation z (Figure 1). Evolution of the perturbation z can be determined by linearization of Equation (1):
d z d t = f x + z ,   t f x ,   t = f x x ,   t z = U x ,   t z
where U x ,   t = f x x ,   t is the Jacobi matrix obtained by differentiation of f with respect to x . If the Jacobi matrix was constant, then the evolution of the perturbation z in directions of subsequent eigenvectors would be specified by corresponding eigenvalues of that matrix. However, as long as the system (1) is nonlinear, the Jacobi matrix varies along the trajectory meaning that the evolution of the perturbation z cannot be directly predicted from properties of the Jacobi matrix. In such a case, Lyapunov exponents are applied to describe an average rate of expansion or contraction of a perturbation. Consequently, Lyapunov exponents can be treated as generalization of eigenvalues [59] of the Jacobi matrix. Moreover, according to [59] during an evolution of the system, eigenvectors connected with the largest eigenvalue spans the linear subspace which tends to align with the direction of the perturbation z t . As such, all the analyses of LLE can be focused on this direction. Following this eigenvalue idea, during numerical integration for each i - step of n integration steps, Equation (2) in the actual z i direction can be presented in the following scalar form:
d z i d t = λ i z i
where λ i tends to the largest eigenvalue of U x ,   t , and its average value is equal to LLE.
Equation (3) can be expressed in the form:
d z i z i = λ i d t
In the case where the perturbation is normalized before each integration step, z i = 1 :
d z i = λ i d t
For n numerical integration steps, from Equation (5), averaged perturbation:
i = 1 n d z i n = i = 1 n λ i d t n i = 1 n d z i n d t = i = 1 n λ i n = LLE
Finally,
1 n d z n t = LLE
From formula (7), one can see that LE can be treated as a dimensionless perturbation change averaged per time unit. It constitutes the basis for the first of the new methods (M1). As perturbation change d z is the scalar obtained from the differences between norms of z before and after each integration step, the method can be applied in the estimation of LLE from any given map. Additionally, it can be also applied for all the systems with any given discontinuities.
Moreover, following Equation (3), when the perturbation is normalized before each integration step:
d z d t = v t = λ
As the value of λ has to be averaged during evolution of the system to obtain LLE, from Equation (8), one can see that LLE equals the averaged speed v ¯ of perturbation changes:
v ¯ = LLE
The above constitutes a basis for the second of the new methods (M2).
Incidentally, both of the methods can be treated as identical. They differ only in the way the computed values are averaged during numerical integration. In the first one (Equation (7)), the values of d z are summed up and then averaged by division by time t of calculations. Additionally, regarding the averaged speed (Equation (9)), the actual speed is computed and summed up and the final value of LLE is obtained by division by number n of times of integration. As n · d t = t , both of the methods are equivalent.

3. Numerical Simulations

3.1. Methodology

All the programs for conducting numerical simulations have been written in C++ by means of the Code: Blocks environment. The Runge–Kutta method of the fourth order (RK4) has been used to solve ordinary differential equations. The integration step has been adjusted for each analyzed system separately, based on its own time scale.
We have studied perturbation change averaged per time unit and averaged speed v ¯ of perturbation change and compared them with three other methods. All of the considered algorithms of the LLE estimation require integration of the system (1) along with the Equation (2) in order to obtain the state vector x t and the perturbation z t in subsequent moments of time. Depending on the method, the vector z t was either normalized before each integration step or normalized only in the case of excessively high or low values of perturbation length z t .
All the programs for estimation of the LLE share the same code for integration of the systems (1), (2). The only difference between these programs is the method of the LLE calculation.

3.1.1. Method 1 (M1)

In the first one, the value given by the Equation (7) is calculated after each integration step. Value dz was obtained from the differences between norms of z before and after each integration step. Vector z was normalized before each integration step.

3.1.2. Method 2 (M2)

In the second case, the value given by Equation (9) is calculated from projection of the vector d z d t onto the direction of normalized vector z according to formula:
λ = d z d t · z z 2 = d z d t · z z · z  
As vector z was normalized before each integration step,
λ = d z d t · z

3.1.3. Method 3 (M3)

The third case involves application of the classical method [59] for vector z normalized in the case of excessively high or low values of perturbation length z t :
λ = 1 t   ln z t z 0

3.1.4. Method 4 (M4)

The fourth case is application of the classical method [59] for vector z normalized before each integration step:
λ = 1 t   ln z t

3.1.5. Method 5 (M5)

The last case is the application of our effective method presented in [58]. In this case, the value given by Equation (9) is calculated from the projection of the vector d z d t onto the direction of not normalized vector z according to:
λ = d z d t · z z · z
In simulation algorithms, conditions for termination of calculations have to be selected. It seems reasonable to finish the estimation procedure if the obtained value of the LLE stabilizes at some fixed value and does not display any relevant fluctuations. In order to measure stabilization of the LLE value, the authors propose to define a buffer of a fixed size. In this research, the buffer capacity equal to 100 was selected. After each calculation step, the current value of the LLE was saved to the buffer. When the buffer was full, the standard deviation of all the LLE values in the buffer was calculated. If the standard deviation related to actual average LLE was below a specified threshold, the value of the LLE was considered as stable and the calculations could be terminated. Failing that, the buffer was cleared and the procedure repeated. The value of the selected threshold corresponded to the desired accuracy of estimation. Lowering the threshold meant higher accuracy, but, consequently, a longer estimation time. Considering the standard deviation threshold, two methods, with relative and not relative deviation value, can be applied. They differ in accuracy of LLE estimation depending on the dynamical state of the system. In the regions of higher absolute LLE values for high accuracy, it proves advantageous to use nonrelative deviation; in the case of quasiperiodic regions, relative value will produce more accurate results. As insignificant differences in values in the periodic and chaotic regions are not of considerable importance, and conversely, detection of the exact bifurcation point is one of the most important considerations in nonlinear systems investigations, relative deviation was applied in our simulations.
As regards the threshold of excessively high or low values of perturbation’s vector z length, the normalization condition was associated with the product of the first two coordinates of vector z . It allowed for introducing a condition, which does not burden simulation procedures much.

3.2. Results of Numerical Simulations

In order to verify the presented methods of the LLE estimation, two typical nonlinear systems have been analyzed. What follows are the results obtained for Duffing and Van der Pol systems with external forcing. Since the details that follow are organized in the same manner, in order to avoid repeating the same description, specification of the graphs is provided only once below.
The first type of the graphs that follow provides the obtained values of the LLE along with computation time lengths for all the investigated methods. Ratios of the program execution times t 1 ,   ,   t 5 for all of the five methods represent the execution time of LLE estimation for the specified bifurcation parameter and method, respectively. In the article, we have associated uniform color schemes and types of curves with respective methods.
Subsequently, efficiency analysis is presented. Special efficiency indicators are introduced. Let T 1 ,   ,   T 5 be sums of t i values, presenting the time measured from the beginning of simulations to the moment a specified bifurcation parameter for each of the five methods has been reached. Let us use these values to introduce efficiency indicators:
η i = T i T 3  
Relations of η i , with respect to bifurcation parameter of the investigated algorithms M1, M2, M4, M5 as compared to the classical method M3 are presented on subsequent charts. The efficiency gain of the four investigated methods in comparison to the commonly applied method M3 is appreciable.
In the following charts, dependence of LLE on bifurcation parameter is presented along with focused analyses presenting the accuracy of LE estimation for three different dynamical states: periodic, quasiperiodic and chaotic.

3.3. Duffing Oscillator

The Duffing oscillator can be described by the following set of differential equations:
x ˙ 1 = x 2 x ˙ 1 = β x 2 α x 1 3 + q c o s ω t  
Based on Equation (2), the Jacobi matrix is necessary to observe evolution of a perturbation. For the Duffing oscillator, the Jacobi matrix is defined as follows:
U = 0 1 3 α x 1 2 β  
The plot of the LLE for different values of the parameter q and graphs depicting computation time ratios are presented in Figure 2.
It is evident that the longest times occur in chaotic regions, and in instances when the system is approaching bifurcation points. This is related to a longer time which is required to stabilize LLE in a chaotic regime in the first case. In the second case, the main reason was given above and is connected with computing relative or not relative standard deviation in the procedure concluding LLE computations. Since minor differences in values in the periodic and chaotic regions are not highly important, and, conversely, the detection of the exact bifurcation point is one of the most important issues in nonlinear systems investigations, relative deviation was applied in our simulations. Obviously, this increased the time needed to satisfy the required LLE value stability condition.
The efficiency analysis of four methods M1, M2, M4, M5 with respect to M3 is presented in Figure 3. From the method of construction of the η i indicators, one can deduce that these values for the specified bifurcation parameter q show the average efficiencies of computations of each method from the beginning of calculations until the parameter q is reached. Therefore, the values η i corresponding to the last values of bifurcation parameter present all the average efficiencies of each of the methods. As is evident from Figure 3, only method M4 has efficiency which is not superior to M3. This is to be expected, as M4 is based on M3 and utilizes normalization of perturbation vector in each integration step, while, in M3, the normalization is carried out only in the cases of excessively high or low values of perturbation’s vector z length. The final efficiency η 4 is equal to 0.997. Both of the new presented methods, M1 and M2, offer better efficiency than M3. Method M1, which has the potential to be applied in non-continuous systems, offers the final efficiency η 1 equal to 0.927. Therefore, on average, M1 saves about 7% of the computation time. The effect will be even more pronounced when applied for the maps. Method M2 has the final η 2 equal to 0.853, so it saves on average about 14% of the computation time. Finally, method M5 has the best average efficiency η 2 equal to 0.780, so M5 saves on average about 22% of the computation time. The results for M1, M2, and M5 will be marginally inferior for more complex systems, as shown in [58]. However, they will be invariably superior to M3.
Accuracy comparison of LLE estimation is presented in Figure 3, where the LLE dependence on bifurcation parameter q on a low scale is shown. It can be seen that there exists correspondence between the results of all five of the methods. Higher scale results for the three different dynamical states of the system are presented in the upper part of Figure 4. There is good agreement for M2…M5 methods and only minor differences exist for the M1 method in the periodic and quasiperiodic regions. As these are fourth order level differences, they do not disqualify the M1 method, especially that its efficiency will be considerably higher when applied in LLE estimation from maps.

3.4. The Van der Pol Oscillator

The Van der Pol oscillator can be described by the following set of differential equations:
x ˙ 1 = x 2 x ˙ 2 = μ x 1 2 + 1 x 2 + x 1 + q c o s ω t  
Jacobi matrix was used to simulate evolution of a perturbation according to Equation (2). For the Van der Pol oscillator, the Jacobi matrix is defined as follows:
U = 0 1 2 μ x 1 x 2 + 1 μ x 1 2 + 1  
The plot of the LLE for different values of the parameter μ , together with computation times, is presented in Figure 5, whereas graphs depicting computation time ratios are presented in Figure 6. For the same reasons as in the case of Duffing system, the longest computations appear in chaotic regions and when the system is approaching bifurcation points. These effects connected with the application of the relative deviation can be also observed in the regions with high absolute values of negative LLE. As the demanded accuracy in these regions decreases together with the values of LLE, one can see short computation times in these areas. What is important and also evident from the lower part of Figure 5, the influence of such variable accuracy on the estimated values of LLE is negligible–no significant noise caused by this effect can be observed on the LLE graph.
It can also be seen in Figure 5 that, even in the Van der Pol system, its divergence varies in time of oscillations—as its dumping is nonlinear, less time is needed to stabilize LLE in chaotic and qusiperiodic regions than in the case of the Duffing system. Maximal values of time for Van der Pol are approximately 120 [s], while, for the Duffing system, they are approximately 160 [s]. As divergence of the system is equal to the Lyapunov exponents’ sum, it would appear that the varying divergence could disrupt the stabilization process of LLE values. Apparently, not only does it disrupt the process, but it speeds it up.
The effects of variable divergence can be also observed in the values of LLE in periodic regions. In the case of the Duffing system, the values of LLE are constant and equal to half of the divergence (the second Lyapunov exponent is equal to LLE). In the case of Van der Pol, there exist no regions of constant LLE.
Efficiency analysis of the four methods M1, M2, M4, M5 with respect to M3 is presented in Figure 6. For the same reasons as in the case of the Duffing system, it is only method M4 that has no superior efficiency compared to M3. Both of the new presented methods, M1 and M2, have better efficiency than M3. Method M1, which has the potential to be applied in non-continuous systems, has the final efficiency η 1 equal to 0.928 (Duffing 0.927). Therefore, on average, M1 saves about 7% of the computation time. The effect will be even more pronounced when applied for maps. Method M2 has the final η 2 equal to 0.867 (Duffing 0.853), so it saves on average about 14% of the computation time. Finally, method M5 has the best average efficiency η 2 Equal to 0.837 (Duffing 0.780), which translates into an average of approximately 16% savings of the computation time. These results confirm conclusions for the Duffing system.
Accuracy comparison of LLE estimation is presented in Figure 7. In the bottom section, one can see LLE dependence on bifurcation parameter q on a low scale. Similarly to the Duffing system, there exists a good agreement between the results of all five methods.
Higher scale results for the three different dynamical states of the system are presented in the upper part of Figure 4. It can be appreciated that, unlike for the Duffing system, the results instead merge and cannot be accurately determined.

4. Largest Lyapunov Exponent (LLE) from Maps

As it was proved in [60], with the use of our method, perturbation behavior of a perturbation can be reconstructed based on the time series of the dynamical system, without reconstruction of the Jacobi matrix. It can be combined with the approach presented above and then applied for dynamical maps.
The first approach comes directly from application of the method M1. In this case, the value of the sum of perturbations was averaged while the trajectory x t crossed the hyperplane π —see Figure 8. A time series comparison with the method M1 is shown in Figure 9. It can be seen that estimation error is within the same range as for the method M1.
The second approach requires an extended analysis. In Figure 8, a trajectory x t of a dynamical system and the perturbed system trajectory y t can be seen. While these trajectories cross the hyperplane π , one obtains perturbation z and then next perturbation z1 from the next points of crossing trajectories through the hyperplane π . After projection of the difference of the vectors z 1 z on to the direction of perturbation z , one obtains a differential d z . It allows for substituting the lengths z and d z into Equation (4) to find λ value. Alternatively, d z can be calculated from the difference of the norms of vectors z 1 z . However, in this case, the estimation error is expected to be higher. During the evolution of the system, obtained values have to be averaged and then recalculated according to the error correction analysis presented below.

Error Correction Analysis

Between the trajectory crossing the hyperplane π, there were calculated i steps of numerical integration. During numerical calculation of LLE, in each integration step, values λ i are obtained, and then averaged in time in order to obtain LLE. Following reasoning that justified scalar notation of Equation (3), we can continue in the same vein in the case of maps. Then, the value of the proposed indicator for a map is:
λ m a p = z i z 0 z  
Consider
z i + 1 = z i + d z i   a n d   d z i = λ i z i d t  
λ m a p d t = λ 0 + λ 1 + + λ i 1 + λ i + δ  
where δ is LLE estimation error.
While the final value of LLE is an average of λ i :
λ m a p i d t = λ m a p T = LLE + δ i  
where T is the time from one to the next crossing the map. To simplify the analysis of the correction error, let us start with i = 5. Then,
δ = d t ( λ 0 λ 1 + λ 0 λ 2 + λ 0 λ 3 + λ 0 λ 4 + λ 1 λ 2 + λ 1 λ 3 + λ 1 λ 4 + λ 2 λ 3 + λ 2 λ 4 + λ 3 λ 4 ) + + d t 2 λ 0 λ 1 λ 2 + λ 0 λ 1 λ 3 + λ 0 λ 1 λ 4 + λ 0 λ 2 λ 3 + λ 0 λ 2 λ 4 + λ 0 λ 3 λ 4 + λ 1 λ 2 λ 3 + λ 1 λ 2 λ 4 + λ 1 λ 3 λ 4 + λ 2 λ 3 λ 4 + d t 3 λ 0 λ 1 λ 2 λ 3 + λ 0 λ 1 λ 2 λ 4 + λ 0 λ 1 λ 3 λ 4 + λ 0 λ 2 λ 3 λ 4 + λ 1 λ 2 λ 3 λ 4   + d t 4 λ 0 λ 1 λ 2 λ 3 λ 4  
Finally, n-th power of dt is connected with i n 1 combinations of products   of   n + 1   of   λ j   , where: j = 0…i−1. Obviously, λ j   are unknown while calculating λ m a p . In order to estimate the value of the correction error, we have assumed that λ j   equals the average value λ a v . Then:
δ = d t i 2 λ a v   2 + d t 2 i 3 λ a v   3 + d t 3 i 4 λ a v   4 + d t 4 i 5 λ a v   5 + + d t i 2 i i 1 λ a v   i 1 + d t i 1 i i λ a v   i i
As λ a v   = LLE and for nondimensional T = 1 i = 1 d t , we obtain the final correction error CE:
CE = j = 2 i d t j i j LLE j  
Finally, LLE can be estimated from the following dependence:
λ m a p T =   LLE   +   CE
The presented approach was applied to estimate LLE of the Duffing system (Equations (16) and (17)). Time series of LLE obtained from numerical simulations can be seen in Figure 10. As is evident, the estimated value of LLE = −0.0237. As the dumping coefficient b = −0.05 and the system remains within the range of the periodic regime LLE = b /2 = −0.025. Thus, the error of the estimated value is 0.0013. From Figure 11, it can be seen that the correction error is within the range of 0.0025. After correction of the obtained LLE value, finally LLE = −0.0262. It means that the error of the presented LLE estimation is about 5%. However, as the value of LLE is computed only while the trajectory intersects the map, the method is expected to be much faster than the continuous ones. In our next article, we will present an extended study of the presented method.

5. Conclusions

The present article introduces new methods of LLE estimation for continuous systems and maps.
We have proved that the sum of dimensionless perturbations, averaged per time unit of measuring the evolution of the system, constitutes the value of the LLE. We have shown that this approach works also in the case of dynamical maps. Additionally, we have proved that LLE can be also equated to the average speed of perturbation change.
The basic background of the methods was presented. The results were compared with other methods. Investigations were carried out for two typical nonlinear systems. We have shown a good agreement of the results obtained with the use of the new approaches with respect to the other methods.
In the case of continuous systems, we have also compared efficiencies of algorithms based on these methods. We have shown that the new presented methods have better efficiency than the commonly applied M3. We have shown that M1 can save about 7% of the computation time. Method M2 is faster and saves on average about 14% of the computation time. We have also shown that the fastest method, M5, saves on average about 16–22% of the computation time.
We have also discussed basic aspects of the application of the presented methods in estimation of LLE from maps and for noncontinuous systems and showed the initial results of our approach. An extended study of this section of the article will be presented in the next publication.

Author Contributions

Conceptualization, A.D. and M.B.; methodology, A.D. and A.S.; software, A.D. and T.S.; validation, A.D., S.Z., V.D. and M.B.; formal analysis, A.D. and M.B.; investigation, A.D.; resources, A.D. and A.S.; data curation, A.D., S.Z. and V.D.; writing—original draft preparation, A.D.; writing—review and editing, M.B., S.Z. and V.D.; visualization, A.D., S.Z. and V.D.; supervision, A.S.; project administration, A.S. and S.Z.; funding acquisition, A.S. and S.Z. All authors have read and agreed to the published version of the manuscript.

Funding

The authors acknowledge Lodz University of Technology for APC funding. This study has been supported by the National Science Centre, Poland, PRELUDIUM Programme (Project No. 2020/37/N/ST8/03448). This study has been supported by the National Science Centre, Poland under project No. 2017/27/B/ST8/01619.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgments

This article has been completed while the fifth author, Sandra Zarychta, was in the Interdisciplinary Doctoral School at the Lodz University of Technology, Poland.

Conflicts of Interest

The authors declare that they have no conflict of interest.

References

  1. Hénon, M.; Heiles, C. The applicability of the third integral of the motion: Some numerical results. Astron. J. 1964, 69, 73–79. [Google Scholar] [CrossRef] [Green Version]
  2. Liao, H. Novel gradient calculation method for the largest Lyapunov exponent of chaotic systems. Nonlinear Dyn. 2016, 85, 1377–1392. [Google Scholar] [CrossRef]
  3. Peixoto, M.L.; Nepomuceno, E.G.; Martins, S.A.; Lacerda, M.J. Computation of the largest positive Lyapunov exponent using rounding mode and recursive least square algorithm. Chaos Solitons Fractals 2018, 112, 36–43. [Google Scholar] [CrossRef]
  4. Zhou, S.; Wang, X. Identifying the linear region based on machine learning to calculate the largest Lyapunov exponent from chaotic time series. Chaos Interdiscip. J. Nonlinear Sci. 2018, 28, 123118. [Google Scholar] [CrossRef]
  5. Zhou, S.; Wang, X.; Wang, Z.; Zhang, C. A novel method based on the pseudo-orbits to calculate the largest Lyapunov exponent from chaotic equations. Chaos Interdiscip. J. Nonlinear Sci. 2019, 29, 033125. [Google Scholar] [CrossRef] [PubMed]
  6. Shuang, Z.; Yong, F.; Wen-Yuan, W.; Wei-Hua, W. A novel method based on the fuzzy C-means clustering to calculate the maximal Lyapunov exponent from small data. Acta Phys. Sin. 2016, 65, 020502. [Google Scholar]
  7. Liao, H.; Wu, W. The Reduced Space Shooting Method for Calculating the Peak Periodic Solutions of Nonlinear Systems. J. Comput. Nonlinear Dyn. 2018, 13, 061001. [Google Scholar] [CrossRef]
  8. Zhang, Y.; Chen, D.; Guo, D.; Liao, B.; Wang, Y. On exponential convergence of nonlinear gradient dynamics system with application to square root finding. Nonlinear Dyn. 2014, 79, 983–1003. [Google Scholar] [CrossRef]
  9. Balcerzak, M.; Dabrowski, A.; Blazejczyk–Okolewska, B.; Stefanski, A. Determining Lyapunov exponents of non-smooth systems: Perturbation vectors approach. Mech. Syst. Signal Process. 2020, 141, 106734. [Google Scholar] [CrossRef]
  10. Liao, H. Stability Analysis of Duffing Oscillator with Time Delayed and/or Fractional Derivatives. Mech. Based Des. Struct. Mach. 2015, 44, 283–305. [Google Scholar] [CrossRef]
  11. Śmiechowicz, W.; Loup, T.; Olejnik, P. Lyapunov Exponents of Early Stage Dynamics of Parametric Mutations of a Rigid Pendulum with Harmonic Excitation. Math. Comput. Appl. 2019, 24, 90. [Google Scholar] [CrossRef] [Green Version]
  12. Liao, H. Nonlinear Dynamics of Duffing Oscillator with Time Delayed Term. Comput. Modeling Eng. Sci. 2014, 103, 155–187. [Google Scholar]
  13. Stefanski, A.; Pikunov, D.; Balcerzak, M.; Dabrowski, A. Synchronized chaotic swinging of parametrically driven pendulums. Int. J. Mech. Sci. 2020, 173, 105454. [Google Scholar] [CrossRef]
  14. Kengne, J.; Njitacke Tabekoueng, Z.; Kamdoum Tamba, V.; Nguomkam Negou, A. Periodicity, chaos, and multiple attractors in a memristor-based Shinriki’s circuit. Chaos Interdiscip. J. Nonlinear Sci. 2015, 25, 103126. [Google Scholar] [CrossRef]
  15. Wadduwage, D.P.; Wu, C.Q.; Annakkage, U.D. Power system transient stability analysis via the concept of Lyapunov Exponents. Electr. Power Syst. Res. 2013, 104, 183–192. [Google Scholar] [CrossRef]
  16. Rajagopal, K.; Jafari, S.; Karthikeyan, A.; Srinivasan, A.; Ayele, B. Hyperchaotic Memcapacitor Oscillator with Infinite Equilibria and Coexisting Attractors. Circuits Syst. Signal Process. 2018, 37, 3702–3724. [Google Scholar] [CrossRef]
  17. Dingwell, J.; Cusumano, J.P. Nonlinear time series analysis of normal and pathological human walking. Chaos Interdiscip. J. Nonlinear Sci. 2000, 10, 848–863. [Google Scholar] [CrossRef]
  18. Eteme, A.S.; Tabi, C.B.; Ateba, J.F.B.; Ekobena, H.P.F.; Mohamadou, A.; Kofane, T.C.; Fouda, H.P.E. Neuronal firing and DNA dynamics in a neural network. J. Phys. Commun. 2018, 2, 125004. [Google Scholar] [CrossRef]
  19. Cao, J.; Lu, J. Adaptive synchronization of neural networks with or without time-varying delay. Chaos Interdiscip. J. Nonlinear Sci. 2006, 16, 013133. [Google Scholar] [CrossRef]
  20. Zhou, S.; Feng, Y.; Wu, W.-Y.; Li, Y.; Liu, J. Low-dimensional chaos and fractal properties of long-term sunspot activity. Res. Astron. Astrophys. 2013, 14, 104–112. [Google Scholar] [CrossRef]
  21. Shuang, Z.; Yong, F.; Wen-Yuan, W. Chaos and fractal properties of solar activity phenomena at the high and low latitudes. Acta Phys. Sin. 2015, 64, 249601. [Google Scholar]
  22. Iwaniec, J.; Uhl, T.; Staszewski, W.J.; Klepka, A. Detection of changes in cracked aluminium plate determinism by recurrence analysis. Nonlinear Dyn. 2012, 70, 125–140. [Google Scholar] [CrossRef]
  23. Panahi, S.; Jafari, S.; Khalaf, A.J.M.; Rajagopal, K.; Pham, V.; Alsaadi, F.E. Complete dynamical analysis of a neuron under magnetic flow effect. Chin. J. Phys. 2018, 56, 2254–2264. [Google Scholar] [CrossRef]
  24. Etémé, A.S.; Tabi, C.B.; Mohamadou, A.; Kofané, T.C. Long-range memory effects in a magnetized Hindmarsh-Rose neural network. Commun. Nonlinear Sci. Numer. Simul. 2020, 84, 105208. [Google Scholar] [CrossRef]
  25. Etémé, A.S.; Tabi, C.B.; Mohamadou, A. Synchronized nonlinear patterns in electrically coupled Hindmarsh–Rose neural networks with long-range diffusive interactions. Chaos Solitons Fractals 2017, 104, 813–826. [Google Scholar] [CrossRef]
  26. Ornelas-Tellez, F.; Sanchez, E.N.; Loukianov, A.G.; Navarro-Lopez, E.M. Speed-gradient inverse optimal control for discrete-time nonlinear systems. In Proceedings of the 50th IEEE Conference on Decision and Control and European Control Conference, Orlando, FL, USA, 12–15 December 2011; IEEE: New York, NY, USA, 2011; pp. 290–295. [Google Scholar] [CrossRef]
  27. Hegger, R.; Kantz, H.; Schreiber, T. Practical implementation of nonlinear time series methods: The TISEAN package. Chaos Interdiscip. J. Nonlinear Sci. 1999, 9, 413–435. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  28. Balcerzak, M.; Sagan, T.; Dabrowski, A.; Stefanski, A. Fast and simple Lyapunov Exponents estimation in discontinuous systems. Eur. Phys. J. Spec. Top. 2020, 229, 2167–2181. [Google Scholar] [CrossRef]
  29. Pikunov, D.; Stefanski, A. Numerical analysis of the friction-induced oscillator of Duffing’s type with modified LuGre friction model. J. Sound Vib. 2019, 440, 23–33. [Google Scholar] [CrossRef]
  30. Balcerzak, M.; Dabrowski, A.; Stefański, A.; Wojewoda, J. Spectrum of Lyapunov exponents in non-smooth systems evaluated using orthogonal perturbation vectors. In Proceedings of the MATEC Web of Conferences, International Conference on Engineering Vibration, Sofia, Bulgaria, 4–7 September 2017; EDP Sciences: Les Ulis, France, 2018; Volume 148, p. 10003. [Google Scholar]
  31. Fuhg, J.N.; Fau, A. Surrogate model approach for investigating the stability of a friction-induced oscillator of Duffing’s type. Nonlinear Dyn. 2019, 98, 1709–1729. [Google Scholar] [CrossRef] [Green Version]
  32. Prakash, P.; Rajagopal, K.; Singh, J.P.; Roy, B.K. Megastability, Multistability in a Periodically Forced Conservative and Dissipative System with Signum Nonlinearity. Int. J. Bifurc. Chaos 2018, 28, 1830030. [Google Scholar] [CrossRef]
  33. Rajagopal, K.; Duraisamy, P.; Weldegiorgis, R.; Karthikeyan, A. Multistability in Horizontal Platform System with and without Time Delays. Shock. Vib. 2018, 2018, 1092812. [Google Scholar] [CrossRef] [Green Version]
  34. Rajagopal, K.; Pham, V.-T.; Tahir, F.R.; Akgül, A.; Abdolmohammadi, H.R.; Jafari, S. A chaotic jerk system with non-hyperbolic equilibrium: Dynamics, effect of time delay and circuit realisation. Pramana 2018, 90, 52. [Google Scholar] [CrossRef]
  35. Rajagopal, K.; Jafari, S.; Laarem, G. Time-delayed chameleon: Analysis, synchronization and FPGA implementation. Pramana J. Phys. 2017, 89, 92. [Google Scholar] [CrossRef]
  36. Rajagopal, K.; Jafari, S.; Akgul, A.; Karthikeyan, A.; Çiçek, S.; Shekofteh, Y. A Simple Snap Oscillator with Coexisting Attractors, Its Time-Delayed Form, Physical Realization, and Communication Designs. Zeitschrift für Naturforschung A 2018, 73, 385–398. [Google Scholar] [CrossRef]
  37. Stefański, A.; Kapitaniak, T.; Dąbrowski, A. The largest Lyapunov exponent of dynamical systems with time delay. In Proceedings of the IUTAM Symposium on Chaotic Dynamics and Control of Systems and Processes in Mechanics, Rome, Italy, 8–13 June 2003; Springer: Dordrecht, The Netherlands, 2005; pp. 493–500. [Google Scholar]
  38. Ren, S.; Panahi, S.; Rajagopal, K.; Akgul, A.; Pham, V.-T.; Jafari, S. A New Chaotic Flow with Hidden Attractor: The First Hyperjerk System with No Equilibrium. Zeitschrift für Naturforschung A 2018, 73, 239–249. [Google Scholar] [CrossRef]
  39. Rajagopal, K.; Jafari, S.; Akgul, A.; Karthikeyan, A. Modified jerk system with self-exciting and hidden flows and the effect of time delays on existence of multi-stability. Nonlinear Dyn. 2018, 93, 1087–1108. [Google Scholar] [CrossRef]
  40. Liao, H. Optimization analysis of Duffing oscillator with fractional derivatives. Nonlinear Dyn. 2014, 79, 1311–1328. [Google Scholar] [CrossRef]
  41. Huang, S.; Wang, B. Stability and stabilization of a class of fractional-order nonlinear systems for Stability and stabilization of a class of fractional-order nonlinear systems for 0 < α < 2. Nonlinear Dyn. 2017, 88, 973–984. [Google Scholar] [CrossRef]
  42. Gomes, C.; Karalis, P.; Navarro-López, E.M.; Vangheluwe, H. Approximated Stability Analysis of Bi-modal Hybrid Co-simulation Scenarios. In Software Engineering and Formal Methods; Lecture Notes in Computer Science; Cerone, A., Roveri, M., Eds.; Springer: London, UK, 2018; Volume 10729, pp. 345–360. [Google Scholar]
  43. Navarro-López, E.M.; Laila, D.S. Group and Total Dissipativity and Stability of Multi-Equilibria Hybrid Automata. IEEE Trans. Autom. 2013, 58, 3196–3202. [Google Scholar] [CrossRef]
  44. Pecora, L.M.; Carroll, T.; Johnson, G.A.; Mar, D.J.; Heagy, J.F. Fundamentals of synchronization in chaotic systems, concepts, and applications. Chaos Interdiscip. J. Nonlinear Sci. 1997, 7, 520–543. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Rajagopal, K.; Naseradinmousavi, P.; Khalaf, A.J.M.; Jafari, S.; Karthikeyan, A.; Çiçek, S. A novel parametrically controlled multi-scroll chaotic attractor along with electronic circuit design. Eur. Phys. J. Plus 2018, 133, 354. [Google Scholar] [CrossRef]
  46. Juan, S.; Xiao-Xia, L.; Jin-Hao, Z.; Yu-Zhuo, S.; Yan-Yu, L. Synchronizability and eigenvalues of multilayer star networks through unidirectionally coupling. Acta Phys. Sin. 2017, 66, 188901. [Google Scholar]
  47. La Guardia, G.G.; Miranda, P.J. Lyapunov exponent for Lipschitz maps. Nonlinear Dyn. 2018, 92, 1217–1224. [Google Scholar] [CrossRef] [Green Version]
  48. Balcerzak, M.; Chudzik, A.; Stefanski, A. Properties of generalized synchronization in smooth and non-smooth identical oscillators. Eur. Phys. J. Spec. Top. 2020, 229, 2151–2165. [Google Scholar] [CrossRef]
  49. Dabrowski, A.; Glogowski, M.; Kubiak, P. Improving the efficiency of four-stroke engine with use of the pneumatic energy accumulator-simulations and examination. Int. J. Automot. Technol. 2016, 17, 581–590. [Google Scholar] [CrossRef]
  50. Etémé, A.S.; Tabi, C.B.; Mohamadou, A. Firing and synchronization modes in neural network under magnetic stimulation. Commun. Nonlinear Sci. Numer. Simul. 2019, 72, 432–440. [Google Scholar] [CrossRef]
  51. Dabrowski, A.; Balcerzak, M.; Pikunov, D. Applications of the New Method of the Lyapunov Exponents Estimation. Proceedings of 23rd International Conference on Engineering Mechanics, Svratka, Czech Republic, 15–18 May 2017; pp. 266–269. [Google Scholar]
  52. Balcerzak, M.; Dąbrowski, A.; Pikunov, D. Tuning the control system of a nonlinear inverted pendulum by means of the new method of Lyapunov exponents estimation. AIP Conf. Proc. 2018, 1922, 100016. [Google Scholar] [CrossRef]
  53. Dąbrowski, A.; Jach, A.; Kapitaniak, T. Application of artificial neural networks in parametrical investigations of the energy flow and synchronization. J. Theor. Appl. Mech. 2010, 48, 871–896. [Google Scholar]
  54. Dabrowski, A. Energy–vector method in mechanical oscillations. Chaos Solitons Fractals 2009, 39, 1684–1697. [Google Scholar] [CrossRef]
  55. Dąbrowski, A.; Kapitaniak, T. Using chaos to reduce oscillations: Experimental results. Chaos Solitons Fractals 2009, 39, 1677–1683. [Google Scholar] [CrossRef]
  56. Dabrowski, A. The construction of the energy space. Chaos Solitons Fractals 2005, 26, 1277–1292. [Google Scholar] [CrossRef]
  57. Dabrowski, A. New design of the impact damper. Mech. Mech. Eng. 2000, 4, 191–196. [Google Scholar]
  58. Dabrowski, A.; Balcerzak, M.; Pikunov, D.; Stefanski, A. Improving efficiency of the largest Lyapunov exponent’s estimation by its determination from the vector field properties. Nonlinear Dyn. 2020, 102, 1869–1880. [Google Scholar] [CrossRef]
  59. Parker, T.S.; Chua, L. Practical Numerical Algorithms for Chaotic Systems; Springer: New York, NY, USA, 1989. [Google Scholar]
  60. Dabrowski, A. The largest transversal Lyapunov exponent and master stability function from the perturbation vector and its derivative dot product (TLEVDP). Nonlinear Dyn. 2012, 69, 1225–1235. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Graphic illustration of the method.
Figure 1. Graphic illustration of the method.
Materials 14 07197 g001
Figure 2. Diagram of the largest Lyapunov exponent of the Duffing system and computation times t i [s] graphs. α = 1, β = 0.05, ω = 0.47.
Figure 2. Diagram of the largest Lyapunov exponent of the Duffing system and computation times t i [s] graphs. α = 1, β = 0.05, ω = 0.47.
Materials 14 07197 g002
Figure 3. Diagram of efficiencies η 1 ,   η 2 ,   η 4 ,   η 5 , of LLE computations of the Duffing system. α = 1, β = 0.05, ω = 0.47.
Figure 3. Diagram of efficiencies η 1 ,   η 2 ,   η 4 ,   η 5 , of LLE computations of the Duffing system. α = 1, β = 0.05, ω = 0.47.
Materials 14 07197 g003
Figure 4. Diagram of accuracy of LLE computations, of the Duffing system. α = 1, β = 0.05, ω = 0.47.
Figure 4. Diagram of accuracy of LLE computations, of the Duffing system. α = 1, β = 0.05, ω = 0.47.
Materials 14 07197 g004
Figure 5. Diagram of the largest Lyapunov exponent of the Van der Pol system and computation times t i [s] graphs. q = 12.95, ω = 4.64.
Figure 5. Diagram of the largest Lyapunov exponent of the Van der Pol system and computation times t i [s] graphs. q = 12.95, ω = 4.64.
Materials 14 07197 g005
Figure 6. Diagram of efficiencies η 1 ,   η 2 ,   η 4 ,   η 5 of LLE computations of the Van der Pol system. q = 12.95, ω = 4.64.
Figure 6. Diagram of efficiencies η 1 ,   η 2 ,   η 4 ,   η 5 of LLE computations of the Van der Pol system. q = 12.95, ω = 4.64.
Materials 14 07197 g006
Figure 7. Diagram of the accuracy of LLE computations of the Van der Pol system. q = 12.95, ω = 4.64.
Figure 7. Diagram of the accuracy of LLE computations of the Van der Pol system. q = 12.95, ω = 4.64.
Materials 14 07197 g007
Figure 8. Graphical illustration of the method for maps.
Figure 8. Graphical illustration of the method for maps.
Materials 14 07197 g008
Figure 9. Time series comparison.
Figure 9. Time series comparison.
Materials 14 07197 g009
Figure 10. Time series of the largest Lyapunov exponent of the Duffing system. α = 1, β = 0.05, ω = 0.47.
Figure 10. Time series of the largest Lyapunov exponent of the Duffing system. α = 1, β = 0.05, ω = 0.47.
Materials 14 07197 g010
Figure 11. Correction error dependent on number of components (Equation (26)).
Figure 11. Correction error dependent on number of components (Equation (26)).
Materials 14 07197 g011
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Dabrowski, A.; Sagan, T.; Denysenko, V.; Balcerzak, M.; Zarychta, S.; Stefanski, A. Alternative Methods of the Largest Lyapunov Exponent Estimation with Applications to the Stability Analyses Based on the Dynamical Maps—Introduction to the Method. Materials 2021, 14, 7197. https://doi.org/10.3390/ma14237197

AMA Style

Dabrowski A, Sagan T, Denysenko V, Balcerzak M, Zarychta S, Stefanski A. Alternative Methods of the Largest Lyapunov Exponent Estimation with Applications to the Stability Analyses Based on the Dynamical Maps—Introduction to the Method. Materials. 2021; 14(23):7197. https://doi.org/10.3390/ma14237197

Chicago/Turabian Style

Dabrowski, Artur, Tomasz Sagan, Volodymyr Denysenko, Marek Balcerzak, Sandra Zarychta, and Andrzej Stefanski. 2021. "Alternative Methods of the Largest Lyapunov Exponent Estimation with Applications to the Stability Analyses Based on the Dynamical Maps—Introduction to the Method" Materials 14, no. 23: 7197. https://doi.org/10.3390/ma14237197

APA Style

Dabrowski, A., Sagan, T., Denysenko, V., Balcerzak, M., Zarychta, S., & Stefanski, A. (2021). Alternative Methods of the Largest Lyapunov Exponent Estimation with Applications to the Stability Analyses Based on the Dynamical Maps—Introduction to the Method. Materials, 14(23), 7197. https://doi.org/10.3390/ma14237197

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop