Set-oriented numerical computation of rotation sets

We establish a set-oriented algorithm for the numerical approximation of the rotation set of homeomorphisms of the two-torus homotopic to the identity. A theoretical background is given by the concept of {\epsilon}-rotation sets. These are obtained by replacing orbits with {\epsilon}-pseudo-orbits in the definition of the Misiurewicz-Ziemian rotation set and are shown to converge to the latter as {\epsilon} decreases to zero. Based on this result, we prove the convergence of the numerical approximations as precision and iteration time tend to infinity. Further, we provide analytic error estimates for the algorithm under an additional boundedness assumption, which is known to hold in many relevant cases and in particular for non-empty interior rotation sets.


Introduction
Rotation theory for orientation-preserving homeomorphisms on the circle was established by H. Poincaré in 1885 [1], who showed that the long-term behaviour of such maps can be classified by the rotation number. This topological invariant provides a dichotomy for the dynamics depending on whether it is a rational or an irrational number, corresponding to periodic or quasiperiodic motion, respectively.
For homeomorphisms of higher dimensional tori, it is well-known that a unique rotation vector does not need to exist anymore. In [2], Misiurewicz and Ziemian therefore introduce the rotation set of a torus homeomorphism, which collects all possible asymptotic rotation vectors and carries strong information about the system's asymptotic behaviour. For homeomorphisms on the two-torus, this compact non-empty set is always convex [2]. If it has non-empty interior, then the system has positive topological entropy [3] and all rational points in the interior of the rotation set are realised by periodic orbits [4]. Apart from these nowadays classical results, considerable progress has been made in the last years on the rotation theory of surface homeomorphisms (e.g. [5]- [10]), including partial results on the well-known Franks-Misiurewicz Conjecture [11]- [17] that excludes the existence of certain line segments as rotation sets.
Concerning the shape of the rotation set, it is further known that genericallyin the C 0 -topology -it is a polygon with rational vertices [18]. For generic areapreserving torus homeomorphisms, this polygon is non-degenerate (has non-empty interior) [19]. Moreover, all rational polygons can be realised as the rotation set of a torus homeomorphism [20]. For a particular parameter family of such maps, derived, through some elaborate inverse limit construction, from symbolic systems related to beta expansions, bifurcations of the rotation set involving new types of convex sets have been described in [21] (see also [22]). Apart from that, however, there is still little knowledge concerning the question which compact convex subsets of the plane may appear. Likewise, there is still not much insight into the behaviour of the rotation set in 'natural' parameter families of torus diffeomorphisms, such as those studied by theoretical physicists in [23].
One serious obstruction for further progress in this direction is the fact that, except for some particular cases, it is usually not possible to analytically determine the rotation set. Moreover, due to the inherent nonlinearity of the problem, the numerical computation has proven to be difficult as well. In [19], pointwise approaches to the numerical approximation of the rotation set are discussed. Based on the detection of periodic orbits of the system on the one hand, and of a discretisation of the system by the projection to a finite lattice on the other hand, the complexity of this numerical task as well as the restrictions of the proposed approaches are discussed. It turns out that an accurate computation is already difficult in the case where the rotation set is still a rational polygon, but its vertices correspond to periodic orbits of larger periods. When the rotation set is not a rational polygon at all, it has to be expected that the situation is even worse.
Thus, our aim here is to establish a more reliable algorithm for the numerical computation of rotation sets based on set-oriented methods. On the theoretical level, this corresponds to considering ε-pseudo-orbits reflecting the inaccuracies of numerical calculations. This approach leads to the definition of an ε-rotation set in Section 3. In the two-dimensional case, we prove that as the perturbation size ε > 0 decreases to zero the respective ε-rotation sets converge to the original rotation set (Theorem 7).
In Section 4, we formulate our set-oriented algorithm for the numerical approximation of the rotation set of a homeomorphism f : T 2 −→ T 2 homotopic to the identity, which is defined via a lift F : R 2 −→ R 2 of f . Based on the fact that the rotation set is approximated with arbitrary precision in the Hausdorff metric by the normalised iterates 1 n F n ([0, 1] 2 ) of the unit square (Corollary 10), our algorithm aims to provide a good visualisation of the latter sets for large n ∈ N. This is backed up by rigorous error estimates and a result on the convergence of the approximations to the true rotation set as the precision and iteration time tend to infinity (Theorem 11). If the considered system has a shadowing property, this leads to a further considerable improvement of the results on the convergence rate (Theorem 15). By a systematic overestimation of the system's orbits, we further ensure that the numerical algorithm always yields a super-set of the actual rotation set. This is particularly important since, as mentioned, previous approaches tend to 'miss out' the vertices of polygonal rotation sets if these are realised by periodic orbits of large period [19].
Numerical results for some parameter families inspired by [2,19,23] are then presented in Section 5. It turns out that in most cases the algorithm performs much better than predicted by the numerical error estimates. The reason for this is possibly the fact that the considered systems possess a shadowing property. As mentioned, this leads to a faster convergence of the approximations.

Preliminaries
By Conv(A) we denote the closed convex hull of a set A ⊆ R m . The open εneighbourhood of points or sets in R m will be denoted by B ε (x), respectively B ε (A). The closure of a set A is denoted by A. By d H (A, B) we denote the Hausdorff distance between two subsets A, B of a metric space and also write A n → H A if a sequence (A n ) n∈N converges to A in the Hausdorff topology as n → ∞.
For m ∈ N, we let T m = R m /Z m denote the m-dimensional torus. The set of all lifts F : R m → R m of continuous maps f : T m −→ T m homotopic to the identity is denoted by C m . The subset of C m consisting of lifts of torus homeomorphisms is denoted by H m . Note that C m consists of all continuous functions F : R m −→ R m that satisfy for all x ∈ R m , k ∈ Z m , and H m is the set of all such F which are in addition homeomorphisms of the plane. In their seminal paper [2] Misiurewicz and Ziemian introduced the rotation set as Thus, (F ) can be viewed as the collection of all possible rotation vectors of the system. By writing for k ∈ N, we alternatively obtain the rotation set as the upper Hausdorff limit of the sequence (K n (F )) n∈N , that is, Moreover, in the two-dimensional case, the approximating sets K n (F ) converge to the rotation set in the Hausdorff distance. ). Let F ∈ C m . Then (F ) ⊆ Conv(K n (F )) for all n ∈ N. A generalisation of Lemma 3 to ε-rotation sets is given by Lemma 8 below. An important technical estimate in [2] is the following.
). Applied to G = F n and taking into account that this immediately implies the following consequence. Corollary 5. If F ∈ H 2 and n ∈ N, then Conv(K n (F )) ⊆ B 3 √ 2 n (K n (F )) .

The ε-rotation set
In order to model numerical approximations of rotation sets, we take into account the fact that inaccuracies are inherent to computer simulations. While the definition of the rotation set is based on proper orbits of the underlying torus map, we introduce an alternative by allowing perturbations of the system's orbits. This corresponds to the well-known concept of ε-pseudo-orbits and leads to the notion of ε-rotation sets. Since the start of this project, these have also been described independently by Guiheneuf and Koropecki in [24].
Let (X, d) be a metric space and F : X −→ X a self-map. For n ∈ N and ε ≥ 0, an (n + 1)-tuple (ξ j ) n j=0 of points ξ j ∈ X is called an for all j ∈ {0, . . . , n − 1}. In the same way infinite sequences (ξ j ) ∞ j=0 or (ξ j ) ∞ j=−∞ satisfying (3.1) for all j ∈ N 0 or j ∈ Z, respectively, are called (infinite) ε-pseudoorbits. Obviously, every proper orbit is an ε-pseudo-orbit for all ε > 0. Note that we slightly deviate from the standard definition by not requiring strictness of the inequality in (3.1). This will be very convenient later on for technical reasons, but has no significance on a conceptual level.
For F ∈ C m and ε ≥ 0, we now define the ε-rotation set to be the set of all accumulation points of sequences of the form for k ∈ N and can alternatively define the ε-rotation set by Note that since F commutes with integer translations, one could also replace ξ 0 ∈ [0, 1] m by ξ 0 ∈ R m in the definition of K ε k (F ). For every k ∈ N, the inclusion K k (F ) ⊆ K ε k (F ) is apparent, so that due to (2.3) and (3.2), we obtain for all ε ≥ 0. We omit the elementary proof of the following lemma.
Then (i) for each k ∈ N the set K ε k (F ) is non-empty and compact, and thus the same holds for ε (F ); (ii) (F ) ⊆ B M (0) and ε (F ) ⊆ B M +ε (0). We now aim to show convergence of the ε-rotation sets as ε decreases to zero. This presents the main result of this section. An alternative proof can be found in [24]. However, we include the proof both for the convenience of the reader and because the employed arguments and estimates will become crucial again in Section 4.
Since for (F ) ⊆ ε (F ) for all ε > 0 by (3.3), it remains to show that for every δ > 0 there is an ε > 0 such that We split the proof into Lemmas 8 and 9, beginning with a formulation of Lemma 3 in terms of ε-rotation sets.
since the tuples (ξ in , ξ in+1 , . . . , ξ (i+1)n ) are ε-pseudo-orbits of F of length n + 1 and Thus, the norm of the first summand in (3.6) is bounded by r k k (M + ε). For the second summand, we have For a fixed value of n, the inclusion (3.7) is valid for all k ≥ n. Therefore, for each n ∈ N, the representation (3.2) of the ε-rotation set yields due to the convergence r k k → 0 as k → ∞.
As δ > 0 was chosen arbitrarily and n ∈ N is fixed, this proves the convergence Proof of Theorem 7. Recall that in order to prove Theorem 7, we have to show the inclusion (3.4). Let δ > 0. By Lemma 1, the rotation set (F ) can be approximated up to an arbitrary precision with respect to the Hausdorff distance by the sets K n (F ). Therefore, we choose n 0 ∈ N sufficiently large, so that Furthermore, by Lemma 9 there exists an ε 0 > 0 sufficiently small, such that for all 0 < ε ≤ ε 0 . Applying Lemma 8, together with (3.8) and (3.9) we obtain is convex itself, since the rotation set is convex by Theorem 2. As δ > 0 was arbitrary, this shows the asserted convergence ε (F ) → H (F ) as ε → 0.

Set-oriented computation of rotation sets
4.1. Why set-oriented numerics -shortcomings of the direct approach.
Before we turn to our algorithm for the set-oriented computation of rotation sets, we briefly want to comment on the problems that arise with a more conventional approach.
The easiest way to compute the rotation set of a torus homeomorphism numerically would be to fix a standard grid Γ ⊆ [0, 1] 2 of N × N points, to compute the normalised displacement vectors 1 for all x ∈ Γ and some large n ∈ N and to plot the collection of these vectors in order to obtain an approximation of the rotation set. However, one may now consider the following situation, which actually turns out to be generic (with respect to C 0 -topology, see [19]): Suppose that an area-preserving torus homeomorphism f has a rotation set with non-empty interior, and that the Lebesgue measure Leb T 2 on T 2 is ergodic with respect to f . Then there exists a well-defined rotation vector x is interpreted as a function on the torus, and we have for Lebesgue-a.e. x ∈ T 2 . Hence, if n above is chosen too large, the rotation vectors coming from starting points in the grid Γ will almost surely be arbitrarily close to f (Leb T 2 ). In this case, the numerical approximation of the rotation set will show the singleton f (Leb T 2 ), or something very close to it, whereas the true rotation set is much bigger due to its non-empty interior. Figure 4.1 illustrates this numerical phenomenon for a standard example f 1,1 (see Section 5) of a torus diffeomorphism given by the lift which exhibits f 1,1 (Leb T 2 ) = (0, 0) and is known to have the rotation set (F 1,1 ) = [−1, 1] 2 (see Lemma 22); as the number of iterations increases the majority of the approximate rotation vectors based on the iteration of grid points tend towards the origin and only the fixed points at the vertices are detected. It is effects like this which prevent a direct approach from producing reasonable results [19], and frequently lead to an underestimation of the actual rotation set. In order to counter this, one would have to choose the grid constant N exponentially large with respect to n (of magnitude N ∼ L n , where L is a Lipschitz constant of f ), which is not feasible for efficient computations. Moreover, the above discussion shows that in the direct approach there is a critical dependency between the parameters N (grid size) and n (number of iterates), which cannot be chosen independent of each other. In contrast to this, it is surely desirable to have an algorithm whose precision increases in both parameters independently, so that computation capacities can be pushed to their limit without having to worry about the precise relation of the parameters. This is exactly what the set-oriented method will allow us to do.

4.2.
The set-oriented approach -description of the algorithm. We start with a basic observation. By Lemma 1, the sets K n (F ) converge to the rotation set (F ), and moreover we have the elementary estimate Thus, we obtain This allows to focus on the sets 1 n F n ([0, 1] 2 ) in order to compute the rotation sets, which is quite convenient from a practical viewpoint. Hence, we aim for an accurate numerical approximation of the sets 1 n F n ([0, 1] 2 ) for some large n ∈ N. Inspired by the software package GAIO (Global Analysis of Invariant Objects), we apply the concept of box coverings to formulate our algorithm. This library, developed by Dellnitz, Froyland and Junge [25], provides numerical methods for the analysis of dynamical systems by set-oriented techniques. As mentioned above, on the theoretical level this corresponds to considering ε-pseudo-orbits.  Note that B 0 can be considered as a covering of T 2 , whose lift to R 2 is then given by The elements of the collections B 0 and B will be referred to as boxes. Further, given η > 0, for each box B ∈ B, let Γ B ⊆ B be a finite set of points which is η-dense in B and chosen such that Γ B+t = Γ B + t for every integer vector t ∈ Z 2 . Finally, fix R > 0. Then, for F ∈ H 2 , we now generate a sequence of sets (Q * n ) n∈N (approximations of the rotation set) according to the following algorithm. Normalisation: The starting point of our algorithm is the unit square [0, 1] 2 . Instead of iterating single test points, we consider collections of boxes B ∈ B and always collect those which are hit by an image of one of the previous boxes. By choosing the parameter R for the numerical calculation of the box images sufficiently large, we can ensure that no boxes are missed out and the images are always overestimated (see Lemma 12). This leads to a systematic overestimation of the sets 1 n F n [0, 1] 2 , which will eventually lead to the fact that the algorithm always yields a superset of the actual rotation set (up to a small shift of order 2 √ 2 /n, see Lemma 20). Hence, the effect of underestimation described for the direct approach in the previous subsection can be excluded (see Theorem 11 below).
Note also that for the implementation of the algorithm, the box images only have to be computed once of each box B ∈ B 0 at the beginning. This immediately provides the box images for the respective integer translates as well, and the information can then be used throughout the whole iterative procedure.

4.3.
Convergence of the approximations -qualitative results and error bounds. Relations between the normalised approximations Q * n and the rotation set (F ) are established by the following two results. Theorem 11 is qualitative in nature and ensures convergence, whereas Theorem 14 uses an additional boundedness assumption to provide error estimates. Theorem 11. Suppose that F ∈ H 2 is Lipschitz continuous with Lipschitz constant L > 1 and for each ε > 0 the constants η, R > 0 are chosen such that Lη ≤ R ≤ ε. Then in the sense that for all δ > 0 there exist ε 0 > 0 and n 0 ∈ N such that if ε in (4.2) satisfies ε ∈ (0, ε 0 ] and n ≥ n 0 , then d H (Q * n , (F )) < δ. Remark 4.1. For theoretical purposes, one may also want to ignore the fact that the images of the boxes B ∈ B can only be approximated via test points, and define alternative sequences Q n n∈N , Q * n n∈N by using the precise images viã Q 0 = Q 0 , This corresponds to setting the parameters η and R to zero. Then the above convergence result is still valid for the new sequence, that is, limn→∞ ε→0 d H (Q * n , (F )) = 0. Moreover, the assumption of Lipschitz continuity is not needed in this case, but it would still be required in the respective analogue of the error estimates given in Theorem 14 below. For the proof of Theorem 11, we will need the following lemma. Lemma 12. Suppose that F ∈ H 2 is Lipschitz continuous with Lipschitz constant L > 1 and for each ε > 0 the constants η, R > 0 are chosen such that Lη ≤ R ≤ ε.
Proof. First, suppose that v = 1 n F n (x 0 ) ∈ 1 n F n [0, 1] 2 and let x j = F j (x 0 ). Choose boxes B j ∈ B such that x j ∈ B j . We claim that B j ⊆ Q j for j = 0, . . . , n, so that finally x n = nv ∈ B n ⊆ Q n and hence v ∈ Q * n . For j = 0 the claim is obvious. Therefore, suppose that B j ⊆ Q j . Then there exists a test point y ∈ Γ B j η-close to x j , so that d(F (y), x j+1 ) ≤ Lη ≤ R. Since x j+1 ∈ B j+1 , this implies B j+1 ∈ I(B) and thus B j+1 ∈ Q j+1 . This shows the first inclusion 1 n F n [0, 1] 2 ⊆ Q * n . In order to show the second inclusion, suppose that v ∈ Q * n . Then, by definition of Q n , there exists a sequence of boxes B 0 , . . . , B n such that nv ∈ B n and B j+1 ∈ I(B j ) for all j = 0, . . . , n − 1. By definition of the box images, there exists a sequence of test points ξ j ∈ Γ B j , for all j = 0, . . . , n − 2, and this remains true for j = n − 1 if we let ξ n = nv ∈ B n . This means that (ξ j ) n j=0 is a 2ε-pseudo-orbit and thus yields the required second inclusion.
Corollary 13. In the situation of Lemma 12, we have Proof of Theorem 11. First, by Theorem 7 and the convexity of the rotation set, there exists ε 0 > 0 such that Further, by the definition of 2ε 0 (F ), there exists n 0 ∈ N such that for all n ≥ n 0 we have K 2ε 0 n (F ) ⊆ Bδ /3 2ε 0 (F ) and hence Conv 2ε 0 (F ) .
As for all n ∈ N we have (F ) ⊆ Conv (K n (F )) ⊆ Conv K 2ε 0 n (F ) (see Lemma 3), we obtain that At the same time, Lemma 1 implies that we can choose n 0 such that d H (K n (F ), (F )) < 2δ 3 (4.9) for all n ≥ n 0 . If n 0 is chosen such that √ 2 /n 0 < δ /3 and n ≥ n 0 , then using (4.9) together with the first inclusion in Corollary 13, we obtain

Conversely, (4.8) together with the second inclusion in Corollary 13 yield
Together with the fact that K 2ε n (F ) ⊆ K 2ε 0 n (F ) for all ε < ε 0 , this means that d H (Q * n , (F )) < δ whenever n ≥ n 0 and ε ∈ (0, ε 0 ], as required. In general, it is not possible to give quantitative error estimates for the numerical computation of rotation sets. The reason is that there are no general apriori bounds for the convergence of the sets K n (F ) to (F ). However, it turns out that in many situations, and in particular whenever the rotation set has non-empty interior, there exists a positive constant c > 0 such that This fact has been proven for diffeomorphisms in [26] and the result was later generalised to homeomorphisms in [12]. This is also referred to as bounded deviation property. In our context, (BD) together with the existence of a Lipschitz constant allows to provide the following quantitative estimates. Theorem 14. Suppose F ∈ H 2 is Lipschitz continuous with Lipschitz constant L > 1, let ε > 0, ηL ≤ R ≤ ε and n ∈ N. Further, assume that (BD) holds for c > 0. Then where γ ε,n = 2r n n (M + ε) + 1 − r n n min Here M = max x∈[0,1] 2 F (x) − x , k n is the number between 1 and n for which the minimum on the right is attained and r n = n mod k n .
The proof is given in Section 4.5 below. Remark 4.2.
(a) It should be noted that the above estimate is rather of theoretical than practical interest. This is exemplified in part (b) of this remark below. We include it nevertheless, since on the one hand it demonstrates what is possible on the analytic side, and on the other hand the proof reflects and demonstrates very well how and why the nonlinearity of the dynamics complicates the efficient computation of rotation sets. (b) In order to see why the above estimates are hardly relevant for the numerical implementation, suppose that the constant c in (BD) is known (which is usually not the case) and relatively small, say, equal to 1. Even in this case, in order to achieve an apriori error bound of order 10 −2 , the integer k in the term 1 k c + 2ε L k −1 L−1 in (4.11) would have to be at least 100 -otherwise c /k > 100 -but then at the same time ε would need to be smaller than Hence, even if L is only 2, this would require to work with a box diameter of ε 2 −100 , which is hardly possible with contemporary computers. (c) Fortunately, it turns out that in the actual implementation the convergence is usually much faster than indicated by the above error estimates. One possible reason is the fact that f may possess a shadowing property. This leads to improved error bounds, where essentially the exponential term ε L k −1 L−1 can be dropped altogether. We discuss this in detail in the following subsection. 4.4. Implications of shadowing. Given a self map g : X → X of some metric space X, we say an orbit (x n ) n∈N of g δ-shadows a sequence (ξ n ) n∈N if d(x n , ξ n ) < δ for all n ∈ N. Given δ, ε > 0, we say g has the δ-shadowing property with constant ε if all ε-pseudo-orbits of g are δ-shadowed by some orbit of g. If such a constant ε = ε(δ) exists for all δ > 0, we simply say that g has the shadowing property.
Proof of Theorem 15. Under the assumptions of the theorem, any finite ε-pseudoorbit (ξ j ) n j=0 of F is δ-shadowed by some F -orbit (x j ) n j=0 , so that ξ n −ξ 0 Therefore d H (K n (F ), K ε n (F )) ≤ 2δ /n, and due to the definition of the sets K n (F ), K ε n (F ) and ε (F ) and the convergence K n (F ) → H (F ) as n → ∞ by Lemma 1, this implies ε (F ) = (F ).
As mentioned before, the shadowing property leads to improved error estimates for the set-oriented computation of the rotation set, and in particular allows to eliminate the exponential term γ ε,n in Theorem 14. Theorem 18. Suppose that f satisfies the assumptions of Theorem 14 and the bounded deviations hypothesis (BD) with constant c > 0. Let F ∈ H 2 be a lift of f and Q * n be defined by (4.5). Then The proof is postponed until the end of the next subsection.

4.5.
Quantitative estimates -proofs of Theorems 14 and 18. Throughout this section, we assume that F ∈ H 2 and Q n , Q * n are defined as in the preceding section. Let P ε n (F ) = ξ n ∈ R 2 | (ξ j ) n j=0 is an ε-pseudo-orbit of F with ξ 0 ∈ [0, 1] 2 and note that P 0 n (F ) = F n [0, 1] 2 . Then the statement of Lemma 12 can be rewritten as This leads to the following initial estimate. Lemma 19. Suppose that F ∈ H 2 is Lipschitz continuous with Lipschitz constant L > 1 and for each ε > 0 the constants η, R > 0, are chosen such that Lη ≤ R ≤ ε. Then, for every n ∈ N, we have Proof. We have Lemma 20. Suppose that F ∈ H 2 is Lipschitz continuous with Lipschitz constant L > 1 and ε > 0. Then, for all n ∈ N, we have Proof. Let (ξ j ) n j=0 be an 2ε-pseudo-orbit of F withξ 0 ∈ [0, 1] 2 . Using the fact that for all j = 1, . . . , n − 1, we recursively obtain the estimate Thus, for any v = 1 n (ξ n − ξ 0 ) ∈ K 2ε n (F ) we have Since conversely we always have K n (F ) ⊆ K 2ε n (F ), this proves the statement.
Proof. Let k ∈ N. Applying the estimate for the Hausdorff distance between the sets K 2ε k (F ) and K k (F ), denoted by κ ε,k in Lemma 20, and the assumption (BD), we obtain Conv K 2ε ⊆ Conv B c k +κ ε,k ( (F )) = B c k +κ ε,k ( (F )) .

(4.13)
For the last equality, note that (F ) is convex. Let k n ∈ {1, . . . , n} be the natural number for which the minimum in the definition of γ ε,n in (4.11) is attained. Further, let m n ∈ N, r n ∈ {0, . . . , k n − 1} be such that n = m n k n + r n . By the inclusion (3.7) in the proof of Lemma 8, since n ≥ k n , we know that (4.14) Combining (4.13) and (4.14) and by the choice of k n and the definition (4.11) of γ ε,n , we obtain Proof of Theorem 14. In Lemma 19 we deduced that Conversely, from (4.12) and Lemma 21 we obtain n +γε,n ( (F )) . Altogether, we obtain the error estimate (4.10).
Proof of Theorem 18. On the one hand, we have by Lemma 19. On the other hand, we have shown in the proof of Theorem 15 that K 2ε n (F ) ⊆ B2δ /n (K n (F )) ⊆ B1 /n (K n (F )) (note that δ < 1 /2 by assumption) and thus obtain This shows the required estimate.
See [23,6] for previous numerical studies and [2] for structurally similar examples. For specific parameter values, the rotation set of F α,β can easily be determined analytically, which allows to test the numerical algorithm in a controlled setting.
Proof. This follows directly from two elementary observations. First, we have the general estimate ( For the initialisation of the algorithm in Section 4.2 , we choose B 0 to be the standard covering of [0, 1] 2 by k 2 squares of side length 1 /k, k ∈ N. Note that we can thus choose ε = √ 2 /k in (4.2). We fix a Lipschitz constant L of F α,β (for example, L = 1 + 4π 2 works for all (α, β) ∈ [0, 1] 2 ) and set R = ε in (4.3). Moreover, for each B ∈ B 0 we choose Γ B as a standard grid of m 2 points in B, so that Γ B is √ 2 /(k(m−1))-dense in B (see Fig. 4.2). Thereby, we choose m = m(k) such that η = √ 2 /(k(m−1)) < ε /L.
In order to keep the dependence on k explicit, we will from now on write Q * k,n , instead of Q * n , for the approximation defined in (4.6). Then the assumptions of Theorem 14 with ε = √ 2 /k are satisfied and we obtain that lim k,n→∞ Q * k,n = (F α,β ), with error bound provided by (4.10) (and by Theorem 18 if F α,β has the shadowing property). Zooming in (Fig. 5.2) on the boundary of Q * 8,100 and Q * 8,200 of Figure 5.1 reveals the difference that exists between these approximations and (F 1,1 ) = [0, 1] 2 , which is of a magnitude of 10 −2 (and hence significantly smaller than the theoretical error bound in (4.10)). By Lemma 19, the 2 √ 2 /n-neighbourhood of Q * k,n covers the rotation set.
In the examples F 1,1 and F1 /2, 1 /2 above, the vertices of the polygonal rotation set are realised by periodic orbits of very low period (1, respectively 2). As discussed in [19], such rotation sets can still be accurately predicted by conventional direct approaches, but these tend to fail if the vertices correspond to periodic points of higher periods. For this reason, we next consider the map G = g 3 • g 2 • g 1 : R 2 −→ R 2 with g 1 (x, y) = x, y + 1 8 sin(5 · 2πx) , g 2 (x, y) = x + 2 5 sin(8 · 2πy), y and g 3 (x, y) = x − 1 5 , y +    Finally, we consider a slightly perturbed versionF = R • F of the above examples by introducing a slight additional rotation R : R 2 → R 2 , (x, y) → (x + r 1 , y + r 2 ), r 1 , r 2 ∈ R .
In Table 1 we collect the specific parameter values for both k and n and the perturbations r 1 and r 2 , on which we base our approximations of the related rotation sets of the perturbed maps.
Although it is difficult to check rigorously, we expect that for small perturbations these modifications should not alter the rotation sets of the above examples due to the generic structural stability of the dynamics [18]. Moreover, we expect that the vertices of the rotation sets are still realised by periodic orbits that lie close to the original ones. This fact could in principle be checked by a qualitative index argument. However, we refrain from going into detail and content ourselves with the numerical confirmation of the stability provided by Figure 5.5.  Table 1. Parameter values for the approximations shown in Figure 5.5. Approximations Q * k,n for the rotation sets of the perturbed mapsF1 /2, 1 /2 ,Ḡ,F 1, 1 /4 ,F3 /5, 3 /5 ,F3 /4,1 andF 1,1 (from top left to bottom right) according to Table 1.

Conclusion
In conclusion, the set-oriented approach to the computation of rotation sets provides better and more stable results than conventional direct approaches. Moreover, it can at least partially be backed up by rigorous convergence results, even if the theoretical error estimates are not useful in practice. The much better performance of the algorithm for specific examples finds a possible explanation in the likely presence of a shadowing property, which can again be backed up by rigorous results.
What remains is to use this new numerical method in order to perform a systematic and detailed study of the behaviour and bifurcations of rotation sets in standard parameter families, as the one given by (5.1). Of course, it is highly likely that the performance of the algorithm becomes increasingly worse as bifurcation parameters are approached, at which the rotation set changes and structural stability and shadowing therefore have to break down. Therefore, it seems feasible to carry out such investigations in collaboration with experts on scientific computing and access to high-performance computing facilities, so that at least the limits of contemporary computing capacities can be exhausted to partially counter these effects. We leave this as a task for future research.
Matlab codes will be made available on the authors' homepages.