**2. Definition of errors**

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>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> � � 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, respectively, where *η* is a unit vector. We consider the normalized displacement *η<sup>n</sup>* at iteration *n* defined by

$$\eta\_n = \lim\_{\varepsilon \to 0} \frac{\mathbf{y}\_n - \mathbf{x}\_n}{\varepsilon} = \lim\_{\varepsilon \to 0} \frac{M(\mathbf{y}\_{n-1}) - M(\mathbf{x}\_{n-1})}{\varepsilon} \tag{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 *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

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

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-

Indeed the anomalies in the behavior of FLI [16], due to the choice of the initial

(MEGNO) [17] allows to filter the oscillations which are still present in LE. The RE

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

vector, are not met. The use of exponential growth factor of nearby orbits

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

though with larger fluctuations.

*Progress in Relativity*

ing LE, which we denote as *eL*

totic growth.

**162**

map and therefore on the choice of coordinates.

<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 cost is quadratic rather than linear in the iteration order.

which satisfies the linear recurrence

$$
\eta\_n = DM(\mathbf{x}\_{n-1})\eta\_{n-1} \qquad \qquad (DM)\_{ij} = \frac{\partial M\_i}{\partial \mathbf{x}\_j} \tag{2}
$$

The global stochastic displacement satisfies a linear nonhomogeneous equation

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

*<sup>n</sup>* <sup>¼</sup> h i *<sup>Ξ</sup><sup>n</sup>* � *<sup>Ξ</sup><sup>n</sup>* <sup>1</sup>*=*<sup>2</sup> <sup>¼</sup> Tr <sup>Σ</sup>2*<sup>F</sup>*

*<sup>n</sup>* <sup>¼</sup> *<sup>Ξ</sup><sup>n</sup> <sup>Ξ</sup><sup>T</sup>*

*n*

*n*

<sup>B</sup>*<sup>k</sup> <sup>ξ</sup>n*�*<sup>k</sup>* <sup>B</sup>*<sup>k</sup>* <sup>¼</sup> *DMk*

B*<sup>k</sup>* ¼ B*<sup>k</sup>*�<sup>1</sup>*DM*ð Þ **x***<sup>n</sup>*�*<sup>k</sup>* B0 ¼ I*:* (12)

*Ξ<sup>n</sup>* ¼ *DM*ð Þ **x***<sup>n</sup>*�<sup>1</sup> *Ξn*�<sup>1</sup> þ *ξ<sup>n</sup>* (9)

� � be the covariance matrix

� � � � <sup>1</sup>*=*<sup>2</sup> (10)

ð Þ **x***<sup>n</sup>*�*<sup>k</sup>* (11)

*:* (13)

, with initial point

ð Þ **<sup>x</sup>***<sup>n</sup>*�*k*þ<sup>1</sup> *<sup>Ξ</sup>n,*�*k*þ<sup>1</sup> <sup>þ</sup> *<sup>ξ</sup>*�*<sup>k</sup>* (15)

� �*ξ*�*<sup>j</sup>* (16)

>0. We

� � <sup>þ</sup> <sup>ϵ</sup>*ξ*�*<sup>k</sup> <sup>k</sup>* <sup>¼</sup> <sup>1</sup>*,* …*, n* (14)

*DM*�ð Þ *<sup>k</sup>*�*<sup>j</sup>* **<sup>x</sup>***<sup>n</sup>*�*<sup>j</sup>*

and is defined by

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

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

the forward error is defined by

*<sup>Ξ</sup><sup>n</sup>* <sup>¼</sup> <sup>X</sup>*<sup>n</sup>*

**2.3 Reversibility error**

*<sup>Ξ</sup>n,*�*<sup>k</sup>* <sup>¼</sup> lim<sup>ϵ</sup>!<sup>0</sup>

**165**

*k*¼1

with initial condition *<sup>Ξ</sup>*<sup>0</sup> <sup>¼</sup> 0. Letting <sup>Σ</sup><sup>2</sup> *<sup>F</sup>*

*e F*

The explicit solution for *Ξ<sup>n</sup>* is given by

where B*<sup>k</sup>* can be evaluated recursively as

*DMn*�*<sup>k</sup>*

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

ð Þ **<sup>x</sup>***<sup>k</sup> <sup>ξ</sup><sup>k</sup>* <sup>¼</sup> <sup>X</sup>*n*�<sup>1</sup>

The expression for the forward error finally reads

*e F <sup>n</sup>* <sup>¼</sup> <sup>X</sup>*<sup>n</sup>*�<sup>1</sup>

backward evolution **<sup>y</sup>***n,*�*<sup>k</sup>*, given by the inverse map *<sup>M</sup>*�<sup>1</sup>

random displacement of the same amplitude at each step:

**<sup>y</sup>***n,*�*<sup>k</sup>* <sup>¼</sup> *<sup>M</sup>*�<sup>1</sup> **<sup>y</sup>***n,*�*k*þ<sup>1</sup>

consider then the stochastic process *Ξn,*�*<sup>k</sup>* defined by

*<sup>Ξ</sup>n,* �*<sup>k</sup>* <sup>¼</sup> *DM*�*<sup>k</sup>*

**<sup>y</sup>***n,*�*<sup>k</sup>* � **<sup>x</sup>***<sup>n</sup>*�*<sup>k</sup>* ϵ

*k*¼0

*k*¼0

map, but the storage of the tangent map along the orbit up to *n* is required.

Tr B*<sup>T</sup> <sup>k</sup>* B*<sup>k</sup>* � � !<sup>1</sup>*=*<sup>2</sup>

The computation cost of B*<sup>n</sup>* is negligible, once we have evaluated the tangent

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

**y***n,*<sup>0</sup> ¼ **y***n*. The point **y***<sup>n</sup>* is reached by iterating *n* times the map *M* with a random displacement of amplitude ϵ at each step, starting from **y**<sup>0</sup> ¼ **x**<sup>0</sup> (see the previous subsection). The orbit **<sup>y</sup>***n,*�*<sup>k</sup>* is obtained by iterating *<sup>k</sup>* times the map *<sup>M</sup>*�<sup>1</sup> with a

The random backward displacements *ξ*�*<sup>k</sup>* are independent from the forward

*<sup>Ξ</sup>n,*�*<sup>k</sup>* <sup>¼</sup> *DM*�<sup>1</sup>

*k*

*j*¼1

displacements *<sup>ξ</sup>k*<sup>0</sup>, namely, *<sup>ξ</sup>*�*<sup>k</sup> <sup>ξ</sup><sup>k</sup>* h i0 <sup>¼</sup> 0 and *<sup>ξ</sup>*�*<sup>k</sup> <sup>ξ</sup>*�*<sup>k</sup>* h i0 <sup>¼</sup> <sup>I</sup> *<sup>δ</sup>k k*<sup>0</sup>, for any *k, k*<sup>0</sup>

with initial condition *Ξn,* <sup>0</sup> ¼ *Ξn*. The solution of the recurrence reads

ð Þ **<sup>x</sup>***<sup>n</sup> <sup>Ξ</sup><sup>n</sup>* <sup>þ</sup><sup>X</sup>

where *DM* is the tangent map. For any finite <sup>ϵ</sup>, we have **<sup>y</sup>***<sup>n</sup>* <sup>¼</sup> **<sup>x</sup>***<sup>n</sup>* <sup>þ</sup> <sup>ϵ</sup>*η<sup>n</sup>* <sup>þ</sup> *<sup>O</sup>* <sup>ϵ</sup><sup>2</sup> ð Þ. We might define the error as the norm of *η<sup>n</sup>* which is closely related to the fast Lyapunov indicator FLI (see [3]) as

$$e\_{\pi}(\eta\_0) = \|\eta\_n\| = \|DM^n(\mathbf{x}\_0)\eta\_0\| \qquad \qquad (FLI)\_n = \max\_{1 \le k \le n} \log e\_k(\eta\_0) \tag{3}$$

and to its variants such as OFLI [4]. The mean exponential growth factor of nearby orbits, MEGNO [17], denoted by *Yn* ¼ *Y e*ð Þ*<sup>n</sup>* is the double average of the slope, and we denote it as Δ*e*<sup>2</sup> *<sup>n</sup>*. When *n* is a continuous variable, then Δ*e*<sup>2</sup> *<sup>n</sup>* <sup>¼</sup> *<sup>d</sup>* log *<sup>e</sup>*<sup>2</sup> *<sup>n</sup>=d* log *n*. When *n* is an integer, the standard definition is

$$
\Delta e\_n^2 = n \left( \log e\_n^2 - \log e\_{n-1}^2 \right) \qquad \qquad Y\_n = \left\langle \left\langle \Delta e\_n^2 \right\rangle \right\rangle. \tag{4}
$$

#### **2.1 Lyapunov error**

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

$$\begin{aligned} \mathbf{c}\_{n}^{L} &= \left(\text{Tr}\left(\mathbf{A}\_{n}^{T}\mathbf{A}\_{n}\right)\right)^{1/2} \\ \mathbf{A}\_{n} &= DM^{n}(\mathbf{x}\_{0}) = DM(\mathbf{x}\_{n-1})\mathbf{A}\_{n-1} \qquad &\mathbf{A}\_{0} = \mathbf{I} \end{aligned} \tag{5}$$

It is immediate to verify that given an orthonormal basis *η*0*<sup>k</sup>*, we have

$$e\_n^L = \left(\sum\_{k=1}^{2d} e\_n^2(\eta\_{0k})\right)^{1/2} \tag{6}$$

and obviously the result does not depend on the choice of the basis. The computational cost of *en <sup>η</sup>*<sup>0</sup> ð Þ is 2*<sup>d</sup>* times higher with respect to *<sup>e</sup><sup>L</sup> <sup>n</sup>*, but this difference is negligible with respect to the computational cost of the matrix *DM*ð Þ **x***<sup>n</sup>*�<sup>1</sup> , which recursively gives *η<sup>n</sup>* and A*n*. A similar definition is proposed in the case of Hamiltonian flows (see the last Subsection 2.6 and [12] for more details). An advantage of the proposed definition is that it takes into account the error growth on all possible directions of the initial displacement vector. As a consequence, no spurious effects due to the choice of the initial vector have to be faced (see [16]).

#### **2.2 Forward error**

When an additive noise of amplitude ϵ is introduced, the reference orbit **x***<sup>n</sup>* is replaced by the noisy one (**y***n*) having the same initial condition:

$$\mathbf{y}\_n = M(\mathbf{y}\_{n-1}) + \epsilon \mathbf{f}\_n \qquad \qquad \mathbf{y}\_0 = \mathbf{x}\_0 \tag{7}$$

where *ξ<sup>n</sup>* are independent random vectors satisfying

$$
\langle \mathfrak{g}\_n \rangle = \mathbf{0} \qquad \qquad \left\langle \mathfrak{g}\_n \mathfrak{g}\_m^T \right\rangle = \mathbf{I} \,\delta\_{nm} \tag{8}
$$

*Fast Indicators for Orbital Stability: A Survey on Lyapunov and Reversibility Errors DOI: http://dx.doi.org/10.5772/intechopen.88085*

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

$$\Xi\_n = \lim\_{\varepsilon \to 0} \frac{\mathbf{y}\_n - \mathbf{x}\_n}{\varepsilon} \qquad \qquad \Xi\_n = DM(\mathbf{x}\_{n-1})\Xi\_{n-1} + \mathfrak{g}\_n \tag{9}$$

with initial condition *<sup>Ξ</sup>*<sup>0</sup> <sup>¼</sup> 0. Letting <sup>Σ</sup><sup>2</sup> *<sup>F</sup> <sup>n</sup>* <sup>¼</sup> *<sup>Ξ</sup><sup>n</sup> <sup>Ξ</sup><sup>T</sup> n* � � be the covariance matrix the forward error is defined by

$$
\varepsilon\_n^F = \langle \Xi\_n \cdot \Xi\_n \rangle^{1/2} = \left( \text{Tr} \left( \Sigma\_n^{2F} \right) \right)^{1/2} \tag{10}
$$

The explicit solution for *Ξ<sup>n</sup>* is given by

$$\mathfrak{S}\_n = \sum\_{k=1}^n DM^{n-k}(\mathbf{x}\_k)\mathfrak{f}\_k = \sum\_{k=0}^{n-1} \mathbf{B}\_k \mathfrak{f}\_{n-k} \qquad \qquad \mathbf{B}\_k = DM^k(\mathbf{x}\_{n-k}) \tag{11}$$

where B*<sup>k</sup>* can be evaluated recursively as

$$\mathbf{B}\_{k} = \mathbf{B}\_{k-1} \mathbf{D} \mathbf{M}(\mathbf{x}\_{n-k}) \qquad \qquad \mathbf{B}\_{0} = \mathbf{I}. \tag{12}$$

The expression for the forward error finally reads

$$\boldsymbol{\varepsilon}\_{n}^{F} = \left(\sum\_{k=0}^{n-1} \operatorname{Tr}\left(\mathbf{B}\_{k}^{T}\mathbf{B}\_{k}\right)\right)^{1/2}.\tag{13}$$

The computation cost of B*<sup>n</sup>* is negligible, once we have evaluated the tangent map, but the storage of the tangent map along the orbit up to *n* is required.

#### **2.3 Reversibility error**

which satisfies the linear recurrence

Lyapunov indicator FLI (see [3]) as

slope, and we denote it as Δ*e*<sup>2</sup>

Δ*e* 2

*eL*

*<sup>n</sup>* ¼ *n* log *e*

*<sup>n</sup>* <sup>¼</sup> Tr A*<sup>T</sup>*

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

*<sup>n</sup>* A*<sup>n</sup>* � � � � <sup>1</sup>*=*<sup>2</sup>

> *e L <sup>n</sup>* <sup>¼</sup> <sup>X</sup> 2*d*

putational cost of *en <sup>η</sup>*<sup>0</sup> ð Þ is 2*<sup>d</sup>* times higher with respect to *<sup>e</sup><sup>L</sup>*

due to the choice of the initial vector have to be faced (see [16]).

replaced by the noisy one (**y***n*) having the same initial condition:

*ξ<sup>n</sup>* h i ¼ 0 *ξ<sup>n</sup> ξ*

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

where *ξ<sup>n</sup>* are independent random vectors satisfying

Δ*e*<sup>2</sup>

*<sup>n</sup>* <sup>¼</sup> *<sup>d</sup>* log *<sup>e</sup>*<sup>2</sup>

*Progress in Relativity*

**2.1 Lyapunov error**

**2.2 Forward error**

**164**

choice of the initial vector:

*<sup>η</sup><sup>n</sup>* <sup>¼</sup> *DM*ð Þ **<sup>x</sup>***<sup>n</sup>*�<sup>1</sup> *<sup>η</sup>n*�<sup>1</sup> ð Þ *DM ij* <sup>¼</sup> *<sup>∂</sup>Mi*

where *DM* is the tangent map. For any finite <sup>ϵ</sup>, we have **<sup>y</sup>***<sup>n</sup>* <sup>¼</sup> **<sup>x</sup>***<sup>n</sup>* <sup>þ</sup> <sup>ϵ</sup>*η<sup>n</sup>* <sup>þ</sup> *<sup>O</sup>* <sup>ϵ</sup><sup>2</sup> ð Þ. We might define the error as the norm of *η<sup>n</sup>* which is closely related to the fast

and to its variants such as OFLI [4]. The mean exponential growth factor of nearby orbits, MEGNO [17], denoted by *Yn* ¼ *Y e*ð Þ*<sup>n</sup>* is the double average of the

*<sup>n</sup>=d* log *n*. When *n* is an integer, the standard definition is

We propose a definition of the Lyapunov error which is independent from the

<sup>A</sup>*<sup>n</sup>* <sup>¼</sup> *DM<sup>n</sup>*ð Þ¼ **<sup>x</sup>**<sup>0</sup> *DM*ð Þ **<sup>x</sup>***<sup>n</sup>*�<sup>1</sup> <sup>A</sup>*<sup>n</sup>*�<sup>1</sup> A0 <sup>¼</sup> <sup>I</sup>

It is immediate to verify that given an orthonormal basis *η*0*<sup>k</sup>*, we have

*k*¼1 *e* 2 *<sup>n</sup> η*0*<sup>k</sup>* ð Þ !<sup>1</sup>*=*<sup>2</sup>

and obviously the result does not depend on the choice of the basis. The com-

When an additive noise of amplitude ϵ is introduced, the reference orbit **x***<sup>n</sup>* is

*T m*

� � <sup>þ</sup> <sup>ϵ</sup>*ξ<sup>n</sup>* **<sup>y</sup>**<sup>0</sup> <sup>¼</sup> **<sup>x</sup>**<sup>0</sup> (7)

� � <sup>¼</sup> <sup>I</sup> *<sup>δ</sup>nm* (8)

negligible with respect to the computational cost of the matrix *DM*ð Þ **x***<sup>n</sup>*�<sup>1</sup> , which recursively gives *η<sup>n</sup>* and A*n*. A similar definition is proposed in the case of Hamiltonian flows (see the last Subsection 2.6 and [12] for more details). An advantage of the proposed definition is that it takes into account the error growth on all possible directions of the initial displacement vector. As a consequence, no spurious effects

2 *n*�1 � � *Yn* <sup>¼</sup> <sup>Δ</sup>*<sup>e</sup>*

*<sup>n</sup>*. When *n* is a continuous variable, then

*en <sup>η</sup>*<sup>0</sup> ð Þ¼ <sup>∥</sup>*ηn*<sup>∥</sup> <sup>¼</sup> <sup>∥</sup>*DMn*ð Þ **<sup>x</sup>**<sup>0</sup> *<sup>η</sup>*0<sup>∥</sup> ð Þ *FLI <sup>n</sup>* <sup>¼</sup> max

*∂xj*

<sup>1</sup>≤*k*≤*<sup>n</sup>* log *ek <sup>η</sup>*<sup>0</sup> ð Þ (3)

� � � � *:* (4)

2 *n* (2)

(5)

(6)

*<sup>n</sup>*, but this difference is

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 **<sup>y</sup>***n,*�*<sup>k</sup>*, given by the inverse map *<sup>M</sup>*�<sup>1</sup> , with initial point **y***n,*<sup>0</sup> ¼ **y***n*. The point **y***<sup>n</sup>* is reached by iterating *n* times the map *M* with a random displacement of amplitude ϵ at each step, starting from **y**<sup>0</sup> ¼ **x**<sup>0</sup> (see the previous subsection). The orbit **<sup>y</sup>***n,*�*<sup>k</sup>* is obtained by iterating *<sup>k</sup>* times the map *<sup>M</sup>*�<sup>1</sup> with a random displacement of the same amplitude at each step:

$$\mathbf{y}\_{n,-k} = M^{-1} \begin{pmatrix} \mathbf{y}\_{n,-k+1} \end{pmatrix} + \epsilon \mathbf{f}\_{-k} \qquad \qquad \qquad k = 1, \ldots, n \tag{14}$$

The random backward displacements *ξ*�*<sup>k</sup>* are independent from the forward displacements *<sup>ξ</sup>k*<sup>0</sup>, namely, *<sup>ξ</sup>*�*<sup>k</sup> <sup>ξ</sup><sup>k</sup>* h i0 <sup>¼</sup> 0 and *<sup>ξ</sup>*�*<sup>k</sup> <sup>ξ</sup>*�*<sup>k</sup>* h i0 <sup>¼</sup> <sup>I</sup> *<sup>δ</sup>k k*<sup>0</sup>, for any *k, k*<sup>0</sup> >0. We consider then the stochastic process *Ξn,*�*<sup>k</sup>* defined by

$$\Xi\_{n,-k} = \lim\_{\epsilon \to 0} \frac{\mathbf{y}\_{n,-k} - \mathbf{x}\_{n-k}}{\epsilon} \qquad \qquad \Xi\_{n,-k} = DM^{-1}(\mathbf{x}\_{n-k+1})\Xi\_{n,-k+1} + \mathbf{f}\_{-k} \tag{15}$$

with initial condition *Ξn,* <sup>0</sup> ¼ *Ξn*. The solution of the recurrence reads

$$\Xi\_{n,\ \ \ \lambda} = DM^{-k}(\mathbf{x}\_n)\,\Xi\_n + \sum\_{j=1}^k DM^{-(k-j)}\left(\mathbf{x}\_{n-j}\right)\xi\_{-j} \tag{16}$$

For *<sup>k</sup>* <sup>¼</sup> *<sup>n</sup>*, we obtain the global normalized displacement *<sup>Ξ</sup><sup>R</sup> <sup>n</sup>* ¼ *Ξn,*�*<sup>n</sup>* after *n* forward and *n* backward iterations with noise of vanishing amplitude:

$$\boldsymbol{\Xi}\_{n}^{R} = DM^{-n}(\mathbf{x}\_{n})\,\boldsymbol{\Xi}\_{n} + \sum\_{j=1}^{n} DM^{-(n-j)} \left(\mathbf{x}\_{n-j}\right)\,\boldsymbol{\xi}\_{-j} \tag{17}$$

Tr L�<sup>1</sup> � �*<sup>T</sup>* <sup>L</sup>�<sup>1</sup> � � <sup>¼</sup> Tr L*T*<sup>L</sup> � �*:* (23)

� � h i <sup>¼</sup> *<sup>e</sup>*

!

ð Þ **x**<sup>0</sup>

*L k* � �<sup>2</sup>

*:* (24)

*:* (25)

ð Þ **x**<sup>0</sup> � �*<sup>T</sup> DMk*

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

*k*¼1

This relation clearly shows how the error due to random kicks along the orbit is

The reversibility error method (REM) is a very simple procedure based on *n* iterations of the map *M* followed by *n* iterations of the inverse map. The distance from the initial point normalized by the roundoff amplitude ϵ defines the REM

<sup>ϵ</sup> ∘ *M<sup>n</sup>*

ϵ

where ϵ is the roundoff amplitude. For the eight-byte representation of reals, we choose <sup>ϵ</sup> <sup>¼</sup> <sup>10</sup>�17. If the map has a sufficiently high computational complexity, the displacement *ξ* defined by *M*ϵð Þ� **x** *M*ð Þ¼ **x** ϵ*ξ* is almost random, but a unique realization is available. (For a discussion on the roundoff error, see [22]). As a

<sup>ϵ</sup> ð Þ� **x**<sup>0</sup> **x**0∥

� � (26)

*eL*ðÞ¼ *<sup>t</sup>* Tr L*<sup>T</sup>*ð Þ*<sup>t</sup>* <sup>L</sup>ð Þ*<sup>t</sup>* � � � � <sup>1</sup>*=*<sup>2</sup> (27)

ð Þ*s* ð Þ *ξ*ðÞ�*s ξ*ð Þ *s* � *t ds:* (28)

*<sup>H</sup>=∂xi∂xj*, and the initial condition for the matrix Lð Þ*<sup>t</sup>* is

*<sup>n</sup>* has a smooth dependence

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

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

ð Þ **x***<sup>k</sup>* � � h i <sup>¼</sup> Tr *DM<sup>k</sup>*

Finally, the relation between LE and RE is given by

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

related to the error due to initial orthogonal kicks.

error. Denoting with *M*<sup>ϵ</sup> the map with roundoff, we have

*<sup>n</sup>* <sup>¼</sup> <sup>∥</sup>*M*�*<sup>n</sup>*

*<sup>n</sup>* has large fluctuations, whereas *e<sup>R</sup>*

on *n* since it is defined by an average overall possible realizations of the stochastic

For Hamiltonian flows, we define the Lyapunov error *eL*ð Þ*t* according to

L 0ð Þ¼ I. The relation with the standard fast indicators is the same as for the symplectic maps. Let *ΞR*ð Þ*t* be the stochastic displacement from **x**<sup>0</sup> after a noisy evolution up to time *t* and backward to *t* ¼ 0, divided by the noise amplitude *ε* in the limit *ε* ! 0. It has been proven [12] that *ΞR*ð Þ*t* satisfies a linear Langevin

*<sup>Ξ</sup><sup>R</sup>*ðÞ¼ *<sup>t</sup>*

ð*t* 0 L�<sup>1</sup>

where Lð Þ*t* is the fundamental matrix for the tangent flow, which satisfies the linear equation *d*L*=dt* ¼ JHL, H denotes the Hessian of the Hamiltonian computed

*e REM*

Tr *DM*�*<sup>k</sup>*

*e R n* � �<sup>2</sup> <sup>¼</sup> <sup>X</sup>*<sup>n</sup>*

ð Þ **x***<sup>k</sup>* � �*<sup>T</sup> DM*�*<sup>k</sup>*

*k*¼1

consequence, the *eREM*

along the orbit H*ij* <sup>¼</sup> *<sup>∂</sup>*<sup>2</sup>

**167**

equation whose solution reads

*e L k* � �<sup>2</sup>

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

**2.5 Roundoff-induced reversibility error**

displacements occurring at each iteration.

**2.6 Errors for Hamiltonian flows**

Letting Σ2*<sup>R</sup> <sup>n</sup>* <sup>¼</sup> *<sup>Ξ</sup><sup>R</sup> <sup>n</sup> Ξ<sup>R</sup> n* � �*<sup>T</sup>* D E be the covariance matrix, the reversibility error (RE) is defined by

$$\left\langle e\_n^R = \left\langle \boldsymbol{\Xi}\_n^R \cdot \boldsymbol{\Xi}\_n^R \right\rangle^{1/2} = \left( \operatorname{Tr} \boldsymbol{\Sigma}\_n^{2R} \right)^2 \tag{18}$$

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

$$\mathfrak{S}\_n^R = \sum\_{k=1}^n DM^{-n}(\mathbf{x}\_n) DM^{n-k}(\mathbf{x}\_k) \mathfrak{f}\_k + \sum\_{k=1}^n DM^{-(n-k)}(\mathbf{x}\_{n-k}) \mathfrak{f}\_{-k}.\tag{19}$$

Taking into account the independence of *<sup>ξ</sup><sup>k</sup>* and *<sup>ξ</sup>*�*k*<sup>0</sup>, the expression for the reversibility error *e<sup>R</sup> <sup>n</sup>* <sup>¼</sup> *<sup>Ξ</sup><sup>R</sup> <sup>n</sup>* � *<sup>Ξ</sup><sup>R</sup> n* � �<sup>1</sup>*=*<sup>2</sup> is immediately obtained.

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

$$DM^{-n}(\mathbf{x}\_{\mathfrak{n}}) DM^{n-k}(\mathbf{x}\_{k}) = DM^{-k}(\mathbf{x}\_{k}) \tag{20}$$

We prove this relation by writing *<sup>M</sup>*�*<sup>n</sup> <sup>M</sup><sup>n</sup>*�*<sup>k</sup>*ð Þ **<sup>x</sup>** � � <sup>¼</sup> *<sup>M</sup>*�*<sup>k</sup>*ð Þ **<sup>x</sup>** , computing the tangent map *DM*�*<sup>n</sup> <sup>M</sup><sup>n</sup>*�*<sup>k</sup>*ð Þ **<sup>x</sup>** � �*DM<sup>n</sup>*�*<sup>k</sup>*ð Þ¼ **<sup>x</sup>** *DM*�*<sup>k</sup>*ð Þ **<sup>x</sup>** , and evaluating it for **<sup>x</sup>** <sup>¼</sup> **<sup>x</sup>***k*. As a consequence the expression for the reversibility error becomes

$$\begin{split} \left(\boldsymbol{e}\_{n}^{R}\right)^{2} &= \operatorname{Tr}\left\langle \boldsymbol{\Xi}\_{n}^{R} \left(\boldsymbol{\Xi}\_{n}^{R}\right)^{T} \right\rangle = \sum\_{k=1}^{n} \left( \operatorname{Tr}\left[ \left(\boldsymbol{\mathcal{D}}\boldsymbol{M}^{-k}(\mathbf{x}\_{k})\right)^{T} \left(\boldsymbol{\mathcal{D}}\boldsymbol{M}^{-k}(\mathbf{x}\_{k})\right) \right] + \\ &\quad + \operatorname{Tr}\left[ \left(\boldsymbol{\mathcal{D}}\boldsymbol{M}^{-(n-k)}(\mathbf{x}\_{n-k})\right)^{T} \boldsymbol{\mathcal{D}}\boldsymbol{M}^{-(n-k)}(\mathbf{x}\_{n-k}) \right] \right) = \\ &= 2\sum\_{k=1}^{n-1} \operatorname{Tr}\left[ \left(\boldsymbol{\mathcal{D}}\boldsymbol{M}^{-k}(\mathbf{x}\_{k})\right)^{T} \left(\boldsymbol{\mathcal{D}}\boldsymbol{M}^{-k}(\mathbf{x}\_{k})\right) \right] + \\ &\quad + \operatorname{Tr}\left[ \left(\boldsymbol{\mathcal{D}}\boldsymbol{M}^{-n}(\mathbf{x}\_{n})\right)^{T} \boldsymbol{\mathcal{D}}\boldsymbol{M}^{-n}(\mathbf{x}\_{n}) \right] + \operatorname{Tr}\left(\mathbf{I} \right) \end{split} \tag{21}$$

Starting from *M*�*<sup>k</sup> M<sup>k</sup>* ð Þ **<sup>x</sup>** � � <sup>¼</sup> **<sup>x</sup>**, computing the tangent map, and evaluating it at **x** ¼ **x**0, it follows that

$$DM^{-k}(\mathbf{x}\_k) = \left( DM^k(\mathbf{x}\_0) \right)^{-1} \tag{22}$$

Given any symplectic matrix L<sup>3</sup> , we can prove that

<sup>3</sup> A symplectic matrix L is defined by LJL*<sup>T</sup>* <sup>¼</sup> J where J is antisymmetric and J2 ¼ �I. As a consequence <sup>L</sup>�<sup>1</sup> ¼ �JL*<sup>T</sup>*J, and L�<sup>1</sup> � �*<sup>T</sup>* ¼ �JLJ so that Tr L�1*<sup>T</sup>* <sup>L</sup>�<sup>1</sup> � � <sup>¼</sup> Tr JLJ2 <sup>L</sup>*<sup>T</sup>*<sup>J</sup> � � <sup>¼</sup> Tr LL*<sup>T</sup>* � �.

*Fast Indicators for Orbital Stability: A Survey on Lyapunov and Reversibility Errors DOI: http://dx.doi.org/10.5772/intechopen.88085*

$$\operatorname{Tr}\left(\left(\mathbf{L}^{-1}\right)^{T}\mathbf{L}^{-1}\right) = \operatorname{Tr}\left(\mathbf{L}^{T}\mathbf{L}\right).\tag{23}$$

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

$$\operatorname{Tr}\left[\left(\boldsymbol{D}\boldsymbol{M}^{-k}(\mathbf{x}\_{k})\right)^{T}\left(\boldsymbol{D}\boldsymbol{M}^{-k}(\mathbf{x}\_{k})\right)\right] = \operatorname{Tr}\left[\left(\boldsymbol{D}\boldsymbol{M}^{k}(\mathbf{x}\_{0})\right)^{T}\left(\boldsymbol{D}\boldsymbol{M}^{k}(\mathbf{x}\_{0})\right)\right] = \left(\boldsymbol{\varepsilon}\_{k}^{L}\right)^{2}.\tag{24}$$

Finally, the relation between LE and RE is given by

$$\left(\left(e\_n^R\right)^2 = \sum\_{k=1}^n \left(\left(e\_k^L\right)^2 + \left(e\_{n-k}^L\right)^2\right) = 2\left(\sum\_{k=1}^{n-1} \left(e\_k^L\right)^2 + \frac{1}{2}\left(e\_0^L\right)^2 + \frac{1}{2}\left(e\_n^L\right)^2\right). \tag{25}$$

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

For *<sup>k</sup>* <sup>¼</sup> *<sup>n</sup>*, we obtain the global normalized displacement *<sup>Ξ</sup><sup>R</sup>*

*<sup>n</sup>* <sup>¼</sup> *DM*�*n*ð Þ **<sup>x</sup>***<sup>n</sup> <sup>Ξ</sup><sup>n</sup>* <sup>þ</sup>X*<sup>n</sup>*

maps is obtained. Indeed the global stochastic displacement reads

*Ξ<sup>R</sup>*

*<sup>n</sup> Ξ<sup>R</sup> n* � �*<sup>T</sup>* D E

> *e R <sup>n</sup>* <sup>¼</sup> *<sup>Ξ</sup><sup>R</sup>*

*DM*�*<sup>n</sup>*ð Þ **<sup>x</sup>***<sup>n</sup> DMn*�*<sup>k</sup>*

*<sup>n</sup>* � *<sup>Ξ</sup><sup>R</sup> n*

**2.4 Analytical relation between RE and LE indicators**

*DM*�*<sup>n</sup>*ð Þ **<sup>x</sup>***<sup>n</sup> DMn*�*<sup>k</sup>*

As a consequence the expression for the reversibility error becomes

<sup>¼</sup> <sup>X</sup>*<sup>n</sup> k*¼1

ð Þ **x***<sup>k</sup>* � �*<sup>T</sup> DM*�*<sup>k</sup>*ð Þ **<sup>x</sup>***<sup>k</sup>* � � h i

�

� ��

ð Þ¼ **<sup>x</sup>***<sup>k</sup> DM<sup>k</sup>*

<sup>3</sup> A symplectic matrix L is defined by LJL*<sup>T</sup>* <sup>¼</sup> J where J is antisymmetric and J2 ¼ �I. As a consequence

, we can prove that

*<sup>n</sup>* <sup>¼</sup> *<sup>Ξ</sup><sup>R</sup>*

*<sup>n</sup> Ξ<sup>R</sup> n* � �*<sup>T</sup>* D E

<sup>þ</sup> Tr *DM*�ð Þ *<sup>n</sup>*�*<sup>k</sup>* ð Þ **<sup>x</sup>***<sup>n</sup>*�*<sup>k</sup>* � �*<sup>T</sup>*

Tr *DM*�*<sup>k</sup>*

<sup>þ</sup> Tr *DM*�*<sup>n</sup>* ð Þ ð Þ **<sup>x</sup>***<sup>n</sup> <sup>T</sup> DM*�*<sup>n</sup>*ð Þ **<sup>x</sup>***<sup>n</sup>* h i

*DM*�*<sup>k</sup>*

*<sup>n</sup>* <sup>¼</sup> *<sup>Ξ</sup><sup>R</sup>*

Letting Σ2*<sup>R</sup>*

*Progress in Relativity*

*Ξ<sup>R</sup> <sup>n</sup>* <sup>¼</sup> <sup>X</sup>*<sup>n</sup> k*¼1

reversibility error *e<sup>R</sup>*

*e R n*

� �<sup>2</sup> <sup>¼</sup> Tr *<sup>Ξ</sup><sup>R</sup>*

¼ 2 X*n*�1 *k*¼1

Starting from *M*�*<sup>k</sup> M<sup>k</sup>*

Given any symplectic matrix L<sup>3</sup>

<sup>L</sup>�<sup>1</sup> ¼ �JL*<sup>T</sup>*J, and L�<sup>1</sup> � �*<sup>T</sup>* ¼ �JLJ so that Tr L�1*<sup>T</sup>* <sup>L</sup>�<sup>1</sup> � � <sup>¼</sup> Tr JLJ2

**x** ¼ **x**0, it follows that

**166**

is defined by

forward and *n* backward iterations with noise of vanishing amplitude:

*<sup>n</sup>* � *<sup>Ξ</sup><sup>R</sup> n* � �1*=*<sup>2</sup> <sup>¼</sup> TrΣ2*<sup>R</sup>*

*j*¼1

and using Eqs. (17) and (8), an explicit expression involving only the tangent

ð Þ **<sup>x</sup>***<sup>k</sup> <sup>ξ</sup><sup>k</sup>* <sup>þ</sup>X*<sup>n</sup>*

Taking into account the independence of *<sup>ξ</sup><sup>k</sup>* and *<sup>ξ</sup>*�*k*<sup>0</sup>, the expression for the

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

We prove this relation by writing *<sup>M</sup>*�*<sup>n</sup> <sup>M</sup><sup>n</sup>*�*<sup>k</sup>*ð Þ **<sup>x</sup>** � � <sup>¼</sup> *<sup>M</sup>*�*<sup>k</sup>*ð Þ **<sup>x</sup>** , computing the tangent map *DM*�*<sup>n</sup> <sup>M</sup><sup>n</sup>*�*<sup>k</sup>*ð Þ **<sup>x</sup>** � �*DM<sup>n</sup>*�*<sup>k</sup>*ð Þ¼ **<sup>x</sup>** *DM*�*<sup>k</sup>*ð Þ **<sup>x</sup>** , and evaluating it for **<sup>x</sup>** <sup>¼</sup> **<sup>x</sup>***k*.

Tr *DM*�*<sup>k</sup>*

� �<sup>1</sup>*=*<sup>2</sup> is immediately obtained.

*k*¼1

ð Þ¼ **<sup>x</sup>***<sup>k</sup> DM*�*<sup>k</sup>*

ð Þ **x***<sup>k</sup>* � �*<sup>T</sup> DM*�*<sup>k</sup>*

*DM*�ð Þ *<sup>n</sup>*�*<sup>k</sup>* ð Þ **<sup>x</sup>***<sup>n</sup>*�*<sup>k</sup>*

þ Tr Ið Þ

ð Þ **x**<sup>0</sup>

� � h i

þ

ð Þ **<sup>x</sup>** � � <sup>¼</sup> **<sup>x</sup>**, computing the tangent map, and evaluating it at

*DM*�ð Þ *<sup>n</sup>*�*<sup>j</sup>* **<sup>x</sup>***<sup>n</sup>*�*<sup>j</sup>*

*n*

be the covariance matrix, the reversibility error (RE)

*<sup>n</sup>* ¼ *Ξn,*�*<sup>n</sup>* after *n*

� �*ξ*�*<sup>j</sup>* (17)

� �<sup>2</sup> (18)

*DM*�ð Þ *<sup>n</sup>*�*<sup>k</sup>* ð Þ **<sup>x</sup>***<sup>n</sup>*�*<sup>k</sup> <sup>ξ</sup>*�*<sup>k</sup>:* (19)

ð Þ **x***<sup>k</sup>* (20)

ð Þ **x***<sup>k</sup>*

¼

� ��<sup>1</sup> (22)

<sup>L</sup>*<sup>T</sup>*<sup>J</sup> � � <sup>¼</sup> Tr LL*<sup>T</sup>* � �.

þ

(21)

The reversibility error method (REM) is a very simple procedure based on *n* iterations of the map *M* followed by *n* iterations of the inverse map. The distance from the initial point normalized by the roundoff amplitude ϵ defines the REM error. Denoting with *M*<sup>ϵ</sup> the map with roundoff, we have

$$\sigma\_n^{\rm REM} = \left(\frac{||\mathbf{M}\_e^{-n} \circ M\_e^n(\mathbf{x}\_0) - \mathbf{x}\_0||}{\epsilon}\right) \tag{26}$$

where ϵ is the roundoff amplitude. For the eight-byte representation of reals, we choose <sup>ϵ</sup> <sup>¼</sup> <sup>10</sup>�17. If the map has a sufficiently high computational complexity, the displacement *ξ* defined by *M*ϵð Þ� **x** *M*ð Þ¼ **x** ϵ*ξ* is almost random, but a unique realization is available. (For a discussion on the roundoff error, see [22]). As a consequence, the *eREM <sup>n</sup>* has large fluctuations, whereas *e<sup>R</sup> <sup>n</sup>* has a smooth dependence on *n* since it is defined by an average overall possible realizations of the stochastic displacements occurring at each iteration.

#### **2.6 Errors for Hamiltonian flows**

For Hamiltonian flows, we define the Lyapunov error *eL*ð Þ*t* according to

$$e\_L(t) = \left(\text{Tr}\left(\mathbf{L}^T(t)\mathbf{L}(t)\right)\right)^{1/2} \tag{27}$$

where Lð Þ*t* is the fundamental matrix for the tangent flow, which satisfies the linear equation *d*L*=dt* ¼ JHL, H denotes the Hessian of the Hamiltonian computed along the orbit H*ij* <sup>¼</sup> *<sup>∂</sup>*<sup>2</sup> *<sup>H</sup>=∂xi∂xj*, and the initial condition for the matrix Lð Þ*<sup>t</sup>* is L 0ð Þ¼ I. The relation with the standard fast indicators is the same as for the symplectic maps. Let *ΞR*ð Þ*t* be the stochastic displacement from **x**<sup>0</sup> after a noisy evolution up to time *t* and backward to *t* ¼ 0, divided by the noise amplitude *ε* in the limit *ε* ! 0. It has been proven [12] that *ΞR*ð Þ*t* satisfies a linear Langevin equation whose solution reads

$$\boldsymbol{\Xi}^{\mathbb{R}}(t) = \int\_{0}^{t} \boldsymbol{\mathcal{L}}^{-1}(s) \left( \boldsymbol{\xi}(s) - \boldsymbol{\xi}(s - t) \right) ds. \tag{28}$$

The reversibility error in this case is defined by the mean square deviation of the random displacement *<sup>e</sup>R*ðÞ¼ *<sup>t</sup> <sup>Ξ</sup>R*ðÞ� *<sup>t</sup> <sup>Ξ</sup>R*ð Þ*<sup>t</sup>* � �1*=*<sup>2</sup> . As shown in [12] and from Eq. (28), we immediately obtain

$$e^{R}(t) = \left(2\int\_{0}^{t} \left(e^{L}(s)\right)^{2} ds\right)^{1/2}.\tag{29}$$

Notice that V is a positive-defined matrix and that its determinant is equal to 1 if

*b c* � � <sup>V</sup>�<sup>1</sup> <sup>¼</sup> *<sup>c</sup>* �*<sup>b</sup>*

If the given map is linear and two-dimensional and *M*ð Þ¼ **x** L**x** with ∣Tr L∣ < 2,

Letting **<sup>x</sup>***<sup>n</sup>* <sup>¼</sup> <sup>L</sup>*<sup>n</sup>* **<sup>x</sup>**<sup>0</sup> and **<sup>X</sup>***<sup>n</sup>* <sup>¼</sup> <sup>R</sup>*<sup>n</sup>***X**0, the orbits in the coordinate **<sup>x</sup>** and the normal

� �<sup>2</sup> <sup>¼</sup> Tr V�1Rð Þ �*n<sup>ω</sup>* VRð Þ *<sup>n</sup><sup>ω</sup>* � � <sup>¼</sup> 2 cos <sup>2</sup>ð Þþ *<sup>n</sup><sup>ω</sup> <sup>a</sup>*<sup>2</sup> <sup>þ</sup> *<sup>c</sup>*<sup>2</sup> <sup>þ</sup> <sup>2</sup>*b*<sup>2</sup> � � sin <sup>2</sup>ð Þ¼ *<sup>n</sup><sup>ω</sup>*

<sup>¼</sup> ð Þ *<sup>a</sup>* <sup>þ</sup> *<sup>c</sup>*

where V <sup>¼</sup> <sup>T</sup>*<sup>T</sup>* T, and we have used the representation given by Eq. (34) where *D*Φ ¼ T. The error is constant in normal coordinate **X** and oscillates between 2 and

<sup>2</sup> � <sup>2</sup> <sup>¼</sup> *<sup>a</sup>*<sup>2</sup> <sup>þ</sup> *<sup>c</sup>*<sup>2</sup> <sup>þ</sup> <sup>2</sup>*b*<sup>2</sup> in the coordinate **<sup>x</sup>**. The geometric interpretation is obvious since the orbits of the map belong to an ellipse rather than a circle. The

*f n*ð Þ

We shall first consider the dependence of the errors on the iteration order *n* from *n* ¼ 1 up to a maximum value *N*. Then, we shall consider the dependence on the initial condition **x**<sup>0</sup> when it is varied on a one-dimensional grid crossing the origin for the value *N* of the iteration number. We choose the linear map L which depends

The rotation Rð Þ *ω* **x** is the linear part for the Hénon map, whereas L**x** is the linear part of the standard map that will be discussed in the next sections. The behavior of

<sup>ð</sup> cos 2ð Þþ *<sup>k</sup><sup>ω</sup>* cos 2ð Þ ð Þ *<sup>n</sup>* � *<sup>k</sup> <sup>ω</sup>* Þ ¼ cos 2ð Þþ *<sup>n</sup><sup>ω</sup>* cos 2ð Þ� ð Þ *<sup>n</sup>* � <sup>1</sup> *<sup>ω</sup>* cos 2ð Þ *<sup>n</sup><sup>ω</sup>*

2 ¼

ffiffi *λ* p

2

2

on a single parameter *λ*, and its relation with the rotation frequency is

� � sin *<sup>ω</sup>*

!

2

<sup>2</sup> <sup>þ</sup> <sup>2</sup> � ð Þ *<sup>a</sup>* <sup>þ</sup> *<sup>c</sup>*

<sup>L</sup> <sup>¼</sup> TRð Þ *<sup>ω</sup>* <sup>T</sup>�<sup>1</sup> <sup>R</sup>ð Þ¼ *<sup>ω</sup>* cosð Þ *<sup>ω</sup>* sin*<sup>ω</sup>*

�*b a* � � (34)

� sin ð Þ *ω* cosð Þ *ω*

� � (35)

2

cos 2ð Þ *nω*

1 � cos 2ð Þ *ω*

<sup>2</sup> <sup>0</sup>≤*λ*≤4 (38)

(36)

(37)

2

!

Φ is symplectic. For a two-dimensional map, we can write

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

<sup>V</sup> � *a b*

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

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

coordinate **<sup>X</sup>** <sup>¼</sup> <sup>T</sup>�<sup>1</sup> **<sup>x</sup>** and the Lyapunov errors are given by

then the map is conjugated to a rotation Rð Þ *ω* :

� �<sup>2</sup> <sup>¼</sup> Tr R<sup>½</sup> ð Þ �*n<sup>ω</sup>* <sup>R</sup>ð Þ *<sup>n</sup><sup>ω</sup>* � ¼ <sup>2</sup>

result for the reversibility error is given by

2

<sup>L</sup> <sup>¼</sup> <sup>1</sup> � *<sup>λ</sup>* <sup>1</sup> �*λ* 1

<sup>2</sup> <sup>þ</sup> <sup>2</sup> � ð Þ *<sup>a</sup>* <sup>þ</sup> *<sup>c</sup>*

where *ac* � *<sup>b</sup>*<sup>2</sup> <sup>¼</sup> 1.

*eL <sup>n</sup>*ð Þ **X**<sup>0</sup>

*eL <sup>n</sup>*ð Þ **x**<sup>0</sup>

ð Þ *a* þ *c*

*eR <sup>n</sup>* ð Þ **X**<sup>0</sup> � �<sup>2</sup> <sup>¼</sup> <sup>4</sup>*<sup>n</sup>*

*eR <sup>n</sup>* ð Þ **x**<sup>0</sup>

**169**

*f n*ð Þ¼ <sup>X</sup>*<sup>n</sup>*

*k*¼1

� �<sup>2</sup> <sup>¼</sup> <sup>2</sup>*<sup>n</sup>* ð Þ *<sup>a</sup>* <sup>þ</sup> *<sup>c</sup>*

If the continuous time *t* is replaced by an integer *n* and we approximate the integral with the trapezoidal rule, we recover the relation in Eq. (25) obtained for a symplectic map.

## **3. Integrable maps**

We evaluate the errors for integrable maps with an elliptic fixed point at the origin, whose normal form is a rotation Rð Þ Ω with a frequency Ω depending on the distance from the origin. The LE asymptotic growth is linear, and oscillations are present unless the coordinates are normal. The RE asymptotic growth follows a power law *<sup>n</sup><sup>α</sup>* with exponent *<sup>α</sup>* <sup>¼</sup> <sup>3</sup>*=*2. If the map is linear, its asymptotic growth follows a power law with *α* ¼ 0 for LE and *α* ¼ 1*=*2 for RE. The oscillations reflect the loss of rotational symmetry when generic coordinates are used. The roundoff induced reversibility error REM is also sensitive to the choice of coordinates, and a comparison between RE and REM is presented in the next sections.

#### **3.1 Change of coordinate system**

In generic coordinates an integrable map *M*ð Þ **x** is conjugated to its normal form *N*ð Þ **X** by a symplectic coordinate transformation **x** ¼ Φð Þ **X** ; as a consequence the conjugation equation and its iterates read

$$M(\mathbf{x}) = \Phi \circ N \circ \Phi^{-1}(\mathbf{x}) \qquad \qquad M^\pi(\mathbf{x}) = \Phi \circ N^\pi \circ \Phi^{-1}(\mathbf{x}) \tag{30}$$

which imply that the orbits **<sup>x</sup>***<sup>n</sup>* <sup>¼</sup> *<sup>M</sup><sup>n</sup>*ð Þ **<sup>x</sup>**<sup>0</sup> and **<sup>X</sup>***<sup>n</sup>* <sup>¼</sup> *<sup>N</sup><sup>n</sup>*ð Þ **<sup>X</sup>**<sup>0</sup> are related by **x***<sup>n</sup>* ¼ Φð Þ **X***<sup>n</sup>* . The tangent maps are given by

$$D\mathbf{M}^{n}(\mathbf{x}\_{0}) = D\Phi(\mathbf{X}\_{\pi})D\mathbf{N}^{n}(\mathbf{X}\_{0})D\Phi^{-1}(\mathbf{x}\_{0}) = D\Phi(\mathbf{X}\_{\pi})D\mathbf{N}^{n}(\mathbf{X}\_{0})(D\Phi(\mathbf{X}\_{0}))^{-1} \tag{31}$$

where we used *<sup>D</sup>*Φð Þ **<sup>X</sup>** *<sup>D</sup>*Φ�<sup>1</sup> ð Þ¼ **x** I, a relation which is proved to hold by computing the Jacobian of Φ ∘ Φ�<sup>1</sup> ð Þ¼ **x x**. As a consequence, the expression for the Lyapunov error in both coordinate systems is

$$\begin{aligned} \left(e\_{n}^{L}(\mathbf{X}\_{0})\right)^{2} &= \operatorname{Tr}\left[\left(DN^{n}(\mathbf{X}\_{0})\right)^{T}\left(DN^{n}(\mathbf{X}\_{0})\right)\right] \\ \left(e\_{n}^{L}(\mathbf{x}\_{0})\right)^{2} &= \operatorname{Tr}\left[\left(DN^{n}(\mathbf{x}\_{0})\right)^{T}\left(DN^{n}(\mathbf{x}\_{0})\right)\right] \end{aligned} \tag{32}$$

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

$$\begin{aligned} \left(e\_n^L(\mathbf{x}\_0)\right)^2 &= \text{Tr}\left[\mathbf{V}^{-1}(\mathbf{X}\_0)(D\mathbf{N}^n(\mathbf{X}\_0))^T \mathbf{V}(\mathbf{X}\_n) D\mathbf{N}^n(\mathbf{X}\_0)\right] \\ \mathbf{V}(\mathbf{X}) &= \left(D\Phi(\mathbf{X})\right)^T D\Phi(\mathbf{X}) \end{aligned} \tag{33}$$

*Fast Indicators for Orbital Stability: A Survey on Lyapunov and Reversibility Errors DOI: http://dx.doi.org/10.5772/intechopen.88085*

Notice that V is a positive-defined matrix and that its determinant is equal to 1 if Φ is symplectic. For a two-dimensional map, we can write

$$\mathbf{V} \equiv \begin{pmatrix} a & b \\ b & c \end{pmatrix} \qquad \qquad \mathbf{V}^{-1} = \begin{pmatrix} c & -b \\ -b & a \end{pmatrix} \tag{34}$$

where *ac* � *<sup>b</sup>*<sup>2</sup> <sup>¼</sup> 1.

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

� �1*=*<sup>2</sup>

*ds*

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

If the continuous time *t* is replaced by an integer *n* and we approximate the integral with the trapezoidal rule, we recover the relation in Eq. (25) obtained for a

We evaluate the errors for integrable maps with an elliptic fixed point at the origin, whose normal form is a rotation Rð Þ Ω with a frequency Ω depending on the distance from the origin. The LE asymptotic growth is linear, and oscillations are present unless the coordinates are normal. The RE asymptotic growth follows a power law *<sup>n</sup><sup>α</sup>* with exponent *<sup>α</sup>* <sup>¼</sup> <sup>3</sup>*=*2. If the map is linear, its asymptotic growth follows a power law with *α* ¼ 0 for LE and *α* ¼ 1*=*2 for RE. The oscillations reflect the loss of rotational symmetry when generic coordinates are used. The roundoff induced reversibility error REM is also sensitive to the choice of coordinates, and a

In generic coordinates an integrable map *M*ð Þ **x** is conjugated to its normal form *N*ð Þ **X** by a symplectic coordinate transformation **x** ¼ Φð Þ **X** ; as a consequence the

which imply that the orbits **<sup>x</sup>***<sup>n</sup>* <sup>¼</sup> *<sup>M</sup><sup>n</sup>*ð Þ **<sup>x</sup>**<sup>0</sup> and **<sup>X</sup>***<sup>n</sup>* <sup>¼</sup> *<sup>N</sup><sup>n</sup>*ð Þ **<sup>X</sup>**<sup>0</sup> are related by

� �<sup>2</sup> <sup>¼</sup> Tr *DN<sup>n</sup>* ð Þ ð Þ **<sup>X</sup>**<sup>0</sup>

� �<sup>2</sup> <sup>¼</sup> Tr *DM<sup>n</sup>* ð Þ ð Þ **<sup>x</sup>**<sup>0</sup>

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

<sup>V</sup>ð Þ¼ **<sup>X</sup>** ð Þ *<sup>D</sup>*Φð Þ **<sup>X</sup>** *<sup>T</sup> <sup>D</sup>*Φð Þ **<sup>X</sup>**

ð Þ **<sup>X</sup>**<sup>0</sup> *DN<sup>n</sup>* ð Þ ð Þ **<sup>X</sup>**<sup>0</sup>

ð Þ **<sup>x</sup>** *<sup>M</sup><sup>n</sup>*ð Þ¼ **<sup>x</sup>** <sup>Φ</sup> <sup>∘</sup> *<sup>N</sup><sup>n</sup>* <sup>∘</sup> <sup>Φ</sup>�<sup>1</sup>

ð Þ¼ **<sup>x</sup>**<sup>0</sup> *<sup>D</sup>*Φð Þ **<sup>X</sup>***<sup>n</sup> DN<sup>n</sup>*ð Þ **<sup>X</sup>**<sup>0</sup> ð Þ *<sup>D</sup>*Φð Þ **<sup>X</sup>**<sup>0</sup> �<sup>1</sup> (31)

ð Þ¼ **x** I, a relation which is proved to hold by com-

ð Þ¼ **x x**. As a consequence, the expression for the

*<sup>T</sup>* <sup>ð</sup>*DN<sup>n</sup>*ð Þ **<sup>X</sup>**<sup>0</sup>

*<sup>T</sup>* <sup>ð</sup>*DM<sup>n</sup>*ð Þ **<sup>x</sup>**<sup>0</sup>

*<sup>T</sup>* <sup>V</sup>ð Þ **<sup>X</sup>***<sup>n</sup> DN<sup>n</sup>*ð Þ **<sup>X</sup>**<sup>0</sup>

h i

h i

h i

. As shown in [12] and from

*:* (29)

ð Þ **x** (30)

(32)

(33)

random displacement *<sup>e</sup>R*ðÞ¼ *<sup>t</sup> <sup>Ξ</sup>R*ðÞ� *<sup>t</sup> <sup>Ξ</sup>R*ð Þ*<sup>t</sup>* � �1*=*<sup>2</sup>

*e*

*<sup>R</sup>*ðÞ¼ *<sup>t</sup>* <sup>2</sup>

comparison between RE and REM is presented in the next sections.

Eq. (28), we immediately obtain

symplectic map.

*Progress in Relativity*

**3. Integrable maps**

**3.1 Change of coordinate system**

conjugation equation and its iterates read

*<sup>M</sup>*ð Þ¼ **<sup>x</sup>** <sup>Φ</sup> <sup>∘</sup> *<sup>N</sup>* <sup>∘</sup> <sup>Φ</sup>�<sup>1</sup>

*DM<sup>n</sup>*ð Þ¼ **<sup>x</sup>**<sup>0</sup> *<sup>D</sup>*Φð Þ **<sup>X</sup>***<sup>n</sup> DN<sup>n</sup>*ð Þ **<sup>X</sup>**<sup>0</sup> *<sup>D</sup>*Φ�<sup>1</sup>

where we used *<sup>D</sup>*Φð Þ **<sup>X</sup>** *<sup>D</sup>*Φ�<sup>1</sup>

puting the Jacobian of Φ ∘ Φ�<sup>1</sup>

*e L <sup>n</sup>*ð Þ **x**<sup>0</sup>

**168**

**x***<sup>n</sup>* ¼ Φð Þ **X***<sup>n</sup>* . The tangent maps are given by

Lyapunov error in both coordinate systems is

*e L <sup>n</sup>*ð Þ **X**<sup>0</sup>

*e L <sup>n</sup>*ð Þ **x**<sup>0</sup>

� �<sup>2</sup> <sup>¼</sup> Tr V�<sup>1</sup>

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

If the given map is linear and two-dimensional and *M*ð Þ¼ **x** L**x** with ∣Tr L∣ < 2, then the map is conjugated to a rotation Rð Þ *ω* :

$$\mathbf{L} = \mathbf{T} \mathbf{R}(\boldsymbol{\phi}) \mathbf{T}^{-1} \qquad \qquad \mathbf{R}(\boldsymbol{\phi}) = \begin{pmatrix} \cos \left( \boldsymbol{\phi} \right) & \sin \boldsymbol{\omega} \\ -\sin \left( \boldsymbol{\omega} \right) & \cos \left( \boldsymbol{\omega} \right) \end{pmatrix} \tag{35}$$

Letting **<sup>x</sup>***<sup>n</sup>* <sup>¼</sup> <sup>L</sup>*<sup>n</sup>* **<sup>x</sup>**<sup>0</sup> and **<sup>X</sup>***<sup>n</sup>* <sup>¼</sup> <sup>R</sup>*<sup>n</sup>***X**0, the orbits in the coordinate **<sup>x</sup>** and the normal coordinate **<sup>X</sup>** <sup>¼</sup> <sup>T</sup>�<sup>1</sup> **<sup>x</sup>** and the Lyapunov errors are given by

$$\begin{aligned} \left(\epsilon\_{\pi}^{L}(\mathbf{X}\_{0})\right)^{2} &= \text{Tr}\left[\mathbf{R}(-n\omega)\mathbf{R}(n\omega)\right] = 2\\ \left(\epsilon\_{\pi}^{L}(\mathbf{x}\_{0})\right)^{2} &= \text{Tr}\left[\mathbf{V}^{-1}\mathbf{R}(-n\omega)\mathbf{V}\mathbf{R}(n\omega)\right] \\ &= \frac{(a+c)^{2}}{2} + \left(2 - \frac{(a+c)^{2}}{2}\right)\cos\left(2n\omega\right) \end{aligned} \tag{36}$$

where V <sup>¼</sup> <sup>T</sup>*<sup>T</sup>* T, and we have used the representation given by Eq. (34) where *D*Φ ¼ T. The error is constant in normal coordinate **X** and oscillates between 2 and ð Þ *a* þ *c* <sup>2</sup> � <sup>2</sup> <sup>¼</sup> *<sup>a</sup>*<sup>2</sup> <sup>þ</sup> *<sup>c</sup>*<sup>2</sup> <sup>þ</sup> <sup>2</sup>*b*<sup>2</sup> in the coordinate **<sup>x</sup>**. The geometric interpretation is obvious since the orbits of the map belong to an ellipse rather than a circle. The result for the reversibility error is given by

$$\begin{aligned} \left(e\_n^R(\mathbf{X}\_0)\right)^2 &= 4n\\ \left(e\_n^R(\mathbf{x}\_0)\right)^2 &= 2n\frac{(a+c)^2}{2} + \left(2 - \frac{(a+c)^2}{2}\right)f(n) \\\ f(n) &= \sum\_{k=1}^n \left(\cos\left(2k\omega\right) + \cos\left(2(n-k)\omega\right)\right) = \cos\left(2n\omega\right) + \frac{\cos\left(2(n-1)a\right) - \cos\left(2n\omega\right)}{1 - \cos\left(2\omega\right)} \end{aligned} \tag{37}$$

We shall first consider the dependence of the errors on the iteration order *n* from *n* ¼ 1 up to a maximum value *N*. Then, we shall consider the dependence on the initial condition **x**<sup>0</sup> when it is varied on a one-dimensional grid crossing the origin for the value *N* of the iteration number. We choose the linear map L which depends on a single parameter *λ*, and its relation with the rotation frequency is

$$\mathbf{L} = \begin{pmatrix} \mathbf{1} - \lambda & \mathbf{1} \\ -\lambda & \mathbf{1} \end{pmatrix} \qquad \qquad \sin \frac{\boldsymbol{\alpha}}{2} \quad = \frac{\sqrt{\lambda}}{2} \qquad \mathbf{0} \le \lambda \le 4 \tag{38}$$

The rotation Rð Þ *ω* **x** is the linear part for the Hénon map, whereas L**x** is the linear part of the standard map that will be discussed in the next sections. The behavior of LE and RE for these maps is provided by Eqs. (36) and (37). The error growth follows a power law with exponent *α* ¼ 0 for LE and *α* ¼ 1*=*2 for RE. Oscillations are present when the coordinates are not normal.

For a generic map such as L defined by Eq. (38), the growth of REM follows a power law exponent *α* ¼ 1*=*2 as RE, for almost any value of *λ*, as shown by **Figure 1**, right panel, where the plot of MEGNO corresponding to 2*α* is shown. The result for the map Rð Þ *ω* in normal coordinates is shown in the left panel of the same figure, and the exponent is *α* ¼ 1 for almost all the values of *ω*.

Letting **<sup>X</sup>** <sup>¼</sup> ð Þ *<sup>X</sup>; <sup>P</sup> <sup>T</sup>* and ð Þ *<sup>ϕ</sup>; <sup>J</sup>* be the action angle coordinates defined by *<sup>X</sup>* <sup>¼</sup> ð Þ <sup>2</sup>*<sup>J</sup>* <sup>1</sup>*=*<sup>2</sup> cos *<sup>ϕ</sup>* and *<sup>P</sup>* ¼ �ð Þ <sup>2</sup>*<sup>J</sup>* <sup>1</sup>*=*<sup>2</sup> sin *<sup>ϕ</sup>*, the rotation in the **<sup>X</sup>** plane becomes a translation on the cylinder:

$$
\phi\_n = \phi\_{n-1} + a \mod 2\pi \qquad \qquad J\_n = J\_{n-1} \tag{39}
$$

the rotation Rð Þ *ω* and the linear map L. The linear dependence is evident in both

An integrable map *M* in normal coordinates and the tangent map *DM<sup>n</sup>* read

*=*2 is the action. The square of the Lyapunov error<sup>4</sup> reads

<sup>¼</sup> <sup>4</sup>*<sup>n</sup>* <sup>þ</sup> <sup>2</sup>*J*Ω<sup>0</sup> ð Þ<sup>2</sup> <sup>2</sup> *<sup>n</sup>*<sup>3</sup>

<sup>¼</sup> <sup>14</sup>*:*<sup>5</sup> *x*2 <sup>0</sup> <sup>þ</sup> *<sup>p</sup>*<sup>2</sup> 0 � �∣Ω<sup>0</sup>

*Yn* � *Y e*ð Þ*<sup>n</sup>* . The range of variation is [1, 3]. One can prove that for a given initial

In **Figure 3**, we show the variation with *n* ∈½ � 1*;* 1000 of *Yn* computed for RE given by Eq. (42), corresponding to the map presented in Eq. (40), where Ωð Þ*J* is a linear function of *J*, and find a perfect agreement with the analytical estimate of the value of *n* for which the value *Yn* ¼ 2 is reached. In **Figure 4**, we show the variation

*<sup>N</sup>* and the corresponding MEGNO filter *YN* with the initial condition chosen on a one-dimensional grid crossing the origin for *N* ¼ 1000. The integrable map is given by Eq. (40), where Ω<sup>0</sup> is constant. The error reaches a minimum value at the origin,

origin so that the behavior is similar even though in this case the fluctuations are large. We notice that MEGNO does not eliminate the fluctuations of REM. In order

cost is of the order of *n*2. This can be avoided by storing the sequence **x***ε,m* and

<sup>∥</sup>**x**∥<sup>2</sup> *<sup>η</sup>*<sup>0</sup> ð Þ � **<sup>x</sup>** <sup>2</sup> <sup>þ</sup> <sup>2</sup>*n*Ω<sup>0</sup> *<sup>η</sup>*<sup>0</sup> ð Þ � **<sup>x</sup>** *<sup>η</sup>*<sup>0</sup> ð Þ � <sup>J</sup>**<sup>x</sup>** and J <sup>¼</sup> 0 1

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

gives the Lyapunov error *e<sup>L</sup>*

<sup>4</sup> The standard definition for an initial displacement along the unit vector *η*<sup>0</sup> is *eL*

� �<sup>2</sup> <sup>¼</sup> Tr *DM<sup>n</sup>* ð Þ*<sup>T</sup> DM<sup>n</sup>* � � <sup>¼</sup> <sup>2</sup> <sup>þ</sup> *<sup>n</sup>*<sup>2</sup> <sup>2</sup>*J*Ω<sup>0</sup> ð Þ<sup>2</sup> (41)

*k*¼1

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

*=d* log *n*, and its double average is given by MEGNO

3 þ *n* 6 � �

R0

ð Þ *<sup>n</sup>* � *<sup>k</sup>* <sup>2</sup> <sup>þ</sup>X*<sup>n</sup>*

!

*k*¼1 *k*2

, whose asymptotic value is

<sup>∣</sup> (43)

*<sup>N</sup>* decreases by approaching the

*<sup>n</sup> <sup>η</sup>*<sup>0</sup> ð Þ¼ <sup>∥</sup>*DMnη*0<sup>∥</sup> where

*<sup>m</sup>* for *m* ¼ 1*,* …*, n* whose computational

�1 0 � �. The sum

*<sup>ε</sup>* ð Þ� **x***ε,n* **x***<sup>n</sup>*�*<sup>m</sup>*∥ which turns out to be comparable with

(42)

ð Þ *<sup>n</sup>*<sup>Ω</sup> **xx***<sup>T</sup>* (40)

cases, even though the fluctuations are large for the linear map.

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

*<sup>M</sup>*ð Þ¼ **<sup>x</sup>** <sup>R</sup>ð Þ <sup>Ω</sup>ð Þ*<sup>J</sup>* **<sup>x</sup>** *DMn*ð Þ¼ **<sup>x</sup>** <sup>R</sup>ð Þþ *<sup>n</sup>*<sup>Ω</sup> *<sup>n</sup>*Ω<sup>0</sup>

� �<sup>2</sup> � � <sup>¼</sup> <sup>4</sup>*<sup>n</sup>* <sup>þ</sup> <sup>2</sup>*J*Ω<sup>0</sup> ð Þ<sup>2</sup> <sup>X</sup>*<sup>n</sup>*

and the square of the reversibility error is given by

For a fixed value of the invariant *J*, the slope of *e<sup>R</sup>*

condition, the intermediate value *Yn* ¼ 2 is reached for

and a similar behavior for *YN* is observed. Also *eREM*

to compute *Yn*, one needs the sequence *eREM*

computing ^*emREM* <sup>¼</sup> <sup>∥</sup>*M*�*<sup>m</sup>*

<sup>∥</sup>*DMnη*0∥<sup>2</sup> <sup>¼</sup> <sup>1</sup> <sup>þ</sup> *<sup>n</sup>*Ω<sup>0</sup> ð Þ<sup>2</sup>

*<sup>n</sup>* <sup>¼</sup> <sup>14</sup>*:*<sup>5</sup> 2*J*∣Ω<sup>0</sup> ∣

**3.3 Anisochronous rotations**

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

*e L n*

> þ *e L n*�*k*

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

where *<sup>J</sup>* <sup>¼</sup> <sup>∥</sup>**x**∥<sup>2</sup>

*k*¼1

2*α*, is defined as *d* log *e<sup>R</sup>*

*e L k* � �<sup>2</sup>

*eR n* � �<sup>2</sup> <sup>¼</sup> <sup>X</sup>*<sup>n</sup>*

of *e<sup>R</sup>*

*eREM <sup>m</sup>* .

� � � � *DMn* <sup>1</sup> 0 � �� � � � 2 þ � � � � *DM<sup>n</sup>* <sup>0</sup> 1 � �� � � � 2

**171**

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 *α* varies in the range 0½ � *;* 1 . Unlike RE, we observe that REM depends linearly on the distance of the initial condition **x**<sup>0</sup> from the origin. In **Figure 2**, we plot *eREM <sup>N</sup>* as a function of the initial condition when it varies on a one-dimensional grid issued from the origin for

**Figure 1.**

*Left frame: twice the asymptotic power law exponent provided by the MEGNO filter YN with N* ¼ *1000 applied to REM for a rotation R*ð Þ *ω where ω=2π varies in the interval* ½ � *0; 1=2 . The initial condition is x0* ¼ *0:1, p0* ¼ *0. Right frame: twice the asymptotic power law exponent provided by the MEGNO filter YN with N* ¼ *1000 applied to REM for a linear map L given by Eq. (38) whose parameter λ varies in* ½ � *0; 1 . Initial condition x0* ¼ *0:1, p0* ¼ *0.*

#### **Figure 2.**

*Left frame: reversibility error due to roundoff eREM <sup>N</sup> for a rotation <sup>R</sup>*ð Þ *<sup>ω</sup> with <sup>ω</sup>* <sup>¼</sup> *<sup>2</sup><sup>π</sup>* ffiffiffi *2* <sup>p</sup> � *<sup>1</sup>* � � *and <sup>N</sup>* <sup>¼</sup> *<sup>1000</sup> when the initial condition is varied. We choose x0* ∈ ½ � *0; 0:5 , p0* ¼ *0. The dependence on x0 is evident and a linear fit f x*ð Þ¼ *<sup>0</sup> 5000x0 is shown, purple line. Right frame: computation of the error for the linear map with <sup>λ</sup>* <sup>¼</sup> *4 sin <sup>2</sup>*ð Þ *<sup>ω</sup>=<sup>2</sup> where <sup>ω</sup> has the same value. The linear fit is the same.*

the rotation Rð Þ *ω* and the linear map L. The linear dependence is evident in both cases, even though the fluctuations are large for the linear map.

## **3.3 Anisochronous rotations**

LE and RE for these maps is provided by Eqs. (36) and (37). The error growth follows a power law with exponent *α* ¼ 0 for LE and *α* ¼ 1*=*2 for RE. Oscillations are

Letting **<sup>X</sup>** <sup>¼</sup> ð Þ *<sup>X</sup>; <sup>P</sup> <sup>T</sup>* and ð Þ *<sup>ϕ</sup>; <sup>J</sup>* be the action angle coordinates defined by

cos *<sup>ϕ</sup>* and *<sup>P</sup>* ¼ �ð Þ <sup>2</sup>*<sup>J</sup>* <sup>1</sup>*=*<sup>2</sup> sin *<sup>ϕ</sup>*, the rotation in the **<sup>X</sup>** plane becomes a

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 *α* varies in the range 0½ � *;* 1 . Unlike RE, we observe that REM depends linearly on the distance of the

initial condition when it varies on a one-dimensional grid issued from the origin for

*Left frame: twice the asymptotic power law exponent provided by the MEGNO filter YN with N* ¼ *1000 applied to REM for a rotation R*ð Þ *ω where ω=2π varies in the interval* ½ � *0; 1=2 . The initial condition is x0* ¼ *0:1, p0* ¼ *0. Right frame: twice the asymptotic power law exponent provided by the MEGNO filter YN with N* ¼ *1000 applied to REM for a linear map L given by Eq. (38) whose parameter λ varies in* ½ � *0; 1 . Initial*

*when the initial condition is varied. We choose x0* ∈ ½ � *0; 0:5 , p0* ¼ *0. The dependence on x0 is evident and a linear fit f x*ð Þ¼ *<sup>0</sup> 5000x0 is shown, purple line. Right frame: computation of the error for the linear map with*

*<sup>N</sup> for a rotation <sup>R</sup>*ð Þ *<sup>ω</sup> with <sup>ω</sup>* <sup>¼</sup> *<sup>2</sup><sup>π</sup>* ffiffiffi

*2*

<sup>p</sup> � *<sup>1</sup>* � � *and <sup>N</sup>* <sup>¼</sup> *<sup>1000</sup>*

*ϕ<sup>n</sup>* ¼ *ϕn*�<sup>1</sup> þ *ω* mod 2*π Jn* ¼ *Jn*�<sup>1</sup> (39)

*<sup>N</sup>* as a function of the

For a generic map such as L defined by Eq. (38), the growth of REM follows a power law exponent *α* ¼ 1*=*2 as RE, for almost any value of *λ*, as shown by **Figure 1**, right panel, where the plot of MEGNO corresponding to 2*α* is shown. The result for the map Rð Þ *ω* in normal coordinates is shown in the left panel of the same figure,

present when the coordinates are not normal.

*<sup>X</sup>* <sup>¼</sup> ð Þ <sup>2</sup>*<sup>J</sup>* <sup>1</sup>*=*<sup>2</sup>

*Progress in Relativity*

**Figure 1.**

**Figure 2.**

**170**

*Left frame: reversibility error due to roundoff eREM*

*<sup>λ</sup>* <sup>¼</sup> *4 sin <sup>2</sup>*ð Þ *<sup>ω</sup>=<sup>2</sup> where <sup>ω</sup> has the same value. The linear fit is the same.*

*condition x0* ¼ *0:1, p0* ¼ *0.*

translation on the cylinder:

and the exponent is *α* ¼ 1 for almost all the values of *ω*.

initial condition **x**<sup>0</sup> from the origin. In **Figure 2**, we plot *eREM*

An integrable map *M* in normal coordinates and the tangent map *DM<sup>n</sup>* read

$$M(\mathbf{x}) = \mathbf{R}(\boldsymbol{\Omega}(f))\mathbf{x} \qquad \qquad DM^n(\mathbf{x}) = \mathbf{R}(n\boldsymbol{\Omega}) + n\boldsymbol{\Omega}^\prime \mathbf{R}^\prime (n\boldsymbol{\Omega})\mathbf{x}\mathbf{x}^T \tag{40}$$

where *<sup>J</sup>* <sup>¼</sup> <sup>∥</sup>**x**∥<sup>2</sup> *=*2 is the action. The square of the Lyapunov error<sup>4</sup> reads

$$\left(\left(\mathbf{e}\_{n}^{L}\right)^{2} = \text{Tr}\left(\left(D\mathbf{M}^{n}\right)^{T}D\mathbf{M}^{n}\right) = 2 + n^{2}\left(\mathbf{2J}\boldsymbol{\Omega}^{\prime}\right)^{2}\tag{41}$$

and the square of the reversibility error is given by

$$\begin{split} \left( \left( \mathbf{e}\_{n}^{R} \right)^{2} = \sum\_{k=1}^{n} \left( \left( \mathbf{e}\_{k}^{L} \right)^{2} + \left( \mathbf{e}\_{n-k}^{L} \right)^{2} \right) &= 4n + \left( 2J\Omega' \right)^{2} \left( \sum\_{k=1}^{n} (n-k)^{2} + \sum\_{k=1}^{n} k^{2} \right) \\ &= 4n + \left( 2J\Omega' \right)^{2} 2 \left( \frac{n^{3}}{3} + \frac{n}{6} \right) \end{split} \tag{42}$$

For a fixed value of the invariant *J*, the slope of *e<sup>R</sup> n* � �<sup>2</sup> , whose asymptotic value is 2*α*, is defined as *d* log *e<sup>R</sup> n* � �<sup>2</sup> *=d* log *n*, and its double average is given by MEGNO *Yn* � *Y e*ð Þ*<sup>n</sup>* . The range of variation is [1, 3]. One can prove that for a given initial condition, the intermediate value *Yn* ¼ 2 is reached for

$$m = \frac{14.5}{2|\!/\Omega'|} = \frac{14.5}{\left(\varkappa\_0^2 + p\_0^2\right)|\!/\Omega'|}\tag{43}$$

In **Figure 3**, we show the variation with *n* ∈½ � 1*;* 1000 of *Yn* computed for RE given by Eq. (42), corresponding to the map presented in Eq. (40), where Ωð Þ*J* is a linear function of *J*, and find a perfect agreement with the analytical estimate of the value of *n* for which the value *Yn* ¼ 2 is reached. In **Figure 4**, we show the variation of *e<sup>R</sup> <sup>N</sup>* and the corresponding MEGNO filter *YN* with the initial condition chosen on a one-dimensional grid crossing the origin for *N* ¼ 1000. The integrable map is given by Eq. (40), where Ω<sup>0</sup> is constant. The error reaches a minimum value at the origin, and a similar behavior for *YN* is observed. Also *eREM <sup>N</sup>* decreases by approaching the origin so that the behavior is similar even though in this case the fluctuations are large. We notice that MEGNO does not eliminate the fluctuations of REM. In order to compute *Yn*, one needs the sequence *eREM <sup>m</sup>* for *m* ¼ 1*,* …*, n* whose computational cost is of the order of *n*2. This can be avoided by storing the sequence **x***ε,m* and computing ^*emREM* <sup>¼</sup> <sup>∥</sup>*M*�*<sup>m</sup> <sup>ε</sup>* ð Þ� **x***ε,n* **x***<sup>n</sup>*�*<sup>m</sup>*∥ which turns out to be comparable with *eREM <sup>m</sup>* .

<sup>4</sup> The standard definition for an initial displacement along the unit vector *η*<sup>0</sup> is *eL <sup>n</sup> <sup>η</sup>*<sup>0</sup> ð Þ¼ <sup>∥</sup>*DMnη*0<sup>∥</sup> where <sup>∥</sup>*DMnη*0∥<sup>2</sup> <sup>¼</sup> <sup>1</sup> <sup>þ</sup> *<sup>n</sup>*Ω<sup>0</sup> ð Þ<sup>2</sup> <sup>∥</sup>**x**∥<sup>2</sup> *<sup>η</sup>*<sup>0</sup> ð Þ � **<sup>x</sup>** <sup>2</sup> <sup>þ</sup> <sup>2</sup>*n*Ω<sup>0</sup> *<sup>η</sup>*<sup>0</sup> ð Þ � **<sup>x</sup>** *<sup>η</sup>*<sup>0</sup> ð Þ � <sup>J</sup>**<sup>x</sup>** and J <sup>¼</sup> 0 1 �1 0 � �. The sum � � � � *DMn* <sup>1</sup> 0 � �� � � � 2 þ � � � � *DM<sup>n</sup>* <sup>0</sup> 1 � �� � � � 2 gives the Lyapunov error *e<sup>L</sup> n* � �<sup>2</sup> .

*xn*þ<sup>1</sup> <sup>¼</sup> ð Þ <sup>1</sup> � *<sup>λ</sup> xn* <sup>þ</sup> *pn pn*þ<sup>1</sup> ¼ �*λxn* <sup>þ</sup> *pn* (45)

ð Þ <sup>2</sup>*<sup>π</sup>* <sup>2</sup> cos 2ð Þ *<sup>π</sup><sup>x</sup>* (46)

<sup>p</sup> (see Eq. (38)) when *<sup>λ</sup>* ! 0. Since the

<sup>p</sup> <sup>≫</sup> 1, the symplectic integrator in

*<sup>π</sup>* cosð Þ *<sup>π</sup><sup>x</sup> :* (47)

� sin*ω* cos*ω*

*λ*

� � (48)

<sup>3</sup> (49)

ð Þ *<sup>ω</sup>* <sup>þ</sup> <sup>2</sup>*<sup>x</sup> <sup>=</sup>*<sup>3</sup> <sup>p</sup> . The

<sup>p</sup> *<sup>=</sup>π*. When *<sup>λ</sup>*

*λ* <sup>p</sup> *<sup>=</sup>*2.

This map is conjugated to a rotation Rð Þ *<sup>ω</sup>* for 0 < *<sup>λ</sup>* < 2 where sin ð Þ¼ *<sup>ω</sup>=*<sup>2</sup> ffiffi

which is the interpolating Hamiltonian of the map when *λ* ! 0. We observe that

*λ*

*λ*

The point *x* ¼ �1*=*2*, p* ¼ 0 is hyperbolic, and for *λ* ≪ 1 the corresponding orbit is

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

<sup>2</sup> � *<sup>λ</sup>*

Eq. (44), obtained for a time step Δ*t* ¼ 1, provides a good approximation to the orbit. Conversely, the Hamiltonian provides a good interpolation to the orbit of the

> ffiffi *λ* p

increases, non-integrable features appear, such as chains of islands corresponding to resonances and a chaotic region near the separatrix due to homoclinic intersections.

Close to the origin, this is just a rotation with frequency *ω*. For *ω* ! 0 this is a

*<sup>H</sup>* <sup>¼</sup> *<sup>ω</sup> <sup>p</sup>*<sup>2</sup> <sup>þ</sup> *<sup>x</sup>*<sup>2</sup>

with time step Δ*t* ¼ 1. The approximation is good since the characteristic time is the period of the linear rotation *T* ¼ 2*π=ω*. The motion is bounded

Birkhoff normal forms provide an integrable approximation to the map and the corresponding interpolating Hamiltonian, from which the errors may be analyti-

We have analyzed the errors *en* for a fixed initial condition by varying *n* up to a maximum value *N*, by varying the initial condition on a one-dimensional grid for *n* ¼ *N* and by choosing a grid in the phase plane for *n* ¼ *N*. The LE shows oscillations with *n*, RE grows without oscillations, and the behavior of RE is similar to RE although with large fluctuations. The results obtained by filtering the errors with

*<sup>x</sup>* <sup>¼</sup> 2 tan <sup>ð</sup>*ω=*2Þ*; <sup>p</sup>* ¼ �2 tan <sup>2</sup> ð Þ ð Þ *<sup>ω</sup>=*<sup>2</sup> which corresponds to the critical point ð Þ *x* ¼ *ω; p* ¼ 0 of the Hamiltonian. The stability boundary is approximated by

*H x*ð Þ¼ *; <sup>p</sup> <sup>ω</sup>*<sup>3</sup>*=*6 whose orbit explicitly reads *<sup>p</sup>* ¼ �ð Þ *<sup>ω</sup>* � *<sup>x</sup>* ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

� � <sup>R</sup>ð Þ¼ *<sup>ω</sup>* cos*<sup>ω</sup>* sin*<sup>ω</sup>*

<sup>2</sup> � *<sup>x</sup>*<sup>3</sup>

*p* ¼ �

As a consequence, for *λ* small the width of the separatrix is 2 ffiffi

approximated by the separatrix of the Hamiltonian:

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

the frequency for small oscillations is *<sup>ω</sup>* <sup>¼</sup> ffiffi

time scale of the Hamiltonian *<sup>H</sup>* is *<sup>T</sup>* <sup>¼</sup> <sup>2</sup>*π<sup>=</sup>* ffiffi

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

<sup>¼</sup> <sup>R</sup>ð Þ *<sup>ω</sup> xn*

good symplectic integrator of the Hamiltonian:

*pn* <sup>þ</sup> *<sup>x</sup>*<sup>2</sup> *n*

by the orbit issued from the hyperbolic fixed point of the map

*xn*þ<sup>1</sup> *pn*þ<sup>1</sup>

cally computed.

**173**

**5. The standard map**

!

map. The equation for the separatrix of *H* is given by

*<sup>H</sup>* <sup>¼</sup> *<sup>p</sup>*<sup>2</sup>

#### **Figure 3.**

*Left frame: plot of MEGNO Yn on the error e<sup>R</sup> <sup>N</sup> for 1*≤*n*≤ *N with N* ¼ *1000 for the integrable map with Ω*<sup>0</sup> ¼ *0:1 and initial condition x0* ¼ *0:5, p0* ¼ *0 (blue line). The green line refers to a modified definition YMn, where n loge2 <sup>n</sup>* � *loge2 n*�*1* � � *is replaced with loge<sup>2</sup> <sup>n</sup>* � *loge2 n*�*1* � �*=*ð Þ *logn* � *log n*ð Þ � *<sup>1</sup> , which, applied to the sequence en* <sup>¼</sup> *<sup>n</sup><sup>α</sup>, gives <sup>2</sup><sup>α</sup> for any <sup>n</sup>. The vertical line gives the theoretical estimate of the value of <sup>n</sup> for which Yn* <sup>¼</sup> *<sup>2</sup> (see Eq. (43)) Right frame: the same for x0* ¼ *1, p0* ¼ *0.*

#### **Figure 4.**

*Left frame: plot or the error e<sup>R</sup> <sup>N</sup> with N* ¼ *1000 for the integrable map as a function of the initial condition x0, p0* ¼ *0 with* Ω<sup>0</sup> ð Þ¼ *J 0:1 (red line) and* Ω<sup>0</sup> ð Þ¼ *<sup>J</sup> <sup>1</sup> (blue line). Center frame: same plot with eREM <sup>N</sup> for <sup>ω</sup>* <sup>¼</sup> *<sup>2</sup><sup>π</sup>* ffiffiffi *2* <sup>p</sup> � *<sup>1</sup>* � � *and* <sup>Ω</sup><sup>0</sup> ð Þ¼ *<sup>J</sup> <sup>0</sup>:1, gray line, compared with <sup>e</sup><sup>R</sup> <sup>n</sup> , red line. Right frame: plot of YN for the integrable map with* Ω<sup>0</sup> ¼ *0:1 (red line) and* Ω<sup>0</sup> ð Þ¼ *J 1(blue line).*

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 1*=*2 in the exponent of the power law *n<sup>α</sup>* occurs close to the origin, if the linear part is a rotation R, as for the Hénon map. If the linear part is L as for the standard map, there is no shift. The better correspondence is not surprising since the computational complexity of the map is higher when the coordinates are not normal.

#### **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 T<sup>2</sup> and reads

$$p\_{n+1} = p\_n - \frac{\lambda}{2\pi} \sin\left(2\pi\mathbf{x}\_n\right) \mod \mathbf{1} \qquad \mathbf{x}\_{n+1} = \mathbf{x}\_n + p\_{n+1} \mod \mathbf{1} \tag{44}$$

where *x, p* belong to the interval ½ � �1*=*2*;* 1*=*2 whose ends are identified. For *λ* ≪ 1 and ∣*p*∣ ≫ ffiffi *λ* <sup>p</sup> , it is just a weakly perturbed rotator, and *x, p* are action angle coordinates. The origin is an elliptic fixed and very close to it; the map is approximated by a linear map

*Fast Indicators for Orbital Stability: A Survey on Lyapunov and Reversibility Errors DOI: http://dx.doi.org/10.5772/intechopen.88085*

$$\boldsymbol{x}\_{n+1} = (1 - \lambda)\boldsymbol{x}\_n + \boldsymbol{p}\_n \qquad \qquad \boldsymbol{p}\_{n+1} = -\lambda \boldsymbol{x}\_n + \boldsymbol{p}\_n \tag{45}$$

This map is conjugated to a rotation Rð Þ *<sup>ω</sup>* for 0 < *<sup>λ</sup>* < 2 where sin ð Þ¼ *<sup>ω</sup>=*<sup>2</sup> ffiffi *λ* <sup>p</sup> *<sup>=</sup>*2. The point *x* ¼ �1*=*2*, p* ¼ 0 is hyperbolic, and for *λ* ≪ 1 the corresponding orbit is approximated by the separatrix of the Hamiltonian:

$$H = \frac{p^2}{2} - \frac{\lambda}{\left(2\pi\right)^2} \cos\left(2\pi\infty\right) \tag{46}$$

which is the interpolating Hamiltonian of the map when *λ* ! 0. We observe that the frequency for small oscillations is *<sup>ω</sup>* <sup>¼</sup> ffiffi *λ* <sup>p</sup> (see Eq. (38)) when *<sup>λ</sup>* ! 0. Since the time scale of the Hamiltonian *<sup>H</sup>* is *<sup>T</sup>* <sup>¼</sup> <sup>2</sup>*π<sup>=</sup>* ffiffi *λ* <sup>p</sup> <sup>≫</sup> 1, the symplectic integrator in Eq. (44), obtained for a time step Δ*t* ¼ 1, provides a good approximation to the orbit. Conversely, the Hamiltonian provides a good interpolation to the orbit of the map. The equation for the separatrix of *H* is given by

$$p = \pm \frac{\sqrt{\lambda}}{\pi} \cos \left( \pi \mathbf{x} \right). \tag{47}$$

As a consequence, for *λ* small the width of the separatrix is 2 ffiffi *λ* <sup>p</sup> *<sup>=</sup>π*. When *<sup>λ</sup>* increases, non-integrable features appear, such as chains of islands corresponding to resonances and a chaotic region near the separatrix due to homoclinic intersections.

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

$$\begin{pmatrix} \mathbf{x}\_{n+1} \\ p\_{n+1} \end{pmatrix} = \mathbf{R}(\boldsymbol{\omega}) \begin{pmatrix} \mathbf{x}\_{n} \\ p\_{n} + \mathbf{x}\_{n}^{2} \end{pmatrix} \qquad \qquad \mathbf{R}(\boldsymbol{\omega}) = \begin{pmatrix} \cos \boldsymbol{\omega} & \sin \boldsymbol{\omega} \\ -\sin \boldsymbol{\omega} & \cos \boldsymbol{\omega} \end{pmatrix} \tag{48}$$

Close to the origin, this is just a rotation with frequency *ω*. For *ω* ! 0 this is a good symplectic integrator of the Hamiltonian:

$$H = \alpha \frac{p^2 + x^2}{2} - \frac{x^3}{3} \tag{49}$$

with time step Δ*t* ¼ 1. The approximation is good since the characteristic time is the period of the linear rotation *T* ¼ 2*π=ω*. The motion is bounded by the orbit issued from the hyperbolic fixed point of the map *<sup>x</sup>* <sup>¼</sup> 2 tan <sup>ð</sup>*ω=*2Þ*; <sup>p</sup>* ¼ �2 tan <sup>2</sup> ð Þ ð Þ *<sup>ω</sup>=*<sup>2</sup> which corresponds to the critical point ð Þ *x* ¼ *ω; p* ¼ 0 of the Hamiltonian. The stability boundary is approximated by *H x*ð Þ¼ *; <sup>p</sup> <sup>ω</sup>*<sup>3</sup>*=*6 whose orbit explicitly reads *<sup>p</sup>* ¼ �ð Þ *<sup>ω</sup>* � *<sup>x</sup>* ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð Þ *<sup>ω</sup>* <sup>þ</sup> <sup>2</sup>*<sup>x</sup> <sup>=</sup>*<sup>3</sup> <sup>p</sup> . The Birkhoff normal forms provide an integrable approximation to the map and the corresponding interpolating Hamiltonian, from which the errors may be analytically computed.

#### **5. The standard map**

We have analyzed the errors *en* for a fixed initial condition by varying *n* up to a maximum value *N*, by varying the initial condition on a one-dimensional grid for *n* ¼ *N* and by choosing a grid in the phase plane for *n* ¼ *N*. The LE shows oscillations with *n*, RE grows without oscillations, and the behavior of RE is similar to RE although with large fluctuations. The results obtained by filtering the errors with

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 1*=*2 in the exponent of the power law *n<sup>α</sup>* occurs close to the origin, if the linear part is a rotation R, as for the Hénon map. If the linear part is L as for the standard map, there is no shift. The better correspondence is not surprising since the computational complexity of the map is

ð Þ¼ *J 1(blue line).*

*Ω*<sup>0</sup> ¼ *0:1 and initial condition x0* ¼ *0:5, p0* ¼ *0 (blue line). The green line refers to a modified definition YMn,*

*en* <sup>¼</sup> *<sup>n</sup><sup>α</sup>, gives <sup>2</sup><sup>α</sup> for any <sup>n</sup>. The vertical line gives the theoretical estimate of the value of <sup>n</sup> for which Yn* <sup>¼</sup> *<sup>2</sup> (see*

*<sup>n</sup>* � *loge2 n*�*1*

*<sup>N</sup> with N* ¼ *1000 for the integrable map as a function of the initial condition*

ð Þ¼ *<sup>J</sup> <sup>1</sup> (blue line). Center frame: same plot with eREM*

*<sup>N</sup> for 1*≤*n*≤ *N with N* ¼ *1000 for the integrable map with*

� �*=*ð Þ *logn* � *log n*ð Þ � *<sup>1</sup> , which, applied to the sequence*

*<sup>n</sup> , red line. Right frame: plot of YN for the*

*<sup>N</sup> for*

We examine here the behavior of the proposed dynamical indicators for two

where *x, p* belong to the interval ½ � �1*=*2*;* 1*=*2 whose ends are identified. For *λ* ≪ 1

nates. The origin is an elliptic fixed and very close to it; the map is approximated by

<sup>p</sup> , it is just a weakly perturbed rotator, and *x, p* are action angle coordi-

<sup>2</sup>*<sup>π</sup>* sin 2ð Þ *<sup>π</sup>xn* mod 1 *xn*þ<sup>1</sup> <sup>¼</sup> *xn* <sup>þ</sup> *pn*þ<sup>1</sup> mod 1 (44)

The **standard map** is defined as a map on the torus T<sup>2</sup> and reads

higher when the coordinates are not normal.

ð Þ¼ *J 0:1 (red line) and* Ω<sup>0</sup>

*integrable map with* Ω<sup>0</sup> ¼ *0:1 (red line) and* Ω<sup>0</sup>

*Left frame: plot of MEGNO Yn on the error e<sup>R</sup>*

*Eq. (43)) Right frame: the same for x0* ¼ *1, p0* ¼ *0.*

*<sup>n</sup>* � *loge2 n*�*1* � � *is replaced with loge<sup>2</sup>*

ð Þ¼ *<sup>J</sup> <sup>0</sup>:1, gray line, compared with <sup>e</sup><sup>R</sup>*

basic models, the standard map and the Hénon map.

**4. Non-integrable maps**

*pn*þ<sup>1</sup> <sup>¼</sup> *pn* � *<sup>λ</sup>*

*λ*

and ∣*p*∣ ≫ ffiffi

**Figure 4.**

**Figure 3.**

*where n loge2*

*Progress in Relativity*

*<sup>ω</sup>* <sup>¼</sup> *<sup>2</sup><sup>π</sup>* ffiffiffi *2* <sup>p</sup> � *<sup>1</sup>* � � *and* <sup>Ω</sup><sup>0</sup>

*x0, p0* ¼ *0 with* Ω<sup>0</sup>

*Left frame: plot or the error e<sup>R</sup>*

a linear map

MEGNO confirm this observation. In **Figure 5**, we plot the errors *en* for *λ* ¼ 0*:*1 and *Y e*ð Þ*<sup>n</sup>* by varying *n*. The fast oscillations of LE and the large fluctuations of REM are clearly visible.

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 1*=*ϵ where ϵ is the machine accuracy. Typically in double precision, the overflow corresponds to 10<sup>300</sup> where <sup>ϵ</sup>�<sup>1</sup> � <sup>10</sup>17. The same limitation is met when the Lyapunov error is computed using the shadow orbit method without renormalization rather than with the variational method. In **Figure 6**, we show the errors for a chaotic orbit when *λ* ¼ 0*:*8. Both LE and RE exhibit an exponential growth after an initial transitory phase. The behavior of REM is very similar until *n*≤ 300. For higher values the saturation to 10<sup>17</sup> is evident, and REM ceases to grow exponentially.

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

**Figure 7** shows the variation of LE, RE, and REM for *λ* ¼ 0*:*1, with the initial condition chosen on a regular grid in the vertical axis *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 *λ* ¼ 0*:*8 at which the dynamical structure is rich due to the presence of many resonances and small chaotic regions. The effectiveness in detecting the resonances

*Left frame: variation for the standard map with λ* ¼ *0:8 of the errors LE (blue line), RE (red line), and REM (gray line), computed at N* ¼ *1000 with the initial condition x0* ¼ *0, p0* ∈½ � *0; 0:5 . Right frame: magnification*

*Left frame: variation for the standard map with λ* ¼ *0:1 of the errors LE (blue line), RE (red line), REM (gray line) computed at N* ¼ *1000 with the initial condition x0* ¼ *0, p0* ∈ ½ � �*0:15; 1; 0:15 . Right frame: the same for*

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

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

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 *Ng* � *Ng* points, where we compute the errors and show the result using a color scale. In order to analyze the details, smaller squares may be chosen eventually increasing the iteration number. In **Figure 9**, we show for *N* ¼ 500 and *Ng* ¼ 200 the color plots for the errors of the standard map with *λ* ¼ 0*:*8 and in **Figure 10** for *λ* ¼ 1*:*5. In the first case, the measure of chaotic orbits is small with respect to the regular ones. We observe that LE has some weak structures within the main regular component surrounding the origin, visible when the figure is sufficiently magnified. Such structures of LE, related to the oscillating growth with *n*, disappear when MEGNO is computed and

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

is evident.

**175**

**Figure 8.**

*in the interval p0* ∈ ½ � *0:25; 0:5 .*

**Figure 7.**

*MEGNO YN.*

#### **Figure 5.**

*Left frame: plot of the errors for the standard map with λ* ¼ *0:1 and initial condition x0* ¼ *0, p0* ¼ *0:075. Lyapunov error eL <sup>n</sup> (blue line), reversibility errors eR <sup>n</sup> (red line), and eREM <sup>n</sup> (gray line), for 1*≤*n*≤*2500. Right frame: plots for the MEGNO filter Yn for the same errors.*

#### **Figure 6.**

*Left frame: plot of the errors for the standard map with λ* ¼ *0:8 and initial condition x0* ¼ *0, p0* ¼ *0:26 corresponding to a chaotic orbit. Lyapunov error e<sup>L</sup> <sup>n</sup> (blue line), reversibility errors e<sup>R</sup> <sup>n</sup> (red line), and eREM <sup>n</sup> (gray line), for n* ≤*500. Right frame: plots for the MEGNO filter Yn for the same errors.*

*Fast Indicators for Orbital Stability: A Survey on Lyapunov and Reversibility Errors DOI: http://dx.doi.org/10.5772/intechopen.88085*

**Figure 7.**

MEGNO confirm this observation. In **Figure 5**, we plot the errors *en* for *λ* ¼ 0*:*1 and *Y e*ð Þ*<sup>n</sup>* by varying *n*. The fast oscillations of LE and the large fluctuations of REM are

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 1*=*ϵ where ϵ is the machine accuracy. Typically in double precision, the overflow corresponds to 10<sup>300</sup> where <sup>ϵ</sup>�<sup>1</sup> � <sup>10</sup>17. The same limitation is met when the Lyapunov error is computed using the shadow orbit method without

renormalization rather than with the variational method. In **Figure 6**, we show the errors for a chaotic orbit when *λ* ¼ 0*:*8. Both LE and RE exhibit an exponential growth after an initial transitory phase. The behavior of REM is very similar until *n*≤ 300. For higher values the saturation to 10<sup>17</sup> is evident, and REM ceases to grow

**Figure 7** shows the variation of LE, RE, and REM for *λ* ¼ 0*:*1, with the initial condition chosen on a regular grid in the vertical axis *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

*Left frame: plot of the errors for the standard map with λ* ¼ *0:1 and initial condition x0* ¼ *0, p0* ¼ *0:075.*

*Left frame: plot of the errors for the standard map with λ* ¼ *0:8 and initial condition x0* ¼ *0, p0* ¼ *0:26*

*line), for n* ≤*500. Right frame: plots for the MEGNO filter Yn for the same errors.*

*<sup>n</sup> (blue line), reversibility errors e<sup>R</sup>*

*<sup>n</sup> (red line), and eREM*

*<sup>n</sup> (gray line), for 1*≤*n*≤*2500. Right*

*<sup>n</sup> (red line), and eREM*

*<sup>n</sup> (gray*

clearly visible.

*Progress in Relativity*

exponentially.

REM still fluctuates.

**Figure 5.**

**Figure 6.**

**174**

*Lyapunov error eL*

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

*<sup>n</sup> (blue line), reversibility errors eR*

*frame: plots for the MEGNO filter Yn for the same errors.*

*corresponding to a chaotic orbit. Lyapunov error e<sup>L</sup>*

*Left frame: variation for the standard map with λ* ¼ *0:1 of the errors LE (blue line), RE (red line), REM (gray line) computed at N* ¼ *1000 with the initial condition x0* ¼ *0, p0* ∈ ½ � �*0:15; 1; 0:15 . Right frame: the same for MEGNO YN.*

**Figure 8.**

*Left frame: variation for the standard map with λ* ¼ *0:8 of the errors LE (blue line), RE (red line), and REM (gray line), computed at N* ¼ *1000 with the initial condition x0* ¼ *0, p0* ∈½ � *0; 0:5 . Right frame: magnification in the interval p0* ∈ ½ � *0:25; 0:5 .*

In **Figure 8**, the same results are shown for a higher value of the parameter *λ* ¼ 0*:*8 at which the dynamical structure is rich due to the presence of many resonances and small chaotic regions. The effectiveness in detecting the resonances is evident.

#### **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 *Ng* � *Ng* points, where we compute the errors and show the result using a color scale. In order to analyze the details, smaller squares may be chosen eventually increasing the iteration number. In **Figure 9**, we show for *N* ¼ 500 and *Ng* ¼ 200 the color plots for the errors of the standard map with *λ* ¼ 0*:*8 and in **Figure 10** for *λ* ¼ 1*:*5. In the first case, the measure of chaotic orbits is small with respect to the regular ones. We observe that LE has some weak structures within the main regular component surrounding the origin, visible when the figure is sufficiently magnified. Such structures of LE, related to the oscillating growth with *n*, disappear when MEGNO is computed and

In normal coordinates, the errors grow as 2*J*∣Ω2∣*n<sup>α</sup>* where *<sup>α</sup>* <sup>¼</sup> 1 for LE and *<sup>α</sup>* <sup>¼</sup> <sup>3</sup>*=*<sup>2</sup> for RE. When the frequency attains a low resonant value, a chain of islands appears. Close to the separatrix *J* ¼ *Js*, the frequency vanishes as Ω � 1*=* log *J*ð Þ *<sup>s</sup>* � *J* and

In **Figure 11**, we show the variation of LE, RE, and REM computed at a fixed order *N* and after filtering them with MEGNO, when the initial conditions are chosen on a one-dimensional grid. The resonance 1*=*5 is met, as shown by the appearance of a large chain of islands, since Ωð Þ*J* is monotonically decreasing. The chaotic layer at the border of the islands chain is very thin so that LE and RE grow

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.

*Left frame: errors for the Hénon with ω* ¼ *0:21* � ð Þ *2π : LE (blue line), RE (red line), and REM (gray line) computed at iteration number <sup>N</sup>* <sup>¼</sup> *<sup>1000</sup> along the line <sup>x</sup>* <sup>¼</sup> *<sup>r</sup>* cos *<sup>α</sup>, p* <sup>¼</sup> *<sup>r</sup>* sin *<sup>α</sup> with <sup>α</sup>* <sup>¼</sup> *14o joining the origin with the center of the first of five islands. Center frame: computation of MEGNO with N* ¼ *1000 for LE (blue line), RE (red line), and REM (gray line). Right frame: phase portrait of the Hénon map. The initial conditions*

*Left frame: Hénon map with ω* ¼ *0:21* � ð Þ *2π color plot of LE in a logarithmic scale for N* ¼ *500 and a grid with Ng* ¼ *200. The white points belong to the unstable region beyond the dynamic aperture. Center frame:*

*color plot of RE in a logarithmic scale. Right frame: color plot of REM in a logarithmic scale.*

by approaching it, as for an integrable map with a separatrix.

*for the errors in the left and right frames are chosen on the red segment.*

The errors diverge by approaching the separatrix even though the power law growth does not change except on the separatrix itself. As a consequence, LE and RE can detect the separatrix. If the remainder in the normal form interpolating Hamiltonian is taken into account, then the separatrix becomes a thin chaotic region where the errors have an exponential growth and MEGNO rises linearly with *n*. The REM behaves as RE neglecting its fluctuations. The Hénon map, we consider here, has a linear frequency *ω=*2*π* ¼ 0*:*21 which is close to the resonance 1*=*5. As a consequence a chain of five islands appears before reaching the dynamic aperture, namely, the boundary of the stability region, beyond which the orbits escape to

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

as *Js* is approached, up to a logarithmic correction.

consequently Ω<sup>0</sup>

infinity.

**Figure 11.**

**Figure 12.**

**177**

ðÞ� *<sup>J</sup> <sup>J</sup>*ð Þ *<sup>s</sup>* � *<sup>J</sup>* �<sup>1</sup>

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

#### **Figure 9.**

*Left frame: standard map with λ* ¼ *0:8 color plot of LE in a logarithmic scale for N* ¼ *500 and a grid with Ng* ¼ *200. Center frame: standard map with λ* ¼ *0:8 color plot of RE in a logarithmic scale for N* ¼ *500 and a grid with Ng* ¼ *200. Right frame: color plot of REM in a logarithmic scale.*

**Figure 10.**

*Left frame: standard map with λ* ¼ *0:8 color plot of LE in a logarithmic scale for N* ¼ *500 and a grid with Ng* ¼ *200. Center frame: standard map with λ* ¼ *1:5 color plot of RE in a logarithmic scale for N* ¼ *500 and a grid with Ng* ¼ *200. Right frame: color plot of REM in a logarithmic scale.*

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<sup>10</sup> for LE and 10<sup>15</sup> for RE and REM. This choice is suggested by the asymptotic behavior *n<sup>α</sup>* of the error for regular orbits where *α* ¼ 1 for LE and *α* ¼ 3*=*2 for RE.

#### **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 *R*ð Þ *ω* . As a consequence the power law exponent of REM varies from 1 to 2, whereas the exponent for RE varies from 1/2 to 3/2. Within the main island, the variation range of the exponent for RE and REM is the same 1½ � *=*2*;* 3*=*2 . The behavior of LE and RE close to the origin is analytically obtained by using the normal forms. The frequency Ωð Þ*J* , from normal forms at the lowest order, reads

$$\begin{aligned} \, \, \Omega \simeq \omega + f \, \Omega\_2 \qquad \qquad \qquad \, \Omega\_2 = -\frac{1}{8} \left( 3 \cot \left( \frac{\alpha \nu}{2} \right) + \cot \left( \frac{3 \omega}{2} \right) \right) \end{aligned} \tag{50}$$

a formula valid for frequencies *ω=*2*π* not approaching the unstable resonances 0 and 1/3 where Ω<sup>2</sup> diverges.

In the normal coordinates ð Þ *X; P* , the behavior of errors is given by Eqs. (41) and (42). In the original coordinates ð Þ *x; p* , the error could be evaluated using Eq. (33).

*Fast Indicators for Orbital Stability: A Survey on Lyapunov and Reversibility Errors DOI: http://dx.doi.org/10.5772/intechopen.88085*

In normal coordinates, the errors grow as 2*J*∣Ω2∣*n<sup>α</sup>* where *<sup>α</sup>* <sup>¼</sup> 1 for LE and *<sup>α</sup>* <sup>¼</sup> <sup>3</sup>*=*<sup>2</sup> for RE. When the frequency attains a low resonant value, a chain of islands appears. Close to the separatrix *J* ¼ *Js*, the frequency vanishes as Ω � 1*=* log *J*ð Þ *<sup>s</sup>* � *J* and consequently Ω<sup>0</sup> ðÞ� *<sup>J</sup> <sup>J</sup>*ð Þ *<sup>s</sup>* � *<sup>J</sup>* �<sup>1</sup> as *Js* is approached, up to a logarithmic correction. The errors diverge by approaching the separatrix even though the power law growth does not change except on the separatrix itself. As a consequence, LE and RE can detect the separatrix. If the remainder in the normal form interpolating Hamiltonian is taken into account, then the separatrix becomes a thin chaotic region where the errors have an exponential growth and MEGNO rises linearly with *n*. The REM behaves as RE neglecting its fluctuations. The Hénon map, we consider here, has a linear frequency *ω=*2*π* ¼ 0*:*21 which is close to the resonance 1*=*5. As a consequence a chain of five islands appears before reaching the dynamic aperture, namely, the boundary of the stability region, beyond which the orbits escape to infinity.

In **Figure 11**, we show the variation of LE, RE, and REM computed at a fixed order *N* and after filtering them with MEGNO, when the initial conditions are chosen on a one-dimensional grid. The resonance 1*=*5 is met, as shown by the appearance of a large chain of islands, since Ωð Þ*J* is monotonically decreasing. The chaotic layer at the border of the islands chain is very thin so that LE and RE grow by approaching it, as for an integrable map with a separatrix.

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.

#### **Figure 11.**

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<sup>10</sup> for LE and 10<sup>15</sup> for RE and REM. This choice is suggested by the asymptotic behavior *n<sup>α</sup>* of the error

*Left frame: standard map with λ* ¼ *0:8 color plot of LE in a logarithmic scale for N* ¼ *500 and a grid with Ng* ¼ *200. Center frame: standard map with λ* ¼ *1:5 color plot of RE in a logarithmic scale for N* ¼ *500 and a*

*Left frame: standard map with λ* ¼ *0:8 color plot of LE in a logarithmic scale for N* ¼ *500 and a grid with Ng* ¼ *200. Center frame: standard map with λ* ¼ *0:8 color plot of RE in a logarithmic scale for N* ¼ *500 and a*

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 *R*ð Þ *ω* . As a consequence the power law exponent of REM varies from 1 to 2, whereas the exponent for RE varies from 1/2 to 3/2. Within the main island, the variation range of the exponent for RE and REM is the same 1½ � *=*2*;* 3*=*2 . The behavior of LE and RE close to the origin is analytically obtained by using the normal forms. The frequency Ωð Þ*J* , from normal forms at the lowest order, reads

<sup>8</sup> 3 cot *<sup>ω</sup>*

a formula valid for frequencies *ω=*2*π* not approaching the unstable resonances 0

In the normal coordinates ð Þ *X; P* , the behavior of errors is given by Eqs. (41) and (42). In the original coordinates ð Þ *x; p* , the error could be evaluated using Eq. (33).

2 

þ cot

3*ω* 2

(50)

for regular orbits where *α* ¼ 1 for LE and *α* ¼ 3*=*2 for RE.

*grid with Ng* ¼ *200. Right frame: color plot of REM in a logarithmic scale.*

*grid with Ng* ¼ *200. Right frame: color plot of REM in a logarithmic scale.*

<sup>Ω</sup> <sup>≃</sup>*<sup>ω</sup>* <sup>þ</sup> *<sup>J</sup>*Ω<sup>2</sup> <sup>Ω</sup><sup>2</sup> ¼ � <sup>1</sup>

**6. The Hénon map**

**Figure 9.**

*Progress in Relativity*

**Figure 10.**

and 1/3 where Ω<sup>2</sup> diverges.

**176**

*Left frame: errors for the Hénon with ω* ¼ *0:21* � ð Þ *2π : LE (blue line), RE (red line), and REM (gray line) computed at iteration number <sup>N</sup>* <sup>¼</sup> *<sup>1000</sup> along the line <sup>x</sup>* <sup>¼</sup> *<sup>r</sup>* cos *<sup>α</sup>, p* <sup>¼</sup> *<sup>r</sup>* sin *<sup>α</sup> with <sup>α</sup>* <sup>¼</sup> *14o joining the origin with the center of the first of five islands. Center frame: computation of MEGNO with N* ¼ *1000 for LE (blue line), RE (red line), and REM (gray line). Right frame: phase portrait of the Hénon map. The initial conditions for the errors in the left and right frames are chosen on the red segment.*

#### **Figure 12.**

*Left frame: Hénon map with ω* ¼ *0:21* � ð Þ *2π color plot of LE in a logarithmic scale for N* ¼ *500 and a grid with Ng* ¼ *200. The white points belong to the unstable region beyond the dynamic aperture. Center frame: color plot of RE in a logarithmic scale. Right frame: color plot of REM in a logarithmic scale.*
