Open access peer-reviewed chapter - ONLINE FIRST

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

By Giorgio Turchetti and Federico Panichi

Submitted: January 27th 2019Reviewed: June 16th 2019Published: July 25th 2019

DOI: 10.5772/intechopen.88085

Downloaded: 89

Abstract

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

Keywords

  • variational principles
  • reversibility error
  • additive noise
  • roundoff

1. Introduction

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

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

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

For an integrable map, the growth of LE and RE follows asymptotically a power law nα, 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λn. 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 α=1with oscillations, whereas RE has an asymptotic power law growth with α=3/2without 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 α=0for LE and α=1/2for RE. For REM the power law exponent αvaries between 0 and 1, its value depending on the computational complexity of the map and therefore on the choice of coordinates.

The definition of LE we propose differs from fast Lyapunov indicator (FLI) [3] or orthogonal fast Lyapunov indicator (OFLI) [4], which are based on the growth along the orbit of the norm of a given initial displacement vector. Indeed, we compute the growth of the vectors of an orthogonal basis, which amounts to defining LE, which we denote as enL, as the square root of TrDMnx0TDMnx0where Mxis the map, DMxdenotes the tangent map, and x0is the initial condition. This definition has the obvious advantage of insuring the correct asymptotic growth.

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

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

2. Definition of errors

Given a symplectic map Mxwhere xR2d, we consider the orbits xn=Mnx0and yn=Mny0for two initial points x0and y0=x0+ϵη0, respectively, where ηis a unit vector. We consider the normalized displacement ηnat iteration ndefined by

ηn=limϵ0ynxnϵ=limϵ0Myn1Mxn1ϵE1

which satisfies the linear recurrence

ηn=DMxn1ηn1DMij=MixjE2

where DMis the tangent map. For any finite ϵ, we have yn=xn+ϵηn+Oϵ2. We might define the error as the norm of ηnwhich is closely related to the fast Lyapunov indicator FLI (see [3]) as

enη0=ηn=DMnx0η0FLIn=max1knlogekη0E3

and to its variants such as OFLI [4]. The mean exponential growth factor of nearby orbits, MEGNO [17], denoted by Yn=Yenis the double average of the slope, and we denote it as Δen2. When nis a continuous variable, then Δen2=dlogen2/dlogn. When nis an integer, the standard definition is

Δen2=nlogen2logen12Yn=Δen2.E4

2.1. Lyapunov error

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

enL=TrAnTAn1/2An=DMnx0=DMxn1An1A0=IE5

It is immediate to verify that given an orthonormal basis η0k, we have

enL=k=12den2η0k1/2E6

and obviously the result does not depend on the choice of the basis. The computational cost of enη0is 2dtimes higher with respect to enL, but this difference is negligible with respect to the computational cost of the matrix DMxn1, which recursively gives ηnand An. 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 xnis replaced by the noisy one (yn) having the same initial condition:

yn=Myn1+ϵξny0=x0E7

where ξnare independent random vectors satisfying

ξn=0ξnξmT=IδnmE8

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

Ξn=limϵ0ynxnϵΞn=DMxn1Ξn1+ξnE9

with initial condition Ξ0=0. Letting Σn2F=ΞnΞnTbe the covariance matrix the forward error is defined by

enF=ΞnΞn1/2=TrΣn2F1/2E10

The explicit solution for Ξnis given by

Ξn=k=1nDMnkxkξk=k=0n1BkξnkBk=DMkxnkE11

where Bkcan be evaluated recursively as

Bk=Bk1DMxnkB0=I.E12

The expression for the forward error finally reads

enF=k=0n1TrBkTBk1/2.E13

The computation cost of Bnis negligible, once we have evaluated the tangent map, but the storage of the tangent map along the orbit up to nis required.

2.3. Reversibility error

We have just defined the forward error, but it will not be used, because it is only an intermediate step toward the definition of the reversibility error. Consider the backward evolution yn,k, given by the inverse map M1, with initial point yn,0=yn. The point ynis reached by iterating ntimes the map Mwith a random displacement of amplitude ϵat each step, starting from y0=x0(see the previous subsection). The orbit yn,kis obtained by iterating ktimes the map M1with a random displacement of the same amplitude at each step:

yn,k=M1yn,k+1+ϵξkk=1,,nE14

The random backward displacements ξkare independent from the forward displacements ξk', namely, ξkξk'=0and ξkξk'=Iδkk', for any k,k'>0. We consider then the stochastic process Ξn,kdefined by

Ξn,k=limϵ0yn,kxnkϵΞn,k=DM1xnk+1Ξn,k+1+ξkE15

with initial condition Ξn,0=Ξn. The solution of the recurrence reads

Ξn,k=DMkxnΞn+j=1kDMkjxnjξjE16

For k=n, we obtain the global normalized displacement ΞnR=Ξn,nafter nforward and nbackward iterations with noise of vanishing amplitude:

ΞnR=DMnxnΞn+j=1nDMnjxnjξjE17

Letting Σn2R=ΞnRΞnRTbe the covariance matrix, the reversibility error (RE) is defined by

enR=ΞnRΞnR1/2=TrΣn2R2E18

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

ΞnR=k=1nDMnxnDMnkxkξk+k=1nDMnkxnkξk.E19

Taking into account the independence of ξkand ξk', the expression for the reversibility error enR=ΞnRΞnR1/2is 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

DMnxnDMnkxk=DMkxkE20

We prove this relation by writing MnMnkx=Mkx, computing the tangent map DMnMnkxDMnkx=DMkx, and evaluating it for x=xk. As a consequence the expression for the reversibility error becomes

enR2=TrΞnRΞnRT=k=1nTrDMkxkTDMkxk++TrDMnkxnkTDMnkxnk==2k=1n1TrDMkxkTDMkxk++TrDMnxnTDMnxn+TrIE21

Starting from MkMkx=x, computing the tangent map, and evaluating it at x=x0, it follows that

DMkxk=DMkx01E22

Given any symplectic matrix L3, we can prove that

TrL1TL1=TrLTL.E23

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

TrDMkxkTDMkxk=TrDMkx0TDMkx0=ekL2.E24

Finally, the relation between LE and RE is given by

enR2=k=1nekL2+enkL2=2k=1n1ekL2+12e0L2+12enL2.E25

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

2.5. Roundoff-induced reversibility error

The reversibility error method (REM) is a very simple procedure based on niterations of the map Mfollowed by niterations of the inverse map. The distance from the initial point normalized by the roundoff amplitude ϵdefines the REM error. Denoting with Mϵthe map with roundoff, we have

enREM=MϵnMϵnx0x0ϵE26

where ϵis the roundoff amplitude. For the eight-byte representation of reals, we choose ϵ=1017. If the map has a sufficiently high computational complexity, the displacement ξdefined by MϵxMx=ϵξis almost random, but a unique realization is available. (For a discussion on the roundoff error, see [22]). As a consequence, the enREMhas large fluctuations, whereas enRhas a smooth dependence on nsince 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 eLtaccording to

eLt=TrLTtLt1/2E27

where Ltis the fundamental matrix for the tangent flow, which satisfies the linear equation dL/dt=JHL, Hdenotes the Hessian of the Hamiltonian computed along the orbit Hij=2H/xixj, and the initial condition for the matrix Ltis L0=I. The relation with the standard fast indicators is the same as for the symplectic maps. Let ΞRtbe the stochastic displacement from x0after a noisy evolution up to time tand backward to t=0, divided by the noise amplitude εin the limit ε0. It has been proven [12] that ΞRtsatisfies a linear Langevin equation whose solution reads

ΞRt=0tL1sξsξstds.E28

The reversibility error in this case is defined by the mean square deviation of the random displacement eRt=ΞRtΞRt1/2. As shown in [12] and from Eq. (28), we immediately obtain

eRt=20teLs2ds1/2.E29

If the continuous time tis replaced by an integer nand 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 nαwith exponent α=3/2. If the map is linear, its asymptotic growth follows a power law with α=0for LE and α=1/2for 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 Mxis conjugated to its normal form NXby a symplectic coordinate transformation x=ΦX; as a consequence the conjugation equation and its iterates read

Mx=ΦNΦ1xMnx=ΦNnΦ1xE30

which imply that the orbits xn=Mnx0and Xn=NnX0are related by xn=ΦXn. The tangent maps are given by

DMnx0=DΦXnDNnX0DΦ1x0=DΦXnDNnX0DΦX01E31

where we used DΦXDΦ1x=I, a relation which is proved to hold by computing the Jacobian of ΦΦ1x=x. As a consequence, the expression for the Lyapunov error in both coordinate systems is

enLX02=TrDNnX0T(DNnX0enLx02=TrDMnx0T(DMnx0E32

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

enLx02=TrV1X0DNnX0TVXnDNnX0VX=DΦXTDΦXE33

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

VabbcV1=cbbaE34

where acb2=1.

3.2. Isochronous rotations: oscillations in LE and RE

If the given map is linear and two-dimensional and Mx=Lxwith TrL<2, then the map is conjugated to a rotation Rω:

L=TRωT1Rω=cosωsinωsinωcosωE35

Letting xn=Lnx0and Xn=RnX0, the orbits in the coordinate xand the normal coordinate X=T1xand the Lyapunov errors are given by

enLX02=TrRR=2enLx02=TrV1RVR=2cos2+a2+c2+2b2sin2==a+c22+2a+c22cos2E36

where V=TTT, and we have used the representation given by Eq. (34) where DΦ=T. The error is constant in normal coordinate Xand oscillates between 2 and a+c22=a2+c2+2b2in the coordinate x. 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

enRX02=4nenRx02=2na+c22+2a+c22fnfn=k=1ncos2+cos2nkω=cos2+cos2n1ωcos21cos2ωE37

We shall first consider the dependence of the errors on the iteration order nfrom n=1up to a maximum value N. Then, we shall consider the dependence on the initial condition x0when it is varied on a one-dimensional grid crossing the origin for the value Nof the iteration number. We choose the linear map Lwhich depends on a single parameter λ, and its relation with the rotation frequency is

L=1λ1λ1sinω2=λ20λ4E38

The rotation Rωxis the linear part for the Hénon map, whereas Lxis 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 α=0for LE and α=1/2for RE. Oscillations are present when the coordinates are not normal.

For a generic map such as Ldefined by Eq. (38), the growth of REM follows a power law exponent α=1/2as 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 α=1for almost all the values of ω.

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 x 0 = 0.1 , p 0 = 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 x 0 = 0.1 , p 0 = 0 .

Letting X=XPTand ϕJbe the action angle coordinates defined by X=2J1/2cosϕand P=2J1/2sinϕ, the rotation in the Xplane becomes a translation on the cylinder:

ϕn=ϕn1+ωmod2πJn=Jn1E39

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 01. Unlike RE, we observe that REM depends linearly on the distance of the initial condition x0from the origin. In Figure 2, we plot eNREMas a function of the initial condition when it varies on a one-dimensional grid issued from the origin for 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.

Figure 2.

Left frame: reversibility error due to roundoff e N REM for a rotation R ω with ω = 2 π 2 − 1 and N = 1000 when the initial condition is varied. We choose x 0 ∈ 0,0.5 , p 0 = 0 . The dependence on x 0 is evident and a linear fit f x 0 = 5000 x 0 is shown, purple line. Right frame: computation of the error for the linear map with λ = 4 sin 2 ω / 2 where ω has the same value. The linear fit is the same.

3.3. Anisochronous rotations

An integrable map Min normal coordinates and the tangent map DMnread

Mx=RΩJxDMnx=RnΩ+nΩ'R'nΩxxTE40

where J=x2/2is the action. The square of the Lyapunov error4 reads

enL2=TrDMnTDMn=2+n22JΩ'2E41

and the square of the reversibility error is given by

enR2=k=1nekL2+enkL2=4n+2JΩ'2k=1nnk2+k=1nk2=4n+2JΩ'22n33+n6E42

For a fixed value of the invariant J, the slope of enR2, whose asymptotic value is 2α, is defined as dlogenR2/dlogn, and its double average is given by MEGNO YnYen. The range of variation is [1, 3]. One can prove that for a given initial condition, the intermediate value Yn=2is reached for

n=14.52JΩ'=14.5x02+p02Ω'E43

In Figure 3, we show the variation with n11000of Yncomputed for RE given by Eq. (42), corresponding to the map presented in Eq. (40), where ΩJis a linear function of J, and find a perfect agreement with the analytical estimate of the value of n for which the value Yn=2is reached. In Figure 4, we show the variation of eNRand the corresponding MEGNO filter YNwith the initial condition chosen on a one-dimensional grid crossing the origin for N=1000. The integrable map is given by Eq. (40), where Ω'is constant. The error reaches a minimum value at the origin, and a similar behavior for YNis observed. Also eNREMdecreases 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 emREMfor m=1,,nwhose computational cost is of the order of n2. This can be avoided by storing the sequence xε,mand computing êmREM=Mεmxε,nxnmwhich turns out to be comparable with emREM.

Figure 3.

Left frame: plot of MEGNO Y n on the error e N R for 1 ≤ n ≤ N with N = 1000 for the integrable map with Ω ' = 0.1 and initial condition x 0 = 0.5 , p 0 = 0 (blue line). The green line refers to a modified definition YM n , where n log e n 2 − log e n − 1 2 is replaced with log e n 2 − log e n − 1 2 / log n − log n − 1 , which, applied to the sequence e n = n α , gives 2 α for any n . The vertical line gives the theoretical estimate of the value of n for which Y n = 2 (see Eq. (43)) Right frame: the same for x 0 = 1 , p 0 = 0 .

Figure 4.

Left frame: plot or the error e N R with N = 1000 for the integrable map as a function of the initial condition x 0 , p 0 = 0 with Ω ' J = 0.1 (red line) and Ω ' J = 1 (blue line). Center frame: same plot with e N REM for ω = 2 π 2 − 1 and Ω ' J = 0.1 , gray line, compared with e n R , red line. Right frame: plot of Y N for the integrable map with Ω ' = 0.1 (red line) and Ω ' 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/2in the exponent of the power law nαoccurs close to the origin, if the linear part is a rotation R, as for the Hénon map. If the linear part is Las 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 T2and reads

pn+1=pnλ2πsin2πxnmod1xn+1=xn+pn+1mod1E44

where x,pbelong to the interval 1/21/2whose ends are identified. For λ1and pλ, it is just a weakly perturbed rotator, and x,pare action angle coordinates. The origin is an elliptic fixed and very close to it; the map is approximated by a linear map

xn+1=1λxn+pnpn+1=λxn+pnE45

This map is conjugated to a rotation Rωfor 0<λ<2where sinω/2=λ/2. The point x=±1/2,p=0is hyperbolic, and for λ1the corresponding orbit is approximated by the separatrix of the Hamiltonian:

H=p22λ2π2cos2πxE46

which is the interpolating Hamiltonian of the map when λ0. We observe that the frequency for small oscillations is ω=λ(see Eq. (38)) when λ0. Since the time scale of the Hamiltonian His T=2π/λ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 His given by

p=±λπcosπx.E47

As a consequence, for λsmall the width of the separatrix is 2λ/π. When λ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

xn+1pn+1=Rωxnpn+xn2Rω=cosωsinωsinωcosωE48

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

H=ωp2+x22x33E49

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 x=2tan(ω/2)p=2tan2ω/2which corresponds to the critical point x=ωp=0of the Hamiltonian. The stability boundary is approximated by Hxp=ω3/6whose orbit explicitly reads p=±ωxω+2x/3. 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 enfor a fixed initial condition by varying nup to a maximum value N, by varying the initial condition on a one-dimensional grid for n=Nand 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 MEGNO confirm this observation. In Figure 5, we plot the errors enfor λ=0.1and Yenby varying n. The fast oscillations of LE and the large fluctuations of REM are clearly visible.

Figure 5.

Left frame: plot of the errors for the standard map with λ = 0.1 and initial condition x 0 = 0 , p 0 = 0.075 . Lyapunov error e n L (blue line), reversibility errors e n R (red line), and e n REM (gray line), for 1 ≤ n ≤ 2500 . Right frame: plots for the MEGNO filter Y n for the same errors.

When the orbit is chaotic, the growth of all errors is exponential. However, LE and RE can grow until the overflow is reached, whereas REM can grow only up to 1/ϵwhere ϵis the machine accuracy. Typically in double precision, the overflow corresponds to 10300where ϵ11017. 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 n300. For higher values the saturation to 1017is evident, and REM ceases to grow exponentially.

Figure 6.

Left frame: plot of the errors for the standard map with λ = 0.8 and initial condition x 0 = 0 , p 0 = 0.26 corresponding to a chaotic orbit. Lyapunov error e n L (blue line), reversibility errors e n R (red line), and e n REM (gray line), for n ≤ 500 . Right frame: plots for the MEGNO filter Y n for the same errors.

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.

Figure 7.

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 x 0 = 0 , p 0 ∈ − 0.15,1,0.15 . Right frame: the same for MEGNO Y N .

In Figure 8, the same results are shown for a higher value of the parameter λ=0.8at 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.

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 x 0 = 0 , p 0 ∈ 0,0.5 . Right frame: magnification in the interval p 0 ∈ 0.25,0.5 .

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×Ngpoints, 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=500and Ng=200the color plots for the errors of the standard map with λ=0.8and 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 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 1010 for LE and 1015 for RE and REM. This choice is suggested by the asymptotic behavior nαof the error for regular orbits where α=1for LE and α=3/2for RE.

Figure 9.

Left frame: standard map with λ = 0.8 color plot of LE in a logarithmic scale for N = 500 and a grid with N g = 200 . Center frame: standard map with λ = 0.8 color plot of RE in a logarithmic scale for N = 500 and a grid with N g = 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 N g = 200 . Center frame: standard map with λ = 1.5 color plot of RE in a logarithmic scale for N = 500 and a grid with N g = 200 . Right frame: color plot of REM in a logarithmic scale.

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/23/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

Ωω+JΩ2Ω2=183cotω2+cot3ω2E50

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

In the normal coordinates XP, the behavior of errors is given by Eqs. (41) and (42). In the original coordinates xp, the error could be evaluated using Eq. (33). In normal coordinates, the errors grow as 2JΩ2nαwhere α=1for LE and α=3/2for 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/logJsJand consequently Ω'JJsJ1as Jsis 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.21which 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 Nand after filtering them with MEGNO, when the initial conditions are chosen on a one-dimensional grid. The resonance 1/5is met, as shown by the appearance of a large chain of islands, since ΩJis 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.

Figure 11.

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 N = 1000 along the line x = r cos α , p = r sin α with α = 14 o 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.

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 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 N g = 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.

7. Conclusions

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

We have first examined the behavior of LE and RE for linear maps and for integrable maps. If the fixed point is elliptic, then the asymptotic growth follows a power law nα, and the exponent does not depend on the chosen coordinates for LE and RE. Conversely, the presence of oscillations and their amplitude depends on the choice of coordinates. The growth of REM also follows a power law, but the choice of coordinates affects the exponent itself.

For a generic map which has regular and chaotic components, the error growth follows a power law and an exponential law, respectively. For the standard map and the Hénon map, the behavior of LE, RE, and REM has been compared first by varying the iteration order nup to a some value N, for a selected initial condition. Then the errors for n=Nhave been compared when the initial point moves on a line. The theoretical predictions concerning the power law growth in the regular regions and the exponential growth in the chaotic ones are confirmed. For two-dimensional maps, the error plots for initial conditions in a rectangular domain of phase space are very similar, and the correspondence with the phase space portraits is excellent. Moreover, the different plots allow a quantitative comparison of the orbital sensitivity to initial displacements, noise, and roundoff. For maps of dimension 4 or higher, the proposed error plots on selected phase planes allow to inspect the orbital stability. Hamiltonian flows must be approximated with a high accuracy by symplectic maps, with algorithms which provide simultaneously the corresponding tangent maps [20, 21], in order to compute the errors discussed so far. A special care is required in comparing RE with REM when the chosen phase plane is invariant. Indeed given an initial point in the invariant plane, the noise brings the orbit out of it, whereas the roundoff usually does not. In this case a random kick before reversing the orbit is sufficient to bring the orbit out of the invariant plane and to restore the correspondence between REM and RE. The satisfactory results obtained so far, not only in the simple models presented here but also in high dimensional models of celestial mechanics, prove that the method we propose has a wide range of applicability.

Notes

  • The forward error (FE) due to additive noise in the forward evolution of a map can be defined and written in terms of the tangent map. However, RE is very simply related to LE, whereas FE is not. In addition the error due to roundoff in the forward evolution requires in principle the evaluation of the exact orbit or, in practice, its evaluation with a much higher accuracy.
  • The application of MEGNO to RE is not necessary due to the absence of oscillations, whereas its application to REM is not recommended because the fluctuations are not filtered and the computational cost is quadratic rather than linear in the iteration order.
  • A symplectic matrix L is defined by LJL T = J where J is antisymmetric and J 2 = − I . As a consequence L − 1 = − JL T J , and L − 1 T = − JLJ so that Tr L − 1 T L − 1 = Tr J LJ 2 L T J = Tr L L T .
  • The standard definition for an initial displacement along the unit vector η 0 is e n L η 0 = ∥ DM n η 0 ∥ where ∥ DM n η 0 ∥ 2 = 1 + n Ω ' 2 ∥ x ∥ 2 η 0 ⋅ x 2 + 2 n Ω ' η 0 ⋅ x η 0 ⋅ J x and J = 0 1 − 1 0 . The sum ‖ DM n 1 0 ‖ 2 + ‖ DM n 0 1 ‖ 2 gives the Lyapunov error e n L 2 .

Download

chapter PDF

© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Giorgio Turchetti and Federico Panichi (July 25th 2019). Fast Indicators for Orbital Stability: A Survey on Lyapunov and Reversibility Errors [Online First], IntechOpen, DOI: 10.5772/intechopen.88085. Available from:

chapter statistics

89total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us