## 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

In the framework of the variational methods, we have proposed two indicators [10, 11, 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.^{1}

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

For an integrable map, the growth of LE and RE follows asymptotically a power law

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 defining LE, which we denote as

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

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^{2}. 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

## 2. Definition of errors

Given a symplectic map

which satisfies the linear recurrence

where

and to its variants such as OFLI [4]. The mean exponential growth factor of nearby orbits, MEGNO [17], denoted by

### 2.1 Lyapunov error

We propose a definition of the Lyapunov error which is independent from the choice of the initial vector:

It is immediate to verify that given an orthonormal basis

and obviously the result does not depend on the choice of the basis. The computational cost of

### 2.2 Forward error

When an additive noise of amplitude

where

The global stochastic displacement satisfies a linear nonhomogeneous equation and is defined by

with initial condition

The explicit solution for

where

The expression for the forward error finally reads

The computation cost of

### 2.3 Reversibility error

We have just defined the forward error, but it will not be used, because it is only an intermediate step toward the definition of the reversibility error. Consider the backward evolution

The random backward displacements

with initial condition

For

Letting

and using Eqs. (17) and (8), an explicit expression involving only the tangent maps is obtained. Indeed the global stochastic displacement reads

Taking into account the independence of

### 2.4 Analytical relation between RE and LE indicators

The RE can be obtained from LE in a very simple way. We first notice that

We prove this relation by writing

Starting from

Given any symplectic matrix ^{3}, we can prove that

As a consequence in Eq. (21), we can use the following relation:

Finally, the relation between LE and RE is given by

This relation clearly shows how the error due to random kicks along the orbit is related to the error due to initial orthogonal kicks.

### 2.5 Roundoff-induced reversibility error

The reversibility error method (REM) is a very simple procedure based on

where

### 2.6 Errors for Hamiltonian flows

For Hamiltonian flows, we define the Lyapunov error

where

The reversibility error in this case is defined by the mean square deviation of the random displacement

If the continuous time

## 3. Integrable maps

We evaluate the errors for integrable maps with an elliptic fixed point at the origin, whose normal form is a rotation

### 3.1 Change of coordinate system

In generic coordinates an integrable map

which imply that the orbits

where we used

Taking Eq. (31) into account, the last equation can be written as

Notice that

where

### 3.2 Isochronous rotations: oscillations in LE and RE

If the given map is linear and two-dimensional and

Letting

where

We shall first consider the dependence of the errors on the iteration order

The rotation

For a generic map such as

Letting

and in this case REM vanishes. These results show that REM strongly depends on the computational complexity of the map. The error growth always follows a power law, but, depending on the choice of the coordinates, the exponent

### 3.3 Anisochronous rotations

An integrable map

where ^{4} reads

and the square of the reversibility error is given by

For a fixed value of the invariant *J*, the slope of

In Figure 3, we show the variation with *n* for which the value

If the coordinates are not normal, which is usually the case for a quasi integrable map, the correspondence between RE and REM is better, and it is confirmed by comparing the results for MEGNO. Just a shift of

## 4. Non-integrable maps

We examine here the behavior of the proposed dynamical indicators for two basic models, the standard map and the Hénon map.

The **standard map** is defined as a map on the torus

where

This map is conjugated to a rotation

which is the interpolating Hamiltonian of the map when

As a consequence, for

The **Hénon map** is defined by

Close to the origin, this is just a rotation with frequency

with time step

## 5. The standard map

We have analyzed the errors

When the orbit is chaotic, the growth of all errors is exponential. However, LE and RE can grow until the overflow is reached, whereas REM can grow only up to

### 5.1 Initial conditions on a one-dimensional grid

Figure 7 shows the variation of LE, RE, and REM for *p* for a fixed order *N*. The LE oscillates when the initial condition varies, RE does not oscillate, and REM fluctuates. When the MEGNO filter is applied, LE and RE are equally smooth, whereas REM still fluctuates.

In Figure 8, the same results are shown for a higher value of the parameter

### 5.2 Initial conditions on a two-dimensional domain

We compare here LE, RE, and REM when the initial conditions are chosen in a two-dimensional phase space domain and the iteration number has a fixed value *N*. The most effective way of analyzing the results is to plot the errors using a logarithmic, color scale. Following the conclusions of our previous section, we show LE, RE, and REM, in a logarithmic color scale. We choose a regular two-dimensional grid in a square (or rectangular) domain of phase space with *n*, disappear when MEGNO is computed and are not present in the RE and REM plots. The spurious structures observed in FLI, which depend on the choice of the initial vector, are not present in LE, because in our definition the error does not depend on the choice of an initial displacement vector. Notice that the chosen scales have maximum equal to 10^{10} for LE and 10^{15} for RE and REM. This choice is suggested by the asymptotic behavior

## 6. The Hénon map

We briefly report in this section the numerical results on the errors computed on domains of dimensions 1 and 2 in phase space. Close to the origin, the linear map in this case is a rotation

a formula valid for frequencies

In the normal coordinates

In Figure 11, we show the variation of LE, RE, and REM computed at a fixed order

In Figure 12, we show the color plots of LE, RE, and REM in a square domain of phase space. The weakly chaotic separatrix is detectable in LE and is more clearly visible in RE. The REM plot differs from RE for the up-shit 1/2 of the power law exponent before the thin chaotic separatrix and for the presence of fluctuations.

## 7. Conclusions

We have presented a detailed analysis of the stability indicators LE, RE, and REM recently proposed. Defining the square of LE as the trace of the tangent map times, its transpose renders this indicator independent from the choice of an initial vector, which can introduce spurious structures. The RE is the reversibility error due to additive random noise, whereas REM is the reversibility error due to the roundoff. A very simple relation is found between RE and LE. The oscillations, which affect the fast Lyapunov indicator, can be filtered with MEGNO. Since RE has a smooth behavior and does not exhibit oscillations, filtering it by MEGNO is not necessary. The asymptotic behavior of REM is similar to RE even though it exhibits large fluctuations. The displacements caused by roundoff are almost random vectors, if the map has a high computational complexity, but since just a single realization of the process is available, the fluctuations cannot be averaged.

We have first examined the behavior of LE and RE for linear maps and for integrable maps. If the fixed point is elliptic, then the asymptotic growth follows a power law

For a generic map which has regular and chaotic components, the error growth follows a power law and an exponential law, respectively. For the standard map and the Hénon map, the behavior of LE, RE, and REM has been compared first by varying the iteration order

## Notes

- 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.
- 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 cost is quadratic rather than linear in the iteration order.
- A symplectic matrix L is defined by LJLT=J where J is antisymmetric and J2=−I. As a consequence L−1=−JLTJ, and L−1T=−JLJ so that TrL−1TL−1=TrJLJ2LTJ=TrLLT.
- The standard definition for an initial displacement along the unit vector η0 is enLη0=∥DMnη0∥ where ∥DMnη0∥2=1+nΩ'2∥x∥2η0⋅x2+2nΩ'η0⋅xη0⋅Jx and J=01−10. The sum ‖DMn10‖2+‖DMn01‖2 gives the Lyapunov error enL2.