Fast Indicators for Orbital Stability: A Survey on Lyapunov and Reversibility Errors

*Giorgio Turchetti and Federico Panichi*

#### **Abstract**

We present a survey on the recently introduced fast indicators for Hamiltonian systems, which measure the sensitivity of orbits to small initial displacements, Lyapunov error (LE), and to a small additive noise, reversibility error (RE). The LE and RE are based on variational methods and require the computation of the tangent flow or map. The modified reversibility error method (REM) measures the effect of roundoff and is computed by iterating a symplectic map forward and backward the same number of times. The smoothest indicator is RE since it damps the oscillations of LE. It can be proven that LE and RE grow following a power law for regular orbits and an exponential law for chaotic orbits. There is a numerical evidence that the growth of RE and REM follows the same law. The application to models of celestial and beam dynamics has shown the reliability of these indicators.

**Keywords:** variational principles, reversibility error, additive noise, roundoff

#### **1. Introduction**

The global stability properties of Hamiltonian systems and symplectic maps have a solid theoretical foundation [1, 2]. Nevertheless, the determination of the orbital stability by computing the maximum Lyapunov exponent is a procedure difficult to implement numerically, because of the *t* ! ∞ limit. For this reason a variety of fast indicators has been developed during the last two decades [3–7]. The variational methods mentioned above measure the sensitivity to initial conditions of the orbit computed for finite times. The spectral methods [8, 9] relate the stability to the behavior of the Fourier spectrum of the orbit computed for finite times.

In the framework of the variational methods, we have proposed two indicators [10–12] the Lyapunov error (LE) and the reversibility error (RE) introducing also the modified reversibility error method (REM). The LE is due to a small displacement of the initial condition, the RE is due to an additive noise, and REM is due to roundoff. The reversibility error due to the roundoff or noise is more convenient with respect to the error occurring in the forward evolution of the map.<sup>1</sup>

<sup>1</sup> The forward error (FE) due to additive noise in the forward evolution of a map can be defined and written in terms of the tangent map. However, RE is very simply related to LE, whereas FE is not. In addition the error due to roundoff in the forward evolution requires in principle the evaluation of the exact orbit or, in practice, its evaluation with a much higher accuracy.

In the limit of a vanishing amplitude of the initial displacement or of the random displacement, the LE and RE are defined by using the tangent map along the orbit. Furthermore, RE is related to LE by a very simple formula. A reversibility error is always present in numerical computations due to roundoff even when no additive noise is introduced. We compute REM by iterating *n* times the map *M*, then its inverse *n* times, and dividing the norm of the displacement from the initial point, by the roundoff amplitude. The procedure is extremely simple and does not require the knowledge of the tangent map. Though the effect of roundoff on a single iteration is not equivalent to a random displacement, after many iterations the cumulative result is comparable if the computational complexity of the map is sufficiently high. The main difference is that for an additive noise, the error is defined as the root mean square deviation of the noisy orbit with respect to the exact one, obtained by averaging over all possible realizations of the noise, whereas for the roundoff a unique realization is available. As a consequence REM fluctuates with the iteration number *n*, whereas RE does not. A statistical analysis of roundoff compared to a random noise was previously performed using the fidelity method [13, 14], and a comparison of REM with other fast indicators was initially carried out for the standard map [15]. The growth of errors, for REM-, RE-, and Lyapunov-based indicators, is governed by the tangent map. For LE a small initial displacement is propagated and amplified along the orbit. For RE or REM, a random or pseudorandom displacement is introduced at any (forward or backward) iteration of the map and is propagated and amplified along the orbit. The final random displacement is the sum of the global displacements triggered by the local displacements (due to noise or roundoff) occurring at any iteration. Therefore, it is not surprising that the square of RE is twice the sum of the squares of LE computed along the orbit and that all the numerical experiments suggest that REM exhibits a similar behavior even though with larger fluctuations.

is obtained from the covariance matrix which is computed from the tangent map.

*Fast Indicators for Orbital Stability: A Survey on Lyapunov and Reversibility Errors*

*n*�1 � �<sup>2</sup> <sup>þ</sup> *<sup>e</sup><sup>L</sup>*

linear maps to explore the behavior of REM. A systematic comparison of LE, RE, and REM is presented for two basic models: the standard map and the Hénon map. The asymptotic power law exponents are computed by using the MEGNO filter. For nonlinear two-dimensional maps, the behavior of the errors has been compared moving along a one-dimensional grid in the phase plane: crossing of islands has a clear signature, the chaotic regions are very neatly distinguished, and good agree-

A rectangular region of phase plane has been examined by choosing a grid and using a color, logarithmic scale for the errors at each point. Also in this case, a good correspondence with the phase space portrait is found. On the basis of the analysis presented here and the experience gained in investigating more complex models from celestial mechanics [18] and beam dynamics [19], we suggest to compare RE and REM with LE, possibly, filtered with MEGNO, to damp the oscillations<sup>2</sup>

maps of dimension 4 or higher, a direct geometric inspection of the orbits is not possible since the Poincaré section requires an interpolation Hamiltonian. As a consequence the use of fast indicators is the only practical approach to analyze the orbital stability. Hamiltonian systems have a continuous time flow, and the errors LE and RE denoted by *<sup>e</sup><sup>L</sup>*ð Þ*<sup>t</sup>* and *<sup>e</sup><sup>R</sup>*ð Þ*<sup>t</sup>* , respectively, are computed by using the fundamental matrix Lð Þ*<sup>t</sup>* of the tangent flow. In this case *eL*ðÞ¼ *<sup>t</sup>* Tr L*<sup>T</sup>*ð Þ*<sup>t</sup>* <sup>L</sup>ð Þ*<sup>t</sup>* � � � � <sup>1</sup>*=*<sup>2</sup>

approximation gives the relation found for the maps [12]. Standard procedures allow to approximate the orbit **<sup>x</sup>**ð Þ*<sup>t</sup>* by the iterates *<sup>M</sup><sup>n</sup>*ð Þ **<sup>x</sup>**<sup>0</sup> of a symplectic integrator map *<sup>M</sup>* (see [20]) and the fundamental matrix Lð Þ*<sup>t</sup>* by *DM<sup>n</sup>*ð Þ **<sup>x</sup>**<sup>0</sup> (see [21]). The paper has the following structure. In Section 2, we recall the definitions for LE and RE and obtain their mutual relation. In Section 3, we present the analytical results on LE and RE together with the numerical results on REM for integrable maps. In Section 4, the key features of two prototype models, the standard map and the Hénon map, are summarized. In Section 5, we present a detailed numerical analysis

of LE, RE, and REM for the standard map. In Section 6, the same analysis is presented for the Hénon map. In Section 7, the summary and conclusions are

Given a symplectic map *<sup>M</sup>*ð Þ **<sup>x</sup>** where **<sup>x</sup>**<sup>∈</sup> <sup>R</sup><sup>2</sup>*<sup>d</sup>*, we consider the orbits

<sup>ϵ</sup> <sup>¼</sup> lim<sup>ϵ</sup>!<sup>0</sup>

<sup>2</sup> The application of MEGNO to RE is not necessary due to the absence of oscillations, whereas its application to REM is not recommended because the fluctuations are not filtered and the computational

**y***<sup>n</sup>* � **x***<sup>n</sup>*

tively, where *η* is a unit vector. We consider the normalized displacement *η<sup>n</sup>* at

*<sup>n</sup>* , which has a very simple relation with LE given by the

. We first analyze the case of

, whose trapezoidal rule

. For

(1)

*n* � �<sup>2</sup>

<sup>0</sup> *ds e<sup>L</sup>*ð Þ*<sup>s</sup>* � �<sup>2</sup>

� � for two initial points **<sup>x</sup>**<sup>0</sup> and **<sup>y</sup>**<sup>0</sup> <sup>¼</sup> **<sup>x</sup>**<sup>0</sup> <sup>þ</sup> <sup>ϵ</sup>*η*0, respec-

� � � *<sup>M</sup>*ð Þ **<sup>x</sup>***<sup>n</sup>*�<sup>1</sup> ϵ

*<sup>M</sup>* **<sup>y</sup>***<sup>n</sup>*�<sup>1</sup>

We denote this error by *e<sup>R</sup>*

0 � �<sup>2</sup> <sup>þ</sup> <sup>2</sup> *eL*

*DOI: http://dx.doi.org/10.5772/intechopen.88085*

1

ment with the theoretical predictions is found.

and *eR*ð Þ*<sup>t</sup>* are given by the square root of 2 <sup>Ð</sup>*<sup>t</sup>*

� �<sup>2</sup> <sup>þ</sup> … <sup>þ</sup> <sup>2</sup> *<sup>e</sup><sup>L</sup>*

square root of *e<sup>L</sup>*

presented.

**163**

**2. Definition of errors**

iteration *n* defined by

**<sup>x</sup>***<sup>n</sup>* <sup>¼</sup> *<sup>M</sup><sup>n</sup>*ð Þ **<sup>x</sup>**<sup>0</sup> and **<sup>y</sup>***<sup>n</sup>* <sup>¼</sup> *<sup>M</sup><sup>n</sup>* **<sup>y</sup>**<sup>0</sup>

*<sup>η</sup><sup>n</sup>* <sup>¼</sup> lim<sup>ϵ</sup>!<sup>0</sup>

cost is quadratic rather than linear in the iteration order.

For an integrable map, the growth of LE and RE follows asymptotically a power law *n<sup>α</sup>*, and the exact analytical result is known. This result can be extended to quasi integrable maps by using the normal forms theory. For uniformly chaotic maps (hyperbolic automorphisms of the torus), the LE and RE have an exponential growth *eλ<sup>n</sup>*. For generic maps, the asymptotic growth of LE and RE follows a power law in the regions of regular motion and an exponential law in the regions of chaotic motion, and the same behavior is observed for REM. For an integrable or quasi integrable map, LE has an asymptotic linear growth *α* ¼ 1 with oscillations, whereas RE has an asymptotic power law growth with *α* ¼ 3*=*2 without oscillations, since they are rapidly damped. The oscillations of LE disappear when the map is written in normal coordinates. For a linear map conjugated to a rotation, the power law exponents are *α* ¼ 0 for LE and *α* ¼ 1*=*2 for RE. For REM the power law exponent *α* varies between 0 and 1, its value depending on the computational complexity of the map and therefore on the choice of coordinates.

The definition of LE we propose differs from fast Lyapunov indicator (FLI) [3] or orthogonal fast Lyapunov indicator (OFLI) [4], which are based on the growth along the orbit of the norm of a given initial displacement vector. Indeed, we compute the growth of the vectors of an orthogonal basis, which amounts to defin-

ing LE, which we denote as *eL <sup>n</sup>*, as the square root of Tr *DM<sup>n</sup>* ð Þ ð Þ **<sup>x</sup>**<sup>0</sup> *<sup>T</sup> DMn*ð Þ **<sup>x</sup>**<sup>0</sup> 

where *M*ð Þ **x** is the map, *DM*ð Þ **x** denotes the tangent map, and **x**<sup>0</sup> is the initial condition. This definition has the obvious advantage of insuring the correct asymptotic growth.

Indeed the anomalies in the behavior of FLI [16], due to the choice of the initial vector, are not met. The use of exponential growth factor of nearby orbits (MEGNO) [17] allows to filter the oscillations which are still present in LE. The RE is obtained from the covariance matrix which is computed from the tangent map. We denote this error by *e<sup>R</sup> <sup>n</sup>* , which has a very simple relation with LE given by the square root of *e<sup>L</sup>* 0 � �<sup>2</sup> <sup>þ</sup> <sup>2</sup> *eL* 1 � �<sup>2</sup> <sup>þ</sup> … <sup>þ</sup> <sup>2</sup> *<sup>e</sup><sup>L</sup> n*�1 � �<sup>2</sup> <sup>þ</sup> *<sup>e</sup><sup>L</sup> n* � �<sup>2</sup> . We first analyze the case of linear maps to explore the behavior of REM. A systematic comparison of LE, RE, and REM is presented for two basic models: the standard map and the Hénon map. The asymptotic power law exponents are computed by using the MEGNO filter. For nonlinear two-dimensional maps, the behavior of the errors has been compared moving along a one-dimensional grid in the phase plane: crossing of islands has a clear signature, the chaotic regions are very neatly distinguished, and good agreement with the theoretical predictions is found.

A rectangular region of phase plane has been examined by choosing a grid and using a color, logarithmic scale for the errors at each point. Also in this case, a good correspondence with the phase space portrait is found. On the basis of the analysis presented here and the experience gained in investigating more complex models from celestial mechanics [18] and beam dynamics [19], we suggest to compare RE and REM with LE, possibly, filtered with MEGNO, to damp the oscillations<sup>2</sup> . For maps of dimension 4 or higher, a direct geometric inspection of the orbits is not possible since the Poincaré section requires an interpolation Hamiltonian. As a consequence the use of fast indicators is the only practical approach to analyze the orbital stability. Hamiltonian systems have a continuous time flow, and the errors LE and RE denoted by *<sup>e</sup><sup>L</sup>*ð Þ*<sup>t</sup>* and *<sup>e</sup><sup>R</sup>*ð Þ*<sup>t</sup>* , respectively, are computed by using the fundamental matrix Lð Þ*<sup>t</sup>* of the tangent flow. In this case *eL*ðÞ¼ *<sup>t</sup>* Tr L*<sup>T</sup>*ð Þ*<sup>t</sup>* <sup>L</sup>ð Þ*<sup>t</sup>* � � � � <sup>1</sup>*=*<sup>2</sup> and *eR*ð Þ*<sup>t</sup>* are given by the square root of 2 <sup>Ð</sup>*<sup>t</sup>* <sup>0</sup> *ds e<sup>L</sup>*ð Þ*<sup>s</sup>* � �<sup>2</sup> , whose trapezoidal rule approximation gives the relation found for the maps [12]. Standard procedures allow to approximate the orbit **<sup>x</sup>**ð Þ*<sup>t</sup>* by the iterates *<sup>M</sup><sup>n</sup>*ð Þ **<sup>x</sup>**<sup>0</sup> of a symplectic integrator map *<sup>M</sup>* (see [20]) and the fundamental matrix Lð Þ*<sup>t</sup>* by *DM<sup>n</sup>*ð Þ **<sup>x</sup>**<sup>0</sup> (see [21]). The paper has the following structure. In Section 2, we recall the definitions for LE and RE and obtain their mutual relation. In Section 3, we present the analytical results on LE and RE together with the numerical results on REM for integrable maps. In Section 4, the key features of two prototype models, the standard map and the Hénon map, are summarized. In Section 5, we present a detailed numerical analysis of LE, RE, and REM for the standard map. In Section 6, the same analysis is presented for the Hénon map. In Section 7, the summary and conclusions are presented.
