Rigorous enclosures of rotation numbers by interval methods

We apply set-valued numerical methods to compute an accurate enclosure of the rotation number. The described algorithm is supplemented with a method of proving the existence of periodic points, which is used to check the rationality of the rotation number. A few numerical experiments are presented to show that the implementation of interval methods produces a good enclosure of the rotation number of a circle map.


Introduction
Given an orientation-preserving homeomorphism f : S 1 → S 1 of the circle, the rotation number ρ is a topological invariant of f , defined as ρ = ρ(f ) = ρ(F )(mod 1), where Here the function F -a lift of f -is a homeomorphism of the real line, satisfying f • π(x) = π • F (x), where π : R → S 1 is the natural projection given by π(x) = x(mod 1). By a classical result of Poincaré [1] the limit (1)is welldefined up to an integer, and independent of the choice of x. The circle map f has a periodic point if and only if its rotation number is rational: if ρ = p/q, where p and q are co-prime integers, then f has a periodic point of prime period q. If, on the other hand, the rotation number is irrational, then there are two possibilities: either all orbits are dense in S 1 and f is topologically conjugate to a rigid rotation or there exists a Cantor set C ⊂ S 1 which is invariant under f and both forward and backward orbits of all points converge to C, and then f is semi-conjugate to a rigid rotation. Apart from a rich flora of theorems concerning rotation numbers for generic classes of circle maps, there exists several results concerning the estimation of the rotation number ρ of a given circle map f , see e.g. [2], [3], [4]. All such existing methods require the computation of very long trajectories of f ; this is apparent from the basic definition (1) of the rotation number. When the dynamics of f displays sensitive dependence (due to regions of expansion), it is well-known that numerical computations of even moderately long trajectories can be very inaccurate. As such it may be warranted to question the reliability of numerically computed estimates of a rotation number. We address this issue by presenting a method that produces explicit, computable bounds on ρ for a given f . As a result, we can e.g. guarantee the correctness of the n first decimals of our computed estimate of ρ.
In order to guarantee a mathematical rigour in our numerical work, we base all computations on set-valued mathematics. This approach is based on interval arithmetic [5,6,7], and is very well suited for controlling the propagation of numerical errors throughout long iterative processes. This paper has the following structure: we begin by recalling some algorithms that -theoretically, and under the right circumstances -can be used to produce enclosures of a rotation number. Next, we make the most efficient of these algorithms practically applicable by deriving explicit bounds that are needed in the computations. The end-product is a tight enclosure of the sought rotation number ρ. By adding a sub-algorithm that proves the existence (and uniqueness) of moderately long periodic orbits, we can also handle the case of rational rotation numbers. Finally, we present some numerical results from our implementation of the described algorithm.

Algorithms for computing the rotation number
The most simple-minded algorithm for computing the rotation number is a direct application of the basic definition (1).
This method produces explicit bounds of the rotation number; its main drawback is that the convergence is only linear. So, given an initial point x 0 ∈ S 1 , and an (interval) enclosure x N of F N (x 0 ), we can define ρ N = x N −x0 N , which produces the bounds: Here, we denote the endpoints of an interval according to In [8] the rotation number ρ is constructed via the continued fraction expansion obtained from the topological behaviour of the circle map: The continued fraction coefficients a i depend on increasingly high iterates of the circle map, and -from a numerical point of view -they are not easy to compute. Nevertheless, using this approach, it is possible to construct a sequence of intervals of rapidly decreasing lengths, all containing the rotation number ρ. In general the convergence is quadratic, and when the rotation number satisfies a Diophantine 1 condition, the convergence can be made cubic, see [8]. Note that the set of rotation numbers satisfying a Diophantine condition has full measure in [0, 1]. On the other hand, the complement lies dense in [0,1]. In what follows, we will not assume any Diophantine properties of the rotation number. Let f be an orientation preserving homeomorphism on the circle S 1 = R/Z. Set b 0 = f (0) ∈ (0, 1), and define I 0 to be the smallest of the two intervals [b 0 , 1) and (0, b 0 ]. If b 0 = 1/2 there is a tie, and we set I 0 = (0, 1/2]. Next, let a 1 be the first return time of f (b 0 ) to the interval I 0 , i.e., the smallest positive integer such that f a1 (b 0 ) ∈ I 0 . If there is no such return, we set a 1 = ∞, and terminate the process. Otherwise we define b 1 = f a1 (0), and define I 1 to be the interval [b 1 , 1) or (0, b 1 ] that is disjoint from I 0 . This ends the initial stage of the general procedure for computing the continued fraction of ρ. Continuing with the general inductive stage, we let {I i } be a sequence of half-open intervals of the form [b i , 1) or (0, b i ], and let {φ i }, i ≥ 1 denote the sequence of first return maps of f from I i to (I i I i−1 ). We define a i+1 to be the smallest positive integer such that φ ai+1+1 i As before, if the return map is not well-defined, we set a i+1 = ∞ and terminate the process. In parallell with this process we recursively define positive integers p i and q i as follows 2 : These integer sequences are used in the following theorem which provides explicit bounds on the rotation number ρ.
Theorem 2.2. [8] Let N i be the number of iterates needed to compute q i . For the sequences defined above the following holds: (a) If ρ is irrational, then pi qi converges to ρ as i → ∞. If ρ is rational, then the process terminates (a i+1 = ∞) and the last estimate pi qi = ρ.
We will use a version of case (b). First, let us begin by stating an obvious fact about continued fractions of rotation numbers. Remark. If we knew that the rotation number was irrational, we would get a tighter bound: By using the procedure of Theorem 2.2, it is clear that we need a means to compute very high iterates of the circle map f . This is needed in order to compute the return maps φ i , which in turn are used to determine the return times a i appearing in the recursions (3), and in the continued fraction of ρ. But computing high iterates of a non-linear dynamical system is not straightforward. We will address this issue in the following section.

The interval Newton method and shooting techniques
As already mentioned, the described algorithms all require long pieces of trajectories in order to produce good approximations (and bounds) on the rotation number. Naturally, this is problematic as it is not always clear how to control the propagation of numerical errors in such computations. Therefore we need a reliable and accurate method to compute long trajectories. A common solution to this problem is to use multi-precision numerical software. This will allow for longer accurate trajectories, but the penalty is high: a decrease in speed by a factor 200-1000 is not uncommon. We will illustrate this in Section 5.
Sticking to standard floating point computations, we will use a version of the shooting method described in [9] but modified to suit our particular set-up. This method is based on a combination of the interval Newton method and multiple shooting.
The interval Newton method is a constructive implementation of the bounds required in Kantorovic's theorem in order to guarantee the convergence of Newton's method. As such, it can be used to prove the existence (or non-existence) of zeros of general n-dimensional maps. Let g : R n → R n be a continuously differentiable function, and suppose that we have an interval extension of its derivative, i.e., given an n-dimensional interval variable z ⊂ R n , that is z = (x 1 , · · · , x n ), we can compute the interval image Dg(z) ⊇ {Dg(z) : z ∈ z}.
To verify the existence of zeros of g, we introduce the interval Newton operator where [Dg(z)] is an interval matrix containing all Jacobian matrices of g of the form Dg(z) for z ∈ z. Herež is an arbitrary point from the interval vector z usually chosen to be the middle point of z. The operator (4) possesses the following key properties: 1) If N g (z) ⊂ z then there exists exactly one point z ∈ z such that g(z ) = 0.
2) If N g (z) ∩ z = ∅ then there are no zeros of g in z.
Let z 0 = z be the initial enclosure of a possible zero of g, and define the sequence of intervals z k+1 = N g (z k ) ∩ z k , k = 0, 1, 2, . . . . If a true zero z is contained in z 0 , and if the interval Newton operator is well-defined on this domain, then the operator remains well-defined for all iterations, we have z ∈ z k , and the intervals z k form a nested sequence converging to the zero of g.
By applying the interval Newton operator to F -a high-dimensional variant of f based on the shooting technique -we can efficiently recast the problem of tracking long orbits of f into one about finding zeros of F . To be more precise: let f : S 1 → S 1 be a circle map as defined in Section 2. In order to compute the n iterates x 1 , x 2 , . . . , x n , where x k = f k (x 0 ), we apply the interval Newton operator to the n-dimensional map F (·; x 0 ) : (S 1 ) n → (S 1 ) n defined as follows: where z = (x 1 , · · · , x n ) and x 0 is the initial point of the orbit. Note that the Jacobian matrix of F (·; x 0 ) at z = (x 1 , · · · , x n ) has the following -very sparse -form: To get a rigorous enclosure of a finite length trajectory starting at x 0 , we begin by computing an approximate trajectory z = (x 1 , · · · ,x n ). We then consider a neighbourhood of this piece of trajectory by widening each component into an interval, producing the box z = (x 1 , · · · , x n ). Next, if we can show that the interval Newton operator maps this interval vector into itself N F (z; x 0 ) ⊂ z, then we have a true trajectory enclosed in z. By reapplying the interval Newton operator a few times, we can tighten the enclosure of the trajectory as much as we want. This gives us very high-quality trajectories, starting from only coarse approximate trajectories.
The computation of the interval matrix inverse can cause difficulties for large n. Fortunately, the explicit inverse is not needed; we can simply solve for the correction term h = DF (z; x 0 ) −1 F (ž; x 0 ). This equation can be recast as We define all h 2 , · · · , h n recursively in the following way On completion of these computation, we can evaluate the interval Newton operator as N F (z; x 0 ) =ž − h. This concludes the description of how to compute long, yet highly accurate (in fact rigorous), iterates of the circle map f .

Verification of the rationality of the rotation number
Recall from part (a) of Theorem 2.2 that if ρ is rational, then the process terminates (a i+1 = ∞) and the last estimate pi qi = ρ. In fact, this means that the circle map f has a stable period-q i point, and all forward trajectories tend to the corresponding periodic orbit 3 . With this in mind, we always search for a period-q i point of the circle map before attempting to determine the next return time a i+1 .
Again we use the interval Newton method to verify the existence of a periodic point. The line of argument is very similar to that of the previous section. We first try to locate a good candidate by iterating a random point many times, producing the finite orbit x 0 , x 1 , . . . , x N . Here N is taken much larger than the sought period q i . If there is a period orbit, these iterates should accumulate on it, which means that x N +1 , x N +2 , . . . , x N +qi should form a good approximation of the periodic orbit. Widening this finite set of points into a q i -dimensional box, and verifying that this box is mapped into itself under the interval Newton map, established the existence of a period-q i point. This -in turn -means that a i+1 = ∞, and that the rotation number is rational: ρ = pi qi .

Numerical experiments
In this section we present numerical results from our implementations of both described algorithms: the linear one (based on Theorem 2.1) and the quadratic one (based on Theorem 2.2). We base our code on the CAPD interval library [10]. For all computations performed, we list the timings in seconds using a single thread on an Intel(R) Core(TM) i7 CPU 920 @ 2.67GHz.

The Arnold family
We consider the two-parameter Arnold family of circle maps given by The map f α, (x) with 2π| | ≤ 1 is an orientation-preserving analytic circle diffeomorphism for every pair of parameters (α, ). The set of parameters (α, ) corresponding to a fixed irrational ρ forms a continuous curve known to be analytic when ρ is a Diophantine number; for the rational ρ the set has an nonempty interior, and formes the well-known the Arnold tongue of the rotation number ρ [3,11]. Table 1 displays enclosures of the rotation number when computed according to Theorem 2.1 using set-valued computations for different pairs of parameters of the Arnold map. We list the final enclosure of the rotation number, its radius, and the number of iterations of the map f . In order to save space, and to highlight the significant digits, we use a compact notation for intervals: for example, 0.   Enclosures of the rotation number obtained using our algorithm (Lemma 2.3) are given in Table 3. Here the convergence is at least quadratic, and the use of interval arithmetic allows us to get a much better enclosure already for a small number of iterations.  Using the interval Newton method as described in Section 3 we can often verify if the rotation number is rational. For the cases shown in Table 3 no periodic points were detected, and so the results encloses a (possibly) irrational rotation number. Similar computations were performed for other pairs of parameters (α, ) given in Table 5. These parameters produce rational rotation numbers, and thus we can provide an exact representation of the form ρ = p q , where p and q correspond to the integers p i , q i in Theorem 2.2. Note that here q is the period of the periodic point.
We end this example by computing rotation numbers for many different parameters. This produces the well-known devil's staircase graph, illustrating the frequency locking phenomena. To be concrete, we will take ∈ {0.01, 0.1, 0.159},   and for each value of we compute the enclosure of the rotation number ρ for 1000 different values of α. The results are illustrated in Figure 1.

The delayed logistic map
Similar computations were carried out for the delayed logistic map, defined by The approximation of the rotation number for this map has been considered in [12]. It is believed that for 2 < λ ≤ 2.16 the map φ possesses a smooth invariant curve. In order to prove the existence of an invariant curve for the delayed logistic map, one may adopt the method described in [13], where a topological approach to prove the existence of a normally hyperbolic manifold for maps is introduced. As this is not the focus of our work, we will simply assume that the invariant curve exists for the parameters we are considering.  Table 6: Computational results for the delayed logistic map, using the algorithm with quadratic convergence. The short computational times are due to the simple form of the map (no trigonometric functions).
To compute an approximate enclosure of the rotation number, we first iterate the delayed logistic map several times in order to avoid transient behaviour. This leaves us with a point very close to the invariant curve (see Figure 3), and we will use this point as our initial point. To compute intervals |I i | from Theorem 2.2 at each iteration we consider an angle defined using arctangent function of two arguments. We present our results for different values of the parameter λ in Table 6. For comparison, the second column ρ V was computed using the method presented in [12].

Discussion
We have presented a method for computing an enclosure of the rotation number of a given circle map. This method is based upon set-valued computations combined with explicit bounds that can be obtained whilst computing the continued fraction of the sought rotation number. We have demonstrated our method by applying it to the classical two-parameter Arnold family of circle maps, and to the delayed logistic equation. In the latter case, we compare our results with those appearing in [12].
Of course, our main goal so far has been to obtain mathematical rigour in the numerical computations necessary. Nevertheless, by utilizing the method of multiple shooting together with the interval Newton method, we can perform all computation in double precision. This makes our method quite fast in comparison to approximative methods that utilize multi-precision.