Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 Contents lists available at ScienceDirect Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns Research paper An algorithm to compute rotation intervals of circle mapsR Lluís Alsedà a,b, Salvador Borrós-Cullell c,∗ a Departament de Matemàtiques, Edifici C, Universitat Autònoma de Barcelona, Bellaterra 08193, Barcelona, Spain b Centre de Recerca Matemàtica, Campus de Bellaterra, Edifici C, Universitat Autònoma de Barcelona, Bellaterra 08193, Barcelona, Spain c Departament de Matemàtiques, Edifici Cc, Universitat Autònoma de Barcelona, Bellaterra 08193, Barcelona, Spain a r t i c l e i n f o a b s t r a c t Article history: In this article we present an efficient algorithm to compute rotation intervals of circle Received 11 December 2020 maps of degree one. It is based on the computation of the rotation number of a mono- Revised 19 May 2021 tone circle map of degree one with a constant section. The main strength of this algorithm Accepted 31 May 2021 is that it computes exactly the rotation interval of a natural subclass of the continuous Available online 5 June 2021 non-invertible degree one circle maps. 2020 MSC: We also compare our algorithm with other existing ones by plotting the Devil’s Staircase Primary 37E10 of a one-parameter family of maps and the Arnold Tongues and rotation intervals of some Secondary 37E45 special non-differentiable families, most of which were out of the reach of the existing algorithms that were centred around differentiable maps. Keywords: Rotation number © 2021 The Authors. Published by Elsevier B.V. Circle maps This is an open access article under the CC BY-NC-ND license Nondecreasing degree one lifting (http://creativecommons.org/licenses/by-nc-nd/4.0/) Algorithm 1. Introduction The rotation interval plays an important role in combinatorial dynamics. For example Misiurewicz’s Theorem [1] links the set of periods of a continuous lifting F of degree one to the set M := {n ∈ N : nk ∈ Rot (F ) for some integer k}, where Rot (F ) denotes the rotation interval of F . Moreover, it is natural to compute lower bounds of the topological entropy depending on the rotation interval [2]. In any case, the knowledge of the rotation interval of circle maps of degree one is of theoretical importance. The rotation number was introduced by H. PoincarȨ to study the movement of celestial bodies [3], and since then has been found to model a wide variety of physical and sociological processes. In the physical sense, it has been recently applied to climate science [4]. In the sociological one, the application to voting theory [5,6] is specially surprising in this context. The computation of the rotation number for invertible maps of degree one from S1 onto itself is well studied, and many very efficient algorithms exist for its computation [7–10]. However, there is a lack of an efficient algorithm for the non- invertible and non-differentiable case. In this article, we propose a method that allows us to compute the rotation interval for the non-invertible case. Our algo- rithm is based on the fact that we can compute exactly the rotation number of a natural subclass of the class of continuous non-decreasing degree one circle maps that have a constant section and a rational rotation number. From this algorithm we R Supported by the Spain’s “Agencial Estatal de Investigación” (AEI) grants MTM2017-86795-C3-1-P and MDM-2014-0445 within the “María de Maeztu” Program. ∗ Corresponding author. E-mail addresses: alseda@mat.uab.cat, lalseda@crm.cat (L. Alsedà), sborros@mat.uab.cat (S. Borrós-Cullell). https://doi.org/10.1016/j.cnsns.2021.105915 1007-5704/© 2021 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/) L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915  Fig. 1. An example of a map from L1 which can be considered as a toy model for the elements of that class. The picture shows F  , and F is globally    [0,1] defined as F (x ) = F  {{x}} + x. [0,1] get an efficient way to compute exactly the rotation interval of this subset of the continuous non-invertible degree one circle maps by using the so called upper and lower maps, which, when different, always have a constant section. When dealing with maps outside the aforementioned class, the algorithm will return an arbitrarily precise rational aproximation of the rotation number. To check the efficiency of our algorithm we will use it to compute some classical results such as a Devil’s Staircase. When doing so, we will compare the efficiency of our algorithm with the performance of some other algorithms that have been traditionally used under the hypothesis of non-invertibility. On the other hand, we will also compute the rotation interval and Arnold tongues for a variety of maps, in the same comparing spirit. These maps include the Standard Map and variants of it but have issues either with the differentiability, or even with the continuity. Of course these variants are not well suited for algorithms that strongly use differentiability. The paper is organised as follows. In Section 2 the theoretical background will be set. In Section 3 the algorithm will be presented, and in Section 4 we will provide the mentioned examples of the use of the algorithm. Finally in Section 5 we will discuss the advantages and disadvantages of the proposed algorithm. 2. A short survey on rotation theory and the computation of rotation numbers We will start by recalling some results from the rotation theory for circle maps. To do this we will follow [11]. The floor function (i.e. the function that returns the greatest integer less than or equal to the variable) will be denoted as ·. Also the decimal part of a real number x ∈ R, defined as x − x ∈ [0, 1 ) will be denoted by {{x}}. In what follows S1 denotes the circle, which is defined as the set of all complex numbers of modulus one. Let e : R → S1 be the natural projection from R to S1 , which is defined by e(x ) := exp(2π ix ). Let f : S1 → S1 be continuous map. A continuous map F : R → R is a lifting of f if and only if e(F (x )) = f (e(x )) for every x ∈ R. Note that the lifting of a circle map is not unique, and that any two liftings F and F  of the same continuous map f : S1 → S1 verify F = F  + k for some k ∈ Z. For every continuous map f : S1 → S1 there exists an integer d such that F (x + 1 ) = F (x ) + d for every lifting F of f and every x ∈ R (that is, the number d is independent of the choice of the lifting and the point x ∈ R). We shall call this number d the degree of f . The degree of a map roughly corresponds to the number of times that the whole image of the map f covers homotopically S1 . In this paper we are interested studying maps of degree 1, since the rotation theory is well defined for the liftings of these maps.  denote the set of all liftings of maps of degree 1 by L1 . Observe We will   that to define a map from L1 it is enough to define F  (see Fig. 1) because F can be globally defined as F (x ) = F  {{x}} + x for every x ∈ R. [0,1] [0,1] Remark 2.1. It is easy to see that, for every F ∈ L1 , F n (x + k ) = F n (x ) + k for every n ∈ N, x ∈ R and k ∈ Z. Consequently, F n ∈ L1 for every n ∈ N. 2 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 Fig. 2. An example of a map F ∈ L1 with its lower map Fl in red and its upper map Fu in blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Definition 2.2. Let F ∈ L1 , and let x ∈ R. We define the rotation number of x as F n (x ) − x ρF (x ) := lim sup . n→∞ n Observe (Remark 2.1) that, ρF (x ) = ρF (x + k ) for every k ∈ Z. The rotation set of F is defined as: Rot (F ) = {ρF (x ) : x ∈ R} = {ρF (x ) : x ∈ [0, 1]}. Ito [12], proved that the rotation set is a closed interval of the real line. So, henceforth the set Rot (F ) will be called the rotation interval of F . Proposition 2.3 (Proposition 3.7.11 in [11]). Let F ∈ L1 be non-decreasing. Then, for every x ∈ R the limit F n (x ) − x lim n→∞ n exists and is independent of x. For a non-decreasing map F ∈ L1 , the number ρF (x ) = limn→∞ F (xn)−x will be called the rotation number of F , and will n be denoted by ρF . Now, by using the notation from [11], we will introduce the notion of upper and lower functions, that will be crucial to compute the rotation interval. Definition 2.4. Given F ∈ L1 we define the F -upper map Fu as Fu (x ) := sup{F (y ) : y ≤ x}. Similarly we will define the F -lower map as Fl (x ) := inf{F (y ) : y ≥ x}. An example of such functions is shown in Fig. 2. It is easy to see that Fl , Fu ∈ L1 are non decreasing, and Fl (x ) ≤ F (x ) ≤ Fu (x ) for every x ∈ R. The rationale behind introducing the upper and lower functions comes from the following result, stating that the rotation interval of a function F ∈ L1 is given by the rotation number of its upper and lower functions. Theorem 2.5 (Theorem 3.7.20 in [11]). Let F ∈ L1 . Then,   Rot (F ) = ρFl , ρFu . Note that this theorem makes indeed sense, since the upper and lower functions are non-decreasing and by Proposition 2.3 they have a single well defined rotation number. 3 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 Let f : S1 → S1 and let z ∈ S1 . The f -orbit of z is defined to be the set Orb f (z ) := {z, f (z ), f 2 (z ), . . . , f n (z ), . . . }. We say that z is an n-periodic point of f if Orb f (z ) has cardinality n. Note that this is equivalent to f n (z ) = z and f k (z ) = z for every k < n. In this case the set Orb f (z ) will be called an n-periodic orbit (or, simply, a periodic orbit). If we have a periodic orbit of a circle map, a natural question that might arise is how it behaves at a lifting level. This motivates the introduction of the notion of a lifted cycle. Given a set A ⊂ R and m ∈ Z we will denote A + m := {x + m : x ∈ A}. Analogously, we set A + Z := {x + m : x ∈ A, m ∈ Z}. Definition 2.6. Let f : S1 → S1 be a continuous map and let F be a lifting of f. A set P ⊂ R is called a lifted cycle of F if e(P ) is a periodic orbit of f. Observe that, then P = P + Z. The period of a lifted cycle is, by definition, the period of e(P ). Hence, when e(P ) is an n-periodic orbit of f,P is called an n-lifted cycle, and every point x ∈ P will be called an n-periodic (mod 1 ) point of F . The relation between lifted orbits and rotation numbers is clarified by the next lemma. Lemma 2.7 (Lemmas 3.7.2 and 3.7.3 in [11]). Let F ∈ L1 . Then, x is an n-periodic (mod 1 ) point of F if and only if there exists k ∈ Z such that F n (x ) = x + k but F j (x ) − x ∈ / Z for j = 1, 2, . . . , n − 1. In this case, F m (x ) − x k ρF (x ) = mlim = . →∞ m n Moreover, let P be a lifted n-cycle of F . Every point x ∈ P is an n-periodic (mod 1 ) point of F , and the above number k does not depend on x. Hence, for every x ∈ P we have ρF (P ) := ρF (x ) = nk . Now we can revisit Proposition 2.3: Proposition 2.3 (Proposition 3.7.11 in [11]). Let F ∈ L1 be non-decreasing. Then, for every x ∈ R the limit F n (x ) − x ρF := nlim →∞ n exists and is independent of x. Moreover, ρF is rational if and only if F has a lifted cycle. In the next two subsections we will survey on two known algorithms that have been already used to compute rotation numbers of non-differentiable and non-invertible liftings from L1 . The first one (Algorithm 1) stems automatically from the definition of rotation number (Definition 2.2); the other one (Algorithm 2) is due to Simó et al. [13]. Algorithm 1 Direct Algorithm pseudocode procedure Rotation_Number   (F , error) 1 n ← CEIL error x←0 k←0 for i ← 1, n do x ← F (x ) s ← floor(x ) k←k+s  k = F n ( 0 ) x←x−s  x = {{F n (0 )}} = F n (0 ) − k end for return k+n x end procedure 2.1. Algorithm 1: the numerical algorithm to compute the rotation interval that stems from the definition of rotation number The first algorithm to compute ρF consists in using Proposition 2.3 and the following approximation: F m (x ) − x ρF = mlim →∞ m  F n (x ) − x  ≈ n  x=0 F n (0 ) = . n 4 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 Algorithm 2 Simó et al. ([13]) Algorithm in pseudocode procedure Rotation_Number(F , n) index[] ← x←0 ρmin ← 0 ρmax ← 1 for i ← 0, n do x ← F (x ) ki ← floor(x) αi ← x − k i index[i] ← i end for sort α [index[i]] by rearranging index[] for i ← 0, n − 1 do kindex[i + 1] −kindex[i] ρaux ← index[ i + 1]−index[i] if index[i + 1] > index[i] then ρmin ← max{ρmin , ρaux } else ρmax ← min{ρmax , ρaux } end if end for return ρmin , ρmax end procedure However, a priori we do not know how good the convergence is. In Lemma 2.8 we will show that the error is of order 1/n. The implementation of the computation of this approximation to the rotation number can be found in the side algorithm pseudocode. Since the maps from L1 are defined so that    F (x ) = F  {{x}} + x, [0,1] we need to evaluate the function FLOOR(· ) = · once per iterate. So, for clarity and efficiency, it seems advisable to split F n (0 ) as {{F n (0 )}} + F n (0 ). The next lemma clarifies the computation error as a function of the number of iterates. In particular it explicitly gives the necessary number of iterates, given a fixed tolerance. For every non-decreasing lifting F ∈ L1 , and every n ∈ N we set (see Fig. 3) F (n ) := min F n (x ) − x = min F n (x ) − x. x∈R x∈[0,1] The second equality holds because F has degree 1, and hence F (n ) is well defined. Lemma 2.8. For every non-decreasing lifting F ∈ L1 and n ∈ N we have (a) either F n (z ) = z + F (n ) + 1 for some z ∈ R, or x + F (n ) ≤ F n (x ) < x + F (n ) + 1 for every x ∈ R;  (n )  (n )+1 (b)  F n ≤ ρF ≤  F n ; and  F n (x )−x  (c) ρF − n  < 1 n for every x ∈ R. Proof. We will prove the whole lemma by considering two alternative cases. Assume first that F n (z ) = z + F (n ) + 1 for  (n )+1 some z ∈ R. Then (a) holds trivially, and Proposition 2.3 and Lemma 2.7 imply that ρF = F n . So, Statement (b) also holds in this case. Now observe that from the definition of F (n ) we have F ( n ) ≤ F n ( x ) − x ≤ F n ( x ) − x (1) for every x ∈ R. Moreover, there exists k = k(x ) ∈ Z such that x ∈ [z + k, z + k + 1 ) and, since F is non-decreasing, so is F n . Thus, F n (x ) − x ≤ F n (z + k + 1 ) − x = F n (z ) + k + 1 − x = F ( n ) + 1 + ( z + k + 1 − x ) < F ( n ) + 2, by Remark 2.1. Consequently, 1 F ( n ) F n (x ) − x 1 ρF − = ≤ < ρF + ; n n n n which proves (c) in this case. 5 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 Fig. 3. Plot of x + F (n ) x + F (n ) + 1, and F n (x ) for two arbitrary non-decreasing maps F ∈ L1 that fit in the two cases of the lemma. Now we consider the case F n ( x ) = x + F ( n ) + 1 for every x ∈ R. In view of the definition of F (n ), we cannot have F n ( x ) − x > F ( n ) + 1 for every x ∈ R. Hence, by the continuity of F n (x ) − x and (1), F ( n ) ≤ F n ( x ) − x < F ( n ) + 1 (2) for every x ∈ R. This proves (a). Now we prove (b). We consider the functions: x −→ F (n ) + x,F n , and x −→ F (n ) + 1 + x. They are all non-decreasing and, by Remark 2.1, they belong to L1 . Hence, by Proposition 2.3, [11, Lemma 3.7.19] and (2), F (n ) = ρx→ ≤ ρF n ≤ ρx→ = F ( n ) + 1. F (n )+x F (n )+1+x Consequently, F ( n ) ρ n F ( n ) + 1 ≤ ρF = F ≤ , n n n and (b) holds. Moreover, (2) is equivalent to F ( n ) F n (x ) − x F ( n ) + 1 ≤ ≤ , n n n which proves (c).  2.2. Algorithm 2: the Simó et al. algorithm to compute the rotation interval First of all, it should be noted that even though the authors propose an algorithm to compute the rotation interval for a general map F ∈ L1 , we will only use it for non decreasing maps. A priori this algorithm is radically different from Algorithm 1 and it gives an estimate of ρF by providing and upper and a lower bound of the rotation number (rotation interval in the original paper) of F . Moreover, it is implicitly assumed that ρF ∈ [0, 1] (in particular that F (0 ) ∈ [0, 1 ) — this can be achieved by replacing the lifting F by the lifting G := F − F (0 ), if necessary). The algorithm goes as follows: (Alg. 2-1) Decide the number of iterates n in function of a given tolerance. (Alg. 2-2) For i = 0, 1, 2, . . . , n compute ki = F i (x0 ) and αi = F i (x0 ) − ki (i.e. αi is the fractionary part of F i (x0 )). 6 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 (Alg. 2-3) Sort the values of ki and αi so that αi0 < αi1 < . . . < αin (this can be achieved efficiently with the help of an index vector). (Alg. 2-4) Initialise ρmin = 0 and ρmax = 1. ki −ki (Alg. 2-5) For j = 0, 1, 2, . . . , n − 1 set ρ j = j+1 i j+1 −i j j , and if i j+1 > i j set ρmin = max{ρmin , ρ j }; otherwise, • • if i j+1 < i j set ρmax = min{ρmax , ρ j }. (Alg. 2-6) Return ρmax and ρmin as upper and lower bounds of the rotation number of F , respectively. The  real  issue in this algorithm consists in dealing with the error. If the rotation number ρF satisfies a Diophantine condition ρ − p  ≤ cq−ν , with c > 0 and ν ≥ 2, then the error verifies F q 1 ε< . (cnν ) ν −1 1 Note that this error depends strongly on the chosen number n of iterates, and that n must be chosen before knowing what the rotation number could possibly be. Hence Algorithm 2 it is not well suited to compute unknown rotation numbers of L1 maps. However, it is excellent in continuation methods where the current rotation number gives a good estimate of the next one. Remark 2.9. Note that the original aim of the algorithm to determine the existence of closed invariant curves on dynamical systems on the plane rather than the computation of rotation numbers of a given map of the circle. The rationale of the algorithm is that if, after computing ρmin and ρmax , we find that ρmin > ρmax then the computed orbit cannot lay on a closed invariant curve. This explains most of the limitations we have encountered, such as the lack of an a priori estimate of the error, or the fact that the algorithm is suited only for rotation numbers ρ ∈ [0, 1]. 3. An algorithm to compute rotation numbers of non-decreasing maps with a constant section The diameter of an interval K which, by definition is equal to the absolute value of the difference between their endpoints, will be denoted as diam(K ). A constant section of a lifting F of a circle map is a closed non-degenerate (i.e. different from a point or, equivalently, with non-empty interior, or such that diam(K ) > 0) subinterval K of R such that F  is constant. In the special case when K F ∈ L1 , we have that F (x + 1 ) = F (x ) + 1 = F (x ) for every x ∈ R. Hence, diam(K ) < 1. The algorithm we propose is based on Lemma 2.8 but, specially, on the following simple proposition which allows us to compute exactly the rotation number of a non-decreasing lifting from L1 that has a constant section, provided that F n (K ) ∩ (K + Z ) = ∅. In this sense, Proposition 3.1 has a completely different strategical aim than Algorithm 1 and Lemma 2.8, which try to (costly) estimate the rotation number. Proposition 3.1. Let F ∈ L1 be non-decreasing and have a constant section K. Assume that there exists n ∈ N such that F n (K ) ∩ (K + Z ) = ∅, and that n is minimal with this property. Then, there exists ξ ∈ R such that F n (K ) = {ξ } ⊂ K + m with m = ξ − min K  ∈ Z,ξ is an n-periodic (mod 1 ) point of F , and ρF = m n. Proof. Since K is a constant section of F ,F (K ) contains a unique point, and hence there exists ξ ∈ R such that F n (K ) = {ξ }. Then, the fact that F n (K ) ∩ (K + Z ) = ∅ implies that ξ ∈ K + m with m = ξ − min K  ∈ Z.   Set ξ := ξ − m ∈ K. Then, F n ξ = F n (K ) = ξ + m . Moreover, the minimality of n implies that F j ξ − ξ ∈ / Z for j = 1, 2, . . . , n − 1. So, Lemma 2.7 tells us that ξ (and hence ξ ) is an n-periodic (mod 1 ) point of F . Thus, ρF = m n by Proposition 2.3.  As already said, Proposition 3.1 is a tool to compute exactly the rotation numbers of non-decreasing liftings F ∈ L1 which have a constant section and have a lifted cycle intersecting the constant section (and hence having rational rotation number). In the next subsection we shall investigate how restrictive are these conditions, when dealing with computation of rotation numbers. 3.1. On the genericity of Proposition 3.1 First observe that the fact that Proposition 3.1 only allows the computation of rotation numbers of non-decreasing liftings F ∈ L1 which have a constant section is not restrictive at all. Indeed, if we want to compute rotation intervals of non- invertible continuous circle maps of degree one, Theorem 2.5 tells us that this is exactly what we want. Clearly, one of the real restrictions that cannot be overcome in the above method to compute exact rotation numbers is that it only works for maps having a rational rotation number. Hence, for maps with non-rational rotation number we can only hope to get a rational approximation like the one given by Algorithm 1, which can be archieved with arbitrary precision. However, as stated by Theorem 4.3, it is not easy to find maps with irrational rotation number, even harder to do so with floating point aproximation. 7 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 Fig. 4. Example of a non-decreasing lifting in L1 with a constant section and rational rotation number which does not verify the assumptions of Proposition 3.1. On the other hand, we also have the formal restriction that Proposition 3.1 requires that the map F has a lifted cycle intersecting the constant section (indeed this is a consequence of the condition F n (K ) ∩ (K + Z ) = ∅). A natural question is whether this restriction is just formal or it is a real one. In the next example we will see that the restriction is not superfluous since there exist maps which do not satisfy it. Consequently, Proposition 3.1 is useless in computing the rotation numbers of non-decreasing liftings in L1 which have a constant section and either irrational rotation number or rational rotation number but do not have any lifted cycle inter- secting the constant section. The only reasonable solution to these problems is to use an iterative algorithm to estimate the rotation number with a prescribed error, such as Algorithm 1, Algorithm 2 or others. Example 3.2. There exist non-decreasing liftings in L1 which have a constant section and rational rotation  number  but do not have any lifted cycle intersecting the constant section: Let F ∈ L1 be as in Figure 4 such that F (x ) = F  {{x}} + x for [0,1] every x ∈ R. 1 The map F is a non-decreasing lifting from L1 , having a constant section K = [0.8, 1] and rotation number 3 given by the 3-lifted cycle P = {0.1, 0.3, 0.4} + Z (c.f. Lemma 2.7 and Proposition 2.3). Now let us see that F does not have any lifted cycle intersecting the constant section. First, observe that F 3 (K ) = F (F (F (K ))) = F (F ({1.2} )) = F ({1.35} ) = {1.75} ⊂ K + Z. Hence, there is no lifted cycle of period 3 intersecting K. On the other hand, again by Lemma 2.7, we have that if x is an n-periodic (mod 1 ) point of F then there exists k ∈ Z such that F n (x ) = x + k and 1 F m (x ) − x k = ρF = lim = ρF ( x ) = . 3 m→∞ m n Moreover, since F is non-decreasing, we know by [11, Corollary 3.7.6] that n and k must be relatively prime. Thus, any lifted cycle of F has period 3, and from above this implies that there is no lifted cycle intersecting K. 3.2. Algorithm 3: A constant section based algorithm arising from Proposition 3.1 From the last paragraph of the previous subsection it becomes evident that Proposition 3.1 does not give a complete algorithm to compute rotation numbers of non-decreasing liftings in L1 which have a constant section. Such an algorithm must rather be a mix-up of Proposition 3.1, and Algorithm 1 to be used when we are not able to determine whether we are in the assumptions of that proposition. As we did for Algorithm 1, we will split F n (0 ) as {{F n (0 )}} + F n (0 ). The goal is twofold, on the one hand splitting helps minimizing the truncation errors. On the other hand, thanks to the splitting we can apply Proposition 3.1 more efficiently, since it requires the computation of m as an integer part. Note that here we are denoting the constant section by K and assuming that 0 ∈ K, which will be justified later. Then, observe that the computations to be performed are exactly the same in both cases (meaning when we can use Proposition 3.1, and when al- ternatively we must end up by using Algorithm 1); except for the conditionals that check whether there exists n ≤ max_iter such that F n (K ) ∩ (K + Z ) = ∅ is verified (that is, whether the assumptions of Proposition 3.1 are verified) before exhausting the max_iter iterates determined a priori. In what follows F n (0 ) will denote the computed value of F n (0 ) with rounding errors for n = 1, 2, . . . , max_iter. The algorithm goes as follows (see Algorithm 3 for a full implementation in pseudocode, and see the explanatory com- ments below): (Alg. 3-1) Re-parametrize the lifting F so that it has a maximal (with respect to the inclusion relation) constant section of the form [−tol, β + tol], where tol is the pre-defined rounding error bound. 8 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 Algorithm 3 Constant Section Based Algorithm For a non-decreasing map F ∈ L1 parametrised so that [−tol, β + tol] is a constant section of F define tol ←  Procedure parameter that bounds the rounding errors in thecomputation of F n (0 )  (F , β , error) procedure Rotation_Number  1 max_iter ← CEIL error  Maximum number of iterates al- lowed(to estimate the rotation x←0 number with the prescribed error m←0 when reached) for n ← 1, max_iter do x ← F (x ) s ← floor(x ) m←m+s  m = F n ( 0 ) x←x−s  x= {{F n (0 )}} = F n (0 ) − m if x ≤ β then return m n  Exact rotation number:Proposition 3.1 holds as- suming that therounding error of F n (0 )is smaller end if than tol end for m+x return max_iter  We do not know whether we are in the assumptions ofProposition 3.1.So, we itera- end procedure tively estimate the rotation number as in- Algorithm 1.The error bound is given by Lemma 2.8 (Alg. 3-2) Set the inputs of the algorithm: • β as in step 1, • F , the map from which we want the rotation number, • error, the maximum error we want our approximation to have.   1 (Alg. 3-3) Decide the maximum number of iterates max_iter = CEIL error to perform in the worst case (i.e. when Proposition 3.1 does not work). (Alg. 3-4) Initialize x = 0 and m = 0. (Alg. 3-5) Compute iteratively x = {{Fn (0 )}} and m = F n (0 ) (so that F n (0 ) = x + m) for n ≤ max_iter. (Alg. 3-6) Check whether x ≤ β . On the affirmative we are in the assumptions of Proposition 3.1, and thus, ρF = m n . Then, the algorithm returns this value as the “exact” rotation number. (Alg. 3-7) If we reach the maximum number of iterates (i.e. n = max_iter) without being in the assumptions of Proposition 3.1 (i.e. with x > β for every x) then, by Lemma 2.8, we have        n (0 )   m+x   F   F n (0 )  1 ρF − max_iter  = ρF − max_iter  ≈ ρF − max_iter  < max_iter ,   and the algorithm returns m+x max_iter as an estimate of ρF with 1 max_iter as the estimated error bound. Remark 3.3. The fact  that we can only check whether the assumptions of Proposition 3.1 are verified before exhausting 1 themax_iter = CEIL error iterates determined a priori does not allow to take into account that F may have a lifted cycle intersecting the constant section but of very large period, i.e. with period larger than max_iter. In practice this problem is totally equivalent to the non-existence (or rather invisibility) of a lifted cycle intersecting the constant section, and it can be considered as a new (algorithmic) restriction to Proposition 3.1. It is solved in (Alg. 3-6) in the same manner as the two other problems related with the applicability of Proposition 3.1 that have already been discussed: by estimating the rotation number as in Algorithm 1. In the last part of this subsection we are going to discuss the rationale of (Alg. 3-2) (and, as a consequence of (Alg. 3-5)). The necessity of this tuning of the algorithm comes again from a challenge concerning the application of Proposition 3.1, which turns to be one of the most relevant restrictions in the use of that proposition. We will begin by discussing how we can efficiently check the condition ξ = F n (0 ) ∈ K + Z (or equivalently F n (K ) ∩ (K + Z ) = ∅) by taking into account that the computation of F (x ) is done with rounding errors, and thus we do not know the exact values of F n (0 ) forn = 1, 2, . . . , max_iter. The next example shows the problems arising in this situation. 9 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 3 Fig. 5. The graph of F 5 . It lies below the graph of x −→ x + 2 but very close to it at five F -preimages of x = 4 . Example 3.4. Fn (0 ) ∈ K + Z but F n (K ) ∩ (K + Z ) = ∅, and this leads to a completely wrong estimate of ρ .    F Let F ∈ L1 be the map such that F (x ) = F  {{x}} + x for every x ∈ R, and let [0,1] 4    x + μ if x ∈ 0, 3 , F (x ) := 3 3 4  [0,1] 1+μ if x ∈ 4 ,1 , with μ = 819 3124 − 10−16 .   For this map F we have K = − 34 , 0 and (see Fig. 5) the graph of F 5 lies above the graph of x −→ x + 1 and below the 3 graph of x −→ x + 2, but very close to it at five F -preimages of x = 4. On the other hand, F (0 ) = 1.74999999999999887 · · · ∈ 5 / K+Z but the distance between F 5 (0 ) and K + Z is 47 − F 5 (0 ) ≈ 1.138 · 10−15 . Should the computations be done with rounding errors of this last magnitude, we may have F 5 (0 )  7 , and accept erroneously that F 5 (0 ) ∈ K + Z. This would lead to the 4 conclusion that ρF = 5 but, as it can be checked numerically, ρF ≈ 0.3983 which is far from 25 . 2 At a first glance this seems to be paradoxical but, indeed, it can be viewed in the following way: The graph of F 5 does not intersect the diagonal (modulo 1) x + 2, but there is a map G close (at rounding errors distance) to F such that the graph of G5 intersects that diagonal, and this gives a lifted periodic orbit of period 5 and rotation number 25 for G. On the other hand, nothing is granted about the modulus of continuity of ρF as a function of F (notice that that everything here is continuous including the dependence of the rotation number of F on the parameter μ), and this example explicitly shows that it may be indeed very big. In short, close functions can have very different rotation numbers. The most reasonable solution to the problem pointed out in the previous example consists in restricting the size of K depending of an a priori estimate of the rounding errors in computing F n (0 ) for n = 1, 2, . . . , max_iter. Thus, we denote by 10 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 tol an upper bound of these rounding errors, so that    n n (0 ) ≤ tol F (0 ) − F  holds for n = 1, 2, . . . , max_iter, and, given a maximal (with respect to the inclusion relation) constant section K such that 0 ∈ K we write K := [α − tol, β + tol]. Then observe that the condition Fn (0 ) ∈ [α , β ] + m for some n ∈ N and m ∈ Z implies ξ = F n (0 ) ∈ K + m, and ρ = m F n by Proposition 3.1. In practice, this “rounding errors free” version of the algorithm imposes a new restriction to the applicability of Proposition 3.1 (in the sense that it reduces even more the class of functions for which we can get the “exact rotation number”). However, as before, the rotation numbers of the maps in the assumptions of Proposition 3.1 for which we cannot compute the “exact rotation number” can be estimated as in Algorithm 1. The computational efficiency of the algorithm strongly depends on how we check the condition F n (0 ) ∈ K + Z. Taking into account the above considerations and improvements of the algorithm, this amounts checking whether α +  ≤ F n (0 ) ≤ β +     for some  ∈ Z, and we have to do so by using x = {{F (0 )}} and m = F (0 ) instead of F (0 ) = x + m, which is the n n n algorithmic available information. Checking whether α +  ≤ Fn (0 ) ≤ β +  for some  ∈ Z is problematic since it requires at least two comparisons, and moreover in general  = m (and thus we need some more computational effort to find the right value of ). A very easy solution to this problem is to change the parametrization of F so that α = 0. In this situation we have m = m + α ≤ F n (0 ), m + β < m + 1 because diam(K ) < 1, and m = F n (0 ) . Consequently, α +  ≤ F n (0 ) ≤ β +  for some  ∈ Z is equivalent to =m and x ≤ β . Thus, by “tuning” F so that α = 0 we get that  = m and we manage to determine whether F n (0 ) ∈ [α , β ] + m just with one comparison (x ≤ β ). To see how we can change the parametrization of F (that is step Alg. 3-1) so that α = 0 consider the map G(x ) := F (x + α ) − α . Clearly, F and G are conjugate by the rotation of angle α : x −→ x + α . Then, obviously, G is a non-decreasing map in L1 , has a constant section [−tol, β − α + tol], and ρF = ρG . So, every lifting can be replaced by one of its re- parametrizations with the same rotation number and constant section [−tol, β + tol], where β < 1 − 2tol. 4. Testing the Algorithm In this section we will test the performance of Algorithm 3 by comparing it against Algorithms 1 and 2 when dealing with different usual computations concerning rotation intervals. First we will compare the efficiency of the three algorithms in computing and plotting Devil’s Staircases. Afterwards we will plot rotation intervals and Arnold tongues for two bi- parametric families that mimic the standard map family. In the latter two cases, we will try to compare our algorithm with Algorithms 1 and 2 whenever possible. 4.1. Computing Devil’s staircases In this subsection we will perform the comparison of algorithms by computing and plotting the Devil’s staircase for the parametric family Fμ ⊂ L1 defined as μ∈[0,1] Definition 4.1.    Fμ (x ) = Fμ  {{x}} + x, [0,1] where (see Fig. 1)   4 x + μ if x ≤ 34 Fμ  (x ) = 3 . (3) [0,1] μ + 1 if x > 34 Before doing this we shall remind the notion of a Devil’s Staircase, and why typically exist for such families. To this end we will first recall and survey on the notion of persistence of a rotation interval. Definition 4.2. Given a subclass A of L1 , we say that F ∈ A has an A-persistent rotation interval if there exists a neighbour- hood U of F in A such that Rot (G ) = Rot (F ) for every G ∈ U. We can now state the Persistence Theorem (c.f. [14]): Theorem 4.3 (Persistence Theorem). Let A be a subclass of L1 . Then the following statements hold: 11 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 Fig. 6. The Devil’s Staircase associated to the family (3) computed with Algorithm 3 (upper picture). The lower pictures show the plots of the differences between the value of ρFμ computed with Algorithm 3 and the value of ρFμ computed with Algorithm 1 (left picture), and with the value of ρFμ computed with Algorithm 2(right picture). (a) The set of all maps with A-persistent rotation interval is open and dense in A (in the topology of A). (b) If F has an A-persistent rotation interval, then ρF and ρFu are rational. l Remark 4.4. If we apply Theorem 4.3 to our family Fμ which verifies that the rotation number of Fμ exists for μ∈[0,1] every μ ∈ [0, 1], we have that the set of parameters μ ∈ [0, 1] for which we have irrational rotation number has measure 0. Furthermore, for any κ ∈ Q such that there exists μ with ρFμ = κ , there exists an interval [α , β ]  μ such that for all η ∈ [α , β ],ρFη = κ . The so-called Devil’s staircase is the result of plotting the rotation number as a function of the parameter μ. By Theorem 4.3 we have that this plot will have constant sections for any rational rotation number, hence the ”Staircase” in the name. To test the algorithms, a μ-parametric grid computation of ρFμ with μ ranging from 0 to 1 with a step of 10−5 has been done. For Algorithms 1 and 3 the error has been set to 10−6 . For Algorithm 3 the tolerance has been set to 10−10 . For Algorithm 2 we have arbitrarily set the number of iterates to 10 0 0. In Fig. 6 we show a plot of the Devil’s Staircase computed with Algorithm 3, and the plots of the differences between ρFμ computed with Algorithms 3 and 1, and the differences between ρFμ computed with Algorithms 3 and 2 Table 1. 12 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 Table 1 Performance of the three algorithms studied for a variety of problems. The cells marked with N/A in blue remark that Algorithm 2 does not work in general for ρ ∈/ [0, 1]. The ones marked with N/A in red denote that the computation lasted more than a 100 processor hours and thus was terminated before it ended. Function Time taken by algorithm (s) Problem Family Classic Simó et al. Proposed Devil’s Staircase Fμ (Def. 4.1) 2425.25 210.648 0.1413 Rotation Standard 354.868 N/A 3.2874 Interval PWLSM (Def. 4.7) 110.892 N/A 0.4737 DSM (Def. 4.8) 63.588 N/A 0.2463 Arnol’d Standard N/A N/A 14948.41 Tongues PWLSM N/A N/A 9729.17 DSM N/A N/A 4562.75 shows the times1 taken by each of the three algorithms in computing the whole Devil’s staircase using the three algo- rithms studied. We remark that in the computation of the Devil’s Staircase, Algorithm 3 has been reduced to Algorithm 1 only for μ = 0 and for μ = 1, as one would expect, since these cases follow the pattern of Example 3.4. As a part of the testing of the algorithms we have also considered the inverse problem: Given a value ρ ∈ R\Q and a tol- erance ε > 0 find the value μ = μ(x ) such that ρFμ ∈ [ρ − ε , ρ + ε ]. This problem has turned to be extremely ill-conditioned: √ by choosing ρ to be any irrational number. We have tried to use algebraic numbers such as the golden mean or 1/ 2 and transcendental numbers such as π /4 or e/3. In any studied case, the continuity module of the function μ → ρFμ around μ(ρ ) was estimated to be at least 1025 , making any attempt to solve the problem unfeasable. Note that this is coherent with Theorem 4.3, since the values of μ that give non-rational values are nowhere dense. 4.2. Rotation intervals for standard-like maps In this subsection we test our algorithm by efficiently computing the rotation intervals and some Arnol’d tongues for three bi-parametric families of maps: the standard map family and two piecewise-linear extensions of it; one continuous but not differentiable, and another one which is not even continuous. We emphasize that the usual algorithms such as the ones from [8–10,15] cannot be used for these last two families while the one we propose here it works very well indeed. First we will recall the notion of Arnol’d tongue. Definition 4.5 (Arnol’d Tongue [16]). Let Fa,b be a two-parameter family of maps in L1 for which the rotation in-   (a,b)∈P terval Rot Fa,b is well defined for every possible point (a, b) ∈ P in the parameter set. Given a point ∈ R we define the −Arnold Tongue of Fa,b as (a,b)∈P   T = {(a, b) ∈ P : ∈ Rot Fa,b } ⊂ P. Next we introduce each of the three families that we study and, for each of them we show the results and we explain the performance of the algorithm. Definition 4.6 (Standard Map). S ,a ∈ L1 is defined as (see Fig. 7): a S ,a (x ) := x + − sin(2π x ). (4) 2π To compute   the rotation   intervals of S ,a we will use Theorem 2.5, together with Algorithm 3. To this end, first we will compute S ,a and S ,a (that is, the lower and upper maps of S ,a ), and then we will use Algorithm 3 to compute the l u rotation numbers ρ(S ) and ρ(S ) of these maps. ,a l ,a u     Note that S ,a is non-invertible for a > 1. Hence, in this case, S ,a and S ,a do not coincide and have constant l u sections. However, the characterization of these constants sections is not straightforward, since their endpoints have to be computed numerically. This is the reason why the computations of the rotation intervals and Arnol’d tongues for the standard map have been the slowest ones. In Fig. 8 we show some graphs of the rotation interval and Arnol’d tongues for the Standard Map. The graphs of the rotation intervals are plotted for three different values of as a function of the parameter a. 1 The simulations have been done with an Intel® CoreTM i7-3770 CPU @3.4GHz. 13 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 Fig. 7. The standard map with a = 2π and = 0, with its lower map in blue and its upper map in red. Fig. 8. Graphs of the rotation interval and Arnol’d tongues for the Standard Map S ,a . The graphs of the rotation intervals are plotted as a function of the parameter a. 14 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 Fig. 9. The piecewise-linear standard map T ,a with a = 52π and = 0. The lower map of T ,a is drawn in blue, and the upper map in red. Definition 4.7 (Piecewise-linear standard map). We start by defining a convenience map τ : [0, 1] → [−1, 1] as follows: ⎧  1 ⎨ 4x when x ∈ 0, ,  1 43  τ (x ) = 2 − 4x when x ∈ 4 , 4 , and (5) ⎩ 3  4(x − 1 ) when x ∈ 4 ,1 . Then, the piecewise-linear standard map T ,a ∈ L1 is defined by (see Fig. 9): a   T ,a (x ) = x + − τ {{x}} , (6) 2π which corresponds to the standard map but using the τ wave function instead of the sin(2π x ) function. The upper and lower maps for this family are very easy to compute. Moreover, T ,a is non-increasing for a > π2 and hence, in this case, the upper and lower maps do not coincide and have constant sections. To compute the rotation intervals and Arnol’d Tongues of T ,a we proceed as for the Standard Map by using Theorem 2.5 and Algorithm 3. In Fig. 10 we show some graphs of the rotation interval and Arnol’d tongues for the piecewise-linear standard map. The graphs of the rotation intervals are plotted for three different values of as a function of the parameter a. Definition 4.8 (The Discontinuous Standard Map). D ,a ∈ L1 is defined as (see Fig. 11): a D ,a (x ) := x + + {{x}}. (7) 2π 15 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 Fig. 10. Graphs of the rotation interval and Arnol’d tongues for the piecewise-linear standard map T ,a . The graphs of the rotation intervals are plotted as a function of the parameter a. Fig. 11. The discontinuous standard map with a = 2π and = 0 with its lower map in blue and its upper map in red. The map D ,a , being discontinuous, belongs to the so called class of old heavy maps [17] (the old part of the name stands for degree one lifting — that is, D ,a ∈ L1 ). A map F ∈ L1 is called heavy if for any x ∈ R, lim F (y ) ≤ F (x ) ≤ lim− F (y ) yx + yx (in other words, the map “falls down” at all discontinuities). 16 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 Fig. 12. Graphs of the rotation interval and Arnol’d tongues for the discontinuous standard map D ,a . The graphs of the rotation intervals are plotted as a function of the parameter a. Observe that for the class of old heavy maps the upper and lower maps in the sense of Definition 2.4 are well defined and continuous. Moreover, the whole family of water functions (c.f. [11]) is well defined and continuous. So, the rotation interval of the old heavy maps is well defined [17, Theorem A] and, moreover, Theorem 2.5 together with Algorithm 3 work for this class. Hence, to compute the rotation intervals and Arnol’d Tongues of D ,a we proceed again as for the Standard Map. As for the piecewise-linear standard maps the upper and lower maps are very easy to compute, and have constant sections for a = 0. In Fig. 12 we show some graphs of the rotation interval and Arnol’d tongues for the discontinuous standard map. The graphs of the rotation intervals are plotted for three different values of as a function of the parameter a. The times taken for all the computation related with the rotation intervals and the Arnol’d Tongues for each of the families studied using Algorithms 1, 2 and 3 can be found in Table 1. 5. Conclusions The proposed algorithm clearly outperforms all the other tested algorithms, both in precision and speed even though the “exact” (and quick) part of the algorithm does not work for all the non-decreasing liftings in L1 which have a con- stant section (and hence the rotation number of these “bad” cases has to be computed with the much more ineffi- cient classical algorithm). For all natural examples for which it has been tested, the computational speed and precision were unparalleled. Moreover, the set of functions becomes very general when one considers the fact that the upper and lower functions inherently have constant sections for any F that is not strictly increasing. Hence, the algorithm becomes a crucial tool to compute rotation intervals for general functions in L1 and, therefore, to find the set of periods of such maps [11]. Moreover, a deeper study has been done on the dependence of the rotation number on the parameters. Our preliminary results have found that for irrational rotation numbers, the dependence of the parameters around them is extremely sen- sitive, with continuity module being at least 1025 . This agrees with Theorem 4.3, which says that non-persistent functions have measure zero. 17 L. Alsedà and S. Borrós-Cullell Commun Nonlinear Sci Numer Simulat 102 (2021) 105915 Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. CRediT authorship contribution statement Lluís Alsedà: Conceptualization, Methodology, Investigation, Software, Writing - review & editing, Supervision, Funding acquisition, Resources. Salvador Borrós-Cullell: Conceptualization, Methodology, Investigation, Software, Writing - review & editing. References [1] Misiurewicz M. Periodic points of maps of degree one of a circle. Ergodic Theory Dyn Syst 1982;2:221–7. [2] Alsedà L, Llibre J, Mañosas F, Misiurewicz M. Lower bounds of the topological entropy for continuous maps of the circle of degree one. Nonlinearity 1988;1:463–79. [3] Poincaré H. Sur les courbes définies par les équations différentielles (III). Journal de mathématiques pures et appliquées 4e série 1885;1:167–244. [4] Marangio L, Sedro J, Galatolo S, et al. Arnold maps with noise: Differentiability and non-monotonicity of the rotation number. J Stat Phys 2020(179):15941624. [5] Janson S, Öberg A. A piecewise contractive dynamical system and election methods. Bulletin de la SociȨtȨ MathȨmatique de France 2019;147(3):395–411. [6] Mora X, Oliver M. Eleccions mitjanȺant el vot d’aprovació. el mètode de phragmén i algunes variants. Butlletȡ de la Societat Catalana de Matemàtiques 2015;30(1):57–101. [7] Herman M. Sur la conjugaison différentiable des difféomorphismes du cercle à des rotations. Publications Mathématiques de l’IHÉS 1979;49:5–233. http://www.numdam.org/item/PMIHES_1979__49__5_0 [8] Pavani R. A numerical approximation of the rotation number. Appl Math Comput 1995(73):191–201. [9] Seara TM, Villanueva J. On the numerical computation of diophantine rotation numbers of analytic circle maps. Physica D 2006(217):107–20. [10] Veldhuizen MV. On the numerical approximation of the rotation number. J Comput Appl Math 1988(21):203–12. [11] Alsedà L, Llibre J, Misiurewicz M. Combinatorial dynamics and entropy in dimension one. World Scientific; 20 0 0. [12] Ito R. Rotation sets are closed. Math Proc CambPhilos Soc 1981;89(1):107111. doi:10.1017/S03050 0410 0 057984. [13] Sánchez J, Net M, Simó C. Computation of invariant tori by Newton–Krylov methods in large-scale dissipative systems. Physica D 2009(239):123–33. [14] Misiurewicz M. Persistent rotation intervals for old maps. Banach Center Publications; 1989. [15] Broer H, Simó C. Hill’s equation with quasi-periodic forcing. Boletim da Sociedade Brasileira de Matemática 1998;29(2):253–93. [16] Boyland PL. Bifurcation of circle maps: Arnol’d tongues, bistability and rotation intervals. Commun Math Phys 1986(106):353–81. [17] Misiurewicz M. Rotation intervals for a class of maps of the real line into itself. Ergodic Theory Dyn Syst 1986;6(1):117132. doi:10.1017/ S014338570 0 0 03321. 18