The largest transversal Lyapunov exponent and master stability function from the perturbation vector and its derivative dot product (TLEVDP)

The behavior of systems of coupled nonlinear oscillators and, connected with it, the synchronization phenomena are of significant interest in many areas of science. One of the most important problems in this field is the stability of the synchronous state. The most often applied tool which allows one to quantify this stability is the largest Transversal Lyapunov Exponent (TLE) and, connected with it, Master Stability Function (MSF) theory (Pecora and Carroll in Phys. Rev. Lett. 80:2109, 1998). Thus there is still need to elaborate fast and simple methods of TLE calculation. The new method of the TLE estimation is presented in this paper. It applies the perturbation vector and its derivative Dot Product to calculate the largest Lyapunov exponent in the direction transversal to the synchronization manifold. The method (TLEVDP) of the TLE calculation is very simple in its background and numerical application. It does not require calculations of the Jacobi matrix eigenvalues and orthonormalization in each integration step. Thus it is fast and simple. Used methodology applies the method of Master Stability Function proposed by Pecora and Carroll (in Phys. Rev. Lett. 80:2109, 1998) which is undependable on the system type. Moreover on the base of the one MSF calculation one can predict complete synchronization of the sets of oscillators coupled in any coupling scheme. The theoretical improvement of the method is introduced. Results of the numerical simulations are shown and compared with other publications. Investigations of the system of Duffing oscillators coupled in ring scheme as well as the coupled Van der Pol oscillators made with use of the (TLEVDP) method are presented. Fast stabilization of the TLE value was shown.


Introduction
Dynamics of oscillating systems can be characterized different types of invariant, depending on the dynamical system type and kind of the information that is useful in its investigations. 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 the 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 Benettin 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 were 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. Thus only the largest Lyapunov exponent (LLE), which determines the predictability of the dynamical system, is frequently referred [18][19][20][21][22][23][24]. That is because the presence of at least one positive Lyapunov exponent, by definition, is the most important evidence for chaos [18].
In the special types of the coupled systems LLE can be used for determining stability of the synchronous state [25][26][27][28][29]. In the cases of the more complicated coupling schemes conditional Lyapunov exponents can be calculated. For the coupled identical systems and the complete synchronization stability, one can use Lyapunov exponent measured in direction transversal to the synchronization manifold. It allows obtaining the MSF which is the base for searching the CS in different coupling schemes. Although the method of Master Stability Function assumes that each node system is of the same type with no parameter mismatch, experiments show that it works for the real systems which of course cannot have the same parameters values. In that case synchronization is almost complete [38]. Nowadays TLE and MSF are employed in many different areas of the scientific research [30][31][32][33][34][35][36]. Thus there is still need to elaborate fast and simple methods of TLE calculation. The new method of the TLE estimation is presented in this paper. The method (TLEVDP) applies the perturbation vector and its derivative dot product to calculate the largest Lyapunov exponent in the direction transversal to the synchronization manifold, which is the TLE.

The method
In the article the new method of the (TLE) calculation is presented. The method was applied to obtain (MSF) and then synchronization study of the chosen coupled oscillators were made.

Master stability function
Consider the network of the coupled n identical oscillatorṡ described in the block forṁ where: G-connectivity matrix; H-engaged equations matrix; α-overall coupling coefficient; ⊗-Kronecker product. Examples of such a notation application for the specific oscillators connections are given in the next paragraphs.
After linearization of (1b) one can obtain variational equation in the following form: The idea of the synchronization study with use of the MSF [30][31][32][33] is based on the analysis of the matrix G eigenvalues (λ 1 , . . . , λ n−1 ) considered in transversal directions to the synchronization manifold. The eigenvalue λ n = 0 corresponds to the eigenvector tangent to the synchronization manifold and has no influence on the stability of the synchronous state.
After substituting αλ = γ +iη, where γ = α Re(λ), η = α Im(λ), into (2) one obtains the generic variational equation: where z symbolizes an arbitrary transverse mode. The connectivity matrix G satisfies condition zero row sum so that the synchronization manifold x 1 = x 2 = · · · = x n is invariant and all the real parts of the eigenvalues (λ 1 , . . . , λ n−1 ) associated with transversal modes are negative. Now, one can define the MSF as a surface representing the largest transversal Lyapunov exponent (TLE), calculated for the generic variational equation, over the complex number plane (γ, η). If all the eigenmodes corresponding to the eigenvalues αλ = γ + iη can be found in the range of negative TLE, then the synchronous state is stable for the considered configuration of couplings. In the case of so called real coupling [32,33], when coupling between each two nodes is symmetrical and causes mutual reaction, all the eigenvalues are real. MSF is reduced then to the dependence of the largest TLE on the γ . Taking into account the relation γ = αλ one can rescale the MSF and obtain the curves representing each eigenvalue separately, and conclude from its proceedings about stability of the synchronization manifold for the real coupling values.

The method (TLEVDP) of the largest transversal Lyapunov exponent (TLE) determination
The method for determination of the TLE employs a concept of the vectors dot product. Generally, the essence of this approach is based on the analysis of the vector z(t) obtained from the variational equation (2) representing evolution of the synchronization manifold disturbance. Te graphical interpretation one can see in (Fig. 1). Vector z(t) determines difference between solution y(t) along the synchronization manifold and the disturbed solution x(t): After linearization one can obtain variational equation describing approximate disturbance behavior or: Analyzing transformation of the vector z in the direction of the chosen transversal eigenvector w * of the DF transformation it can be written After transformation: and integration for z * (0) = z * 0 one obtains where: z * -component of the vector z in direction of the eigenvector w * ; λ * -eigenvalue in direction of the eigenvector w * . Equation (8) describes the mean transformation of the vector z * after time t. For t → ∞λ * is a value of Lyapunov exponent in direction of w * . For the n-dimensional transversal eigenvectors set (w * 1 , w * 2 , . . . , w * n ) and respective eigenvalues (λ * 1 , λ * 2 , . . . , λ * n ) one can obtain the whole disturbance z transformation in the following form: where a 1 , a 2 , . . . , a n depend on the initial conditions. In the case λ * > λ * 1 , λ * 2 , . . . , λ * n as t → ∞ terms e λ * 1 t , e λ * 2 t , . . . , e λ * n t become negligible and in the space transversal to the synchronization manifold perturbation evolves as e λ * t what is of the most importance for presented method (TLEVDP) of the largest Transversal Lyapunov Exponent calculation. LLE is of course value of λ * .
In the presented method (5b) is used for the system stability determination. Using vectors dot product properties ( Fig. 1): After transformations: The average valueλ * of the λ * is the parameter allowing one to determine stability of the synchronization manifold in the direction of the disturbance z.
Taking also into account that the direction of the disturbance vector z is convergent to the eigenvector respective to the largest transversal Lyapunov of the transformation DF(x(t)) in (5b) one can conclude that obtained valueλ * is similar to the value of the TLE. However, valuesλ * and TLE can differ a bit. That is because the disturbance vector z is influenced not only by the largest eigenvalue (see (9)). As was mentioned earlier, the direction with respect to the TLE has the main influence on changes of the vector z. Simultaneously during the evolution of the system, the other eigenvalues have influence in that direction too. 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 in it. These are eigenvectors of the DF transformation (5b) treated as the mean transformation along the disturbance z trajectory. Eigenvectors w 1 and w 2 can form the new base of the space, and vector z can be decomposed into directions of these eigenvectors. One gets the component vectors z 1 and z 2 then. Transformation DF acts in that directions by the eigenvalues λ 1 and λ 2 the following way: It means that differences dz 1 and dz 2 depend straightly on the eigenvalues λ 1 , λ 2 and components z 1 and z 2 .
Assuming that Re(λ 1 ) > Re(λ 2 ) > 0 and 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 one can see that vector z rotates and after transformation, the direction of the vector z + dz is closer to the vector w 1 connected with the larger eigenvalue λ 1 . Following further the results of that transformation one can see from Fig. 2 that vector z grows more in direction of w 1 than w 2 . From (12) one can conclude that in the next step of the transformation, growth of the vector z in that direction will be bigger than in the first step. Thus direction of the vector z in each step asymptotically approaches the direction of the vector w 1 connected with the larger eigenvalue λ 1 . Concluding that reasoning the track proposed in the method of the paper of the TLE determination applies the feature of the convergence of the vector z to the direction of the eigenvector connected with the largest eigenvalue of the mean transformation DF which is in our case TLE value.
The most important attributes of the method proposed is its simplicity and connected with it possibility of building the fast acting algorithms. There is no need of the calculating eigenvalues of the Jacobi matrix in each step of the integration, no need of the Gram-Schmidt vector orthonormalization, the vector z can be just normalized, which preserves it from the fast growth.

Numerical simulations
Program code was written in Delphi. The numerical simulations were based on integration of the differential (8), (9) with use of the Runge-Kutta (RK4) method of the fourth order. An integration step was fit to the system behavior as a function of the mean oscillation period. It was important to find the optimal value of the step taking into account the accuracy of simulations on the one side and time of the numerical integration on the second. Especially in some cases TLE stabilized after thousands of the oscillations period.
Depending on the system type and parameters the proper value of the integration step can differ much, thus an integration step was estimated for each parameter case. The value of the required steps for the aver-age period T was obtained by comparing results obtained with use of the procedure of integration with control of the integration error (RungeKutta4system [37]) and standard procedure (RK4). The proper step value was referred to the average oscillations period T and fit to each system parameters set with use of the required steps per period T . Then it was applied into (4) with constant integration step, and in both the equations sets, the system variables x and the z, were integrated in each step of RK4.
To obtain the value λ * at first dot product of the vectors z and dz/dt was calculated. The vector z coordinates were estimated from (4) with use of the procedure (RK4). dz/dt vector was obtained with use of (4) too, and got from within the RK4 procedure during estimation of z coordinates. Vector z was normalized before integration in the RK4 procedure, which preserved it from the fast growth.
Then standard dot product in the form m-the dimension of the node system space; was applied according (8) to calculate the actual value λ * : The value λ * was calculated at each integration step. To obtain its average valueλ * actual values λ * were summed up in time of all the computations, and after a fixed number of the average periods T the average valueλ * of the λ * was calculated. Computations were continued till the increment δλ * of theλ * in the number of the consecutiveλ * calculations periods was smaller than the matched ε value.
The last value attainedλ * was accepted for TLE value.

Synchronization of the four systems coupled in the ring scheme
In this section synchronization of four Duffing systems and then four Van der Pol ones coupled in the ring Fig. 3 Ring coupling scheme scheme (Fig. 3) will be analyzed with use of the presented method. The prepared results concern synchronization of the node systems working in the periodic and chaotic regimes. The mathematical model of such a set of Duffing oscillators can be described by eight differential equations of the first order: The mathematical model of such a set of Van der Pol oscillators can be described by eight differential equations of the first order: Such the networks of the four coupled oscillators can be described in the following block form: where: G-connectivity matrix; H-engaged equations matrix; α-overall coupling coefficient; ⊗-Kronecker product.
Due to the fact that eigenvalues of the H matrix can achieve only values of zero, or one, only α and eigenvalues of the G matrix have influence on the general coupling coefficient. In the considered case the eigenvalues of the G matrix are as follows: and the corresponding eigenvectors: respectively. As was mentioned in Sect. 2.1 the last eigenvalue λ 4 = 0 corresponds to the eigenvector tangent to the Fig. 6 Time dependence of theλ * for the synchronization and desynchronization states Fig. 7 Comparison of the MSF obtained with MSFDP method with results presented in [35] synchronization manifold and has no influence to the stability of the synchronous state.
For comparing results obtained with use of the presented method with other results, simulations were followed with ones presented in [35].
The first analyzed is the ring of the four Duffing oscillators working in the periodic regime. For periodic systems synchronization analysis the following system parameters values were assumed: β = 0.1, η = 1.0, q = 7 (16). Comparison of the MSF obtained with use of the presented method with the one presented in [35] can be seen in Fig. 4. Convergence of the results prove the efficiency of the TLEVDP method. After rescaling MSF with use of the eigenvalues λ 1 = −4, λ 1,2 = −2, with respect to the dependence αλ = γ , Fig. 5, one obtains graphs of the functions dependent on the real coupling coefficient α, which show the stability of the synchronization manifold in the directions of the eigenvectors v 1 , v 2 , v 3 . One can conclude about complete synchronization feature with use of these graphs. For the range of α where both of the rescaled func-  (22) was contained in the presented in Fig. 5, chart. One can see that for the range of α where both of the rescaled functions achieve negative values, synchronization er-ror achieves values close to zero. All these facts prove the efficiency of the presented TLEVDP method.
The time dependent results of theλ * estimation is shown in Fig. 6. The diagrams presented show the behavior of theλ * for negative and positive TLE values. In the case TLE > 0 the value of theλ * stabilizes in almost 100 excitation periods. Computations were continued till the increment δλ * of theλ * in the number of the consecutiveλ * calculations periods was smaller than the matched ε = 1e−4. In the case TLE < 0 value of theλ * stabilization takes a little more time, but after 400 oscillations period differences became very small. The second analyzed is the ring of the four Duffing oscillators working in the chaotic regime. For chaotic systems synchronization analysis the following system parameters values were assumed: β = 0.1, η = 1.0, q = 5.6 (16). Comparison of the MSF obtained with use of the presented method with the one presented in [35] can be seen in Fig. 7. One can see at first that points with zero TLE values agree each other. Convergence of the results prove the efficiency of the TLEVDP method in the synchronization feature study. However, TLE values differ a bit. That is because the disturbance vector z is influenced not only by the largest eigenvalue, what was explained in par.
(1.2). After rescaling MSF with use of the eigenvalues λ 1 = −4, λ 1,2 = −2, with respect to the dependence αλ = γ , Fig. 8, one obtains graphs of the functions dependent on the real coupling coefficient α, which show the stability of the synchronization manifold in directions of the eigenvectors v 1 , v 2 , v 3 . One can conclude to complete synchronization feature with use of these graphs. For the range of α where both of the rescaled functions achieve negative values, complete synchronization in the system of coupled oscillators occurs. For the synchronization proof the bifurcation diagram of the synchronization error was contained in the presented in Fig. 8, chart. One can see that for the range of α where both of the rescaled functions achieve negative values, and the synchronization error achieves values close to zero. All these facts prove the efficiency of the presented TLEVDP method.
The time dependent results ofλ * estimation are shown in Fig. 9. The diagrams presented show the behavior of theλ * for negative and positive TLE values. One can see that in this case stabilization of theλ * values takes much more time for chaotic coupled systems stabilization than for periodic ones. In the case TLE < 0 value of theλ * stabilizes longer and it takes almost 10 000 oscillations periods. As for the periodic systems computations were continued till the increment δλ * of theλ * in the number of the consecutivê λ * calculations periods was smaller than the matched ε = 1e−4. In the case TLE < 0 value of theλ * stabilization takes less time, but there is still long time needed forλ * stabilization.
The third analyzed is the ring of the four Van der Pol oscillators. For systems synchronization analysis the following system parameters values were assumed: β = 0.401, η = 0.366 (17). MSF obtained with use of the presented method is presented in Fig. 10. After rescaling MSF with use of the eigenvalues λ 1 = −4, λ 1,2 = −2, with respect to the dependence αλ = γ , one obtains graphs of the functions dependent on the real coupling coefficient α, which show stability of the synchronization manifold in directions of the eigenvectors v 1 , v 2 , v 3 Fig. 11. One can conclude to complete synchronization feature with use of these graphs. For the range of α where both of the rescaled functions achieve negative values, complete synchronization in the system of coupled oscillators occurs. For the synchronization proof the bifurcation diagram of the synchronization error was contained in the presented chart. One can see that for the range of α where both of the rescaled functions achieve negative values, the synchronization error achieves values close to zero. All these facts prove the efficiency of the presented TLEVDP method.

Conclusions
The new method of the TLE estimation is presented in this paper. It applies the perturbation vector and its derivative Dot Product to calculate the largest Lyapunov exponent in the direction transversal to the synchronization manifold. It was shown that the method (TLEVDP) is very simple in its background and numerical application. Theoretical improvement was introduced. Results of the numerical simulations were shown and compared with other publications [35]. Investigations of the Duffing oscillators coupled in a ring scheme were made with use of the (TLEVDP) method. Predicted regions of the synchronization existence were confirmed with use of the synchronization error diagrams. Times of the TLE stabilization were compared for the cases presented. Fast stabilization of the TLE value was proved.