Numerical Approximation of Lyapunov-Exponents for Quasiperiodic Motions

This paper proposes an approach to approximate the Lyapunov-spectrum of quasiperiodic flows on isolated invariant manifolds numerically. Once the invariant manifold has been determined, integrations over the infinite, one dimensional time interval – as calculating the Lyapunov-spectrum for instance – can be transformed into an integral over a finite, p-dimensional domain, where p is the dimension of the manifold. The application of the proposed approach is demonstrated by calculating the Lyapunov-spectrum of periodic and quasiperiodic motions of a forced van-der-Pol equation. The results are compared to results from a classical time integration based method using a continuous Gram-Schmidt orthonormalization.


Introduction & Prerequisites
Dynamical systems excited by mechanisms with one single frequency usually exhibit periodic oscillations. As soon as a second mechanism excites the dynamical system simultaneously, oscillations of different physical origin interfere which may lead to quasiperiodic motions, unless synchronization is existent. Due to the infinite period length the description of quasiperiodic motions using standard methods (which are mainly tailored for periodic processes) is usually either not possible or very disadvantageous. The alternative approach proposed in this contribution is based on isolated (normally hyperbolic) invariant manifolds of the flow, which can be determined by solving a system of partial differential equations (invariance equations) [3,[6][7][8].
Due to the invariance property, solutions on the manifold will never leave it. Moreover, it is assumed that the manifold is isolated, i.e. neighboring solutions are either attracted or repelled. Finally, it is assumed that the flow fills the manifold densely as t → ∞. Since on the other hand the Lyapunov-exponents converge for time approaching infinity (Oseledec's Theorem), the description using an invariant p-dimensional manifold for quasiperiodic motions basically includes all required informations; here, p is the number of internal (base) frequencies. Now, the transition of the infinite time interval to a finite pdimensional domain is performed by applying the meanvalue theorem [6]. A prerequisite for the application of the mean-value theorem is the presence of a parallel flow on the manifold. Since the method described in [3] -which is based on a physical amplitude-phase-parametrizationdoes not directly provide parallel flows, a transformation from the physical to a parallel parametrization has been introduced. By calculating the corresponding diffeomorphism, the mean-value theorem becomes applicable.

Calculation of the Lyapunov-exponents
The Lyapunov-exponents are an established quantity to assess the stability of arbitrary motions. In general they describe the local exponential growth or decrease rate of perturbations applied to a known solution. Consider an autonomous dynamical systeṁ which is linearized about a reference solution z 0 . This solution has n Lyapunov-exponents of first order σ 1 ≥ σ 2 ≥ ... ≥ σ n which describe the averaged local exponential decay rates near this solution. Together they constitute the Lyapunov-spectrum. If no exponents is positive the corresponding motion is considered as stable. Beyond these exponents of first order one may define m th. order exponents σ (m) (m = 1, ..., n), which describe the growth or decrease of m-dimensional volumes in the state-space near z 0 and are calculated by all possible sums of m different Lyapunov-exponents of first order [1]. In order to evaluate the Lyapunov-exponents over a finitedimensional invariant manifold, a connection between time parametrization and the invariant manifold describing the quasiperiodic motion is required.

The mean-value theorem
In [6] a possible way to connect quantities over a time parametrization and an invariant manifold describing quasiperiodic motions is given through where g is an arbitrary function, z is a quasiperiodic solution andθ = [θ 1 , ...,θ p ] is the parametrization (torus coordinates) of an invariant manifold exhibiting a parallel flow.
In the following variables with a bar( .) refer to parallel flows. Since the assumption was made, that an invariant manifold is calculated with the method described in [3], the flow on the manifold is not parallel because of the usage of physical coordinates. But since the parallel flow is a canonical form of the physical flow, there has to be a diffeomorphism connecting both parametrization, hence where θ = [θ 1 , ..., θ p ] is the physical parametrization (torus coordinates) of an invariant manifold and φ(θ) =θ is a diffeomorphism. It is important to note, that if an invariant manifold is calculated with the methods described in [7,8] a parallel flow is given and therefore section 2.2 can be skipped.

Invariance equation of the diffeomorphism
In [3] equation (1) is first transformed to amplitudes A(θ) and phases with a subsequent partitioning of the phases in dependent phase angles ψ(θ) and torus coordinates (parametrization) θ resulting iṅ Using these for solving the invariance equation the quantities h(z), Ω(z) and Ψ(z) are known. By assuming a parallel flow, the equations of motion read: It is important to note, that by this assumption the time derivation of the parametrization is constantω. It is easily comprehensible by looking at figure 1. The flow on the invariant manifold represents the solution of the dynamical system (1) in time. By choosing an arbitrary point (θ * 1 ,θ * 2 ) on the invariant manifold, the change ofθ 1 andθ 2 in time has always to be the same because otherwise there would be no parallel flow and therefore hasθ to be constant. In section 2.1 is the diffeomorphism introduced as a function describing the deformation to a parallel flowθ = φ(θ), whereby a substitution in equation (5) is possible. For further considerations it is sufficient to investigate the resulting equations for the parametrization Applying the chain rule and identifying the time derivation of the parametrization of the physical flow in (4) lead to the invariance equation of the diffeomorphism where Ω(θ) is known through the calculation of the invariant manifold. Unknown in this equations are the internal frequenciesω and the wanted diffeomorphism φ(θ). Equation (8) is just as the invariance equation for the manifold 2π-periodic in θ. Since the diffeomorphism is only in the derivatives existent and does not have any restrictions concerning the parametrization, the solution of equation (8) is not unique. Through additional boundary conditions, e.g. φ(0) = 0 which corresponds to equal origins of the parametrization with θ andθ, the equation gets solvable.
From this an equation to calculate the Lyapunov-exponent of n th. order can be derived [1] Applying the mean-value theorem with respect to the physical flow (3) a criterion for a known manifold results as (11) This equation can be evaluated once the invariance equation of the diffeomorphism (8) has been solved. Note that this criterion is also applicable on manifolds with a parallel flow, since the Jacobian determinant is then one.

The Lyapunov-spectrum
In contrast to the Lyapunov-exponent of n th. order, the exponents of first order are associated with characteristic directions. This results in numerical difficulties which can be solved by different approaches. An overview can be found in [5]. An established method for the calculation of the whole spectrum or the m largest ones bases upon a differential formulation of the Gram-Schmidt process [4], which is enhanced by a stability parameter in [2]. The equations of motion reaḋ where J is the Jacobian matrix, E = [e 1 , ..., e k ] is a set of time-dependent orthonormal vectors and G m = J mm + β (e m · e m − 1) with the stabilization parameter β > −σ k , which requires some a priori knowledge. L lm = J lm + J ml + 2β e l · e m and J lm = e l · (J · e m ). The Lyapunov-exponents of first order can be calculated by whereby they can be seen as the temporal mean value. These values converge almost surly for arbitrary initial values z 0 and E 0 (Oseledec's theorem). By integration ofΛ m under consideration of the initial values and substituting it in equation (13) results. Considering the mean-value theorem lead to where the basis vectors e m are unknown. Since they depend only on time, it is possible to define them on the manifold. Applying the chain rule, a system of partial differential equations for calculating the basis vectors results. As long as the spectrum is not degenerated, i.e. when two or more Lyapunov-exponents are equal or from a numerical point of view very close to each other, equation (16) has periodic boundary conditions in θ and has unique e m . Inserting the results from equations (8) and (16) in (15) a criterion for the Lyapunov-spectrum is defined.

Application to the forced van-der-Pol equation
To demonstrate and analyze the proposed approach, it is applied to the forced van-der-Pol equation where ε = 0.3 and Γ = 0.6. The excruciation frequency is chosen as η ∈ [0.3, 1.5]. Due to interfering mechanisms of excitation quasiperiodic motions as well as periodic motions (due to synchronization) can be observed. In order to verify the numerical results for the Lyapunovexponents, equation (12) is solved with a time integration for τ = [0, 20000]. Results calculated with a time integration are identified by τ (.). For a given invariant manifold the calculation of the diffeomorphism to assure for parallel flows is the first task, since both the criterion for σ (n) and the one for the Lyapunov-spectrum require a parallel flow. If σ (n) is of interest, equation (11) can be evaluated with no further calculations. Is the calculation of the spectrum of interest, the invariant basis E on the manifold has to be determined. Solving equation (16) the basis is known and the criterion described by equation (15) can be evaluated. In figure   3 3

A [−]
A Time Sim. Stable periodic Unstable periodic Stable quasiperiodic 2 are basis at random points depicted. It is noteworthy that the tangential space of the invariant manifolds at the given points is visible and therefore the directions of the Lyapunov-exponents equal to zero. Knowing the invariant basis at every point, the Lyapunov-spectrum can be calculated. Figure 3 shows the resulting σ (n) for invariant manifolds describing periodic and quasiperiodic motions. It is conspicuous that the quasiperiodic solution can not be calculated to the bifurcation point. The reason for that is the parametrization of the underlying method of the manifold approximation. By a modification of the method this behavior can be disabled but should not be discussed here further. The comparison with the time simulation reveals, that both methods exhibit the same results. In the case of the van-der-Pol equation the calculation of the sum is sufficient to state stability of a quasiperiodic motions since a T 2 in R 3 has two Lyapunov-exponents equal to zero. In the case of higher dimensional spaces R n with embedded T p the necessity of knowing the Lyapunov-spectrum, or at least the largest exponent, is essential for identifying stability of motions, since the absolute value of the sum of the negative exponents can be larger than the one of the positive ones. The outcome is a negative sum containing positive Lyapunov-exponents resulting in an unstable motion. In contrast to that states a positive sum always an unstable motion since then at least one exponent is positive. To demonstrate the identification of the spectrum, the proposed approach is used to calculate it at 1 and 2 .  Comparing both methods shows that the results are almost equal. Only σ 3 differs from the time simulation: however, this difference stems from accumulated numerical errors of the underling methods. The advantages resulting from the proposed approach is on the one hand that especially for quasiperiodicity the Lyapunov-spectrum is known without concerns about convergence and on the other hand can the approach be applied to unstable solution which is not possible with time simulations. Reversing the time is not expedient since both the unstable and stable exponents revers their behavior resulting in an unstable motion. Concluding it has to be mentioned that for large systems the calculation of the spectrum is costly, because for equation (16) n 2 equations have to be solved. However the evaluation of equation (11) itself is fast once the invariant manifold has been determined with the methods proposed in [3,7,8], for instance.