Set-oriented numerical computation of rotation sets

TL;DR

This paper introduces a set-oriented numerical algorithm for approximating the rotation set of certain homeomorphisms on the two-torus, with proven convergence and error estimates.

Contribution

It develops a novel set-oriented method for computing rotation sets, including theoretical convergence proofs and error bounds, advancing numerical analysis in dynamical systems.

Findings

Algorithm converges as precision and iteration increase

Error estimates are provided under boundedness assumptions

Method is applicable to systems with non-empty interior rotation sets

Abstract

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.

Tables

Table 1. Table 1. Parameter values for the approximations shown in Figure 5.5 .

	${\overline{F}}_{\frac{1}{2},\frac{1}{2}}$	$\overline{G}$	${\overline{F}}_{1,\frac{1}{4}}$	${\overline{F}}_{\frac{3}{5},\frac{3}{5}}$	${\overline{F}}_{\frac{3}{4},1}$	${\overline{F}}_{1,1}$
$k$	$50$	$60$	$16$	$50$	$45$	$8$
$n$	$130$	$130$	$140$	$100$	$80$	$100$
$r_1$	$0.012$	$0.008$	$0.012$	$0.01$	$0.002$	$0.022$
$r_2$	$0.014$	$0.001$	$0.002$	$0.011$	$0.013$	$0.015$

Equations

$$F(x+k)=F(x)+k$$

$$\displaystyle\varrho\left(F\right)\ =\ \left\{v\in\mathbb{R}^{m}\;\Big{|}\;\exists n_{i}\nearrow\infty,x_{i}\in\mathbb{R}^{m}:\ \lim_{i\rightarrow\infty}\frac{F^{n_{i}}(x_{i})-x_{i}}{n_{i}}\to v\right\}.$$

$$\displaystyle K_{k}(F)\ =\ \left\{\frac{F^{k}(x)-x}{k}\;\bigg{|}\;x\in[0,1]^{m}\right\}$$

$$\varrho \left(F\right)=\bigcap _{n\ge 1}\overline{\bigcup _{k\ge n}{K}_k(F)}.$$

$${\mathrm{d}}_{\mathrm{H}}\left(K_n(F),\frac{F^n([0,1{]}^2)}{n}\right)<\frac{\sqrt{2}}{n},$$

$$\mathrm{Conv}(K_n(F))\subseteq \overline{B_{\frac{3\sqrt{2}}{n}}(K_n(F))}.$$

$$\mathrm{d}\left(F({\xi }_j),{\xi }_{j+1}\right)\le \epsilon$$

$${\left(\frac{{\xi }_{n_i}^i-{\xi }_0^i}{n_i}\right)}_{i\in \mathbb{N}}$$

$$\displaystyle K_{k}^{\varepsilon}(F)\ =\ \left\{\frac{\xi_{k}-\xi_{0}}{k}\;\bigg{|}\;\left(\xi_{j}\right)_{j=0}^{k}\text{ is an }\varepsilon\text{-pseudo-orbit of }F\text{ with }\xi_{0}\in[0,1]^{m}\right\}$$

$${\varrho }^{\epsilon }(F)=\bigcap _{n\ge 1}\overline{\bigcup _{k\ge n}{K}_k^{\epsilon }(F)}.$$

$$\varrho \left(F\right)\subseteq {\varrho }^{\epsilon }(F)$$

$${\varrho }^{\epsilon }(F)\subseteq B_{\delta }(\varrho \left(F\right)).$$

$$\frac{1}{kn}\left({\xi }_{kn}-{\xi }_0\right)=\frac{1}{k}\sum _{i=0}^{k-1}\frac{1}{n}\left({\xi }_{(i+1)n}-{\xi }_{in}\right)\in \mathrm{Conv}\left(K_n^{\epsilon }(F)\right),$$

$$K_{kn}^{\epsilon }(F)\subseteq \mathrm{Conv}(K_n^{\epsilon }(F)).$$

$$K_k^{\epsilon }(F)\ni \frac{1}{k}\left({\xi }_k-{\xi }_0\right)=\frac{1}{k}\left({\xi }_{m_{k}n+r_k}-{\xi }_{m_{k}n}\right)+\frac{1}{k}\left({\xi }_{m_{k}n}-{\xi }_0\right).$$

$$\Vert {\xi }_{m_{k}n+r_k}-{\xi }_{m_{k}n}\Vert =\Vert \sum _{i=0}^{r_k-1}{\xi }_{m_{k}n+i+1}-{\xi }_{m_{k}n+i}\Vert$$

$$\frac{1}{k}\left({\xi }_{m_k}-{\xi }_0\right)=\frac{m_{k}n}{k}\cdot \frac{1}{m_{k}n}\left({\xi }_{m_{k}n}-{\xi }_0\right)\in \frac{m_{k}n}{k}\cdot K_{m_{k}n}^{\epsilon }(F).$$

$$K_k^{\epsilon }(F)$$

$$\subseteq \overline{B_{\frac{r_k}{k}(M+\epsilon )}\left(\left(1-\frac{r_k}{k}\right)\mathrm{Conv}\left(K_n^{\epsilon }(F)\right)\right)}.$$

$${\varrho }^{\epsilon }(F)$$

$$\subseteq \bigcap _{l\ge n}\overline{\bigcup _{k\ge l}{B}_{\frac{r_k}{k}(M+\epsilon )}\left(\left(1-\frac{r_k}{k}\right)\mathrm{Conv}\left(K_n^{\epsilon }(F)\right)\right)}=\mathrm{Conv}\left(K_n^{\epsilon }(F)\right),$$

$$\frac{{\xi }_n-{\xi }_0}{n}=\frac{F^n({\xi }_0)-{\xi }_0}{n}+\frac{{\xi }_n-F^n({\xi }_0)}{n}\in B_{\delta }(K_n(F)).$$

$$K_{n_0}(F)\subseteq \overline{B_{\frac{\delta }{2}}(\varrho \left(F\right))}.$$

$$K_{n_0}^{\epsilon }(F)\subseteq \overline{B_{\frac{\delta }{2}}(K_{n_0}(F))}$$

$${\varrho }^{\epsilon }(F)$$

$$\subseteq \mathrm{Conv}\left(B_{\delta }(\varrho \left(F\right))\right)=\overline{B_{\delta }(\varrho \left(F\right))}$$

$${\varrho }_f({\mathrm{Leb}}_{{\mathbb{T}}^2})={\int }_{{\mathbb{T}}^2}\phi (x)dx,$$

$$\underset{n\to \infty }{\lim }\frac{F^n(x)-x}{n}=\frac{1}{n}\sum _{i=0}^{n-1}(\phi \circ f^i)(x)={\varrho }_f({\mathrm{Leb}}_{{\mathbb{T}}^2})$$

$$F_{1,1}:\mathbb{R}^{2}\to\mathbb{R}^{2}\ ,\quad(x,y)\mapsto\big{(}x+\sin(2\pi(y+\sin(2\pi x))),\ y+\sin(2\pi x)\big{)},$$

$${\mathrm{d}}_{\mathrm{H}}\left(\frac{1}{n}{F}^n([0,1{]}^2),K_n(F)\right)\le \frac{\sqrt{2}}{n}.$$

$$\underset{B\in {\mathcal{B}}_0}{\sup }\mathrm{diam}(B)\le \epsilon \text{ and }\bigcup _{B\in {\mathcal{B}}_0}B=[0,1{]}^2.$$

$$\displaystyle\mathcal{B}\ =\ \left\{B+t\;\big{|}\;B\in\mathcal{B}_{0},\ t\in\mathbb{Z}^{2}\right\}\ .$$

Full text

*2010 MSC numbers. Primary: 65P99, Secondary: 37M25, 37E45 *

1. Introduction

Rotation theory for orientation-preserving homeomorphisms on the circle was established by H. Poincaré in 1885 [2], 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 [3], 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 [3]. If it has non-empty interior, then the system has positive topological entropy [4] and all rational points in the interior of the rotation set are realised by periodic orbits [5]. Apart from these nowadays classical results, considerable progress has been made in the last years on the rotation theory of surface homeomorphisms (e.g. [6]–[11]), including partial results on the well-known Franks-Misiurewicz Conjecture [12]–[18] that excludes the existence of certain line segments as rotation sets.

Concerning the shape of the rotation set, it is further known that generically – in the $\mathcal{C}^{0}$-topology – it is a polygon with rational vertices [19]. For generic area-preserving torus homeomorphisms, this polygon is non-degenerate (has non-empty interior) [20]. Moreover, all rational polygons can be realised as the rotation set of a torus homeomorphism [21]. 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 [22] (see also [23]). 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 [24].

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 [20], 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 $\varepsilon$-pseudo-orbits reflecting the inaccuracies of numerical calculations. This approach leads to the definition of an $\varepsilon$-rotation set in Section 3. In the two-dimensional case, we prove that as the perturbation size $\varepsilon>0$ decreases to zero the respective $\varepsilon$-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:\mathbb{T}^{2}\longrightarrow\mathbb{T}^{2}$ homotopic to the identity, which is defined via a lift $F:\mathbb{R}^{2}\longrightarrow\mathbb{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 $\frac{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\in\mathbb{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 [20].

Numerical results for some parameter families inspired by [3, 20, 24] 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.

Acknowledgements. TJ has been supported by a Heisenberg professorship of the German Research Council (grant OE 538/6-1). TJ and KPG acknowledge funding from EU Marie-Skłodowska-Curie Innovative Training Network Critical Transitions in Complex Systems (H2020-MSCA-2014-ITN 643073 CRITICS).

2. Preliminaries

By $\mathrm{Conv}(A)$ we denote the closed convex hull of a set $A\subseteq\mathbb{R}^{m}$. The open $\varepsilon$-neighbourhood of points or sets in $\mathbb{R}^{m}$ will be denoted by $B_{\varepsilon}(x)$, respectively $B_{\varepsilon}(A)$. The closure of a set $A$ is denoted by $\overline{A}$. By $\mathrm{d}_{\mathrm{H}}(A,B)$ we denote the Hausdorff distance between two subsets $A,B$ of a metric space and also write $A_{n}\rightarrow_{\mathrm{H}}A$ if a sequence $(A_{n})_{n\in\mathbb{N}}$ converges to $A$ in the Hausdorff topology as $n\to\infty$.

For $m\in\mathbb{N}$, we let $\mathbb{T}^{m}=\mathbb{R}^{m}/\mathbb{Z}^{m}$ denote the $m$-dimensional torus. The set of all lifts $F:\mathbb{R}^{m}\to\mathbb{R}^{m}$ of continuous maps $f:\mathbb{T}^{m}\longrightarrow\mathbb{T}^{m}$ homotopic to the identity is denoted by $\mathcal{C}_{m}$. The subset of $\mathcal{C}_{m}$ consisting of lifts of torus homeomorphisms is denoted by $\mathcal{H}_{m}$. Note that $\mathcal{C}_{m}$ consists of all continuous functions $F:\mathbb{R}^{m}\longrightarrow\mathbb{R}^{m}$ that satisfy

[TABLE]

for all $x\in\mathbb{R}^{m},k\in\mathbb{Z}^{m}$, and $\mathcal{H}_{m}$ is the set of all such $F$ which are in addition homeomorphisms of the plane. In their seminal paper [3] Misiurewicz and Ziemian introduced the rotation set as

[TABLE]

Thus, $\varrho\left(F\right)$ can be viewed as the collection of all possible rotation vectors of the system. By writing

[TABLE]

for $k\in\mathbb{N}$, we alternatively obtain the rotation set as the upper Hausdorff limit of the sequence $\left(K_{n}(F)\right)_{n\in\mathbb{N}}$, that is,

[TABLE]

Moreover, in the two-dimensional case, the approximating sets $K_{n}(F)$ converge to the rotation set in the Hausdorff distance.

Lemma 1** ([3]).**

Let $F\in\mathcal{H}_{2}$. Then $K_{n}(F)\rightarrow_{\mathrm{H}}\varrho\left(F\right)\text{ as }n\to\infty.$

Furthermore, we have

Theorem 2** ([3]).**

Let $F\in\mathcal{H}_{2}$. Then the rotation set $\varrho\left(F\right)$ is convex.

Lemma 3** ([3]).**

Let $F\in\mathcal{C}_{m}$. Then $\varrho\left(F\right)\subseteq\mathrm{Conv}(K_{n}(F))$ for all $n\in\mathbb{N}$.

A generalisation of Lemma 3 to $\varepsilon$-rotation sets is given by Lemma 8 below. An important technical estimate in [3] is the following.

Lemma 4** ([3]).**

Let $G\in\mathcal{H}_{2}$. Then $\mathrm{Conv}\left(G([0,1]^{2})\right)\subseteq\overline{B_{\sqrt{2}}(G([0,1]^{2}))}$.

Applied to $G=F^{n}$ and taking into account that

[TABLE]

this immediately implies the following consequence.

Corollary 5**.**

If $F\in\mathcal{H}_{2}$ and $n\in\mathbb{N}$, then

[TABLE]

3. The $\boldsymbol{\varepsilon}$-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 $\varepsilon$-pseudo-orbits and leads to the notion of $\varepsilon$-rotation sets. Since the start of this project, these have also been described independently by Guiheneuf and Koropecki in [25].

Let $(X,\mathrm{d})$ be a metric space and $F:X\longrightarrow X$ a self-map. For $n\in\mathbb{N}$ and $\varepsilon\geq 0$, an $(n+1)$-tuple $(\xi_{j})_{j=0}^{n}$ of points $\xi_{j}\in X$ is called an $\varepsilon$-pseudo-orbit of length $n\!+\!1$ if

[TABLE]

for all $j\in\{0,\ldots,n-1\}$. In the same way infinite sequences $(\xi_{j})_{j=0}^{\infty}$ or $(\xi_{j})_{j=-\infty}^{\infty}$ satisfying (3.1) for all $j\in\mathbb{N}_{0}$ or $j\in\mathbb{Z}$, respectively, are called (infinite) $\varepsilon$-pseudo-orbits. Obviously, every proper orbit is an $\varepsilon$-pseudo-orbit for all $\varepsilon>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\in\mathcal{C}_{m}$ and $\varepsilon\geq 0$, we now define the $\varepsilon$-rotation set to be the set of all accumulation points of sequences of the form

[TABLE]

where $n_{i}\to\infty$ as $i\to\infty$ and $(\xi^{i}_{j})_{j=0}^{n_{i}}$ is an $\varepsilon$-pseudo-orbit of $F$ of length $n_{i}+1$ for each $i\in\mathbb{N}$. Analogous to (2.3), we let

[TABLE]

for $k\in\mathbb{N}$ and can alternatively define the $\varepsilon$-rotation set by

[TABLE]

Note that since $F$ commutes with integer translations, one could also replace $\xi_{0}\in[0,1]^{m}$ by $\xi_{0}\in\mathbb{R}^{m}$ in the definition of $K_{k}^{\varepsilon}(F)$. For every $k\in\mathbb{N}$, the inclusion $K_{k}(F)\subseteq K_{k}^{\varepsilon}(F)$ is apparent, so that due to (2.3) and (3.2), we obtain

[TABLE]

for all $\varepsilon\geq 0$. We omit the elementary proof of the following lemma.

Lemma 6**.**

Let $F\in\mathcal{C}_{m},\varepsilon>0$ and $M=\max\nolimits_{x\in[0,1]^{m}}\|F(x)-x\|$. Then

• (i)

for each $k\in\mathbb{N}$ the set $K_{k}^{\varepsilon}(F)$ is non-empty and compact, and thus the same holds for $\varrho^{\varepsilon}(F)$;

• (ii)

$\varrho\left(F\right)\subseteq\overline{B_{M}(0)}$* and $\varrho^{\varepsilon}(F)\subseteq\overline{B_{M+\varepsilon}(0)}$.*

We now aim to show convergence of the $\varepsilon$-rotation sets as $\varepsilon$ decreases to zero. This presents the main result of this section. An alternative proof can be found in [25]. 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.

Theorem 7**.**

Let $F\in\mathcal{H}_{2}$. Then $\varrho^{\varepsilon}(F)\rightarrow_{\mathrm{H}}\varrho\left(F\right)\text{ as }\varepsilon\to 0.$

Since for $\varrho\left(F\right)\subseteq\varrho^{\varepsilon}(F)$ for all $\varepsilon>0$ by (3.3), it remains to show that for every $\delta>0$ there is an $\varepsilon>0$ such that

[TABLE]

We split the proof into Lemmas 8 and 9, beginning with a formulation of Lemma 3 in terms of $\varepsilon$-rotation sets.

Lemma 8**.**

Let $F\in\mathcal{C}_{m}$ and $\varepsilon\geq 0$. Then $\varrho^{\varepsilon}(F)\subseteq\mathrm{Conv}(K_{n}^{\varepsilon}(F))$ for all $n\in\mathbb{N}$.

Proof.

Let $k,n\in\mathbb{N}$. For an $\varepsilon$-pseudo-orbit $(\xi_{j})_{j=0}^{kn}$ with $\frac{1}{kn}\left(\xi_{kn}-\xi_{0}\right)\in K_{kn}^{\varepsilon}(F)$ we have

[TABLE]

since the tuples $(\xi_{in},\xi_{in+1},\ldots,\xi_{(i+1)n})$ are $\varepsilon$-pseudo-orbits of $F$ of length $n+1$ and hence $\frac{1}{n}(\xi_{(i+1)n}-\xi_{in})\in K_{n}^{\varepsilon}(F)$ for each $i\in\{0,\ldots,k-1\}$. Consequently, for all $k,n\in\mathbb{N}$ we have

[TABLE]

Now let $k,n\in\mathbb{N}$ with $k\geq n$. Each number $k$ can uniquely be split into $k=m_{k}n+r_{k}$, where $m_{k}\in\mathbb{N}$ is the integer part of $\frac{k}{n}$ and $r_{k}\in\{0,\ldots,n-1\}$. If $(\xi_{j})_{j=0}^{k}$ is an $\varepsilon$-pseudo-orbit of length $k+1$, we find

[TABLE]

Since the tuple $(\xi_{m_{k}n+j})_{j=0}^{r_{k}}$ is an $\varepsilon$-pseudo-orbit of $F$ of length $r_{k}+1$, we obtain that

[TABLE]

Thus, the norm of the first summand in (3.6) is bounded by $\frac{r_{k}}{k}(M+\varepsilon)$. For the second summand, we have

[TABLE]

By (3.5), $K_{m_{k}n}^{\varepsilon}(F)$ is included in $\mathrm{Conv}\left(K_{n}^{\varepsilon}(F)\right)$. Altogether, we obtain

[TABLE]

For a fixed value of $n$, the inclusion (3.7) is valid for all $k\geq n$. Therefore, for each $n\in\mathbb{N}$, the representation (3.2) of the $\varepsilon$-rotation set yields

[TABLE]

due to the convergence $\frac{r_{k}}{k}\to 0$ as $k\to\infty$. ∎

Lemma 9**.**

Let $F\in\mathcal{C}_{m}$ and $n\in\mathbb{N}$. Then $K_{n}^{\varepsilon}(F)\rightarrow_{\mathrm{H}}K_{n}(F)\text{ as }\varepsilon\to 0.$

Proof.

The inclusion $K_{n}(F)\subseteq K_{n}^{\varepsilon}(F)$ is evident for every $\varepsilon>0$. Let $\delta>0$. Due to the continuity of $F$, there exists an $\varepsilon_{0}>0$ such that $\|F^{n}(\xi_{0})-\xi_{n}\|<n\delta$ for every $\varepsilon$-pseudo-orbit $(\xi_{j})_{j=0}^{n}$ of $F$ with $0<\varepsilon\leq\varepsilon_{0}$. Therefore, we obtain

[TABLE]

As $\delta>0$ was chosen arbitrarily and $n\in\mathbb{N}$ is fixed, this proves the convergence $K_{n}^{\varepsilon}(F)\rightarrow_{\mathrm{H}}K_{n}(F)$ as $\varepsilon\to 0$. ∎

Proof of Theorem 7. Recall that in order to prove Theorem 7, we have to show the inclusion (3.4). Let $\delta>0$. By Lemma 1, the rotation set $\varrho\left(F\right)$ 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}\in\mathbb{N}$ sufficiently large, so that

[TABLE]

Furthermore, by Lemma 9 there exists an $\varepsilon_{0}>0$ sufficiently small, such that

[TABLE]

for all $0<\varepsilon\leq\varepsilon_{0}$. Applying Lemma 8, together with (3.8) and (3.9) we obtain

[TABLE]

for every $0<\varepsilon\leq\varepsilon_{0}$. $\overline{B_{\delta}(\varrho\left(F\right))}$ is convex itself, since the rotation set is convex by Theorem 2. As $\delta>0$ was arbitrary, this shows the asserted convergence $\varrho^{\varepsilon}(F)\rightarrow_{\mathrm{H}}\varrho\left(F\right)$ as $\varepsilon\to 0$. ∎

4. 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 $\Gamma\subseteq[0,1]^{2}$ of $N\times N$ points, to compute the normalised displacement vectors $\frac{1}{n}(F^{n}(x)-x)$ for all $x\in\Gamma$ and some large $n\in\mathbb{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 $\mathcal{C}^{0}$-topology, see [20]): Suppose that an area-preserving torus homeomorphism $f$ has a rotation set with non-empty interior, and that the Lebesgue measure $\mathrm{Leb}_{\mathbb{T}^{2}}$ on $\mathbb{T}^{2}$ is ergodic with respect to $f$. Then there exists a well-defined rotation vector

[TABLE]

where $\varphi(x)=F(x)-x$ is interpreted as a function on the torus, and we have

[TABLE]

for Lebesgue-a.e. $x\in\mathbb{T}^{2}$. Hence, if $n$ above is chosen too large, the rotation vectors coming from starting points in the grid $\Gamma$ will almost surely be arbitrarily close to $\varrho_{f}(\mathrm{Leb}_{\mathbb{T}^{2}})$. In this case, the numerical approximation of the rotation set will show the singleton $\varrho_{f}(\mathrm{Leb}_{\mathbb{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

[TABLE]

which exhibits $\varrho_{f_{1,1}}(\mathrm{Leb}_{\mathbb{T}^{2}})=(0,0)$ and is known to have the rotation set $\varrho\left(F_{1,1}\right)=[-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 [20], 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\sim 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 $\varrho\left(F\right)$, and moreover we have the elementary estimate

[TABLE]

Thus, we obtain

Corollary 10**.**

Let $F\in\mathcal{H}_{2}$. Then $\frac{1}{n}F^{n}([0,1]^{2})\rightarrow_{\mathrm{H}}\varrho\left(F\right)$ as $n\to\infty$.

This allows to focus on the sets $\frac{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 $\frac{1}{n}F^{n}([0,1]^{2})$ for some large $n\in\mathbb{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 [26], provides numerical methods for the analysis of dynamical systems by set-oriented techniques. As mentioned above, on the theoretical level this corresponds to considering $\varepsilon$-pseudo-orbits.

For a given maximal diameter $\varepsilon>0$ let $\mathcal{B}_{0}$ be a collection of compact sets $B\subset[0,1]^{2}$ with

[TABLE]

Note that $\mathcal{B}_{0}$ can be considered as a covering of $\mathbb{T}^{2}$, whose lift to $\mathbb{R}^{2}$ is then given by

[TABLE]

The elements of the collections $\mathcal{B}_{0}$ and $\mathcal{B}$ will be referred to as boxes. Further, given $\eta>0$, for each box $B\in\mathcal{B}$, let $\Gamma_{B}\subseteq B$ be a finite set of points which is $\eta$-dense in $B$ and chosen such that $\Gamma_{B+t}=\Gamma_{B}+t$ for every integer vector $t\in\mathbb{Z}^{2}$. Finally, fix $R>0$. Then, for $F\in\mathcal{H}_{2}$, we now generate a sequence of sets $\left(Q_{n}^{\ast}\right)_{n\in\mathbb{N}}$ (approximations of the rotation set) according to the following algorithm.

Initialisation:

$Q_{0}$

$\ =\ [0,1]^{2}$

Box images:

Given $B\in\mathcal{B}$, let the box image of $B$ be defined as

$\mathcal{I}(B)$

$\ =\ \{B^{\prime}\in\mathcal{B}\mid\exists x\in\Gamma_{B}:d(F(x),B^{\prime})\leq R\}$

[TABLE]

(see Fig. 4.2).

Iteration:

For $k=0,\ldots,n-1$ generate the box coverings

$\mathcal{B}_{k+1}$

$\ =\ \left\{B^{\prime}\in\mathcal{B}\mid\exists B\in\mathcal{B}_{k}:B^{\prime}\in\mathcal{I}(B)\right\},$

[TABLE]

$Q_{k+1}$

$\ =\ \bigcup\nolimits_{B\in\mathcal{B}_{k+1}}\!\!B.$

[TABLE]

Normalisation:

$Q_{n}^{\ast}$

$\ =\ \frac{1}{n}Q_{n},\ n\in\mathbb{N}$

[TABLE]

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\in\mathcal{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 $\frac{1}{n}F^{n}\left([0,1]^{2}\right)$, 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 $\nicefrac{{2\sqrt{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\in\mathcal{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}^{\ast}$ and the rotation set $\varrho\left(F\right)$ 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\in\mathcal{H}_{2}$ is Lipschitz continuous with Lipschitz constant $L>1$ and for each $\varepsilon>0$ the constants $\eta,R>0$ are chosen such that $L\eta\leq R\leq\varepsilon$. Then

[TABLE]

in the sense that for all $\delta>0$ there exist $\varepsilon_{0}>0$ and $n_{0}\in\mathbb{N}$ such that if $\varepsilon$ in (4.2) satisfies $\varepsilon\in(0,\varepsilon_{0}]$ and $n\geq n_{0}$, then $\mathrm{d}_{\mathrm{H}}(Q_{n}^{\ast},\varrho\left(F\right))<\delta$.

Remark 4.1**.**

For theoretical purposes, one may also want to ignore the fact that the images of the boxes $B\in\mathcal{B}$ can only be approximated via test points, and define alternative sequences \big{(}\tilde{Q}_{n}\big{)}_{n\in\mathbb{N}},\big{(}\tilde{Q}_{n}^{*}\big{)}_{n\in\mathbb{N}} by using the precise images via

[TABLE]

This corresponds to setting the parameters $\eta$ and $R$ to zero. Then the above convergence result is still valid for the new sequence, that is, $\lim_{\begin{subarray}{c}n\to\infty\\ \varepsilon\to 0\end{subarray}}\ \mathrm{d}_{\mathrm{H}}(\tilde{Q}_{n}^{\ast},\varrho\left(F\right))=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\in\mathcal{H}_{2}$ is Lipschitz continuous with Lipschitz constant $L>1$ and for each $\varepsilon>0$ the constants $\eta,R>0$ are chosen such that $L\eta\leq R\leq\varepsilon$. Then *

[TABLE]

Proof.

First, suppose that $v=\frac{1}{n}F^{n}(x_{0})\in\frac{1}{n}F^{n}\left([0,1]^{2}\right)$ and let $x_{j}=F^{j}(x_{0})$. Choose boxes $B_{j}\in\mathcal{B}$ such that $x_{j}\in B_{j}$. We claim that $B_{j}\subseteq Q_{j}$ for $j=0,\ldots,n$, so that finally $x_{n}=nv\in B_{n}\subseteq Q_{n}$ and hence $v\in Q_{n}^{*}$.

For $j=0$ the claim is obvious. Therefore, suppose that $B_{j}\subseteq Q_{j}$. Then there exists a test point $y\in\Gamma_{B_{j}}$ $\eta$-close to $x_{j}$, so that $d(F(y),x_{j+1})\leq L\eta\leq R$. Since $x_{j+1}\in B_{j+1}$, this implies $B_{j+1}\in\mathcal{I}(B)$ and thus $B_{j+1}\in Q_{j+1}$.

This shows the first inclusion $\frac{1}{n}F^{n}\left([0,1]^{2}\right)\subseteq Q^{*}_{n}$. In order to show the second inclusion, suppose that $v\in Q^{*}_{n}$. Then, by definition of $Q_{n}$, there exists a sequence of boxes $B_{0},\ldots,B_{n}$ such that $nv\in B_{n}$ and $B_{j+1}\in\mathcal{I}(B_{j})$ for all $j=0,\ldots,n-1$. By definition of the box images, there exists a sequence of test points $\xi_{j}\in\Gamma_{B_{j}}$, $j=0,\ldots,n-1$, such that $d(F(\xi_{j}),B_{j+1})\leq R$. Hence

[TABLE]

for all $j=0,\ldots,n-2$, and this remains true for $j=n-1$ if we let $\xi_{n}=nv\in B_{n}$. This means that $(\xi_{j})_{j=0}^{n}$ is a $2\varepsilon$-pseudo-orbit and thus yields the required second inclusion. ∎

Corollary 13**.**

In the situation of Lemma 12, we have

[TABLE]

Proof of Theorem 11.

First, by Theorem 7 and the convexity of the rotation set, there exists $\varepsilon_{0}>0$ such that

[TABLE]

Further, by the definition of $\varrho^{2\varepsilon_{0}}(F)$, there exists $n_{0}\in\mathbb{N}$ such that for all $n\geq n_{0}$ we have $K_{n}^{2\varepsilon_{0}}(F)\subseteq B_{\nicefrac{{\delta}}{{3}}}\left(\varrho^{2\varepsilon_{0}}(F)\right)$ and hence

[TABLE]

As for all $n\in\mathbb{N}$ we have $\varrho\left(F\right)\subseteq\mathrm{Conv}\left(K_{n}(F)\right)\subseteq\mathrm{Conv}\left(K_{n}^{2\varepsilon_{0}}(F)\right)$ (see Lemma 3), we obtain that

[TABLE]

At the same time, Lemma 1 implies that we can choose $n_{0}$ such that

[TABLE]

for all $n\geq n_{0}$. If $n_{0}$ is chosen such that $\nicefrac{{\sqrt{2}}}{{n_{0}}}<\nicefrac{{\delta}}{{3}}$ and $n\geq n_{0}$, then using (4.9) together with the first inclusion in Corollary 13, we obtain

[TABLE]

Conversely, (4.8) together with the second inclusion in Corollary 13 yield

[TABLE]

Together with the fact that $K_{n}^{2\varepsilon}(F)\subseteq K_{n}^{2\varepsilon_{0}}(F)$ for all $\varepsilon<\varepsilon_{0}$, this means that

[TABLE]

whenever $n\geq n_{0}$ and $\varepsilon\in(0,\varepsilon_{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 $\varrho\left(F\right)$. 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

[TABLE]

This fact has been proven for diffeomorphisms in [27] and the result was later generalised to homeomorphisms in [13]. 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\in\mathcal{H}_{2}$ is Lipschitz continuous with Lipschitz constant $L>1$, let $\varepsilon>0$, $\eta L\leq R\leq\varepsilon$ and $n\in\mathbb{N}$. Further, assume that (BD) holds for $c>0$. Then

[TABLE]

where

[TABLE]

Here $M=\max_{x\in[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\bmod 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 \frac{1}{k}\big{(}c+2\varepsilon\frac{L^{k}-1}{L-1}\big{)} in (4.11) would have to be at least $100$ – otherwise $\nicefrac{{c}}{{k}}>100$ – but then at the same time $\varepsilon$ would need to be smaller than \big{(}\frac{2(L^{k}-1)}{L-1}\big{)}^{{}_{-1}}. Hence, even if $L$ is only $2$, this would require to work with a box diameter of $\varepsilon\simeq 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 $\varepsilon\frac{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\to X$ of some metric space $X$, we say an orbit $(x_{n})_{n\in\mathbb{N}}$ of $g$ $\delta$-shadows a sequence $(\xi_{n})_{n\in\mathbb{N}}$ if $d(x_{n},\xi_{n})<\delta$ for all $n\in\mathbb{N}$. Given $\delta,\varepsilon>0$, we say $g$ has the $\delta$-shadowing property with constant $\varepsilon$ if all $\varepsilon$-pseudo-orbits of $g$ are $\delta$-shadowed by some orbit of $g$. If such a constant $\varepsilon=\varepsilon(\delta)$ exists for all $\delta>0$, we simply say that $g$ has the shadowing property.

Theorem 15**.**

Suppose $f$ is a torus homeomorphism homotopic to the identity with lift $F\in\mathcal{H}_{2}$ and $\delta,\gamma\in(0,\nicefrac{{1}}{{2}})$ are such that $d(x,y)<\delta$ implies $d(f(x),f(y))<\nicefrac{{1}}{{2}}-\gamma$ for all $x,y\in\mathbb{T}^{2}$. Further, assume that $f$ has the $\delta$-shadowing property with constant $\varepsilon\in(0,\gamma)$. Then

[TABLE]

Corollary 16**.**

If a torus homeomorphism $f$ homotopic to the identity has the shadowing property, then there exists $\varepsilon_{0}>0$ such that $\varrho^{\varepsilon}(F)=\varrho\left(F\right)$ for any $\varepsilon\in(0,\varepsilon_{0}]$ and any lift $F$ of $f$.

For the proof of Theorem 15, we need the following elementary statement.

Lemma 17**.**

Suppose $f,F,\delta,\gamma,\varepsilon$ are chosen as in Theorem 15. Then $F$ has the $\delta$-shadowing property with constant $\varepsilon$.

Proof.

Suppose $(\hat{\xi}_{n})_{n\in\mathbb{N}}$ is an $\varepsilon$-pseudo-orbit for $F$. Then $\xi_{n}=\pi(\hat{\xi}_{n})$ defines an $\varepsilon$-pseudo-orbit of $f$, and we can therefore find some $x_{0}\in\mathbb{T}^{2}$ such that $d\left(f^{n}(x_{0}),\xi_{n}\right)<\delta$ for all $n\geq 0$. Let $\hat{x}_{0}$ be the unique lift of $x_{0}$ in B_{\delta}\big{(}\hat{\xi}_{0}\big{)}. We claim that $d(F^{n}(\hat{x}_{0}),\hat{\xi}_{n})<\delta$ for all $n\geq 0$, so that the orbit of $\hat{x}_{0}$ is the required $\delta$-shadowing orbit.

For $n=0$ there is nothing to prove. If $d(F^{n}(\hat{x}_{0}),\hat{\xi}_{n})<\delta$ for some $n\geq 0$, then F^{n+1}(\hat{x}_{0})\in B_{1/2-\gamma}\big{(}F(\hat{\xi}_{n})\big{)}\subseteq B_{1/2}\big{(}\hat{\xi}_{n+1}\big{)}. However, as B_{1/2}\big{(}\hat{\xi}_{n+1}\big{)} projects injectively to $\mathbb{T}^{2}$ and $d(\pi(F^{n+1}(\hat{x}_{0}),\pi(\hat{\xi}_{n+1}))=d(f^{n+1}(x_{0}),\xi_{n+1})<\delta$, we also obtain $d(F^{n+1}(\hat{x}_{0}),\hat{\xi}_{n+1}))<\delta$ as required. ∎

Proof of Theorem 15.

Under the assumptions of the theorem, any finite $\varepsilon$-pseudo-orbit $(\hat{\xi}_{j})_{j=0}^{n}$ of $F$ is $\delta$-shadowed by some $F$-orbit $(x_{j})_{j=0}^{n}$, so that

[TABLE]

Therefore $\mathrm{d}_{\mathrm{H}}(K_{n}(F),K^{\varepsilon}_{n}(F))\leq\nicefrac{{2\delta}}{{n}}$, and due to the definition of the sets $K_{n}(F)$, $K^{\varepsilon}_{n}(F)$ and $\varrho^{\varepsilon}(F)$ and the convergence $K_{n}(F)\rightarrow_{\mathrm{H}}\varrho\left(F\right)$ as $n\to\infty$ by Lemma 1, this implies $\varrho^{\varepsilon}(F)=\varrho\left(F\right)$. ∎

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 $\gamma_{\varepsilon,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\in\mathcal{H}_{2}$ be a lift of $f$ and $Q_{n}^{\ast}$ be defined by (LABEL:E_def_Qn). Then

[TABLE]

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\in\mathcal{H}_{2}$ and $Q_{n},\ Q_{n}^{\ast}$ are defined as in the preceding section. Let

[TABLE]

and note that $P^{0}_{n}(F)=F^{n}\left([0,1]^{2}\right)$. Then the statement of Lemma 12 can be rewritten as

[TABLE]

This leads to the following initial estimate.

Lemma 19**.**

Suppose that $F\in\mathcal{H}_{2}$ is Lipschitz continuous with Lipschitz constant $L>1$ and for each $\varepsilon>0$ the constants $\eta,R>0$, are chosen such that $L\eta\leq R\leq\varepsilon$. Then, for every $n\in\mathbb{N}$, we have

[TABLE]

Proof.

We have

[TABLE]

∎

Lemma 20**.**

Suppose that $F\in\mathcal{H}_{2}$ is Lipschitz continuous with Lipschitz constant $L>1$ and $\varepsilon>0$. Then, for all $n\in\mathbb{N}$, we have

[TABLE]

Proof.

Let $(\hat{\xi}_{j})_{j=0}^{n}$ be an $2\varepsilon$-pseudo-orbit of $F$ with $\hat{\xi}_{0}\in[0,1]^{2}$. Using the fact that $\|\hat{\xi}_{1}-F(\hat{\xi}_{0})\|\leq 2\varepsilon$ and

[TABLE]

for all $j=1,\ldots,n-1$, we recursively obtain the estimate

[TABLE]

Thus, for any $v=\frac{1}{n}(\xi_{n}-\xi_{0})\in K_{n}^{2\varepsilon}(F)$ we have

[TABLE]

so that

[TABLE]

Since conversely we always have $K_{n}(F)\subseteq K^{2\varepsilon}_{n}(F)$, this proves the statement. ∎

Lemma 21**.**

Suppose that $F\in\mathcal{H}_{2}$ is Lipschitz continuous with Lipschitz constant $L>1$ and $\varepsilon>0$. Further, assume that $F$ additionally satisfies (BD) with $c>0$. Then for all $n\geq 0$ we have

[TABLE]

where $\gamma_{\varepsilon,n}$ is defined as in (4.11).

Proof.

Let $k\in\mathbb{N}$. Applying the estimate for the Hausdorff distance between the sets $K_{k}^{2\varepsilon}(F)$ and $K_{k}(F)$, denoted by $\kappa_{\varepsilon,k}$ in Lemma 20, and the assumption (BD), we obtain

[TABLE]

For the last equality, note that $\varrho\left(F\right)$ is convex. Let $k_{n}\in\{1,\ldots,n\}$ be the natural number for which the minimum in the definition of $\gamma_{\varepsilon,n}$ in (4.11) is attained. Further, let $m_{n}\in\mathbb{N},\ r_{n}\in\{0,\ldots,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\geq k_{n}$, we know that

[TABLE]

Combining (4.13) and (4.14) and by the choice of $k_{n}$ and the definition (4.11) of $\gamma_{\varepsilon,n}$, we obtain

[TABLE]

For the inclusion from the second to the third line, note that $\varrho\left(F\right)\subseteq\overline{B_{M}(0)}$, so that $\left(1-\frac{r_{n}}{n}\right)\varrho\left(F\right)\subseteq\overline{B_{\frac{r_{n}}{n}M}\left(\varrho\left(F\right)\right)}$. ∎

Proof of Theorem 14..

In Lemma 19 we deduced that

[TABLE]

Conversely, from (4.12) and Lemma 21 we obtain

[TABLE]

Altogether, we obtain the error estimate (4.10).  ∎

Proof of Theorem 18..

On the one hand, we have

[TABLE]

by Lemma 19. On the other hand, we have shown in the proof of Theorem 15 that $K_{n}^{2\varepsilon}(F)\subseteq B_{\nicefrac{{2\delta}}{{n}}}(K_{n}(F))\subseteq B_{\nicefrac{{1}}{{n}}}(K_{n}(F))$ (note that $\delta<\nicefrac{{1}}{{2}}$ by assumption) and thus obtain

[TABLE]

This shows the required estimate. ∎

5. Numerical implementation and results

In order to implement the above algorithm and to apply it to some specific examples, we consider a standard family of (lifts of) torus diffeomorphisms given by

[TABLE]

where $\alpha,\beta\in\mathbb{R}$. Note that $F_{\alpha,\beta}$ is obtained as the composition of two skew shifts, $F_{\alpha,\beta}=F_{\alpha}\circ F_{\beta}$, where

[TABLE]

See [24, 7] for previous numerical studies and [3] for structurally similar examples. For specific parameter values, the rotation set of $F_{\alpha,\beta}$ can easily be determined analytically, which allows to test the numerical algorithm in a controlled setting.

Lemma 22**.**

$\varrho\left(F_{1,1}\right)=[-1,1]^{2}.$**

Proof.

This follows directly from two elementary observations. First, we have the general estimate $\varrho\left(F_{\alpha,\beta}\right)\subseteq[-\alpha,\alpha]\times[-\beta,\beta]$, as $\alpha$ and $\beta$ are the maximal step sizes in the horizontal and vertical direction. Hence, $\varrho\left(F_{1,1}\right)\subseteq[-1,1]^{2}$.

Conversely, it is easily checked that the rotation vectors $(1,1),(-1,1),(1,-1)$ and $(-1,-1)$ are realised by the fixed points $(\nicefrac{{1}}{{4}},\nicefrac{{1}}{{4}}),(\nicefrac{{3}}{{4}},\nicefrac{{1}}{{4}}),(\nicefrac{{1}}{{4}},\nicefrac{{3}}{{4}})$ and $(\nicefrac{{3}}{{4}},\nicefrac{{3}}{{4}})$. By convexity, this means that $[-1,1]^{2}\subseteq\varrho\left(F_{1,1}\right)$. ∎

For the initialisation of the algorithm in Section 4.2 , we choose $\mathcal{B}_{0}$ to be the standard covering of $[0,1]^{2}$ by $k^{2}$ squares of side length $\nicefrac{{1}}{{k}}$, $k\in\mathbb{N}$. Note that we can thus choose $\varepsilon=\nicefrac{{\sqrt{2}}}{{k}}$ in (4.2). We fix a Lipschitz constant $L$ of $F_{\alpha,\beta}$ (for example, $L=1+4\pi^{2}$ works for all $(\alpha,\beta)\in[0,1]^{2}$) and set $R=\varepsilon$ in (LABEL:e.box_images). Moreover, for each $B\in\mathcal{B}_{0}$ we choose $\Gamma_{B}$ as a standard grid of $m^{2}$ points in $B$, so that $\Gamma_{B}$ is $\nicefrac{{\sqrt{2}}}{{(k(m-1))}}$-dense in $B$ (see Fig. 4.2). Thereby, we choose $m=m(k)$ such that $\eta=\nicefrac{{\sqrt{2}}}{{(k(m-1))}}<\nicefrac{{\varepsilon}}{{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 (LABEL:e.Qnstar). Then the assumptions of Theorem 14 with $\varepsilon=\nicefrac{{\sqrt{2}}}{{k}}$ are satisfied and we obtain that $\lim_{k,n\to\infty}Q_{k,n}^{*}=\varrho\left(F_{\alpha,\beta}\right)$, with error bound provided by (4.10) (and by Theorem 18 if $F_{\alpha,\beta}$ has the shadowing property). Figure 5.1 shows $Q_{k,n}^{*}$ for $F_{1,1}$ for $k=8$ and different values of $n$.

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 $\varrho\left(F_{1,1}\right)=[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 $\nicefrac{{2\sqrt{2}}}{{n}}$-neighbourhood of $Q_{k,n}^{\ast}$ covers the rotation set.

As further examples, we consider the maps $F_{\nicefrac{{1}}{{2}},\nicefrac{{1}}{{2}}},F_{\nicefrac{{3}}{{5}},\nicefrac{{3}}{{5}}}$ and $F_{\nicefrac{{3}}{{4}},1}$. In the first case, we still have an a priori lower bound for the rotation set: $\varrho\left(F_{\nicefrac{{1}}{{2}},\nicefrac{{1}}{{2}}}\right)$ contains the square spanned by the points $(\pm\nicefrac{{1}}{{2}},0)$ and $(0,\pm\nicefrac{{1}}{{2}})$, since these rotation vectors are realised by the two-periodic points $(0,\nicefrac{{1}}{{4}}),(\nicefrac{{1}}{{4}},0)$,$(0,\nicefrac{{3}}{{4}})$ and $(\nicefrac{{3}}{{4}},0)$. The numerical approximation in Figure 5.3 indicates that this square indeed is the rotation set of $F_{\nicefrac{{1}}{{2}},\nicefrac{{1}}{{2}}}$ (recall here that our algorithm never underestimates).

In the examples $F_{1,1}$ and $F_{\nicefrac{{1}}{{2}},\nicefrac{{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 [20], 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}\circ g_{2}\circ g_{1}:\mathbb{R}^{2}\longrightarrow\mathbb{R}^{2}$ with

[TABLE]

The related rotation set can be determined analytically as the rectangle $\varrho\left(G\right)=[-\nicefrac{{3}}{{5}},\nicefrac{{1}}{{5}}]\times[\nicefrac{{1}}{{8}},\nicefrac{{3}}{{8}}]$ since its vertices correspond to elliptic periodic orbits of period $40$ (those of the points $(\nicefrac{{3}}{{20}},\nicefrac{{3}}{{32}}),(\nicefrac{{1}}{{20}},\nicefrac{{3}}{{32}}),(\nicefrac{{1}}{{20}},\nicefrac{{1}}{{32}})$ and $(\nicefrac{{3}}{{20}},\nicefrac{{1}}{{32}})$). Figure 5.4 shows the approximate rotation set $Q_{60,130}^{*}$ for $G$.

Finally, we consider a slightly perturbed version $\bar{F}=R\circ F$ of the above examples by introducing a slight additional rotation

[TABLE]

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 [19]. 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.

6. 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.