Birkhoff Averages and Rotational Invariant Circles for Area-Preserving Maps
BBirkhoff Averages and Rotational Invariant Circles forArea-Preserving Maps
E. Sander and J.D. Meiss ∗ Department of Mathematical Sciences Department of Applied MathematicsGeorge Mason University University of ColoradoFairfax, VA 22030, USA Boulder, CO [email protected] [email protected] 3, 2020
Abstract
Rotational invariant circles of area-preserving maps are an important and well-studied exam-ple of KAM tori. John Greene conjectured that the locally most robust rotational circles haverotation numbers that are noble, i.e., have continued fractions with a tail of ones, and that,of these circles, the most robust has golden mean rotation number. The accurate numericalconfirmation of these conjectures relies on the map having a time reversal symmetry, and thesemethods cannot be applied to more general maps. In this paper, we develop a method based ona weighted Birkhoff average for identifying chaotic orbits, island chains, and rotational invariantcircles that do not rely on these symmetries. We use Chirikov’s standard map as our test case,and also demonstrate that our methods apply to three other, well-studied cases.
The dynamics of an integrable Hamiltonian or volume-preserving system consists of periodic andquasi-periodic motion on invariant tori. When such a system is smoothly perturbed, Kolmogorov-Arnold-Moser (KAM) theory [dlL01] implies that some of these tori persist and some are replaced byisolated periodic orbits, islands, or chaotic regions. On each KAM torus, the dynamics is conjugateto a rigid rotation with some fixed frequency vector. Typically, as the perturbation grows theproportion of chaotic orbits increases and more of the tori are destroyed. Invariant tori can befound numerically by taking limits of periodic orbits [Gre79] and by iterative methods based on theconjugacy to rotation [HdlL06]. In these methods, one fixes a frequency vector and attempts to findinvariant sets on which the dynamics has this frequency. In this paper we explore an alternativetechnique, based on windowed Birkhoff averages [DSSY16], to distinguish between chaotic, resonant, ∗ JDM was supported in part by NSF grant DMS-181248. ES and JDM acknowledge support from NSF grantDMS-140140 while they were at residence at the Mathematical Sciences Research Institute in Berkeley, CA, duringthe Fall 2018 semester. Useful conversations with Xinzhi Rao are gratefully acknowledged. a r X i v : . [ n li n . C D ] D ec nd quasiperiodic dynamics. Since we do not fix the rotation vector in advance, this method permitsone to accurately compute the rotation vector for each initial condition that lies on a regular orbit.As such the method is analogous to Laskar’s frequency analysis [LFC92, BBGT96], which uses awindowed Fourier transform to compute rotation numbers.As an illustrative example, we will primarily study Chirikov’s standard map [Chi79], thoughin the last section we will consider generalizations of this map. Two-dimensional, area-preservingmaps are simplest, nontrivial case of Hamiltonian dynamics (for a review, see [Mei92]). Letting f : M → M , where M = T × R , the cylinder, the standard map can be written as ( x t +1 , y t +1 ) = f ( x t , y t ) with x t +1 = x t + Ω( y t +1 ) mod 1 ,y t +1 = y t + F ( x t ) . (1)For Chirikov’s case, the “frequency map” and “force” are given byΩ( y ) = y, F ( x ) = − k π sin(2 πx ) , respectively. When the parameter k = 0, the action y is constant, and every orbit lies on a rotationalinvariant circle with rotation number ω = Ω( y ). When ω is irrational the orbit is dense on thecircle, and the dynamics is conjugate to the quasiperiodic, rigid rotation θ → θ + ω mod 1 (2)for θ ∈ T , under the trivial conjugacy ( x, y ) = C ( θ ) = ( θ, ω ).As k increases, some of these rotational invariant circles persist, as predicted by KAM theory,but those with rational or “near” rational rotation numbers are destroyed. On each KAM circle,the dynamics is still conjugate to (2), for some irrational ω , under a smooth map C : T → T × R .As an exampel, Fig. 1 depicts the dynamics for the Chirikov map for k = 0 .
7. In the top row,we distinguish between non-chaotic and chaotic dynamics, and in the bottom row we distinguishbetween two types of non-chaotic behavior, namely island chains and rotational invariant circles.The methods for doing this will be discussed in § { ( x t , y t ) : t ∈ Z } has a rotation number ω if the limit ω = lim T →∞ T T − (cid:88) t =0 Ω( y t ) (3)exists. Of course, if k = 0, ω is simply the value of Ω on the conserved action. If an orbit isperiodic, say ( x n , y n ) = ( x , y ) + ( m,
0) for some integers m, n , then ω = mn is rational. Indeed,(1) implies that if we lift x to R , then x T − x = T (cid:88) t =1 Ω( y t )so for the periodic case ω = ( x n − x ) /n . Note that ω , as a rotation number, is measured withrespect to rotation in x . For an invariant circle within an island chain, the effect of the rotation ofthe orbit about the island center will average out, and ω will equal m/n , the value for the periodic2rbit it encloses. This can be seen in the lower left portion of Fig. 1, where each elliptic islandhas a single solid color due to having the same value of ω . In particular, the rotational invariantcircles are the only non-chaotic orbits with the property that ω is irrational. In §
3, we develop anumerical method to determine whether a floating point number is (with high probability) rationalor irrational. With this method, we are able to use the rotation number computed with the weightedBirkhoff average to distinguish between rotational and non-rotational invariant circles.The invariant circles that persist by KAM theory have Diophantine rotation numbers, i.e., thereis a τ ≥ c > | nω − m | > c | n | τ , ∀ n ∈ N , m ∈ Z . (4)Such rotation numbers are hard to approximate by rationals (see § M in which it is the last invariant circle; i.e., it existsfor 0 ≤ k ≤ k cr ( ω ) and k cr is a local maximum. It is known from careful numerical studies thatinvariant circles with “noble” rotation numbers (their continued fractions have an infinite tail ofones) are robust [Gre79, Mac93]. Since these continued fractions are asymptotically periodic, theserotation numbers are quadratic irrationals and satisfy (4) with τ = 1. John Greene discoveredthat the last rotational circle of (1) has rotation number given by the golden mean γ , and thatit is destroyed at k = k cr ( φ ) ≈ . k > k cr no rotational invariant circles areobserved.The average (3) need not exist for orbits that are neither periodic nor quasiperiodic. For exampleif an orbit is heteroclinic between two periodic orbits with different rotation numbers, the forwardand backward time averages of y will be different. Moreover, when k is large enough, y can beunbounded, and the limit (3) need not even converge. However, if an orbit ergodically covers abounded region, then Birkhoff’s ergodic theorem implies that the time average of Ω does exist.More generally a finite-time Birkhoff average on a orbit of a map f beginning at a point z ∈ M for any function h : M → R is given by B T ( h )( z ) = 1 T T − (cid:88) t =0 h ◦ f t ( z ) . (5)This average need not converge rapidly. Even if the orbit lies on a smooth invariant circle withirrational rotation number, the convergence rate of (5) is O ( T − ), due to edge effects at the twoends of the finite orbit segment. By contrast, for the chaotic case, the convergence rate of (5) isobserved to be O ( T − / ), in essence as implied by the central limit theorem [LM10].We can significantly improve the convergence of a Birkhoff sum on a quasiperiodic set by usingthe method of weighted Birkhoff averages developed in [DSSY16, DDS +
16, DSSY17], see § f , the function h , and the quasiperiodic set are C ∞ , and the rotation number is Diophan-tine, this method is superconvergent, meaning that the error decreases faster than any power of T [DSSY17, DY18]. However, the weighted Birkhoff method does not speed up the convergencerate on chaotic sets since these lack smoothness. Therefore weighted Birkhoff averages have two Or any integer shift of this value by a discrete symmetry of (1). Often one thinks of y as diffusing in this case, but it can also grow linearly in time due to “accelerator modes”[Chi79]). ynamics at k = 0.7 ChaosIslands Rotational circles ω Figure 1:
The dynamics of the standard map for k = 0 .
7. Using the weighted Birkhoff method we are able todistinguish chaotic orbits (upper right), islands (lower left), and rotational circles (lower right). The rotation numberof each nonchaotic orbit is color-coded (color bar at right). The computations were performed for an evenly spacedgrid of 1000 points in [0 , with T = 10 , using the distinguishing criteria in (23). distinct uses: (a) to distinguish chaotic from regular dynamics, and (b) to give a high precisioncomputation of the rotation number. This method has also been used to find a high precision com-putation of the Fourier series expansion of the conjugacy map for the invariant circle [DSSY17],but we do not make use of this in the current paper.A finite-time computation of the rotation number [SSC +
13] has been used to define coherentstructures by considering ridges in a finite time sum (3). This method also can distinguish betweentrapped and escaping orbits [SMS +
19] by monitoring the gradient of (3) with respect to initialcondition, and to determine the breakup of circles in nontwist maps [SMS + k cr = 0 . § § § § § §
4, after removing chaotic orbits andisland chains, we are left with the rotational circles. We are able to create the critical functiondiagram, and describe the number theoretic properties of the rotation numbers for rotational circles,showing that their behavior does not match that of randomly chosen irrational numbers. In §
5, weapply our methods to three different generalizations of the standard map, namely a symmetric two-harmonic generalized standard map, a standard non-twist map, and an asymmetric two-harmonicmap. We conclude in § In this section, we introduce the weighted Birkhoff method, and we compare it to two differentmethods for distinguishing chaos from regular dynamics, namely Lyapunov exponents and the 0–1test of Gottwald and Melbourne [GM09].
We now describe in more detail the method of weighted Birkhoff averages [DSSY16, DDS + g ( t ) ≡ (cid:26) e − [ t (1 − t )] − t ∈ (0 , t ≤ t ≥ , be an exponential bump function that converges to zero with infinite smoothness at 0 and 1, i.e., g ( k ) (0) = g ( k ) (1) = 0 for all k ∈ N . To estimate the Birkhoff average of a function h : M → R efficiently and accurately for a length T segment of an orbit, we modify (5) to compute WB T ( h )( z ) = T − (cid:88) t =0 w t,T h ◦ f t ( z ) , (6)where S = T − (cid:88) t =0 g (cid:0) tT (cid:1) , w t,T = 1 S g (cid:0) tT (cid:1) . (7)That is, the weights w are chosen to be normalized and evenly spaced values along the curve g ( t ).For a quasiperiodic orbit, the infinitely smooth convergence of g to the zero function at the edgesof the definition interval preserves the smoothness of the original orbit. Indeed it was shown in[DY18] that given a C ∞ map f , a quasiperiodic orbit { f t ( z ) } with Diophantine rotation number,5nd a C ∞ function h , it follows that (6) is super-convergent: there are constants c n , such that forall n ∈ N (cid:12)(cid:12)(cid:12)(cid:12) WB T ( h )( z ) − lim N →∞ B N ( h )( z ) (cid:12)(cid:12)(cid:12)(cid:12) < c n T − n . (8)Laskar [LFC92] used a similar method to compute frequencies with a sin ( πs ) function instead ofa bump function, but this function is fourth order smooth rather than infinitely smooth at the twoends, implying that the method converges as O ( T − ), see e.g., [DSSY17, Fig. 7]. By contrast, whenan orbit is chaotic (i.e., has positive Lyapunov exponents), then (6) typically converges much moreslowly; in general it converges no more rapidly than the unweighted average of a random signal,i.e., with an error O ( T − / ) [LM10, DSSY17].A graph of the error in W B T for h = cos(2 πx ) as a function the number of iterates T isshown on the left panel of Fig. 2. Here we have chosen 50 orbits of (1) for the parameter k = 0 . x = 0 .
45, and y evenly spaced between 0 and 0 .
5. For orbits that areindependently identified as chaotic (red)
W B T essentially does not decrease with T , however fororbits that lie on rotational (blue) or island (green) invariant circles, the error for all but three hasdecreased to machine precision, 10 − , when T reaches 10 . Further, right panel of Fig. 2 showsthe convergence rate as a function of y , this time for 1000 orbits. Note that there is no evidence ofsuperconvergence in Fig. 2: the convergence rate for (8) has n = 2 − .
5. Indeed, superconvergencewas only observed in [DSSY17] when extended precision computations were done. Nevertheless,there is a clear distinction between chaotic and regular orbits for T ≥ . -20 -16 -12 -8 -4 T e rr o r chaoticislandscircles
0 0.1 0.2 0.3 0.4 0.5 y c onv e r g e n ce r a t e chaoticislandscircles Figure 2:
Convergence of the weighted Birkhoff average (6) for orbits of the standard map at k = 0 . h = cos (2 πx ). On left, the error of the computation is shown as a function of the number of iterates T for50 initial conditions at x = 0 .
45 with y on a grid in [0 , . x and k values, where convergence rate was calculated using the errors before the valuesflattened out due to floating point errors, measured by where they have dropped below 10 − . In each case, thevalues are compared with W B T (Ω) at T = 10 . Using the distinguishing criteria in (23), the red curves are identifiedas chaotic, the green curves as islands, and the remaining blue curves are thus the rotational invariant circles.
6o distinguish chaotic from regular dynamics, we compute (6) for two segments of an orbit,using iterates { , . . . , T } and { T + 1 , . . . , T } . In the limit T → ∞ , these values should be thesame. Therefore we can measure convergence rate by comparing them. In order to distinguishchaotic sets, we compute the number of consistent digits beyond the decimal point in our twoapproximations of WB ( h ), which is given by dig T = − log (cid:12)(cid:12) WB T ( h )( z ) − WB T ( h )( f T ( z )) (cid:12)(cid:12) . (9)If dig T is relatively large, then the convergence is fast, meaning the orbit is regular. If dig T is small,then the convergence is slow, meaning the orbit is chaotic. Three examples are shown in Fig. 3(a)for a set of 1000 initial points on a vertical line segment at x = 0 .
321 for three different values of k .For the smallest parameter, k = 1 .
0, a substantial fraction of the orbits are regular, and these havea distribution of dig centered around 14, nearing the maximum possible for a double precisioncomputation. By contrast, when k = 2 . dig centered around 2. Note that when k = 1 . dig T ∈ [6 , dig T between regular and chaotic orbits, we computeda histogram (not shown) of dig for the Chirikov standard map for 500 ,
000 different startingpoints: a grid of 500 k -values between 0 . .
5, with 1000 distinct initial conditions for each.This histogram has two large peaks, one at around 2 and the second around 15 (corresponding tothe machine epsilon value). As is consistent with the case k = 1 . dig T = 5. In our calculations of chaos, we wish to err on the side of falsepositives of chaos, and thus we use a value of 5 . x = 0 .
0, 0 .
321 and 0 . k ranging from 0 . .
5. When x = 0 .
0, the figure is dominated by the regular region around the fixed point islands surrounding(0 ,
0) and (0 , k = 4 .
0. Other islands can also be seen; for examplefor x = 0 . x = 0 .
5, we can see the period two orbit (0 , ) (cid:55)→ ( , ), which is elliptic until k = 2 .
0, where it period doubles. By contrast the line x = 0 .
321 intersects fewer islands, and thereappear to be no regular orbits when k ≥ .
915 with initial conditions on this line.
In this section we recall two other standard tests for chaos: positive Lyapunov exponents and the0–1 test. The finite-time Lyapunov exponent is defined by λ T ( v ) = 1 T log (cid:0) | Df T ( x , y ) v | (cid:1) , (10)where Df is the Jacobian matrix, and v is a generic deviation vector with | v | = 1. Histograms of λ T for three values of k are shown in Fig. 3(b). As noted by [SLV05], when there are regular andchaotic orbits, these histograms are typically bimodal. For example, we observe that when k = 1 . λ = 0 with width of order 0 .
02. This peak is well separatedfrom the broader peak centered near λ = 0 .
3. The peaks are less well separated for smaller values7igure 3:
Histograms of (a) Weighted Birkhoff accuracy, dig T , and (b) finite-time Lyapunov exponent, λ T , fororbits of the standard map with k = 1 .
0, 1 . .
0. The initial conditions are (0 . , y ) with 1000 values of y ona uniform grid in [0 , dig T (9), for h ( x, y ) = cos(2 πx ) and T = 10 . (b) Histograms of λ T (10),for T = 2(10) , and v = (0 , T . of k ; for example when k = 1 . λ T < .
01, and there a broader peakof presumably chaotic trajectories with λ ∈ [0 . , . λ ≈ .
05. As k grows, the mean value of λ increases and the lower peak of regularorbits disappears.To visualize the dependence of the exponents on k , we chose the same three lines of initialconditions shown in Fig. 4 for the weighted Birkhoff average. The resulting exponent, as a functionof y and the parameter k of (1) is shown in Fig. 5. In this figure, orbits with λ < . − arecolored black: these correspond to the regular orbits. As k grows, the distribution in the chaoticregion is peaked around a growing value that reaches a maximum of λ = 0 .
605 when k = 2 .
5. Notethat each of the panes of this figure is essentially the negative of the corresponding pane in Fig. 4.The fraction of chaotic orbits can be estimated by removing orbits with λ T in the range of thelower peak. Since the value λ = 0 .
05 is a minimum in a histogram (not shown) of λ values over k in [0 . , . k is shown in Fig. 6. This fraction is strongly dependent on theline of initial conditions. For the lines of symmetry (e.g. x = 0 or 0 .
5) of the standard map, thefraction of orbits trapped in regular islands is larger.Another test for chaos is the 0–1 test of Gottwald and Melbourne [GM09]. This test involvescomputing a time series (here we use { sin(2 πx t ) : t ∈ [0 , T ] } ) from which a supplemental time series(called ( p t , q t ) in [GM09]) is constructed and tested for diffusive behavior. This ultimately gives aparameter, K median , that is ideally either 0, when the orbit is quasiperiodic, or 1 when it is chaotic,8igure 4: Using the weighted Birkhoff method, this figure shows the number of digits (9) for orbits of the standardmap (1) with (a) x = 0 and (b) x = 0 .
321 and (c) x = 0 . y ∈ [0 ,
1] and k ∈ [0 . , . dig T is computedby computing the average (6) for the function h ( x, y ) = cos(2 πx ) for T = 2(10) steps. Initial conditions with dig T < . dig T for the regular orbits is indicated in the color bar. Figure 5:
Lyapunov exponent λ for orbits of the standard map (1) with (a) x = 0 and (b) x = 0 .
321 and (c) x = 0 . y ∈ [0 ,
1] and k ∈ [0 . , . v = (0 , T and T = 2(10) . Blackcorresponds to λ < . λ for the chaotic orbits is indicated in the color bar. k C h a o ti c F r ac ti on x = 0.321 x = 0.5 x = 0.0 Figure 6:
Using the Lyapunov method, this figure shows the fraction of orbits that have λ > .
05 for three valuesof x , a uniform grid of 1000 values of y , and T = 20000. When k < ,
0) that decreases the number of chaotic orbits found for x = 0 (red curve), and when k < . , .
5) has a smaller, but similar effect for x = 0 . x = 0 .
321 (green curve), the regular regions around both of these elliptic orbits are not sampled when k ≥ .
4. Thesevariations can be observed in Fig. 4 and Fig. 5. When k ≥ . y values are deemed tonot be chaotic. and we use the cutoff K median > . T = 1000 and 100 random samples gives an algorithm that isabout 200 times slower than computing Lyapunov exponents. The resulting dichotomy betweenregular and chaotic orbits for this test is shown in Fig. 7 for initial conditions at x = 0 .
0. Thisfigure agrees well with those in Fig. 4(a) and Fig. 5(a), though it appears to identify slightly fewerorbits as chaotic than the Lyapunov test: some orbits that were designated chaotic by Lyapunovexponent do not have K median > . In this section we systematically compare the detection of chaos for the Chirikov standard mapusing the three different methods: Lyapunov exponents, 0–1 test, and the weighted Birkhoff method.Figure 8 shows the fraction of orbits identified as chaotic for initial conditions on the line x = 0 . y ∈ [0 , k = 1, 1 . .
3; these values correspond to majorbifurcations in which regular islands and circles are destroyed. Nevertheless, the mean absolute10 .5 1 1.5 2 2.5 k y Figure 7:
Chaotic region of the standard map with x = 0 for y ∈ [0 ,
1] and k ∈ [0 . , .
5] using a time series { sin(2 πx t : 0 < t ≤ } and the 0–1 method of [GM09]. An orbit is deemed to be chaotic (colored white) if the0–1 parameter K median > . deviation between the weighted Birkhoff and 0–1 test results is 2 . TSS = T PT P + F N − F PF P + T N .
Here
T P (“true positive”) is the fraction of initial conditions that are classified correctly as chaoticby the test method over the reference standard,
F P is the fraction classified incorrectly as chaotic,
T N is the fraction that are correctly as non-chaotic, and
F N is the fraction classified incorrectly asnon-chaotic. The
TSS ranges from − TSS is that it does not depend upon the numberof trials, just on the relative accuracy. However, if we are comparing two predictions, the skillstatistic does depend upon which prediction is designated as the “ground truth”: changing thisdesignation is equivalent to exchanging
F P ↔ F N .An alternative measure is to simply count the percentage of correctly classified initial conditions, R = T P + T NT P + F P + T N + F N , the “ratio” of [Woo78]. 11 k C h a o ti c fr ac ti on Figure 8: (b) The fraction of chaotic orbits as a function of k for initial conditions on the line x = 0 .
0. Thered curve shows the fraction of chaotic orbits computed using Lyapunov exponents with λ > .
05, the blue curveshows the fraction with the 0–1 parameter K median > .
5, and the yellow curve shows the fraction of chaotic orbitscomputed using the weighted Birkhoff average with dig < .
5. For each k value, 1000 initial conditions were used forthe Lyapunov and weighted Birkhoff methods. For the Lyapunov method T = 2(10) , and for the weighted Birkhoffmethod T = 10 (which involves 2(10) calculations, so is equivalent to the Lyapunov method). For the 0–1 test, T = 1000 and 500 initial conditions were used for each k -value. For the comparison, in Fig. 8 we computed the Lyapunov exponent and weighted Birkhoffaverage using a total orbit length of 2(10) iterates, for k ∈ [0 . , .
5] on an evenly spaced gridof 500 values. At each parameter value, we chose initial conditions on the line x = 0 .
321 with y ∈ [0 ,
1] on an evenly spaced grid of 1000 points: thus thus there are 250 ,
000 trials. Since the 0–1method is slower, we chose a smaller number of iterates, T = 1000, and only 500 initial conditionsfor each parameter value.The agreement between every pair of the three methods gave R ≈ dig T ,and the best agreement occurs when dig T = 4. However, R is not very sensitive to this choice, andin fact varying dig T between 3 . TSS = 0 .
96 for both weighted Birkhoff and0–1 comparisons. Comparing the 0–1 test to the weighted Birkhoff average as the ground truthgives
TSS = 0 .
96 as well. Indeed, the number of false positives and false negatives are roughlyequal in the each of three comparisons: if they were exactly equal then
T SS would not changeupon the choice of ground truth.To validate the methods, we also compared each method to itself with double the number ofiterates. Comparing T = 10 to 2(10) for the λ T gave a different answer 0 .
8% of the time. Bycontrast for the WB T , the results disagreed 0 .
4% of the time. Thus the 2% difference between12ethods cannot be explained as a result of the number of iterates but is a true difference inidentification of chaotic orbits.The computation time for the Lyapunov exponent and weighted Birkhoff methods are roughlythe same. For example with 10 initial conditions and T = 10 iterates, each method took around6 minutes using Matlab on a Mac laptop. The 0–1 method was significantly slower. The fact thatthe computation time for the weighted Birkhoff average and Lyapunov exponent is roughly thesame is related to the fact that the Jacobian of the map (1) is quite simple. The implication isthat the time needed to do the averaging required for the weighted Birkhoff method is roughlythe same as the time needed to compute derivatives for the Lyapunov exponent. However, if thederivative were computationally more expensive, then the weighted Birkhoff method would have aspeed advantage.The weighted Birkhoff method has another advantage, as we will illustrate in the next section:it gives an accurate calculation of the rotation number ω that we can use to distinguish betweenrotational invariant circles and island chains. The regular orbits of the Chirikov standard map are of two distinct topological types: rotationalinvariant circles and orbits within the island chains. We are primarily interested in studying therotational invariant circles, and thus must look for a way to distinguish and remove orbits withinisland chains.For a twist map Birkhoff’s theorem implies that the rotational invariant circles are graphs, x (cid:55)→ ( x, c ( x )). Generically the dynamics on each such circle is conjugate to an incommensuraterotation, implying that ω in (3) is irrational.By contrast, around each elliptic period- n orbit there is generically a family of trapped orbitsforming a chain of n islands. The regular orbits in these island chains are further partitioned intoorbits that are quasiperiodic and those that are periodic relative to the n th power of the map. Thelatter, if elliptic, can again be the center of chains of islands. This gives rise to the familiar island-around-island structure. Each regular, aperiodic orbit within a period- n island chain is genericallydense on a family of topological circles: these are oscillational invariant circles. Nevertheless, therotation number ω , (3) will average out the internal dynamics, resulting in a rational value that isthe rotation number of the central period- n orbit. Of course if one were to measure the rotationnumber of an oscillational circle relative to the periodic orbit that it encloses, one would genericallyfind it to be irrational as well.In § ω for regular orbits. In § ω values are “rational,” and which are “irrational.”In § We are interested in establishing a numerical method to determine whether a numerically computednumber determined using floating point arithmetic is representative of a rational or an irrational.At the outset, this is not a well-posed question, since floating point representations of numbers13re rational. The question becomes whether a numerical value is—with high probability—theapproximation of a rational or an irrational number. In this section, we concentrate on a closelyrelated question, and in the next section we show how the answer can be applied to establishrationality. Our question is: given a number x , and an interval I δ ( x ) ≡ ( x − δ, x + δ ) , (11)with some tolerance δ , what is the rational number p/q with the smallest denominator in I δ ( x )?If, for a small δ , there is a rational p/q ∈ I δ ( x ) with a sufficiently small denominator q , we wouldexpect that x is—to a good approximation—given by this rational. Whereas if all such rationalshave large denominators, we would expect that x is an approximation of an irrational number.Actually, we will argue that if q is too large , x is more likely an approximation of a rational numberthat just missed being in the interval. We will return to the question of what constitutes small,large, and too large for values of q , but first we discuss the question of how to actually find thevalue p/q in a prescribed interval.We denote the smallest denominator for a rational in an interval I by q min ( I ) ≡ min { q ∈ N : pq ∈ I, p ∈ Z } . (12)The question of finding q min has been considered previously in [BBdA98, For07, CP16], and aclosely related question is considered in [CB09].Given an interval I in R , one would imagine that there are standard algorithms for finding therational number p/q in I with q = q min ( I ). Indeed packages such as Mathematica and Matlab bothhave commands that appear to do this. However these algorithms use truncations of the continuedfraction expansion [HW79], and neither of them work correctly in the sense of finding the smallestdenominator [For07]. Recall that the continued fraction expansion for x ∈ R + is x = a + 1 a + a + ... ≡ [ a ; a , a , . . . ] , a i ∈ N , a ∈ N ∪ { } . (13)Truncation of this path after a finite number of terms gives a rational “convergent” of x : p k q k = [ a ; a , a , . . . , a k ] . (14)Convergents are best approximants in the sense that if (cid:12)(cid:12)(cid:12)(cid:12) pq − x (cid:12)(cid:12)(cid:12)(cid:12) < q , (15)then p/q is a convergent to x [HW79, Theorem 184]. Moreover, at least one of any two successiveconvergents satisfies (15).However, the convergents are not necessarily the rationals with the smallest denominators in agiven interval. As a simple example, the rational with the smallest denominator within δ = 10 − of π is , i.e., q min ( I − ( π )) = 64. However, this rational is not a convergent of the continuedfraction π = [3; 7 , , , , , , ... ]; indeed, the first convergent in the interval is p q = .14 correct algorithm (e.g., that proposed by Forisek [For07]), is easiest to explain based on the Stern-Brocot or Farey tree. Every number in R + has a unique representation as a path on thisbinary tree: x = s s . . . , s i ∈ { L, R } . (16)The tree, whose first levels are sketched in Fig. 9, is constructed beginning with the root values and . Subsequent levels are obtained by taking the mediants of each neighboring pair: p m q m = p l q l ⊕ p r q r ≡ p l + p r q l + q r . (17)Level zero of the tree is the mediant of the roots, ; it is defined to have the null path. If x < ,then its first symbol is L , and if x > , then its first symbol is R . At level (cid:96) of the tree, 2 (cid:96) new rationals are added, the mediants of each consecutive pair. The left and right parents areneighboring rationals that have level less than (cid:96) . Every consecutive pair of rationals at level (cid:96) are neighbors in the sense that p r q l − p l q r = 1 . (18)A consequence is that p m and q m are coprime.For (cid:96) = 1, the new mediants = ⊕ and = ⊕ are added to give the level-two Fareysequence , , , , . Then the 2 mediants of each neighboring pair are added to give 2 level-threeintervals, see Fig. 9. Since the level-three rational > , to the right of its level-two parent, then = LR . Similarly is to the left of its level-two parent = R , so = RL .
01 10[0; 2] = 12 = [0 ; 1, 1][0; 1] = 11 = [1] [1; 1] = 21 = [2] [0; 3] = 13 = [0 ; 2, 1] [0; 1 ] = 23 = [0 ; 1, 2] [2; 1] = 31 = [3] LL RR RR LLLL [0;1, 1, 2] = 35 = [0;1 ][0;1, 1, 3] = 47 = [0;1, 1, 2, 1] [1;3] = 43 = [1;2, 1] [1; 2] = 32 = [1 ; 1, 1] [1;2, 1, 1] = 75 = [1;2, 2] Figure 9:
The continued fraction expansion and some entries on the Stern-Brocot or Farey tree. Each rational hastwo possible finite continued fractions but a unique Farey path. See the appendix for the relationship between thetwo.
The Farey path (16) for any x ∈ R + is the unique path of left and right transitions that leadto x starting at . Every rational has a finite path and every irrational number has an infinitepath [HW79]. Algorithm 1 in the appendix computes the Farey path, up to issues of floating pointaccuracy and a stopping criterion.The Farey expansion allows one to find the rational with smallest denominator in any interval:15 emma 1 (Smallest Rational) . The smallest denominator rational in an interval I ⊂ R + is thefirst rational on the Stern-Brocot tree that falls in I . The proof of this lemma, from [For07], is given in Appendix A. An alternative version of thisresult using continued fractions can be found in [BBdA98].An algorithm for finding q min ( I δ ( x )), based on Lem. 1 is given in the appendix in Algorithm 2.For example, for x = 0 .
12 = = L R = [0; 8 , , , , , , , , , , , . Given δ = 0 .
005 for example, the Farey interval ( , ) ⊃ I δ , and the mediant ∈ I δ since0 . − ≈ . < δ . Thus from the algorithm we obtainSmallDenom(0 . , . , ⇒ q min ( I . (0 . . As noted in [For07], the built-in routines of standard mathematical software do not always computethe smallest rational approximation correctly. For example, the built-in Matlab command “rat”gives rat(0 . , . , x itself, since the second convergent = 0 .
12 + 0 . I . The point is that the intermediate convergents of the Farey path can satisfy the approximationcriterion before the principal convergent of the continued fraction, and this can happen wheneverthe Farey path is not alternating . . . LR . . . or equivalently the continued fraction elements are notall 1’s.To determine the “typical” size of a denominator in an interval I , we show in Fig. 10 a histogramof the minimal denominator computed using Algorithm 2 in the appendix for randomly chosenfloating point numbers in (0 ,
1) with a uniform distribution. For this case, when δ = 10 − , themean minimal denominator appears to be close to 10 = δ − / . The distribution is not log-normal:the data is significantly more concentrated around the mean than a normal distribution with thesame standard deviation. Over the range δ = [10 − , − ], the mean log-denominator obeys therelation (cid:104) log q min (cid:105) = − log δ − . ± . , (19)and in this same range of δ values, the standard deviation is nearly constant, σ = 0 . ± . . (20)Further support for this statement is found in Fig. 11, which shows that for δ = 10 − tol , theprobability that q min is in the range 10 tol/ ± s does not depend on the choice of tol . Indeed, thecurves in this graph were obtained from only 10 random trials: if more values were randomlychosen, it would be impossible to distinguish between these distribution plots.The mean of our observations (19) is consistent with the expectation from (15). Indeed, for any δ , then there is a convergent with | x − p/q | < δ , with a denominator that must satisfy q ≥ (2 δ ) − / .Since the minimum denominator is no more than this, we expect that q min ∼ (2 δ ) − / , and thuslog q min ∼ − log δ − . , which is not far from the observation (19). 16igure 10: Probability density of log ( q min ) computed by appendix Algorithm 2 with δ = 10 − for 10 randomlychosen numbers in (0 ,
1) (black). This distribution has mean 5 . . σ = 0 . . randomly chosen numbers of constant type with A = 10 (green). s P r ob a b ilit y tol = 4 tol = 14Constant TypeNormal Figure 11:
A graph of the probability that q min ∈ tol/ [10 − s , s ], for δ = 10 − tol for 10 randomly chosen x ∈ [0 , tol = 4 and the red for tol = 14 (These curves are nearly indistinguishable).The yellow curve shows the probability for numbers of constant type with A = 10. The probability for a normaldistribution with standard deviation (20) is the dashed curve. A related result was obtained by [Ste13]: for intervals of the form J N = ( i − N , iN ], the mean17mallest denominator in grows asymptotically as (cid:104) q min ( J N ) (cid:105) ∼ CN / , with a coefficient 1 . < C < .
04. Since these intervals are of size N = 2 δ , this giveslog (cid:104) q min (cid:105) ∼ − log δ + K, K ∈ [ − . , . . Note that since the logarithm is convex, Jenson’s inequality implies that (19) is no larger than thisresult. We are not aware, however, of any results in the literature that imply the validity of (19)or (20).As a second numerical experiment, we consider numbers of constant type ; that is numbers thathave bounded continued fraction elements: sup k { a k } = A < ∞ . Such numbers can be thoughtof as “highly irrational” in the sense that they are Diophantine (4), with τ = 1, and c > A +2 .Conversely, if x is Diophantine with constant c then A < c [Sha92]. This class of numbers isespecially important in the context of area-preserving maps: it was conjectured that invariantcircles with constant type rotation numbers are locally robust and that every circle that is isolatedfrom at least one side has constant type [MS92].For the numerical experiment shown in Fig. 10, we chose rational numbers with continuedfractions of length 40, with a i ≤ , i = 1 , . . .
40 chosen as iid random integers. Note that thismeans that every trial x is rational; however, the denominator of these rationals is at least as largeas the case a i = 1, which gives the Fibonacci F ≈ . . The resulting smallest denominatordistribution is the green histogram in Fig. 10. The cumulative distribution of these numbers is alsoshown in Fig. 11, which shows that the probability that P rob ( | log ( q min ) − tol/ | > . x can result in bothextremely small and extremely large values of q min ( I δ ( x )). To demonstrate this, Fig. 12 shows aplot of q min ( I − ( x )) for evenly spaced x values between 0 .
095 and 0 . x -axis are centered at each rational with a denominator q ≤
80; the size of each dot is inverselyproportional to q . Note that in the vicinity of each dot, there is a small region where q min drops tothe corresponding small value of q , but additionally, there is a larger interval in which q min becomesmuch larger than average, with a larger jump near smaller denominators. Dynamically these orbitscorrespond to orbits that are limiting on the separatrices of islands, and hence are chaotic.The main takeaway message from the “typical size” experiments in this section is that thatnumbers outside the main peak of the distribution in Fig. 10 correspond to those “close” to rationals.In the next section, we will discard such rotation numbers to filter for candidates for rotationalinvariant circles. In this section, we use the weighted Birkhoff method to obtain an accurate computation of therotation number ω defined for the Chirikov standard map in (3). Namely, ω ( z ) = WB (Ω)( z ) . (21)Using this, we can distinguish rotational invariant circles from orbits in island chains by determiningwhether the computed value of ω is an approximation of a rational or irrational number as follows.18 x q m i n (I δ ( x )) Figure 12:
A plot of the smallest denominator, q min ( I δ ( x )), in an interval (11) for δ = 10 − and 10 values of x ∈ [0 . , . x -axis indicate the size of the denominator of each rational number in theinterval with a denominator up to 80; larger dots correspond to smaller denominators. There is a spike in denominatorsize immediately outside the interval around small denominator rationals. The mean log-denominator, (19), is shownby the dashed (red) line. Fix tol and let δ = 10 − tol . For a rotation number ω , we find q min ( I δ ( ω )) in (12), the smallestdenominator of a rational number within distance δ of ω . In most of our numerics, we have chosen tol = 8. To distinguish between rationals and irrationals, for each rotation number ω define theabsolute deviation dev ω = | log ( q min ( I δ ( ω ))) − tol/ | . (22)For a fixed cutoff value s , we remove the orbits within island chains as follows. Let z be an initialcondition of a regular orbit with associated rotation number ω . If dev ω > s , then we discard z asa member of an island chain. Note that this is equivalent to saying that q min is outside the range10 tol/ ± s .It remains to choose a cutoff value s . In our numerics, when we wish to be conservative aboutidentifying rotational circles, we have used the cutoff value s = 0 . dig T < . WB T (cos(2 πx )) , Irrationality criterion: dev ω < . tol = 8 . (23)In each case, for each initial condition ( x , y ), we compute an orbit and determine whether the19rbit is chaotic using the above chaos criterion. We also compute the rotation number ω using WB T (Ω) and determine whether the orbit is a rotational invariant circle using the irrationalitycriterion. Using the strategy of § ω for initial conditions ( x , y ) that areidentified to lie on rotational invariant circles using distinguishing criteria (23). ky ω x = 0.321 x = 0.0 x = 0.5 k k Figure 13:
The rotation number computed for rotational invariant circles orbits of the standard map (1) with (a) x = 0 and (b) x = 0 .
321 and (c) x = 0 . y ∈ [0 ,
1] and k ∈ [0 . , . T = 2(10) , using the distinguishing criteria in (23). The panels in Fig. 13 strongly resemble the critical function computed using Greene’s methodfor the standard map [MS91, MS92]. In this method, one typically chooses a set of noble irrationalnumbers, and finds the threshold of instability for a long periodic orbit that is close to each of thesenobles. The periodic orbits used in these computations are those that are symmetric under thereversor for (1); for example, every elliptic, symmetric rotational orbit is observed to have a pointon the line x = 0. An advantage of our current method is that symmetry is not required.It is believed that there are no rotational invariant circles for the standard map above k cr =0 . k > [MP85]. In Fig. 14, we show how the fraction of initial conditionsthat are identified as rotational circles in Fig. 13 varies with k . By k = 0 . .
9% of the circlesare destroyed and the fraction drops to zero at k = 0 . k = 0 . y .As a more precise test of the efficacy of the weighted Birkhoff average to determine k cr , we usedcontinuation to find an orbit on the line (0 . , y ) with the fixed rotation number γ − = ( √ − T = 2(10) . A computation of dig T , (9) can then be used to determine if the orbit is20 .1 0.3 0.5 0.7 0.9 k P r opo r ti on o f r o t a ti on a l c i r c l e s x = 0 x = 0.321 x = 0.5 Figure 14:
Fraction of orbits of (1) from Fig. 13 that are on invariant circles for 1000 initial conditions with y ∈ [0 ,
1] on three vertical lines as shown. The largest fraction occurs when x = 0 .
5, as this line tends to avoid manyof the larger islands. not chaotic. For the computation shown in Fig. 15, dig T = 12 at k = 0 . k = 0 . k cr . As an example, when k = 9697 / ≈ . ω = γ − has y = 0 . k > k cr , iteration shows that it remainslocalized to what appears to be a circle for hundreds of millions of iterations.We now focus on the number theoretic properties of rotation numbers for robust circles. Itis thought that the rotation numbers of the more robust invariant circles should have continuedfraction elements with more elements a i = 1 [Gre79, MS92]. To test this, we plot the distributionof continued fraction elements, a n , for the rotation number of invariant circles in Fig. 16. Theexpected distribution for randomly chosen irrationals is the Gauss-Kuzmin distribution [Sha92], P ( a i = k ) = log (1 + 1 / ( k ( k + 2))). When k is relatively small, the observed distribution followsthe Gauss-Kuzmin distribution closely, at least for a n ≤
10; but for k = 0 .
95, when most circleshave been destroyed, the probability of a n = 1 or 2 is larger than would be predicted for randomirrational numbers, and the probability that a n ≥ The method we have developed to find rotational invariant circles works equally well for otherarea-preserving maps. As a first example, we consider two-harmonic generalized standard map (1)with the force F ( x ) = − k π (sin( ψ ) sin(2 πx ) + cos( ψ ) sin(4 πx )) , (24)21 kdig T k cr Figure 15:
Computation of dig T for the golden mean circle for T = 2(10) for 500 values of k ∈ [0 . , . . , y ) that has ω = γ − . The drop of dig T from 12 to 6indicates that the circle is destroyed for a parameter value in (0 . , . that was first studied in [GJSS87] (see [Sim18] for later references). A phase portrait of this map,analogous to that shown for the standard map in Fig. 1, is shown in Fig. 17 for the value ψ = 0 . k is shown in Fig. 18. This figure is similar to [FM14, Figure 12(b)], where the criticalparameters were computed for a set of 256 noble rotation numbers. In that case the last invariantcircle, with ω ≈ . k ≈ . .
9% of the invariant circles are destroyed when k > . ω = 0 . , , , , , , , , , , . . . ] , with q min ( I − ) = 13153.A similar, well-studied map is the standard nontwist map (see [FM14, SMS +
18] for references).This map is of the form (1) with the standard force, but with the frequency mapΩ( y ) = y − δ. (25)The phase space of the dynamics for k = 1 . δ = 0 .
3. At these parametervalues, there are large chaotic regions around islands with rotation number 0 (colored green) and aband of rotational circles near the minimum of Ω (colored blue). The most robust circles tend to bethe shearless circles ; they cross the line y = 0 where Ω (cid:48) ( y ) = 0. The fraction of chaotic orbits androtational circles as k varies is shown in Fig. 20. For the 1000 initial y values in our experiment,the last detected rotational circle is at k = 2 . x , y ) = (0 . , − . ω = − . −
1; 1 , , , , , , , , . . . ] , -4 -3 -2 -1 P r obab ilit y k=0.300000k=0.950000Gauss-Kuzmin a n Figure 16:
Probability distribution for the occurrence of continued fraction elements of the rotation number forrotational invariant circles of the standard map for two values of k . These were computed using T = 10 iterates,and a grid of 5(10) initial conditions at x = 0 . k = 0 . k = 0 .
95 (red), we found 682. The black curve shows the Gauss-Kuzmindistribution, which is the distribution of elements for a random irrational chosen with uniform probability in [0 , with q min ( I − ) = 7260. We have independently verified that there are no rotational circles for k ≥ .
79 by direct iteration.Finally, we consider an asymmetric two-harmonic map (1) with the force F ( x ) = − k π (sin( ψ ) sin(2 πx ) + cos( ψ ) cos(4 πx )) , (26)studied in [FM14]. This map does not have the usual x (cid:55)→ − x reversor of the standard map (1),and therefore its periodic orbits are not aligned by a symmetry. Phase portraits for k = 0 . ψ = 0 . k in Fig. 22. The weighted Birkhoff average (6) and the distinguishing criteria (23) have been shown to efficientlycategorize orbits as chaotic, trapped in islands, or quasiperiodic on rotational circles. Using only T = 10 iterations, the rotation number of regular orbits is typically known to machine precision, asshown in Fig. 2. By contrast the weighted Birkhoff average of chaotic orbits converges much moreslowly, and this allowed us to identify chaotic trajectories. Orbits trapped in islands have rationalrotation numbers, and we are able to identify these using the distribution, shown in Fig. 10, of theminimal denominator in an interval of size δ defined by q min ( I δ ) in (12).The weighted Birkhoff method has the advantage that it does not rely on the reversing symmetryused to find periodic orbits in Greene’s residue method. Using a total orbit length of 2(10) , we23 Chaos ω Rotational circlesIslands
Dynamics at k = 0.5
Figure 17:
The dynamics of the two-harmonic map, with force (24) for k = 0 . ψ = 0 . initial conditions in [0 , × [ − . , .
75] with T = 10 , using the distinguishing criteria in (23). estimated the break-up parameter for the golden mean invariant circle to 0 .
3% accuracy, as seenin Fig. 15. While this accuracy does not compete with that of Greene’s method, our method canbe applied, as we have seen in § d -dimensional invariant tori. Here (with one exception [FM13]),Greene’s method no longer applies. There are several difficulties in any attempt to extend Greene’smethod; one is that there is no completely satisfactory continued fraction algorithm for multi-dimensional frequency vectors. To generalize our method will require computing the minimaldenominator q min for resonance relations, e.g., finding a minimal ( p, q ) ∈ Z d +1 such that | qω − p | issmall. One possible approach is to use generalized Farey path methods [KO86] that may provide a24 k F r ac ti on chaosrotational circles k ω y Figure 18: (a) Fraction of orbits of the map (1), with two-harmonic force (24) for ψ = 0 . x = 0 .
35. (b) The rotation number of the rotational circlesas a function of initial y and the parameter k . As in Fig. 17, T = 10 with distinguishing criteria in (23). version of Lem. 1 for this case. Appendices
A Farey Paths and the Smallest Denominator
The Farey path (16) for any number x can be computed by the simple method given in Algorithm 1.In a practical calculation, a stopping criterion based on precision must be included. This gives, forexample = RRLRRRR = [2; 1 , , , = LRRLL = [0; 1 , , , ,e = RRLR LRL RLR LRL . . . = [2; 1 , , , , , , , , , , , , , , . . . ] ,π = R L R LR LRLR LR L . . . = [3; 7 , , , , , , , , , , , . . . ] . Note that each element of the continued fraction records the number of repeated Farey symbols.The value of a is nonzero if the Farey path begins with R , otherwise a = 0, and a counts thenumber of leading L ’s in the path. For the rational case there is an additional last element, whichis fixed to be 1.The Stern-Brocot tree gives a method for finding the rational with the smallest denominator q min ( I ) (12) in an interval I . Here we prove Lem. 1 to show that q min is the denominator of thefirst rational on the tree that falls in I : Proof of Lem. 1.
Suppose that for all levels up to (cid:96) on the Stern-Brocot tree no Farey rational is25
Islands
0 0.2 0.4 0.6 0.8-0.6-0.4-0.200.20.40.6 0 0.2 0.4 0.6 0.8
ChaosDynamics at k = 1.5 -0.6-0.4-0.200.20.40.6 Rotational circles ω Figure 19:
The dynamics of the standard nontwist map, with frequency (25) for k = 1 . δ = 0 .
3. The weightedBirkhoff method distinguishes chaotic orbits (upper right), islands (lower left), and rotational circles (lower right).The rotation number of each nonchaotic orbit is color-coded (color bar at right). The computations were performedfor a grid of 1000 initial conditions in [0 , × [ − . , .
75] with T = 10 , using the distinguishing criteria in (23). in I . Since the Farey intervals partition (0 , ∞ ), there must be a Farey interval J = ( p l q l , p r q r ) ⊃ I for neighbors p l q l and p r q r . Note that every number in J and thus every number in I must then be adescendent of these parents. Denote the mediant (17) by p m /q m . Without loss of generality, we canassume that p m /q m ∈ I . Every rational in the level (cid:96) + 1 daughter interval ( p l q l , p m q m ) is a descendentof p m q m and since all of these are formed by the mediant operation all of these denominators are largerthan q m . The same is true for the upper interval ( p m q m , p r q r ). Since I ⊂ ( p l q l , p m q m ) ∪ { p m q m } ∪ ( p m q m , p r q r ),all remaining rationals in I have denominator greater than q m . Consequently q min ( I ) = q m , and,moreover, the rational with minimal denominator is unique.This result is encapsulated in Algorithm 2 to give a computation of the smallest denominatorrational in I δ ( x ) (11). For example, this algorithm gives q min ( I − ( π )) = 32085 , p m q m = [3; 7 , , , ,q min (( I − ( e )) = 154257 , p m q m = [2; 1 , , , , , , , , , , , , , , , , . Neither of these are convergents of the continued fraction expansions. Algorithm 2 ignores issues26 k F r ac ti on chaosrotational circles ky –0.3–0.2–0.10.00.10.2 ω Figure 20: (a) Proportion of orbits of the standard nontwist map with frequency map (25) when δ = 0 . y and theparameter k for orbits with x = 0 .
35. As in Fig. 19, T = 10 with distinguishing criteria in (23). of finite precision arithmetic, and is not efficient if the Farey path has a long string of repeatedsymbols. An algorithm that does not have this deficit is given in [CP16].We can obtain some additional understanding of the smallest denominator for the specific casewhen the bounds of the interval I are arbitrary rationals [Siv16], I = ( p l q l , p r q r ) . (27)To find the smallest denominator rational we expand each of the boundary points in their Fareypaths: a = p l q l = a a a . . . a m , b = p r q r = b b b . . . b n , with a i , b i ∈ { L, R } . Then, as shown by [Siv16, Thm. 1], there are three cases:1. When the boundary points of (27) are Farey neighbors, the smallest rational in I is themediant, so q min = q l + q r .2. If one Farey path is a subsequence of the other but they are not neighbors, then the smallestrational is a daughter of the shorter path and an ancestor of the longer. For example, if a = b b . . . b n a n +1 a n +2 . . . a m = ba n +1 . . . a m , then the smallest rational has the path pq = ba n +1 . . . a k , for some k < m . This is the appropriate daughter of b and ancestor of a . Note that when a < b ,then it must be the case that a n +1 = L . If, for example a = bLRL . . . , then bL < a < bLR < b ,so then we set k = n + 2, and obtain pq = bLR .27 lgorithm 1 Compute the Farey path for x ∈ R + assuming exact arithmetic procedure FareyPath (x) i ← while x (cid:54) = 1 doif x < then s i = Lx ← x − x else s i = Rx ← x − end if i ← i + 1 end whileend procedureAlgorithm 2 Find the smallest rational in the interval I δ ( x ) procedure SmallDenom ( x, δ )( n, d ) = ( p l , q l ) = (0 , p r , q r ) = (1 , while | x − nd | ≥ δ do ( n, d ) = ( p l + p r , q l + q r ) (cid:46) Find the mediant if x < n/d then ( p r , q r ) = ( n, d ) (cid:46) x ∈ ( p l q l , nd ) else ( p l , q l ) = ( n, d ) (cid:46) x ∈ [ nd , p r q r ) end ifend whilereturn ( n, d ) (cid:46) The smallest rational is nd end procedure ω Rotational circlesDynamics at k = 0.2
ChaosIslands
Figure 21:
The dynamics of the asymmetric two harmonic map (26) for k = 0 . ψ = 0 . initial conditions in [0 , with T = 10 with distinguishing criteria in (23).
3. If neither path is a subsequence of the other, then the smallest rational is the unique rationalthat is a common ancestor of both on the tree: the longest Farey path for which they agree.For example, if a i = b i for i = 0 , . . . k < min( m, n ), and a k +1 (cid:54) = b k +1 then pq = a a , . . . , a k isthe smallest rational in I .Finally, for an interval bounded by irrationals we can prove the following lemma. Lemma 2 (Smallest Rational in an Irrational Interval) . If I = ( a, b ) , < a < b , a, b ∈ R \ Q , then q min ( I ) is the denominator of the common Farey ancestor of a and b , if there is one; otherwise q min ( I ) = 1 .Proof. Denote the infinite Farey paths of the irrationals by a = a a . . . and b = b b . . . , where a i , b i ∈ { L, R } , and let (cid:96) ∈ N be chosen so that the common ancestor of a and b is p (cid:96) q (cid:96) = a a . . . a (cid:96) = b b . . . b (cid:96) , a (cid:96) +1 (cid:54) = b (cid:96) +1 . If (cid:96) does not exist, then since a < b , a = L and b = R , which means that ∈ I , so that q min = 1.Now suppose that there is a common ancestor of length (cid:96) ≥
1. Then since a < b , we musthave a (cid:96) +1 = L and b (cid:96) +1 = R and a < p (cid:96) q (cid:96) < b . Denote a “left truncation” of a path as a rational29 k F r ac ti on chaosrotational circles ky ω Figure 22: (a) Proportion of orbits of the asymmetric standard map with force (26) and ψ = 0 . y and theparameter k for orbits with x = 0 .
35. As in Fig. 21, T = 10 with distinguishing criteria in (23). a L = a a . . . a j < a and a “right truncation” as a rational a R = a a . . . a k > a , see Fig. 23. Forexample if a j +1 = R , and a k +1 = L then we know that a a . . . a j < a < a a . . . a k . Note thatsuch truncations always exist for any irrational and any choice of minimal length since the infinitepaths with tails . . . L ∞ and . . . R ∞ are rationals. Now, by item (3) above [Siv16, Thm. 1], for theinterval I outer = ( a L , b R ), the smallest denominator is that of the common Farey ancestor of a L and b R : q min ( I outer ) = q (cid:96) . Thus, whenever these rational truncations are both longer than (cid:96) , then I outer contains the common Farey ancestor p (cid:96) q (cid:96) and this has the smallest denominator. Note thatsince a L < a and b R > b , then I ⊂ I outer . Thus q min ( I ) is no less than q (cid:96) . Moreover since p (cid:96) q (cid:96) ∈ I ,then q min ( I ) is no more than q (cid:96) . Thus q min ( I ) = q (cid:96) . References [AC15] C. V. Abud and I. L. Caldas. On Slater’s criterion for the breakup of invariant curves.
Physica D , 308:34–39, 2015. https://doi.org/10.1016/j.physd.2015.06.005 .[ACP06] E.G. Altmann, G. Cristadoro, and D. Paz. Nontwist non-Hamiltonian systems.
Phys.Rev. E , 73(5):056201, 2006. http://link.aps.org/abstract/PRE/v73/e056201 .[BBdA98] S.J. Beslin, D.J. Baney, and V. de Angelis. Small denominators: No small prob-lem.
Mathematics Magazine , 71(2):132–138, 1998. . 30 )) aa L a L b R b R b ( p ℓ q ℓ R RR RRR L LL LL L I = a a a ...LR Figure 23:
An interval I = ( a, b ) bounded by a pair of irrationals with the outer interval ( a L , b R ), and their commonFarey Ancestor p (cid:96) q (cid:96) [BBGT96] R. Bartolini, A. Bazzani, M. Giovanniozzi, and E. Todesco. Tune evaluation in simula-tions and experiments. Particle Accelerators , 52:147–177, 1996. http://cdsweb.cern.ch/record/292773?ln=en .[CB09] E. Charrier and L. Buzer. Approximating a real number by a rational numberwith a limited denominator: A geometric approach.
Discrete Applied Mathemat-ics , 157(16):3473–3484, August 2009. https://linkinghub.elsevier.com/retrieve/pii/S0166218X09000894 .[Chi79] B.V. Chirikov. A universal instability of many-dimensional oscillator systems.
Phys.Rep. , 52(5):263–379, 1979. https://doi.org/10.1016/0370-1573(79)90023-1 .[CP16] M. Citterio and R. Pavani. A fast computation of the best k -digit rational approx-imation to a real number. Mediterranean Journal of Mathematics , 13(6):4321–4331,December 2016. http://link.springer.com/10.1007/s00009-016-0747-z .[DDS +
16] S. Das, C.B. Dock, Y. Saiki, M. Salgado-Flores, E. Sander, J. Wu, and J.A. Yorke.Measuring quasiperiodicity.
Euro. Phys. Lett. , 114:40005, 2016. https://doi.org/10.1209/0295-5075/114/40005 .[dlL01] R. de la Llave. A tutorial on KAM theory. In
Smooth Ergodic Theory and Its Appli-cations (Seattle, Wa, 1999) , volume 69 of
Proc. Sympos. Pure Math. , pages 175–292.Amer. Math. Soc., Providence, 2001.[DSSY16] S. Das, Y. Saiki, E. Sander, and J.A. Yorke. Quasiperiodicity: Rotation numbers.In C. Skiadas, editor,
The Foundations of Chaos Revisited: From Poincar´e to Recent dvancement , Understanding Complex Systems. Springer, Cham, 2016. https://doi.org/10.1007/978-3-319-29701-9_7 .[DSSY17] S. Das, Y. Saiki, E. Sander, and J.A. Yorke. Quantitative quasiperiodicity. Nonlinearity ,30(11):4111, 2017.[DY18] S. Das and J.A. Yorke. Super convergence of ergodic averages for quasiperiodic orbits.
Nonlinearity , 31(2):491–501, 2018. https://doi.org/10.1088/1361-6544/aa99a0 .[EV01] K. Efstathiou and N. Voglis. A method for accurate computation of the rotation andtwist numbers for invariant tori.
Physica D , 158:151–163, 2001.[FM13] A.M. Fox and J.D. Meiss. Greene’s residue criterion for the breakup of invariant tori ofvolume-preserving maps.
Physica D , 243(1):45–63, 2013. https://doi.org/10.1016/j.physd.2012.09.005 .[FM14] A.M. Fox and J.D. Meiss. Critical invariant circles in asymmetric, multiharmonic gener-alized standard maps.
Comm. Nonlinear Sci. Numer. Simulat. , 19(4):1004–1026, 2014. https://doi.org/10.1016/j.cnsns.2013.07.028 .[For07] M Forisek. Approximating rational numbers by fractions. In P. Crescenzi, G. Prencipe,and G. Pucci, editors,
Fun with Algorithms , volume LNCS 4475, pages 156–165.Springer-Verlag, Berlin, 2007. https://doi.org/10.1007/978-3-540-72914-3 .[GJSS87] J.M. Greene, H. Johannesson, B. Schaub, and H. Suhl. Scaling anomaly at the criticaltransistion of an incommensurate structure.
Phys. Rev. A , 36:5858–5861, 1987. https://doi.org/10.1103/PhysRevA.36.5858 .[GM09] G.A. Gottwald and I. Melbourne. On the implementation of the 01 test for chaos.
SIAMJ. Appl. Dyn. Sys.. , 8(1):129–145, 2009. https://doi.org/10.1137/080718851 .[Gre79] J.M. Greene. A method for determining a stochastic transition.
J. Math. Phys. , 20:1183–1201, 1979. https://doi.org/10.1063/1.524170 .[HdlL06] A. Haro and R. de la Llave. A parameterization method for the computation of invarianttori and their whiskers in quasi-periodic maps: Numerical algorithms.
Disc. Cont. Dyn.Sys. , B6(6):1261–1300, 2006. https://doi.org/10.3934/dcdsb.2006.6.1261 .[HW79] G.H. Hardy and E.M. Wright.
An Introduction to the Theory of Numbers . Oxford Univ.Press, Oxford, 1979.[KO86] S. Kim and S. Ostlund. Simultaneous rational approximations in the study of dynamicalsystems.
Phys. Rev. A , 34:3426–3434, 1986. https://doi.org/10.1103/PhysRevA.34.3426 .[LFC92] J. Laskar, C. Froeschl, and A. Celletti. The measure of chaos by the numerical analysis ofthe fundamental frequencies. Application to the standard mapping.
Physica D , 56:253–269, 1992. 32LM10] Z. Levnaji and I. Mezi. Ergodic theory and visualization. I. Mesochronic plots forvisualization of ergodic partition and invariant sets.
Chaos: An Interdisciplinary Journalof Nonlinear Science , 20(3):033114, September 2010. http://doi:10.1063/1.3458896 .[LV09] A. Luque and J. Villanueva. Numerical computation of rotation numbers ofquasi-periodic planar curves.
Physica D , 238(20):2025–2044, 2009. .[Mac93] R.S. MacKay.
Renormalisation in Area-Preserving Maps , volume 6 of
Adv. Series inNonlinear Dynamics . World Scientific, Singapore, 1993.[May88] D.H. Mayer. On the distribution of recurrence times in nonlinear systems.
Lett. Math.Phys. , 16(2):139–143, 1988.[Mei92] J.D. Meiss. Symplectic maps, variational principles, and transport.
Rev. Mod. Phys. ,64(3):795–848, 1992. https://doi.org/10.1103/RevModPhys.64.795 .[MP85] R.S. MacKay and I.C. Percival. Converse KAM: Theory and practice.
Comm. Math.Phys. , 98:469–512, 1985. https://doi.org/10.1007/BF01209326 .[MS91] S. Marmi and J. Stark. On the standard map critical function.
Nonlinearity , 5(3):743–761, 1991. https://doi.org/10.1088/0951-7715/5/3/007 .[MS92] R.S. MacKay and J. Stark. Locally most robust circles and boundary circles forarea-preserving maps.
Nonlinearity , 5:867–888, 1992. http://iopscience.iop.org/0951-7715/5/4/002 .[Sha92] J. Shallit. Real numbers with bounded partial quotients: A survey.
LEnseignementMathematique , 38:151–187, 1992. .[Sim18] C Simo. Some questions looking for answers in dynamical systems.
Disc. Cont. Dyn.Sys. , 38(12):62156239, 2018. https://doi.org/10.3934/dcds.2018267 .[Siv16] I. Sivignon. A note on the computation of the fraction of smallest denominator inbetween two irreducible fractions.
Discrete Applied Mathematics , 202:197–201, 2016. .[Sla50] N.B. Slater. The distribution of the integers N for which { N θ } < (cid:15) . Proc. CambridgePhilos. Soc. , 46:525–534, 1950. https://doi.org/10.1017/S0305004100026086 .[Sla67] N.B. Slater. Gaps and steps for the sequence nθ mod 1. Proc. Cambridge Philos. Soc. ,63:1115–1123, 1967. https://doi.org/10.1017/S0305004100042195 .[SLV05] J.D. Szezech, S.R. Lopes, and R.L. Viana. Finite-time Lyapunov spectrum for chaoticorbits of non-integrable Hamiltonian systems.
Phys. Lett. A , 335(5-6):394–401, 2005. https://doi.org/10.1016/j.physleta.2004.12.058 .[SMS +
18] M.S. Santos, M. Mugnaine, J.D. Szezech, A/M. Batista, I.L. Caldas, M.S. Baptista,and R.L. Viana. Recurrence-based analysis of barrier breakup in the standard nontwistmap.
Chaos , 28(8):085717, 2018. https://doi.org/10.1063/1.5021544 .33SMS +
19] M.S. Santos, M. Mugnaine, J.D. Szezech, A.M. Batista, I.L. Caldas, and R.L. Viana.Using rotation number to detect sticky orbits in Hamiltonian systems.
Chaos ,29(4):043125, 2019. https://doi.org/10.1063/1.5078533 .[SSC +
13] J.D. Szezech, A.B. Schelin, I.L. Caldas, S R. Lopes, P.J. Morrison, and R.L. Viana.Finite-time rotation number: A fast indicator for chaotic dynamical structures.
Phys.Lett. A , 377:452–456, 2013. https://doi.org/10.1016/j.physleta.2012.12.013 .[Ste13] C. L. Stewart. On the distribution of small denominators in the Farey series of order N . In Advances in combinatorics , pages 275–286. Springer, Heidelberg, 2013.[SV06] T.M. Seara and J. Villanueva. On the numerical computation of Diophantine rotationnumbers of analytic circle maps.
Physica D , 217(2):107–120, 2006.
Monthly Weather Review , 104(10):1209–1214, 1978. https://doi.org/10.1175/1520-0493(1976)104<1209:TEOYFF>2.0.CO;2 .[ZTRK07] Y. Zou, M. Thiel, M.C. Romano, and J. Kurths. Characterization of stickiness bymeans of recurrence.
Chaos , 17:043101, 2007. http://link.aip.org/link/?CHAOEH/17/043101/1http://link.aip.org/link/?CHAOEH/17/043101/1