# Analysis of Round Off Errors With Reversibility Test as a Dynamical Indicator

###### Abstract

We compare the divergence of orbits and the reversibility error for discrete time dynamical systems. These two quantities are used to explore the behavior of the global error induced by round off in the computation of orbits. The similarity of results found for any system we have analysed suggests the use of the reversibility error, whose computation is straightforward since it does not require the knowledge of the exact orbit, as a dynamical indicator.

The statistics of fluctuations induced by round off for an ensemble of initial conditions has been compared with the results obtained in the case of random perturbations. Significant differences are observed in the case of regular orbits due to the correlations of round off error, whereas the results obtained for the chaotic case are nearly the same.

Both the reversibility error and the orbit divergence computed for the same number of iterations on the whole phase space provide an insight on the local dynamical properties with a detail comparable with other dynamical indicators based on variational methods such as the finite time maximum Lyapunov characteristic exponent, the mean exponential growth factor of nearby orbits and the smaller alignment index. For 2D symplectic maps the differentiation between regular and chaotic regions is well full-filled. For 4D symplectic maps the structure of the resonance web as well as the nearby weakly chaotic regions are accurately described.

## 1 Introduction

The discrete time dynamical systems are widely used because
they present a rich structure (regular, chaotic, intermittent and other
types of orbits) in one or two dimensions and because the
numerical evaluation of the orbits is straightforward. A large number of dynamical tools known as indicators of stability allow us
to improve our understanding of dynamical systems: Lyapunov Characteristic Exponents (LCEs) are known from a long time
(Wolf et al., 1985), (Rosenstein et al., 1993), (Skokos, 2010)
as well as Return Times Statistics (Kac, 1934), (Gao, 1999), (Hu et al., 2004) , (Buric et al., 2005). In the recent past many other indicators have been introduced not only to address the
same problem of quantifying the degree of chaoticity of an orbit but also to perform the task fast.
The Smaller Alignment Index (SALI), widely described in Skokos et al. (2002) and Skokos et al. (2004),
allows to discriminate regular from chaotic orbits. Similarly, the Generalized Alignment Index (GALI), introduced
in Skokos et al. (2007), is a family of highly efficient algorithms. Besides, the Mean Exponential Growth factor of
Nearby Orbits (MEGNO) discussed in Cincotta et al. (2003) and Goździewski et al. (2001) is a quantity
that gives a fast identification of the chaoticity of the orbit while its average slope estimates
the maximum LCE (mLCE), see Maffione et al. (2011) for a test of the MEGNO.
Fidelity and correlation decay are also tools that can be successfully used to
characterise stability properties as explained in Vaienti et al. (2007) and Turchetti et al. (2010b) as well as
Frequency Map Analysis (Laskar, 1999), (Robutel & Laskar, 2001).

As it is already known, in any numerical computation of a given trajectory, there is a round off error, and it would
be interesting to study its relation with the chaoticity of the orbit and if possible determine
its effect with the help of some dynamical indicator.

The error between an exact and a numerical orbit is due to the finite
precision used to represent real numbers and to the arithmetics with
round off. This is unavoidable because the length of the binary strings
representing real numbers must not change after arithmetic operations.
The shadowing lemma is often invoked to state the existence of a true orbit
close to a numerical orbit for chaotic systems (Katok & Hasselblatt, 1997), (Hammel et al., 1987), (Chow & Palmer, 1992), (Chow & Palmer, 1991). However it does not provide
information on error growth for a given numerical orbit and the case
of regular numerical orbits is not covered by such a lemma.
The global error between the exact and the numerical orbit
is unknown because the first one is not computable. Nevertheless an
estimate can be provided by replacing the exact orbit with another one
having very high accuracy. If the map is invertible the reversibility
error can be computed without any reference to the
exact orbit.
Both errors have a similar behavior, namely an average linear growth
for regular orbits and an average exponential growth for chaotic orbits and
consequently can be used as dynamical indicators.
Due to the correlation between the single step errors, there is a substantial
difference with respect to the case in which the exact system is perturbed with
random uncorrelated noise.

This difference has been analysed in Turchetti et al. (2010b)
by using the fidelity which measures the deviation of the orbits of a given map and its
perturbation by integrating over all the initial conditions with the appropriate measure.
In the case of regular orbits with random perturbations the decay of fidelity is exponential
whereas with round off errors the decay follows a power law. In the case of chaotic orbits
the asymptotic limit is approached super-exponentially for both situations (random uncorrelated
perturbations and round off). The symplectic maps of physical interest are generally provided by the composition of a linear map with another one whose generating function is the identity plus a function of position or momentum only. In this case the
inversion is immediately obtained in analytic form.
As most of the symplectic maps used in the literature are of this kind,
the condition of invertibility of the map is not too restrictive. In addition, the reversibility error can also be applied to
time-reversible systems of differential equations.

Here we present an analysis of the reversibility error and its
comparison with the divergence of orbits due to round off and other dynamical indicators of stability such as
SALI, MEGNO and a finite time numerical estimation of the mLCE.
Moreover we talk about the similarities and differences between an irreversibility
due to the single precision round off and one due to the application of uncorrelated random noise
in an orbit iterated with double precision.
For two dimensional area preserving maps the reversibility error for a
fixed number of iterations detects the
various regions of phase space with different stability properties quite
effectively as well as other dynamical indicators.
Besides, the reversibility error allows to study the structure of the resonance web of a
four dimensional symplectic map.

## 2 Round off error methods

In a computational device a real number can be represented by a floating-point number , that according to Goldberg (1991) and working with base , can be written as

(1) |

where is called the significand and has binary digits , whose value is or , and where the exponent is an integer that satisfies . The IEEE 754 standard states that for single precision , and while for double precision , and .

Consequently, a floating point number differs from the real number it represents and the relative error , defined by , satisfies , as analysed by Goldberg (1991) and Knuth (1973). Therefore, according to IEEE 754 we have that and for single and double precision which roughly corresponds to 7 and 16 decimal digits, respectively. The arithmetic operations such as sums or multiplications imply a round off, which propagates the error affecting each number. Round off algebraic procedures are hardware dependent as detailed in Knuth (1973). Unlike the case of stochastic perturbations, the error strongly depends on . Suppose we are given a map then the error with respect to the numerical map after the first step is defined by:

(2) |

Analogously, we define the local error produced in the step by where .

The global error

(3) |

accumulates all the local errors and explicit expressions can be written at first order in . In the example of a regular map we take the translation on the torus defined as:

(4) |

so and the global error, which includes also the error to the mod 1 operation, becomes:

(5) |

where is a time average defined as the limit of for and is a bounded fluctuation.

For the chaotic Bernoulli map

(6) |

we have that and that the single step and global error satisfy, respectively, and .

In order to compute the global error we need the knowledge of the exact map, which is usually precluded. A practical way to overcome this difficulty is to replace the exact map with a map computed with an accuracy where . For instance, as corresponds to single precision one might choose corresponding to double precision. If then one might choose and so forth. There are available libraries which allow to compute with any fixed number of bits or significant decimal digits. The computation of the reference orbit is expensive if high precision is used, but there is no other way to evaluate the global error. As a consequence the “exact" orbit is achievable for a definite number of iterations which depends on and the nature of the map. Taking into account what we have just mentioned, in the forthcoming numerical experiments, we will use the divergence of orbits, defined by:

(7) |

where and stand for single and double precision iterations respectively.

If the map is invertible there is another option to overcome the difficulty of not possessing the true map. We define the reversibility error as

(8) |

which is non zero since the numerical inverse of the map
is not exactly the inverse of namely .
Obviously the reversibility error is much easier to compute than the divergence
of orbits (if we know explicitly the inverse map) and the information it provides
is basically the same as the latter.
Both quantities give an average linear growth for a regular map together with an
exponential growth for a chaotic map having positive Lyapounov exponents and strong mixing properties. When computing we will set
in order to compare with .

## 3 Variational methods

In the forthcoming sections we will compare the performance of the indicators presented above with three well known and widely accepted dynamical indicators that are based on the behavior of the solution of the variational equations of the system. These are the finite time mLCE, the cumulative moving time average of MEGNO, and SALI.

Let us briefly state them for discrete time dynamical systems of the form:

(9) |

where is the state vector at time and is a vector valued function.

The concomitant discrete time variational equations, also called tangent map dynamics, associated to a given orbit are the following:

(10) |

where is the Jacobian matrix of the function and is a deviation vector at time .

Skokos (2010) makes a historical review of the definition of the LCEs and its connection with the divergence of nearby orbits. He also states the theorems that guarantee the existence of the spectrum of LCEs and, in particular, the current definition of the mLCE in terms of the solution of the variational equations:

(11) |

with some norm. For a chaotic orbit the mLCE is positive and this implies an exponential divergence of nearby orbits. On the other hand, for regular orbits mLCE is zero.

In order to have a numerically computable quantity we define the finite time mLCE at time as

(12) |

so that equation Eq. (11) can be reformulated as

(13) |

Cincotta et al. (2003) defines a biparametric family of MEGNO indicators:

(14) |

where and are integer numbers. They made experiments with , and and concluded that the last one allows both a fast classification between chaotic and regular orbits, and a clear identification of stable and unstable periodic orbits. Due to this fact, we will use throughout the rest of this article. In order to reduce the fast oscillations that the time evolution of the MEGNO presents, in Cincotta et al. (2003) they use a time average of this quantity, namely:

(15) |

Theoretically, the asymptotic evolution of for any dynamical regime can be put into a single expression:

(16) |

where and for chaotic and regular motion, respectively.

We also compare our results with the SALI, introduced by Skokos (2001), that measures the degree in which a pair of initially linearly independent deviation vectors tend to become aligned. The underlying principle is that, for a chaotic orbit, a deviation vector under the tangent map dynamics changes in order to become aligned with the instantaneous most unstable direction. In other words, for almost every pair of initial deviation vectors (, ), the more chaotic the orbit the faster that the angle between them will reduce to zero. In the case of regular orbits, the behavior depends strongly on the dimensionality of the map: in maps with dimension the deviation vectors generally remain unaligned so the SALI tends to a positive non-zero value, while in 2D maps the two deviation vectors tend to align with a time rate that follows a power law. See (Skokos et al., 2004) and (Skokos et al., 2007) for numerical tests of the SALI and its generalisation, the GALI family. Denoting the Euclidean norm with the SALI is defined as:

(17) |

To control the exponential increase of the norm of the vectors and avoid overflow problem, Skokos et al. (2004) have normalised them, at every time step, keeping their norm equal to 1. Defined like this SALI and SALI = 0 if and only if the two normalised vectors have the same direction, being equal or opposite.

## 4 Error behavior on simple maps

From the linear and exponential bounds on the errors it follows that the numerical orbit remains close to the exact one for a number of iterations proportional to in the regular case and to in the chaotic one.

We consider two types of models where the error grows linearly and exponentially respectively. The first one is the translation on the torus defined by Eq. (4).

This is equivalent to the rotations on the unit circle defined by the map

(18) |

The correspondence between the sequences and is evident after writing . In spite of the rigorous mathematical equivalence between the translations on the torus and the rotations on the unit circle, there is not such similitude in the numerical evaluation of these maps, as will be explained in the following paragraph.

For map (4) the divergence of orbits and the reversibility
error have a linear growth, basically due to the the fact that ,
as shown theoretically for the global error in Eq. (5). However there
are architectures and/or compilers in which the reversibility error may be zero for this simple map
that only requires the computation of sums and modulus operations.
Without knowing the precise implementation of these peculiar operations it is impossible to
establish a priori with which compiler and in which architecture we may obtain this result.
We believe this anomaly is due to the peculiar arithmetic operations involved
and that it does not occur for a generic map. To support this claim we
have checked that the reversibility error never vanishes for the map (18)
which involves multiplications and the evaluation of trigonometric
functions.

In Figure 1-left we compare the divergence and reversibility error for
the torus translation with and .
We can see that both errors approximately satisfy the same expression: in log-log scale it is close to a straight
line with unitary slope. This implies that both quantities are linear functions of time.
We have checked that the behavior is similar for almost every initial condition and frequency but in some cases the fluctuations around the average linear growth have a larger amplitude.
Figure 1-right shows that for the map (18) there is an analogous
behavior. However, in this case the straight line that corresponds to orbit divergence is a bit shifted
upwards with respect to the one of reversibility error. This is telling us that the orbit divergence’s linear
growth has a bigger slope.

The fluctuations have variance computed respect to an ensemble of initial conditions which first grows linearly, but rapidly saturates to a constant value with small oscillations ( see Figure 7 in (Turchetti et al., 2010b) ). In this respect it is quite different with respect to a random perturbation where are independent variables, since in this case we have and the variance is given by for any value of . This different behavior is reflected in the decay of fidelity (Vaienti et al., 2007) which follows a power law for round off errors and an exponential law for random perturbations (Turchetti et al., 2010a) , (Turchetti et al., 2010b). This shows that the round off errors decorrelate very slowly unlike the random errors which are uncorrelated.

For chaotic maps, characterised by an exponential increase in the distance between two nearby orbits, the reversibility error due to round off has the same growth. Bounded two dimensional chaotic maps, like the cat map, exhibit this behaviour during a short time before reaching saturation.

## 5 The case of the standard map

As an example of generic 2D map we have chosen the standard map in a torus:

(19) |

For very low values of the divergence of orbits has an average growth linear with as for the translations on the torus and they depend weakly on the initial condition. As is increased a power law with is observed and the dependence on the initial conditions becomes appreciable. Reversibility error also shows this behavior for .

For approaching one we have coexistence of regular and chaotic regions so that the domain is splitted into several islands of stable orbits and a chaotic sea. In figures 2, 3, 4, 5 and 6 we show, respectively, the value of the reversibility error (8), the orbit divergence (7), the finite time mLCE (12), the time average of the MEGNO (15) and the SALI (17) obtained iterating eq. (19).

The standard map with has been iterated times to compute all the dynamical indicators. Even if might be used to describe the chaotic region, the chosen value allows to highlight differences within regular regions where the growth is very slow.

The pictures in the left side show the value of the dynamical indicators in a chromatic scale for a grid of initial conditions in the two dimensional torus whereas the right pictures show the corresponding graph for a fixed value of the action variable () corresponding to the horizontal line in the left figures. Except for mLCE we have used the natural logarithm of the absolute value for all the dynamical indicators.

Figure 2 and Figure 3 were obtained by taking into account the error only in the action variable. It is evident that both the reversibility error and the orbits divergence (with respect to the exact one) discriminate regular from chaotic orbits. Minor differences exist within the islands of stability: the orbit divergence approaches the minimum value close to the centre, whereas the minimum of the reversibility error occurs on a line crossing the islands possibly due to a mechanism similar to the one reported by Barrio et al. (2009) which analysed spurious errors for variational methods.

The plots obtained for mLCE and the MEGNO (figures 4 and 5 respectively) show no structure within the resonance islands whereas variations and fluctuations appear in the chaotic regions as for the previous indicators.

Figure 6 shows the logarithm of the SALI value. Due to the fact that this indicator converges to zero extremely fast for chaotic orbits, we have used a cut off value at SALI . The presence of this cut-off erases any structure within the strongly chaotic region.

To summarise, all the indicators discriminate regular and chaotic regions but their sensitivity in these regions is different .

Another aspect to consider when comparing the efficiency of chaos indicators is the computational cost. Each variational method needs to iterate both the map and the tangent map forward for steps. The tangent map is the computationally most expensive since it needs the evaluation of the Jacobian matrix at every step. When computing SALI, two deviations vectors must be simultaneously computed. We also use two deviation vectors in evaluating MEGNO, selecting at every step the one that stretches more, in order to reduce the probability of having a vector almost orthogonal to the most unstable direction. In the case of the computation of the mLCE we have used the orthogonalization algorithm developed by Benettin et al. (1980). On the other hand, we notice that the reversibility error method requires only the iteration of the map whereas the orbit divergence method requires the iteration of the single and double precision (or double and higher precision) map. As a consequence, the computationally most economic method is the one based on the reversibility error which does not require any algorithm except the evaluation of the map itself. Typically, the relevant information can be extracted from a computation in single precision.

From now on, within this section, we will focus on some ensemble quantities in order to compare the effect of round off errors with the effect of random perturbations in the standard map. For the standard map with all the orbits are regular and we follow the evolution of an ensemble of 10001 initial conditions randomly chosen in . For each iteration we compute the variance of in action () and angle () variables (see in Figure 7-left) and compare it against the same quantities of a double precision orbit stochastically perturbed with a uniform uncorrelated noise of amplitude (shown in Figure 7-right). In the latter case we have found that and grow according a power law with exponents one and three respectively up to numerical uncertainties.

It is interesting to compare the two pictures in Figure 7 with each other. For large we see that the behaviors of variances due to round off and random perturbations are similar. The presence of a transient for the reversibility case is very likely due to the initial presence of correlation.

In the fully chaotic regime for the standard map the random perturbation produces very similar results to the round off error as shown in Figure 8 for since in presence of a chaotic dynamics the round off error correlations are lost.

The results for the standard map in the small regime can be compared to the variances for the skew map:

(20) |

to which the standard map in equation (19) reduces for . If the skew map is randomly perturbed:

(21) |

where and are random uncorrelated variables the growth of the variances and follows a linear and cubic law respectively (Turchetti et al., 2006). The random perturbation of a standard map with a small value of (shown in Figure 7-right) shows exactly the same growth.

In the case of the round off error the behavior of the skew map with respect to the standard map with a small value of is quite different. Indeed, the round off error affects the skew map just as the translation on the 1-D Torus: it was observed that the global error grows linearly and the variance saturates at a very small value with respect to the size of the torus (see Figure 7 in Turchetti et al. (2010a)). For the standard map, the coupling between action and angle, even for very small , causes a growth of the variance of the fluctuations due to round off as shown in Figure 7-left. As a consequence the effect of round off in a very weakly perturbed standard map is similar to a random perturbation.

## 6 A 4D Map

In this section we show how either the reversibility error or the divergence of orbits can be used to analyze the resonance structure of four dimensional non integrable maps. As an example we consider a symplectic nearly integrable map extensively used in the literature (see Guzzo & Lega (2004)). This map is defined as:

(22) |

where , with and is the perturbation parameter.
We have taken a grid of initial conditions with actions
and a fixed pair of angles, namely and computed , , mLCE(n), and
for .
Associating these values to each initial condition while using a chromatic scale we have performed
figures 9, 10 and 11, where the parameter values
and have been used.
As we did in the previous section, in all except for the mLCE we have used the natural logarithm of the
absolute value of the concomitant indicator and again we have used a cut off value at SALI .
Figure 9 was done taking into account the Euclidean error in only the action plane.
As it happened for the standard map both reversibility error and orbit divergence present the same order
of magnitude. The resonance web appears in every figure and its structure is the same. From a qualitative
viewpoint no substantial differences are found and one might conclude that the dynamical information extracted
from the reversibility error is the same as for the other four dynamical indicators we have considered. In the case of mLCE, the highest values are located on the
diagonal even if they are not well visible in the plot. This is anyway coherent with the fact that on
the diagonal we find the highest instability.

The frequency map analysis was not presented for comparison because it is computationally heavier. The error reaches its highest values in small neighbourhoods of single resonance lines, because it detects the perturbed separatrices of the
pendulum models that one could associate with the resonances, and in relatively large neighbourhoods of the intersections of resonances, where the dynamics are widely non integrable, as shown by the computation
of the interpolating Hamiltonian.
These results are certainly not exhaustive but show that the behavior of the reversibility error is strictly related
to the divergence of orbits and consequently it is very weak in the integrable regions where it does not have a diffusive
character as for a random error.

## 7 Conclusions

We have examined the orbit divergence and the reversibility error in order to determine the effect of round off error for invertible maps. The knowledge of the exact map is not required to compute the reversibility error and the results obtained are about the same with respect to the case in which we compute the orbit divergence.

By choosing an ensemble of initial data we have examined the statistics of the fluctuations due to round off with respect to the time average of the error. There are different behaviours according to the degree of chaoticity that characterises the ensemble. For chaotic orbits the variances have an exponential growth similar to the one observed for random perturbations and correspondingly the decay of fidelity is super-exponential in both cases. For regular or quasi-regular orbits differences are observed between the effects generated by round off errors and random perturbations. For a quasi-integrable map the random perturbations produce a growth of variances linear and cubic for actions and angles, respectively. On the other hand, the round off errors produce an initial transient that possibly lasts the time that the round off errors of the ensemble orbits need to decorrelate. After this initial transient, the behaviour of variances affected by round off and noise are similar to each other. Finally, for a regular map such as the translation of the 1D torus, the round off variance, after a linear growth, saturates before the distribution of errors fills all the torus, unlike in the randomly perturbed case where the growth is always linear until saturation of the full torus. The fidelity has a power law decay for round off whereas it decays exponentially for random perturbations. This is a clear signature of the correlation between the errors due to round off.

To conclude, the reversibility error provides basically the same information as the divergence of orbits and it is easily accessible from a computational view point. This is due to the facts that it does not require the solution of the variational equations and that both the forward and backward orbits can be computed using single precision. When the initial point is varied and the number of iterations is kept fixed the reversibility error due to round off provides an insight on the dynamical structure of the map. For the standard map the reversibility error provides a picture of the dynamical behavior of the map where not only large scale features but also small scale details can be detected. In the case of a 4D symplectic map the reversibility error in action space provides a similar picture where the resonance web and the nearby regions of weakly chaotic motions can be easily highlighted.

Even if no really new information with respect to the standard indicators is provided by the reversibility error, we point out that this type of analysis takes into account not only the dynamics of the map but also the unavoidable effect of finite accuracy due to round off. Moreover, by increasing the accuracy of numerical computations, the time on which the reversibility test is computed can be increased and finer details on the phase space structure can be detected.

## 8 Aknowledgements

We are gratefull to the anonymous reviewers for their important comments. The stay of MM at the Physics Department of the University of Bologna was fully supported by a grant from the Erasmus Mundus External Cooperation Window Lot 16 Programme, EADIC, financed by the European Commission. Besides, MM acknowledges financial support of a grant from the Consejo Nacional de Investigaciones Científicas y Técnicas de la República Argentina (CONICET). Also, he is grateful to A. Bazzani for his helpfull discussions. DF acknowledges the financial support of the EU FP7-ERC project NAMASTE Thermodynamics of the Climate System.

## References

- Barrio et al. (2009) Barrio, R., Borczyk, W. & Breiter, S. [2009] “Spurious structures in chaos indicators maps,” Chaos. Soliton. Fract. 40, 1697–1714.
- Benettin et al. (1980) Benettin, G., Galgani, L., Giorgilli, A. & Strelcyn, J. M. [1980] “Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part 1: Theory,” Meccanica 15, 9–20.
- Buric et al. (2005) Buric, N., Rampioni, A. & Turchetti, G. [2005] “Statistics of Poincaré recurrences for a class of smooth circle maps,” Chaos. Soliton. Fract. 23, 1829–1840.
- Chow & Palmer (1991) Chow, S. N. & Palmer, K. J. [1991] “On the numerical computation of orbits of dynamical systems: the one-dimensional case,” J. Dyn. Differ. Equ. 3, 361–379.
- Chow & Palmer (1992) Chow, S. N. & Palmer, K. J. [1992] “On the numerical computation of orbits of dynamical systems: the higher dimensional case,” J. Complexity 8, 398–423.
- Cincotta et al. (2003) Cincotta, P. M., Giordano, C. M. & Simó, C. [2003] “Phase space structure of multi-dimensional systems by means of the mean exponential growth factor of nearby orbits,” Physica D 182, 151–178.
- Gao (1999) Gao, J. B. [1999] “Recurrence time statistics for chaotic systems and their applications,” Phys. Rev. Lett. 83, 3178–3181.
- Goldberg (1991) Goldberg, D. [1991] “What every computer scientist should know about floating point arithmetic,” ACM Comput. Surv. 23, 5–48.
- Goździewski et al. (2001) Goździewski, K., Bois, E., Maciejewski, A. J. & Kiseleva-Eggleton, L. [2001] “Global dynamics of planetary systems with the MEGNO criterion,” Astron. Astrophys. 378, 569–586.
- Guzzo & Lega (2004) Guzzo, M. & Lega, E. [2004] “First numerical evidence of global Arnold diffusion in quasi–integrable systems,” Discrete Cont. Dyn-B 5, 687–698.
- Hammel et al. (1987) Hammel, S. M., Yorke, J. A. & Grebogi, C. [1987] “Do numerical orbits of chaotic dynamical processes represent true orbits?” J. Complexity 3, 136–145.
- Hu et al. (2004) Hu, H., Rampioni, A., Rossi, L., Turchetti, G. & Vaienti, S. [2004] ‘‘Statistics of Poincaré recurrences for maps with integrable and ergodic components,” Chaos 14, 160.
- Kac (1934) Kac, M. [1934] “On the notion of recurrence in discrete stochastic processes,” Proc. Nat. Acad. Sci. USA 20, 376–379.
- Katok & Hasselblatt (1997) Katok, A. & Hasselblatt, B. [1997] Introduction to the modern theory of dynamical systems (Cambridge Univ Pr).
- Knuth (1973) Knuth, D. E. [1973] The art of computer programming, Vol. 2.
- Laskar (1999) Laskar, J. [1999] “Introduction to frequency map analysis,” Hamiltonian systems with three or more degrees of freedom, pp. 134–150.
- Maffione et al. (2011) Maffione, N. P., Giordano, C. M. & Cincotta, P. M. [2011] “Testing a fast dynamical indicator: The MEGNO,” Int. J. Nonlin. Mech. 46, 23–34.
- Robutel & Laskar (2001) Robutel, P. & Laskar, J. [2001] “Frequency Map and Global Dynamics in the Solar System I:: Short Period Dynamics of Massless Particles,” Icarus 152, 4–28.
- Rosenstein et al. (1993) Rosenstein, M. T., Collins, J. J. & De Luca, C. J. [1993] “A practical method for calculating largest Lyapunov exponents from small data sets,” Physica D 65, 117–134.
- Skokos (2001) Skokos, C. [2001] “Alignment indices: a new, simple method for determining the ordered or chaotic nature of orbits,” J. Phys. A-Math. Gen. 34, 10029–10043.
- Skokos (2010) Skokos, C. [2010] “The lyapunov characteristic exponents and their computation,” Lect. Notes Phys., Berlin Springer Verlag, pp. 63–135.
- Skokos et al. (2002) Skokos, C., Antonopoulos, C., Bountis, T. & Vrahatis, M. [2002] “Smaller alignment index (SALI): Detecting order and chaos in conservative dynamical systems,” Proc. of 4th GRACM Congress on Computational Mechanics (Citeseer).
- Skokos et al. (2004) Skokos, C., Antonopoulos, C., Bountis, T. C. & Vrahatis, M. [2004] “Detecting order and chaos in Hamiltonian systems by the SALI method,” J. Phys. A-Math. Gen. 37, 6269.
- Skokos et al. (2007) Skokos, C., Bountis, T. C. & Antonopoulos, C. [2007] “Geometrical properties of local dynamics in Hamiltonian systems: The Generalized Alignment Index (GALI) method,” Physica D 231, 30–54.
- Turchetti et al. (2006) Turchetti, G., Bassi, G., Bazzani, A., Giorgini, B. & Mais, H. [2006] “Hamiltonian dynamics with a weak noise and the echo effect for the rotator model,” J. Phys. A-Math. Gen. 39, 11417–11440.
- Turchetti et al. (2010a) Turchetti, G., Vaienti, S. & Zanlungo, F. [2010a] “Asymptotic distribution of global errors in the numerical computations of dynamical systems,” Physica A 389, 4994–5006.
- Turchetti et al. (2010b) Turchetti, G., Vaienti, S. & Zanlungo, F. [2010b] “Relaxation to the asymptotic distribution of global errors due to round off,” Europhys. Lett. 89, 40006.
- Vaienti et al. (2007) Vaienti, S., Liverani, C. & Marie, P. [2007] “Random Classical Fidelity,” J. Stat. Phys. 128, 1079.
- Wolf et al. (1985) Wolf, A., Swift, J. B., Swinney, H. L. & Vastano, J. A. [1985] “Determining Lyapunov exponents from a time series,” Physica D 16, 285–317.