Estimation of the largest Lyapunov exponent-like (LLEL) stability measure parameter from the perturbation vector and its derivative dot product (part 2) experiment simulation

Controlling system dynamics with use of the Largest Lyapunov Exponent (LLE) is employed in many different areas of the scientific research. Thus, there is still need to elaborate fast and simple methods of LLE calculation. This article is the second part of the one presented in Dabrowski (Nonlinear Dyn 67:283–291, 2012). It develops method LLEDP of the LLE estimation and shows that from the time series of two identical systems, one can simply extract value of the stability parameter which value can be treated as largest LLE. Unlike the method presented in part, one developed method (LLEDPT) can be applied to the dynamical systems of any type, continuous, with discontinuities, with time delay and others. The theoretical improvement shows simplicity of the method and its obvious physical background. The proofs for the method effectiveness are based on results of the simulations of the experiments for Duffing and Van der Pole oscillators. These results were compared with ones obtained with use of the Stefanski method (Stefanski in Chaos Soliton Fract 11(15):2443–2451, 2000; Chaos Soliton Fract 15:233–244, 2003; Chaos Soliton Fract 23:1651–1659, 2005; J Theor Appl Mech 46(3):665–678, 2008) and LLEDP method. LLEDPT can be used also as the criterion of stability of the control system, where desired behavior of controlled system is explicitly known (Balcerzak et al. in Mech Mech Eng 17(4):325–339, 2013). The next step of development of the method can be considered in direction that allows estimation of LLE from the real time series, systems with discontinuities, with time delay and others.


Introduction
Depending on the dynamical system type and kind of the information that is useful in its investigations, there are applied different types of invariants characterizing the system dynamics. One can use for instance Kolmogorov entropy [1] or correlation dimension [2,3] to determine chaotic level or complexity of the system dynamics [4]. But when there is a need to predict behavior of the real system with possibility of different disturbances existence, Lyapunov exponents are one of the most often applied tools. That is because, these exponents determine the exponential convergence or divergence of trajectories that start close to each other. The existence of such numbers has been proved by Oseledec theorem [5], but the first numerical study of the system behavior using Lyapunov exponents has been done by Henon and Heiles [6], before the Oseledec theorem publication. The most important algorithms for calculating Lyapunov exponents for a continuous systems have been developed by Benettin et al. [7] and Shimada and Nagashima [8], later improved by Benet-tin et al. [9,10] and Wolf [11]. For the system with discontinuities or time delay, one possible approach is the estimation of Lyapunov exponents from the scalar time series basing on Takens procedure [12]. Numerical algorithms for such estimation have been developed by Wolf et al. [13], Sano and Sawada [14], and later improved by Eckmann et al. [15], Rosenstein et al. [16] and Parlitz [17].
The set of Lyapunov exponents contains much physical information characterizing the considered dynamical system, but calculation of the full spectrum demands much time and labor. Hence, only the largest Lyapunov exponent (LLE), which determines the predictability of the dynamical system, is frequently referred. That is because, excluding systems with Perron effect [39], the presence of at least one positive Lyapunov exponent, by definition, is the most important evidence for chaos [18]. The algorithm for calculating the largest Lyapunov exponent was independently presented by Rosentein et al. [16] and Kantz [19] These methods make use of the statistical properties of the local divergence rates of nearby trajectories. An improved algorithm based on Rosentein and Kantz was recently presented by Kim and Choe [20]. The next method of the LLE calculation was introduced by Stefanski [21][22][23][24]. This method based on the synchronization phenomena allows the LLE estimation for both, continuous and not continuous systems, and thus can be applied for system with flow and maps dynamics representation.
Nowadays, LLE is employed in many different areas of the scientific research [25][26][27][28][29][30][31][32][33][34][35][36][37]. Thus, there is still need to elaborate fast and simple methods of LLE calculation. The new method of the LLE estimation is presented in this paper. It is shown that from the behavior of the two identical systems, one can extract stability parameter LLEL which value can be treated as the LLE.

The LLEDPT method
Generally presented method bases on the analysis of the disturbance changes dz(t) in direction of the general disturbance vector z(t) (Fig. 1) and was discussed in [38]. Note that presented here and in [38] theoretical background is just only general way of thinking and cannot be treated as a formal proof. Proofs of the method efficiency one can find in results of the numerical simulations. Vector z(t) determines difference between reference flow y(t) and disturbed flow x(t): The main idea follows the fact that system is stable when the average changes of perturbation dz(t) cause average decrease of z(t).
At each fixed time t, disturbance vector z(t) can obtained from: While in [38] behavior of z(t) was analyzed with use of the linearized variational Eq. (2) where: U(y(t)) is the Jacobi matrix in the point y(t) in present case, we can calculate actual z(t) and dz(t) from the behavior of the two real systems x(t) and y(t).
Thus, U(y(t)) is no longer needed in the procedure of studying behavior of the perturbation z(t). It is one of the advantages of the presented method. One can investigate stability of the dynamical system without any knowledge about analysed systems equations. Analysis of the transformation of the vector z in the direction of the chosen transversal eigenvector w* ( Fig. 1): In each step of the simulations for fixed time t transformation, U can be considered as linear transformation of the vector z. Assuming that λ * is eigenvalue of the transformation U in direction of the eigenvector w * it can be written: After the transformation: (4) and the integration for z * (0) = z * 0 , one obtains: where: |z * |-norm of the projection of vector z on to the direction of the eigenvector w * . λ * -eigenvalue in the direction of the eigenvector w * .
Equation (5) describes a mean transformation of the vector z * after time t. For t → t * , where t * is a time of stabilization of mean λ * , it can be treated as a value of the stability parameter that allows to estimate exponential convergence, or divergence of nearby trajectories measured in the direction of w * . For the m-dimensional transversal eigenvec- , one can obtain the whole disturbance z transformation in the following form: where a 1 , a 2 , . . ., a m depend on the initial conditions. In the case λ * > λ * 1 , λ * 2 , . . ., λ * m as t → t * , the terms e λ * 1 t , e λ * 2 t , . . . , e λ * m t become negligible and the perturbation evolves as e λ * t , which is of the most importance for the presented method. The LLES is of course a mean value of λ * .
In the presented method, Eq. (3) is used for the system stability determination. Using vectors dot product properties ( Fig. 1): After simple transformations: The average valueλ * of real part of λ * is the parameter allowing us to determine stability of the system along the direction of the disturbance z.
Taking also into account that the direction of the disturbance vector z and directions of the initial frame of the orthonormal vectors during the evolution converge to the eigenvector with the largest eigenvalue, of the mean transformation U (y(t)) of all transformations U(y(t)) along the phase trajectory, one can conclude that the obtained valueλ * is similar to the value of the LLE.
However, valuesλ * and the LLE can differ slightly. That is because the disturbance vector z is influenced not only by the largest eigenvalue (see Eq. 6). Directions of the eigenvector with the biggest eigenvalue and main principal axis with the biggest length change are convergent to direction of the disturbance z. But these directions are not the same. As has been mentioned earlier, the direction respective to the LLE has the main influence on vector z changes. Simultaneously, during the evolution of the system, other eigenvalues exert their influence in that direction as well. It causes thatλ * shows general stability of the synchronization manifold. A simplified graphical interpretation of such deliberations is presented in Fig. 2. One can see vectors w 1 and w 2 on it. These are eigenvectors of the U transformation (Eq. 3) treated as the mean transformation along the disturbance z trajectory. The eigenvectors w 1 and w 2 can form a new basis of the space, and the vector z can be decomposed into directions of these eigenvectors. One gets the component vectors z 1 and z 2 then. The transformation U acts in those directions by the eigenvalues λ 1 and λ 2 in the following way: It means that the differences dz 1 and dz 2 depend straightly on the eigenvalues λ1, λ2 and the components z 1 and z 2 . Assuming that Re(λ 1 ) > Re(λ 2 ). > 0 Im(λ 1 ) = Im(λ 2 ). = 0, one can see in Fig. 2 that in the case z 1 = z 2 , it causes dz 1 > dz 2 . From Fig. 2 it follows that the vector z rotates and after some transformations, the direction of the vector z + dz is closer to the vector w 1 connected to the larger eigenvalue λ1. Following further the causes of that transformation, it can be seen from Fig. 2 that the vector z grows more in the direction of w 1 than w 2 . From Eqs. (9), one can conclude that in the next step of the transformation, the growth of the vector z in that direction will be bigger than in the first step. Thus, the direction of the vector z in each step approaches asymptotically the direction of the vector w 1 connected to a larger eigenvalue λ1. Con- Fig. 2 Simplified graphical interpretation of the U transformation influence onto vector z cluding that reasoning, the proposed method of the LLE determination applies the feature of the convergence of the vector z to the direction of the eigenvector connected to the largest eigenvalue of the mean transformation U (y(t)), which is in our case the LLEL value.
In results of the numerical simulations, one can find that LLEL has the same value as LLE. Thus, since that moment LLEL and LLE will be treated and named as the same stability parameter LLE. The most important attribute of the LLEDPT method is its simplicity and a possibility of building fast acting algorithms. There is no need to calculate eigenvalues of the Jacobi matrix in each step of the integration, no need of the Gram-Schmidt vector orthonormalization. The vector z has to be only normalized, which stops it from the growth which introduces inaccuracy into the method. That inaccuracy has its background in the fact, that the method bases on the behavior of the linearized system. Thus, in the case, when z grows too much, Eq. (2) is no longer valid. The LLEDPT is straightly connected to the real stability of the system because of the analysis of the exponential divergence and convergence in the direction of the perturbation vector z. Thus, it could be more efficient in experimental investigations of the real systems.

Numerical simulations
Behavior of two real dynamical systems was simulated with use of the program written in Delphi. The numerical simulations were based on integration of the differential equations (8,9) with use of the Runge Kutta (RK4) method of the fourth order. The step of integration was fit to the system type with use of the procedure of integration with control of the integration error. Then, it was applied into the procedure RK4 with constant integration step. Both the equations sets x(t) and y(t), the system variables x and the perturbation variables z, were integrated in each step of RK4. After each integration step considered in Eq. (6), the value λ * was calculated.
An average value of theλ * was calculated in time of all the computations. The attained value of the LLE was theλ * value for theλ * increment smaller than the matched ε value. Dynamical state of the perturbed system has been reduced into reference one neighborhood on each 15 periods of external excitation.

Van Der Pol oscillator with external excitation
First, analyzed example shows results of LLE extraction from behavior of two Van der Pol oscillators with external excitation (11) In the bifurcation diagram (Fig. 3), one can see different types of the system dynamics and values of LLE obtained with use of the presented method confirming its efficiency. One can observe negative LLE in periodic regions, positive in chaotic ones and zero LLE values that determine period doubling bifurcation points are visible as well. The time of the LLE estimation for each bifurcation parameter value depends on the actualλ * increment. The attained value of the LLE was theλ * value for theλ * increment smaller than the matched ε value. To determine the exact bifurcation parameter values maintaining the fast algorithm operation, ε was changed depending on the actual LLE value.
Comparison of the LLE estimated with use of the LLEDPT method and LLEDP method [38] is shown in Fig. 4. As the results are very similar, one cannot easily see differences between white and black curves. Next, analyzed example shows results of LLE extraction from behavior of two Duffing oscillators with external excitation (12). 3 1 + q sin(ηt) y 2 = −βy 2 − αy 3 1 + q sin(ηt) (12) Analogically to Van der Pol example in Fig. 5 [20] and LLEDPT method [38] is shown in Fig. 6. For the next time, it proves the validity of the taken assumptions.

Conclusions
The next step of development of new LLE estimation method was presented. As the method does not need equations of motion, it can be applied to the dynamical systems of any type, continuous, with discontinuities, with time delay or as the method that allows estimation of LLE from the real time series. Other proposition is to exploit the method as criterion of the control system stability, where desired behavior of controlled system is explicitly known. The theoretical improvement was also introduced. It shows simplicity of the method and its obvious physical background. Presented results are based on simulations of the experiments for Duffing and Van der Pol oscillators. Results of the numerical simulations were compared with ones obtained with use of the Stefanski method [21][22][23][24] and LLEDPT method presented in [38]. The study confirms effectiveness of the presented method of LLE estimation.