Improving efficiency of the largest Lyapunov exponent’s estimation by its determination from the vector field properties

Controlling dynamics of nonlinear systems is one of the most important issues in science and engineering. Thus, there is continuous need to study and develop numerical algorithms of control methods. Among the most frequently applied invariants characterizing different aspects of a systems’ dynamics are Lyapunov exponents, fast Lyapunov index, angles of small deviations, fractal dimension or entropy. There exist many different methods of estimation of these indicators. In this paper, modification of our novel method is presented. We have shown that LLE can be estimated from the vector field properties by means of the most basic mathematical operations. Results of efficiency measurements for typical mechanical, electrical and random systems were discussed. We have proved that discussed modification introduced to our method makes the LLE estimation 17–53% faster than using classical algorithms. In addition, unlike the results presented in our previous publication, an improvement in performance was achieved for each of the analyzed cases. As such, the new approach lends to prospective application of LLE not only in dynamical systems' stability investigations, but also in real-time control of systems since the basic calculations and fast, effective method of LLE estimation can be applied even in simple microcontrollers. Our approach could be also applied in investigations of vector field properties, global stability or basins of attraction analyses, allowing for huge time savings.


Introduction
Dynamics of nonlinear systems can be investigated with emphasis on different aspects, features and systems' properties. In stability analyses, there applied are Lyapunov exponents (LE) [1], angles of small deviations [2], fast Lyapunov index [3], smaller alignment index [4]. When studying topological properties with their applications, the most applicable invariant is fractal dimension [5]. When analyzing the amount of information in different data types, entropy is the superior invariant to be applied [6]. Apparently, there exist correlations between all these invariants. For instance, the Lyapunov exponents of bounded trajectory and the fractal dimension of an attractor are invariant under diffeomorphism of the phase space [7], and interrelated by the Kaplan-Yorke formula [8]. Conversely, according to Pesin's theorem [9], the sum of all positive Lyapunov exponents gives an A. Dabrowski  Moreover, based on the sum of the Lyapunov exponents, one can also have an estimate divergence of the system [10]. However, as far as investigations into the stability of dynamical systems are concerned, application of the largest LE is warranted. Approximately 60% of scientific research utilizes this simpler and faster stability indicator. In view of the above, we decided to extend studies of LLE's properties and present the results of our new investigations.
LE are employed in the whole spectrum of different areas of scientific and engineering research. Investigations of different types of systems, with different types of nonlinear phenomenon apply LE to study many aspects of these systems' dynamics. . Thus, there can be observed a need for continuous development of estimation methods of these indicators [47][48][49][50][51][52][53].
Recently, we have studied different aspects of the nonlinear systems control, introducing different new nonlinear methods. We have investigated stability of continuous systems [54], control system's optimization [55], synchronization phenomena [56,57], chaos-based energy flow control [58][59][60]. We have also investigated the efficiency of our novel method in [61 ] and demonstrated that it allows for huge computation time savings. As we mentioned in the abstract, we develop our method with hope of possible application in realtime control, even in simple microcontrollers. This is why we are focused on the methods which give current value of LLE in each integration step. Moreover, they have to consist the simplest and consequently fastest calculations. They have also may not contain any unnecessary calculations, like, for instance, Gram-Schmidt orthogonalization, which is no longer demanded, when instead of the whole spectrum, only the Largest Lyapunov exponent is estimated.

The method
The method presented in this paper enables simple and efficient estimation of the LLE. The LLE is an indicator which shows how fast an infinitesimal disturbance in a dynamical system expands or shrinks. If k is the LLE, then almost any perturbation in the system under consideration grows, on average, as e kt [62]. Consequently, the LLE describes stability of a solution: a positive value of the LLE means amplification of a disturbance, a negative value-suppression. If a chaotic orbit exists in a dynamical system, then the value of the LLE is positive. Therefore, the LLE is widely applied for stability analysis.
Assume that a dynamical system is described by the set of differential equations in the form: 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. Evolution of the perturbation z can be determined by linearization of Eq. (1): where U x; t ð Þ ¼ of ox 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 [62] of the Jacobi matrix.
The classical method of the LLE estimation originates from one of its basic properties. Due to the fact that almost any perturbation grows on average as e kt , the following approximate expression can be introduced: where z t ð Þ ¼ z t ð Þ j j is the length of the perturbation vector and t 0 is the initial time. Assume that the initial perturbation has been normalized, i.e., z t 0 ð Þ ¼ z t 0 ð Þ j j¼ 1. Then, a simple transformation of Eq. (3) produces: where Dt ¼ t À t 0 . The expression (4) is the foundation of the classical method of the LLE estimation.
However, implementation of the method has to take into consideration numerical issues connected with too high or too low values of z t ð Þ. These problems are solved by normalization of the vector z t ð Þ. Practically, the systems (1), (2) are integrated for a selected period of time Dt and calculate the value (4). Then, the vector z t ð Þ is normalized and integration of Eqs. (1) and (2) is continued with the new, normalized initial perturbation. Subsequently, the value (4) is computed again. Such procedure is repeated in a loop and the final estimation of the LLE is the average of all the values of Formula (4) obtained throughout the calculations. A more detailed presentation of this classical algorithm can be found in [62]. Contrary to the classical method, our novel method does not focus on logarithm of perturbation's norm. Instead, a geometrical reasoning is proposed. Consider two trajectories of the system (1): the reference trajectory (undisturbed) x t ð Þ and the perturbed trajec- Both the trajectories, the perturbation and its derivative are presented in Fig. 1.
The vector of perturbation's time derivative dz=dt can be separated into two components: the one parallel to the perturbation z, which is named dz Ã =dt, and the one orthogonal to the perturbation z. The length of the parallel component can be determined using the dot product: It can be easily demonstrated that dz Ã =dt is the only component of the derivative dz=dt that influences the length z. Indeed, derivative of z with respect to time results in: which is identical to the value in formula (5). Based on that observation, it can be concluded that evolution of perturbation's length is a one-dimensional problem which involves just two parallel vectors: z and dz Ã =dt.
Consequently, a novel method of the LLE estimation can be proposed. Differentiation of the approximate expression (3) yields: Note that k Ã has been substituted for k, the reason of which will be given in the subsequent paragraph. Dividing the expression (7) by z t ð Þ and substitution from Eq. (6) results in: The expression (8) constitutes a novel, alternative method of the LLE estimation. Clearly, calculating the value of (8) once in order to obtain the exact value of the LLE is not sufficient. Instead, the expression (8) needs to be evaluated multiple times during the simulation, preferably at each integration step. The average of all the calculated values k Ã converges to the LLE k.
To prove this fact, consider time derivative of the perturbation's length logarithm: Using the expression (9), the equivalence of the novel method to the classical one (Eq. 4) can be easily demonstrated: where the last expression in Eq. (10) is simply an average of all the values k Ã (Eq. 8) from time t ¼ 0 to t ¼ s. Therefore, the value obtained from Formula (4) is identical to the average of all the values (8). This proves that the novel method will produce the same results as the classical method.
One can see, that, unlike in commonly applied methods, formula (8) can be evaluated before the integration step. Moreover, this value consists of the elements, which are computed inside Runge-Kutta integration procedure, thus can be evaluated inside it. Additionally, possibility of LE calculation before an integration step improves accuracy of the method. Application of this novelty provides improvement of the efficiency of our method presented in [49].
Additionally, Eq. (8) can be presented in the different form: Therefore, as long as for investigated system, behavior of the perturbance can be analyzed using its linearized form, the value of k Ã can be treated as eigenvalue of U x; t ð Þ in direction of eigenvector, to which perturbation vector z aligns. As such, k Ã do not depend on the perturbation vector length.

Methodology
All the programs for conducting numerical simulations have been written in C ? ? . 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 individually, taking into account its own timescale.
Both the considered methods of the LLE estimation, the novel method and the classical method, require integration of the system (1) along with Eq. (2) in order to obtain the state vector x t ð Þ and the perturbation z t ð Þ in subsequent moments of time. In order to present robustness of the method with respect to different initial conditions, random values of initial conditions were applied. Additionally to preclude problems connected with excessively high or low values of perturbation length z t ð Þ, the vector z t ð Þ was normalized after each integration step. Moreover, Formula (8) used in the novel method requires the perturbation derivative dz=dt. However, this derivative is known before the integration step as a property of the vector field in the analyzed point. Nevertheless, it has to be calculated inside the RK4 method in order to integrate Eq. (2). Therefore, instead of calculating dz=dt again, it is saved in the integration procedure and used later for estimation of the LLE. This approach increases efficiency of calculations and leads, unlike in [61], to the improvement of efficiency for all the investigated cases. Programs for estimation of the LLE by means of the classical method (Eq. 4) and the novel one (Eq. 8) share the same code for integration of the systems (1), (2) and normalization of the vector z t ð Þ. The only difference between these two programs is the method of the LLE calculation. In the first one, the value given by Eq. (4) is calculated after each integration step. In the latter, the expression (8)  Lastly, 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 selected 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 was below a specified threshold, the value of the LLE was considered 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.

Results of numerical simulations
In order to verify our novel method of the largest Lyapunov exponent estimation, different systems have been analyzed. What follows are the results obtained for Duffing and Van der Pol systems with external forcing, Lorenz system and two coupled Duffing systems. Moreover, abstract systems of the order 2, 3,., 9 with random Jacobi matrices have also been investigated in order to reduce complexity of the simulation program and to focus on execution time of the LLE estimation procedures.
All the following graphs provide the obtained values of the LLE for all the systems under investigation along with computation time comparison between the classical method and the novel method. 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.
Each figure presents values of the largest Lyapunov exponents and ratios of execution time for both programs, depending on a selected system's parameter. Let t 1 be time of calculations by means of the novel method (Eq. 8) and t 2 -computations time when the classical method (Eq. 4) is used. Then, the relative computation time (RCT) of the novel method is defined as the ratio t 1 =t 2 for each parameter. It is represented by a black point in the graph. Moreover, let T 1 be the total time of calculations for all parameter values when the novel method is applied, whereas T 2 corresponds to the total time for the classical method. The average relative computation time (ARCT) is defined as the ratio T 1 =T 2 and it is marked with gray, horizontal line. The LLE values for the whole parameter range are marked with a black, solid line so as to illustrate the correlations with the computation time.

Duffing oscillator
The Duffing oscillator can be described by the following set of differential equations: According to Eq. (2), Jacobi matrix is necessary to observe evolution of a perturbation. For the Duffing oscillator, the Jacobi matrix is defined as follows: The plot of the LLE for different values of the parameter q and graphs depicting computation time ratios are presented in Fig. 2. Apparently, there is no correlation between the relative efficiency and the value of the LLE. The value of RCT t 1 =t 2 varies from 0.68 to 0.83. The ARCT T 1 =T 2 equals 0.78. It means that, for the Duffing system, estimation of the LLE by means of the novel method saves more than 20% of the computation time.
The values of relative efficiency were expressed with only four values: 0.69, 075, 0.77, 0.83. This is due to the fact that results of time measurement were rounded to seconds. Simulation for each parameter lasted less than 3 min, so the relative efficiency t 1 =t 2 is strongly discretized.

The Van der Pol oscillator
The Van der Pol oscillator can be described by the following set of differential equations: Jacobi matrix was used to simulate evolution of a perturbation according to Eq. (2). For the Van der Pol oscillator, the Jacobi matrix is defined as follows: The LLE plot and time efficiency diagrams for the Van der Pol system (13) are shown in Fig. 3 at the bottom. It can be noticed that the RCT t 1 =t 2 varies from the level of 0.78 to 0.85, but no evident correlation between system's dynamics changes and calculations efficiency can be discerned. The results of the performed tests demonstrate that the estimation of the LLE with the use of the new method is 15-20% faster than by means of the classical algorithm. The ARCT is approximately equal to 0.83.
The reason for which the RCT attains higher values for larger q remains unknown. The authors did their utmost to avoid influence of other applications and processes on the calculation time. Therefore, it is assumed that the change of RTC stems in some way from the dynamics of the system (14).

Coupled Duffing oscillators
The system of two coupled Duffing oscillators can be described by the following system of differential equations: For such a system, the Jacobi matrix (2) is as follows: Fig. 2 Diagram of the largest Lyapunov exponent of the Duffing system and relative computation time graphs. a = 1, b = 0.05,g = 0.33 Fig. 3 Diagram of the largest Lyapunov exponent of the Van der Pol system and relative computation time graphs. l = 8 where u x ¼ À3ax 2 À c and u z ¼ À3az 2 À c. The LLE plot and time efficiency diagrams for the system (15) are shown in Fig. 4 at the bottom. Contrary to previous examples, a clear relationship between system dynamics and the RCT t 1 =t 2 can be observed. In ranges of the parameter q for which a chaotic orbit exists, a large dispersion of RCT values t 1 =t 2 can be observed. In such ranges, RCT varies from the level of 0.2 to above 1. However, the ARC is equal to 0.81, which means that, in total, almost 20% of the computation time can be saved when the novel method is applied.
Reasons for the presented large dispersion of RCT, when the system (15) is unstable, have not been identified.

The Lorenz system
The Lorenz system can be described by the following system of differential equations: For such a system, the Jacobi matrix (2) is as follows: The LLE plot and time efficiency diagrams for the system (17) are shown in Fig. 5 at the bottom. In this case, no correlation between the LLE value and the RCT t 1 =t 2 can be observed. The RCT varies from the level of 0.63 up to 0.77. The mean value ARCT T 1 =T 2 is 0.65, which means that the novel method is on average 35% faster than the traditional one for the Lorenz system (17). Similarly to Duffing system, the RCT values are discretized due to the fact that time measurements were conducted in seconds.

Random Jacobi matrices
The last test case is the most simplified one. In order to focus on efficiency of both methods of the LLE estimation, only the perturbation z was integrated according to Eq. (2). Then, the value of the LLE was estimated by means of the classical method (Eq. 4) and the novel method (Eq. 8). The values of all the terms in the Jacobi matrix U were selected randomly: Fig. 4 Diagram of the largest Lyapunov exponent of the coupled Duffing oscillators system (15) and corresponding relative computation time graphs. a = 1, b = 0.05,g = 0.33, The dependence of the estimated LLE values and the corresponding RCT values t 1 =t 2 for random Jacobi matrices of different orders is shown in Fig. 6 at the bottom. It can be noticed that in such a simplified test, the RCT is practically independent of the LLE value. However, the RCT increases with the order of the Jacobi matrix. The ARCT T 1 =T 2 ranges from 47% for the order 2 up to 79% for the order 9.
Comparing these results, it can be seen that for all the analyzed cases algorithm of LLE estimation with use of the proposed method is at least 20% faster, than commonly applied one. It can be also seen that generally efficiency improvement depends on the complexity of the investigated system. But there also exist other properties of the systems that influence on this improvement. It can be seen while comparing Duffing system of the order 2 and Lorenz system of the order 3. It can be seen that for more complex Lorenz system of the higher order, efficiency was 13% better than for Duffing system. Additionally, for Duffing system of the order 2 and coupled Duffing systems of the order 4, efficiency was within the same range, while the same difference in the orders for the random systems gives about 10% difference in efficiency of LLE estimation.

Conclusions
The present article shows that modifying LLE estimation from the dot product of a perturbation vector and its time derivative, more efficient algorithm can be obtained. One can notice that the novel method is based on very simple computations, involving only basic mathematical operations, such as addition, subtraction, multiplication and division. A simple, geometrical approach leading to the novel method has been presented. Equivalence of the novel method with the classical one has been proven. Efficiency of the novel method has been confirmed by simulations of various dynamical systems of different orders. Dependence of the computation time on the complexity of the system, and the stability of solution have been analyzed. We have shown that our novel method of the LLE estimation is 17-53% faster than the traditionally applied one. It has been also demonstrated that the  Fig. 6 The estimated LLE value vs. the RCT value t 1 =t 2 for random Jacobi matrices of orders from 2 up to 9 difference in efficiency appears to be the largest for simple systems of small order. The amount of time to be saved decreases with the order and the complexity of the system under consideration, but there exist other properties that influence on efficiency of LLE estimation, which can change this dependence. This feature should be further investigated.