Accurate Bounds on Lyapunov Exponents for Expanding Maps of the Interval

In this short note we describe a simple but remarkably effective method for rigorously estimating Lyapunov exponents for expanding maps of the interval. We illustrate the applicability of this method with some standard examples.


Introduction
Lyapunov exponents give a well known characterization of the instability in a dynamical system by quantifying how nearby orbits separate. In particular, a non-zero Lyapunov exponent with respect to an invariant ergodic measure implies that typical nearby orbits separate exponentially quickly. It is therefore useful to have a rigorous and effective estimate of these values, in particular, in the setting of one dimensional expanding maps for an absolutely continuous invariant probability measure. This problem has attracted the attention of many authors who have employed a variety of different methods (see [5,10,13,22]).
In this paper we will consider the nice class N C of expanding piecewise analytic mixing Markov maps of the interval. We recall the definition. Definition 1.1. Let I = [a, b] be a closed interval. We say that a map f belongs to the class N C(I ) if there exists a partition a = x 1 < x 2 < · · · < x n+1 = b such that: 4. The map f is topologically mixing 1 (i.e., for any non-empty open sets U, V ⊂ I there exists n 0 ≥ 1 such the for all n ≥ n 0 we have U ∩ f −n V = ∅).
For every map f ∈ N C(I ) there exists a unique absolutely continuous f -invariant probability measure dμ = ρ(x)dx on I [7]. In particular, the measure μ is ergodic. Furthermore, every map of the class N C(I ) is invertible on each of the intervals (x j , x j+1 ), in particular, there exist analytic maps f jk : (x j , x j+1 ) → (x k , x k+1 ), such that f ( f jk (x)) = x whenever f ((x k , x k+1 )) ∩ (x j , x j+1 ) = ∅. The maps f jk are called inverse branches of f .
A standard approach to constructing the measure μ is to use transfer operators. Let us denote by B the space of analytic functions on the disjoint union n k=1 [x k , x k+1 ]. We can introduce a one-parameter family of linear operators L t : B → B (t ∈ R) called transfer operators defined in terms of inverse branches of f : where χ [x j ,x j+1 ] is the indicator function of the interval [x j , x j+1 ]. In the special case that f is full branched, i.e., f ((x k , x k+1 )) = I for k = 1, . . . , n, we can denote the inverses f k : I → (x k , x k+1 ), i.e., f ( f k (x)) = x for all a ≤ x ≤ b. The transfer operators L t : C ω (I ) → C ω (I ) (t ∈ R) then take the form Most of our examples will be of this type. It is well known [7] that the positive density ρ ∈ C ω (I ) of the measure μ is characterized as a fixed point for the operator L 1 , corresponding to the parameter choice t = 1. Nevertheless, including this operator into a one parameter family will serve us well later.
We can now define the Lyapunov exponent of the system (I, f, μ) which quantifies the sensitivity of typical orbits on initial conditions. Definition 1.2. We define the Lyapunov exponent for the map f and its stationary measure μ by Remark 1.3. This value coincides with the metric entropy h(μ) of the measure μ by the Rokhlin's formula [17].
Since the measure μ is ergodic, applying the Birkhoff ergodic theorem one can see that for μ-almost all x ∈ I we get There are various methods used to estimate the Lyapunov exponents. Probably the most famous are Ulam's method and finite section methods [13]. Another approach is based on periodic points method [11]. Recent work by Wormell [22] is based on the Galerkin spectral method originally developed for PDEs. In this note we present an alternative approach, which starts with the spectral Chebyshev collocation method, also initially developed for PDEs [6, §3]. In dynamical systems it has been used succsefully by Babenko and Yuriev in their solution of the Gauss problem [1] and by Babenko in his computation of the fixed point of the renormalisation operator for the period-doubling map [2]. More recently, Bandtlow and Slipantschuk [4] applied the spectral Chebyshev collocation method in the setting of holomorphic contractions and showed it to be an effective approach to obtaining the spectral data of the transfer operator. In constrast to our approach in [4] the error is estimated directly using the spectral properties of the transfer operator. In our approach, we combine the Chebyshev collocation method with a small amount of thermodynamic formalism (involving the pressure function) and a classical min-max method. The main advantage of this combination of ideas is that it provides an efficient and effective way to estimate Lyapunov exponents and gives rigorous estimates with validated error bounds.
The main results we present in this note are the following. The first theorem gives a method for obtaining rigorous bounds on the Lyapunov exponent. Theorem 1.4. Let f : I → I be an expanding piecewise analytic mixing Markov map of the interval with absolutely continuous probability measure μ. Assume that for some ε > 0 there exists a pair of positive functions 2 p, q : I → R + and a pair of numbers 0 < α < β such that Then the following double inequality holds: Remark 1.5. The idea behind Theorem 1.4 is that for the class of maps we consider for any positive function p the supremum of the ratio L t p p gives an upper bound on the leading eigenvalue of L t .
Note that if the function p is close to the leading eigenfunction of the operator L t then the ratio L t p p is close to a constant function. This observation allows us to estimate the ratios rigorously in practice. Remark 1.6. If we do not assume that f is Markov then the statement of the Theorem remains true, however, in this setting the construction of the functions p and q is more challenging since the eigenfunctions of L t might be non-analytic (but of bounded variation). As we will see later, in practical applications, the interval (α, β) λ( f, μ) depends on the quality of approximation of the leading eigenfunction of L t by polynomials p and q.
The next theorem guarantees that the previous theorem can be used to get bounds on the Lyapunov exponent which are arbitrary accurate. Note that Theorem 1.4 also holds under the weaker assumption that f : I → I is an expanding piecewise C 2 mixing Markov map of the interval, however, in this case it is much harder to compute the functions p and q which will give us good estimates on the Lyapunov exponent. In addition, it is convenient to assume analyticity in order to apply the following theorem.

Examples
In this section we will demonstrate how Theorem 1.4 can be used in practice. To this end we consider four examples, and compare the estimates we obtain for the Lyapunov exponents with previously known results. Theorem 1.4 allows us to obtain rigorous bounds using the built-in MaxValue routine in Mathematica, and the implementation is relatively straightforward. However, some care is required in choosing parameters during the construction of the functions p and q. In Sect. 4.2 we give more details on the practicalities of the implementation.
2.1. Classical example: the Lanford map. We will first illustrate our approach with the standard example of the Lanford map [14] (see also [3]). The original Lanford map and the graph of f L is shown in Fig. 1. Observe that the map is uniformly expanding with T (x) ≥ T (1) = 3 2 for 0 ≤ x ≤ 1. The inverse branches of f L are contractions given by The transfer operator L t therefore takes the form Due to simplicity of the formulae involved we shall attempt to obtain estimates on the Lyapunov exponents of particularly high accuracy, to demonstrate the power of our method. Namely, we shall choose ε = 10 −180 . Then we fix N = 400 and compute the nodes of the Chebyshev polynomial T 400 to 512-digits precision. Subsequently, we want to construct the functions p and q as polynomials of degree 399 using the spectral Chebyshev collocation method. We validate that the polynomials p and q are positive using the method described in Sect which are each presented to 130 significant figures. In particular, with these choices This value has previously been computed by Wormell [22] and her result agrees with the above. See also [10, §10.1]. In the present approach, the simplicity of the functions p and q is the source of the efficiency of the approach. In particular, this estimate was obtained in approximately 2 hours on a personal Macbook pro laptop with 2.8 GHz Quad-Core Intel Core i7 and 16 GB 2133 MHz LPDDR3 using Mathematica.
Remark 2.1. In addition to using the internal MaxValue function, whose code is not available to the public, we can apply a simple Monte-Carlo type method to numerically verify the value we obtained. More precisely, we generate N mc = 1000 pseudo-random points x j , j = 1, . . . , 1000 in the interval [0, 1] and evaluate both ratios at these points to get the values Then we compute a 1 := − log max j y + j and b 1 := log max j y − j . Repeating this procedure a total of t mc = 100 times, we obtain the values that are within a distance of 10 −345 from α and β, respectively. In particular, we see that our estimate agrees with the estimate given by the function MaxValue.
Observe that for the chosen parameter values f c ( and so the map f c is expanding. Then the inverse branches f 1 , f 2 : I → I are contractions defined by Following the formula (1) we obtain the associated transfer operator L t : We next want to compute the Lyapunov exponent λ(c) := λ( f c , μ c ) for forty equally spaced values c = c j = 0.001 + j−1 40 , with j = 1, . . . , 40 with an error of 10 −3 to sketch a graph of λ as a function of c. For this purpose we choose ε = 10 −3 and m = 60 and compute the nodes of the Chebyshev polynomial T 60 with accuracy of 256 digits. We then apply Theorem 1.4 and obtain lower and upper bounds for the Lyapunov exponent. The precision of the MaxValue routine in the computation was set to 128 digits.
Based on this calculation, we sketch the functions α(c) ε (dashed curve) and β(c) ε (solid curve) in Fig. 2. We see that for 0.01 < c < 0.96 the two curves are indistinguishable. However in the interval 0.96 < c < 0.99 they appear to be different. This reflects the fact that f c (1) → 1 as c → 1, i.e. the map f c has weak hyperbolicity for c close to 1. Uniform hyperbolicity is essential for Theorems 1.4 and 1.7 to be applicable.
In addition, we may also calculate the Lyapunov exponent for a selected parameter value c = 1 4 , for example, with high accuracy. To this end, we choose ε = 10 −180 and compute 300 zeros of the Chebyshev polynomial T 300 with accuracy of 400 digits. Then we apply the spectral collocation method to construct polynomials p and q of degree 299. As before, we verify this this functions are positive, and apply MaxValue with working precision 400.
We obtain the following values (for which we give 166 digits): Using the Monte-Carlo method with N mc = 1000 pseudo-random points in the interval [0, 1] and t mc = 100 samples, we can numerically check the output of the routine MaxValue. Namely, taking the maximum of the ratios L 1+ε p p and L 1−ε q q computed at 1000 different points a hundred times, we obtain the values which lie within the distance of 2.0 × 10 −345 from α or β, respectively.

A family of full branch piecewise Möbius maps.
We next consider a family of examples studied by Slipantschuk, Bandtlow and Just in [19] in connection with their study of relation between Lyapunov exponents and mixing rates.
Following [19], for − 1 When c = 0 this reduces to a piecewise linear "tent map" (see Fig. 3 for a plot). In the special case c = 0.11, of particular importance to the authors of [19], they assert that the Lyapunov exponent is λ( f 0.11 , μ 0.11 ) = 0.685 . . ., although the paper does not provide any details as to how this value was computed.
The inverse branches f 1 , and In particular, The associated transfer operator is given by .
We shall recover and improve the estimate of [19]. For this purpose, we choose m = 400 and compute Chebyshev nodes with accuracy of 600 digits. Then we choose ε = 10 −175 and apply Chebyshev collocation method to obtain two polynomials p and q of degree 399. We then verify that they are positive and evaluate x mod 1.
The inverse branches f 1 , f 2 : I → I are defined by and we can associate the transfer operators L t for t ∈ R according to (1). We next want to choose the following parameters for the computation. This is consistent with, and improves on, Froyland's estimate of λ( f, μ) = 0.64946. We see that in this case the accuracy is less than in other examples we have considered so far. One cause is the character of the inverse branches f 1 and f 2 : the formulae implies that providing we know the value of x ∈ (0, 1) with an error of 10 −k , we have the value of f 1 (x) and f 2 (x) with an error of 10 −k/6 .
Another source of complication is the diminished hyperbolicity. A straightforward calculation gives that f (x) ≥ f ( 1 2 ) = 2 + 2 2 3 − √ 6 = 1.1835 . . .. This relatively weak hyperbolicity also suggests an explanation for why the estimates are not as good as in the previous examples. In particular, the maximal eigenfunction for L t may be less regular (e.g., analytic on a relatively small Bernstein ellipse) which make the polynomial approximation used in §4.2 less effective.

Proof of Theorem 1.4
In order to explain the proof of Theorem 1.4 it helps to introduce the following famous function from thermodynamic formalism.

Pressure function.
We begin by introducing the following well known definition.

Definition 3.1.
To any map f ∈ N C(I ) we can associate the pressure function P : R → R defined by where the summation is taken over all fixed points for the iterates f n , n ≥ 1. . 6. The pressure function P(t) and the inequalities in part (4) of the Lemma This is one of many equivalent definitions of the pressure [21]. In fact, the pressure function is well defined in the setting of continuous maps. However, the method we give in the present manuscript requires the pressure function to be differentiable in a neighbourhood of 1. The usefulness of the pressure function to study the Lyapunov exponent is shown by the following simple lemma, the first three parts of which are well-known.

Lemma 3.2.
The pressure function has the following properties: 1. P(1) = 0; 2. P is an analytic convex function; 3. We can write λ( f, μ) = − d P(t) dt | t=1 , where μ is absolutely continuous f -invariant probability measure, and; 4. For any ε > 0 we can write Proof. The first three parts are essentially due to Ruelle [18] (see Corollary 5.27 and Exercise 5 (a) on p.99). The last observation follows easily from the convexity (see Fig. 6).
This leads to the following useful bound on the Lyapunov exponent.

Corollary 3.3. For any ε > 0 the following double inequality holds
Proof. This comes by substituting the identity in part 3 of Lemma 3.2 into the inequality in part 4. . This means that one has to estimate P(1 ± ε) with the double accuracy of the desired estimate for the Lyapunov exponent. Nevertheless it turns out that this approach is quite practical since it is quite easy to estimate the pressure to high precision.

Transfer operator for interval maps.
For definiteness, let us choose coordinates such that I = [−1, 1] (i.e., a = −1 and b = 1 in Definition 1.1) after a simple change of coordinates. In addition, we shall also assume for simplicity that the map f is fullbranch, i.e., f ((x k , x k+1 )) = (−1, 1), with inverse branches f k : (−1, 1) → (x k , x k+1 ), for k = 1, . . . , n, the general case being similar. The approach to estimating the pressure is based on its interpretation in terms of the family of transfer operators introduced in the Introduction. These operators act on the Banach space B of bounded analytic functions on domain U ρ ⊃ I enclosed by the Bernstein ellipse with the foci at −1 and 1 and given by We define the norm on B by f = sup z∈U | f (z)|. In particular, by choosing ρ sufficiently close to 1 we can assume that the inverse branches f j ( j = 1, . . . , n) of the map f ∈ N C(I ) have analytic extensions to U ρ . Therefore the maps f j : U ρ → U ρ are well defined and their derivatives are non-zero. Furthermore since all f j are contractions, we may assume without loss of generality that ∪ j f j U ρ ⊂ U ρ otherwise we can replace f by an iterate f n , where n is chosen sufficiently large that the condition holds with the inverse branches f j for f replaced by the inverse branches of f n . Of course, this would increase the number of contractions we are working with. We formally extend the definition of the transfer operators from (2) as follows: for x ∈ U ρ , and t ∈ R.
Remark 3.6. In this definition the functions | f j | t are real valued and real analytic on [−1, 1]. Thus by a slight abuse of notation we interpret | f j | t as being the complex analytic extension of these functions to U ρ .
We can estimate the pressure values P(t) using the maximal eigenvalue for the transfer operator L t . Lemma 3.7. Let f ∈ N C(I ) and let L t : B → B be the transfer operator defined by (6). Then

The spectral radius of L t is e P(t) . 2. The rest of the spectrum is contained in a disk of radius strictly smaller than e P(t) .
3. For any h ∈ B for which the restriction to I is strictly positive and any x ∈ I we have lim n→+∞ L n t h(x) 1 n = e P(t) .
Proof. Parts 1 and 2 are essentially due to Ruelle [18] (see Proposition 5.13 and 5.24). Part 3 follows directly from Part 1 and Part 2 and the classical spectral radius theorem (cf. [18], Proposition 5.13 and 5.14). Another exposition can be found in [16].
We can use this lemma to estimate the pressure values P(1 ± ε) in Corollary 3.3. In particular, in order to estimate e P(t) for t ∈ R we will use the following simple result.

Lemma 3.8.
Assume that for t ∈ R there exist a function p ∈ B, strictly positive on I , and a constant ρ ∈ R such that sup x∈I Proof. Since L t p(x) ≤ e ρ p(x) for all x ∈ I we can deduce that L n t p(x) ≤ e nρ p(x) for all n ≥ 1. Thus by Part 3 of Lemma 3.7 we have e P(t) = lim n→+∞ L n t p(x) 1 n ≤ e ρ for any x ∈ I .
We now combine Lemma 3.8 and Corollary 3.3 to prove Theorem 1.4.
Remark 3.9. The above arguments extend easily to all maps f ∈ N C(I ), not necessary full branch. In particular, instead of a single domain U ρ , we consider the disjoint union ] each bounded by a Bernstein ellipse with foci x n and x n+1 (k = 1, · · · , n). The Banach space is now taken to be B = ⊕ n k=1 B (k) where B (k) is the space of bounded analytic functions on U (k) ρ . Finally, we use the extension of (1) to define L t : B → B by where h = (h 1 , . . . , h n ) ∈ B. The argument then proceeds as above.

Practical realisation
We next want to explain how to apply Theorem 1.4 in practice. Below we give one way of constructing test functions p and q that we used in order to obtain estimates in the examples we considered. It is based on the spectral Chebyshev collocation method [20]. There are other methods one might consider, such as spline interpolation methods, proposed by Falk and Nussbaum [8], but this approach suffices for our needs.
A particular care is needed in implementation of the algorithm to avoid numerical errors. We implemented the method using the Mathematica built-in functions which use interval arithmetic to control possible errors. The method itself however can be implemented in any programming language where there libraries for interval arithmetic are readily available, e.g. C, Julia, or Python (in alphabetic order).

Constructing test functions p and q.
For notational simplicity, we will describe our construction in the special case of full branch maps. The generalization to the general case of Markov maps is fairly straightforward where I is replaced by the disjoint union of intervals. Definition 4.1. Let x 1 < x 2 < · · · < x n be a collection of n distinct real numbers. The Lagrange polynomials associated to {x j } n j=1 are the polynomials The Lagrange polynomials have the property that p j (x k ) = δ k j for all j = 1, . . . , m and k = 1, . . . , n. In a special case when the points {x j } n j=1 are the roots of a certain polynomial p n , they can be written as We assume below that I = [−1, 1], the general case being similar after a simple change of coordinates, and f ∈ N C(I ). Let us assume that one wishes to compute the Lyapunov exponent with an error of δ, in other words, we assume that one wishes to find an interval (λ 1 , λ 2 ) λ( f, μ) such that |λ 2 − λ 1 | ≤ δ. In order to define the functions p and q, we begin by choosing a natural number m = m(δ). Then we calculate numerically, with help of a computer, the following objects: Here we use the formula (8) to evaluate p j ( f i (x k )), using an inbuilt routine for evaluation of Chebyshev polynomials, which has guaranteed accuracy. 3. The leading left eigenvectors v t corresponding to the maximal eigenvalue of the matrices M t for t = 1 ± ε. They can be efficiently computed using the power method. 4. The polynomials p and q then given by linear combinations of Lagrange polynomials p j with coefficients coming from the eigenvectors: However, the formula (10) is prone to numerical errors. The polynomials p and q can also be written as a linear combination of Chebyshev polynomials, and this has the advantage of being more computationally stable than the more direct expansion in terms of Lagrange polynomials above. More precisely, the following expansion is well known.
where the coefficients are given in terms of the eigenvectors v (1−ε) and v (1+ε) : This allows us to evaluate p and q efficiently. 5. The supremums of the ratios L 1+ε p p and L 1−ε q q over the interval I is computed using internal routine MaxValue with working precision set to D = D(δ) digits.
In addition to exploiting the internal routine MaxValue we can also apply Monte Carlo type method in order to carry out a heuristic check on its output.
In the next subsection we show that in the setting of uniformly expanding piecewise analytic Markov maps the polynomials p and q satisfying the hypothesis of Theorem 1.4 can always be constructed, and moreover the conclusion of Theorem 1.7 holds.
where x j , j = 0, . . . , m are the roots of the Chebyshev polynomial T m+1 , i.e. the Chebyshev nodes and p j ( j = 0, · · · , m) are the Lagrange polynomials on I associated to x j , j = 0, . . . , m defined by (7). In particular, we see that the restriction π m | P m of π m to P m is the identity.
The transfer operator L t : B → B defined by (6) is compact, even nuclear, although this will not be needed. We require an estimate on the operator norm of the difference L t − L t π m defined by Lemma 4.2 (see [4], Theorem 3.3). Let f ∈ N C(I ) and let L t be the associated transfer operator defined by (6). Then there exists C > 0 and 0 < θ < 1 such we can bound that L t − L t π m ≤ C L t θ m for m ≥ 1.
This can be compared to [22, §2.2] where it is shown that L t π m −π m L t π m BV → 0 as m → ∞ using the bounded variation norm.

Remark 4.3.
Although we cannot expect that I − π m B→B → 0 as m → ∞, the composition with the operator L t allows the bound in Lemma 4.2 since for any function f analytic on U ρ for some ρ > 1, there exists ρ > ρ such that the image L t f is analytic on U ρ .
It follows from the properties of L t , Lemma 4.2 and classical analytic perturbation (see the book of Kato [12]) that we have the following: Lemma 4.4. Let f ∈ N C(I ) and let L t be the associated transfer operator defined by (6). Then for δ > 0 sufficiently small and m sufficiently large: 1. L t π m : B → B has a simple maximal eigenvalue λ m with |λ m − e P(t) | < δ; 2. The rest of the spectrum is contained in {z ∈ C : |z| < e P(t) − 2δ}; and 3. The corresponding eigenfunction h m for L t π m has a restriction to I which is strictly positive (i.e, h m (x) > 0 for x ∈ I ).
By perturbation theory the positivity of the restriction of the eigenfunction h associated to e P(t) for L t onto I implies the same for h m (since inf x∈I |h(x) − h m (x)| ≤ h − h m will be arbitrary small for m sufficiently large). The restriction π m L t | P m is a finite rank operator π m L t : P m → P m given by where x j , j = 0, . . . , m are the Chebyshev nodes introduced in Sect. 4.1. Observe that in the basis of Lagrange polynomials { p j } m j=0 given by (7) the operator π m L t is given by the matrix M t defined by (9).
In particular, maximal eigenvalue λ m for L t π m is also an eigenvalue for the matrix M t corresponding to the eigenvector π m (h m ) ∈ P m . This completes the proof of Theorem 1.7.
Given f ∈ N C(I ) and δ > 0 there exist N = N (δ) and D = D(δ) such that the polynomials p and q of degree m ≥ N with coefficients computed to D decimal places lead to estimates with 1 ε log β α < δ. In particular, the exponential convergence in Lemma 4.2 implies that N (δ) = O(log(δ/ε)) and D(δ) = O(− log δ, − log ε). Remark 4.5. (Heuristic estimates on the accuracy of approximation) For small |ε| 1 an O(ε 2 ) approximation to the eigenvalue e P(t) should give an O(ε) estimate on the Lyapunov exponent (since we divide out by ε in the formulae). Furthermore, this error is related to the (uniform) approximation error of the associated eigenfunction f by the interpolating polynomial based on m points, say, which is well known to be bounded by f k ∞ (k + 1)!. Even for very regular (e.g., analytic) functions f one only expects f k ∞ to tend to zero at best exponentially fast. Therefore, we might want ε ∼ √ 1/(k + 1)!. In particular, for degree k = 10 one gets ε = 5.10 −4 , for k = 20 one gets ε = 6.10 −10 , and for k = 100 one gets ε = 1.10 −79 .

Potential generalisations.
In principle the method in this paper can be used in more general settings. It should apply with a few modifications to Markov expanding maps in arbitrary dimensions. However, relaxing the Markov condition may present serious difficulties. Whereas the properties of the pressure in §3.1 persist in the case when log | f | is a function of bounded variation the construction of the test functions p and q in §4.1 in this space may be more problematic.
Similarly, relaxing the hyperbolicity assumption may prove challenging in this case because of potential lack of differentiability of P(t) at t = 0. This happens for instance in the case of phase transition associated to the Manneville-Pomeau and Liverani-Saussol-Vaienti maps [15].
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.