Rigorous Approximation of Diffusion Coefficients for Expanding Maps

We use Ulam’s method to provide rigorous approximation of diffusion coefficients for uniformly expanding maps. An algorithm is provided and its implementation is illustrated using Lanford’s map.


Introduction
The use of computers is essential for predicting and understanding the behaviour of many physical systems. Sensitive dependence on initial conditions is typical in many physical systems. This sensitivity problem raises nontrivial reliability and stability issues regarding any computational approach to such systems. Moreover, it strongly motivates the study of reliable computational methods for understanding statistical properties of physical systems.
In this note we consider the rigorous computation of diffusion coefficients in a class of systems where a central limit theorem holds. Such coefficients are focal in the study of limit theorems and fluctuations for dynamical systems (see [8,12,13,17,23,28] and references therein). Given a piecewise expanding map, an observable, and a pre-specified tolerance on error, we approximate in a certified way the diffusion coefficient up to the per-specified error (see Theorem 2.3).
In [32], following the approach of [18], a Fourier approximation scheme was used to estimate diffusion coefficients for expanding maps. The approach of [32] requires the map to have a Markov partition and to be piecewise analytic. Although the result of [32] provides an order of convergence, it does not compute the constant hiding in the rate of convergence. In our approach, we do not require the map to admit a Markov partition and we only assume it is piecewise C 2 . More importantly, our approximation is rigorous. To give the reader a flavour of what we mean by rigorous, we close this section by providing in part (b) of the following theorem a prototype result of this paper: 1

admits a unique absolutely continuous invariant measure ν and if ψ is a function of bounded variation the Central Limit Theorem holds:
In Sect. 2, we first introduce our framework and the assumptions on it. We then state the problem and introduce the method of approximation. The statement of the general results (Theorems 2.3, 2.5) and an application to expanding maps with a neutral fixed point are also included in Sect. 2. Section 3 contains the proofs and an algorithm. Section 4 contains an example, using Lanford's map, that illustrates the implementation of the algorithm of Sect. 3 and proves part (b) of Theorem 1.1.
operator (Perron-Frobenius) [4] associated with T , P : L 1 → L 1 is defined by duality: for f ∈ L 1 and g ∈ L ∞ Moreover, for f ∈ L 1 we have For f ∈ L 1 , we define We denote by BV the space of functions of bounded variation on I equipped with the norm || · || BV = V (·) + || · || 1 . Further, we introduce the mixed operator norm which will play a key role in our approximation:

Remark 2.1
It is important to remark that the constants α and B 0 in (A1) depend only on the map T and have explicit analytic expressions (see [22]).
The above assumptions imply that T admits a unique absolutely continuous invariant measure ν, such that dν dm := h ∈ BV . Moreover, the system (I, B, ν, T ) is mixing and it enjoys exponential decay of correlations for observables in BV (see [4] for a profound background on this topic). Footnote 3 continued using these assumptions, this work can be extended to the multidimensional case [24] by taking care of the dimension [25] and by working with appropriate observables since the space of functions of bounded variations in higher dimension is not contained in L ∞ . 4 It is well known that the systems under consideration satisfy a Lasota-Yorke inequality. What we are assuming in (A1) is that there is no constant in front of α. Such an assumption is satisfied for instance when inf x |T (x)| > 2 or when T is piecewise onto. When the original map T does not satisfy the assumption (A1), one can find an iterate of T where (A1) is satisfied, and then apply the results of this paper.

The Problem
Let ψ ∈ BV and define (2.1) Under our assumptions the limit in (2.1) exists (see [12]), and by using the summability of the correlation decay and the duality property of P, one can rewrite σ 2 as whereψ := ψ − μ and μ := I ψdν.
The number σ 2 is called the variance, or the diffusion coefficient, of n−1 i=0 ψ(T i x). In particular, for the systems under consideration, it is well known (see [12]) that the Central Limit Theorem holds: The goal of this paper is to provide an algorithm whose output approximates σ 2 with rigorous error bounds. The first step in our approach will be to discretize P as follows:

Ulam's Scheme
k=1 be a partition of [0, 1] into intervals of size λ(I k ) ≤ ε. Let B η be the σ -algebra generated by η and for f ∈ L 1 define the projection P ε , which is called Ulam's approximation of P, is finite rank operator which can be represented by a (row) stochastic matrix acting on vectors in R d(η) by left multiplication. Its entries are given by The following lemma collects well known results on P ε . See for instance [25] for proofs of (1)-(4) of the lemma, and [15,25] and references therein for statement (5) of the lemma.

Statement of the General Result
Defineψ ε := ψ − μ ε and μ ε : Remark 2.4 Theorem 2.3 says that given a pre-specified tolerance on error τ > 0, one finds l * > 0 and ε * > 0 so that σ 2 ε * ,l * approximates σ up to the pre-specified error τ . In Sect. 3.1 we provide an algorithm that can be implemented on a computer to find l * and ε * , and consequently σ 2 ε * ,l * .
To illustrate the issue of the rate of convergence and to elaborate on why we define the approximate diffusion by σ 2 ε,l as a truncated sum, let us define Remark 2.6 Note that σ 2 ε can be written as Since P ε has a matrix representation, and consequently (I − P ε ) −1 is a matrix, one may think that σ 2 ε provides a more sensible formula to approximate σ 2 than σ 2 ε,l . However, from the rigorous computational point of view one has to take into account the errors that arise at the computer level when estimating (I − P ε ) −1 . Indeed (I − P ε ) −1 can be computed rigorously on the computer by estimating it by a finite sum plus an error term coming from estimating the tail of the sum. 5 This is what we do in Theorem 2.3.

Remark 2.7
In [5] an example of a highly regular expanding map (piecewise affine) was presented where the exact rate of Ulam's method for approximating the invariant density h is ε ln ε −1 . In Theorem 2.5 the rate for approximating σ 2 is ε(ln ε −1 ) 2 . This is due to the fact that ||h − h ε || 1 is an essential part in estimating σ 2 and the extra ln ε −1 appears because of the infinite sum in the formula of σ 2 .
Remark 2.8 By using the representation (2.3) of σ 2 ε , it is obvious that the main task in the proof of Theorem 2.5 is to estimate Thus, it would be tempting to use estimate (9) in Theorem 1 of [19], which reads: where θ = ln(r/α) ln(1/α) , r ∈ (α, 1), and c 1 , c 2 are constants that dependent only on α, B 0 and r . On the one hand, this would lead to a shorter proof than the one we present in Sect. 3; however, estimate (2.4) would lead to a convergence rate of order ε θ , where 0 < θ < 1 which is slower than the rate obtained in Theorem 2.5. Naturally, this have led us to opt for using the proofs of Sect. 3.

Approximating the Diffusion Coefficient for Non-uniformly Expanding Maps
We now show that Theorem 2.3 can be used to approximate the diffusion coefficient for nonuniformly expanding maps. We restrict the presentation to the model that was popularized by Liverani-Saussol-Vaienti [27]. Such systems have attracted the attention of both mathematicians [27,37] and physicists because of their importance in the study of intermittent transition to turbulence [33]. Let where the parameter γ ∈ (0, 1). S has a neutral fixed point at x = 0. It is well known that S admits a unique absolutely continuous probability measureν, and the system enjoys polynomial decay of correlation for Hölder observables [37]. For γ ∈ (0, 1 2 ) it is known that the system satisfies the Central Limit Theorem. 6 To study such systems it is often useful to first induce S on a subset of I where the induced map T is uniformly expanding. In particular for the map (2.5), denoting its first branch by S 1 and the second one by S 2 , one can induce S on := [ 1 2 , 1]. For n ≥ 0 we define For n ≥ 1, we define Then we define the induced map T : → by where R Z n is the first return time of Z n to . For x ∈ , we denote by R(x) the first return time of x to . Let f be Hölder with I f dν = 0. Then diffusion coefficient of the system S can be written using the data of the induced map T (see [17]). In particular, for x ∈ , , the diffusion coefficient is given by where h is the unique invariant density of induced map T , P is the Perron-Frobenius operator associated with T , and m is normalized Lebesgue measure on . Thus, for ψ ∈ BV one can use, 7 Theorem 2.3 to approximate σ 2 .

Proofs and an Algorithm
We first prove two lemmas that will be used to prove Theorem 2.3. The explicit estimates of Lemma 3.2 below will also be used in Sect. 3.1 where we present our algorithm to rigorously estimate diffusion coefficients.
7 Although T has countable (infinite) number of branches, one can still implement the approximation on a computer. One way to do so is as follows: first one may perform an intermediate step by considering a map T identical to T on I \ H , such thatT has finite number of branches on I \ H while on H it has, say, one expanding linear branch, with m(H ) ≤ δ and δ τ is sufficiently small. The diffusion coefficients of T andT can be made arbitrarily close using the result of [20], and then one can apply Ulam's method and Theorem 2.3 toT .

Lemma 3.2 For any l
We know estimate (I I ): 3) Finally we estimate (I I I )

:= (I ) + (I I ) + (I I I ).
We start with (I I I ). Since Iψ hdm = 0, there exists a computable constant C * and a computable number 8 ρ * , where α < ρ * < 1, such that Consequently, Thus, choosing l * such that Fix l * as in (3.5). Now using Lemmas 2.2, 3.1 and 3.2, we can find ε * such that This completes the proof of the theorem.
(3) Find ε * = mesh(η) such that (4) Output σ 2 ε * ,l * := Iψ Remark 3.3 Note that the split of τ 2 between items (1) and (2) in Algorithm 3.1 to lead to an error of at most τ can be relaxed in following way. One can compute the error in item (1) to be at most τ k and in item (2) to be k−1 k τ for any integer k ≥ 2. We exploit this fact in the implementation in Sect. 4. We first get an estimate on (I I I ) and (I V ). There exists a computable constant C * and a computable number ρ * , where α < ρ * < 1, such that For (I I ), as in Lemma 3.2, in particular (3.4), and by using Lemma 2.2, we have

Implementation of the Algorithm and Estimating the Diffusion Coefficient for Lanford's Map
Let The map defined in (4.1) is known as Lanford's map [21]. In this section we let ψ = x 2 and compute the diffusion coefficient up to a pre-specified error τ = 0.035. A plot of T on [0, 1] and an approximation of its invariant density computed through Ulam's approximation are plotted in Fig. 1.

Rigorous Projections on the Ulam Basis
To compute the diffusion coefficient rigorously we have to compute rigorously the projection of an observable on the Ulam basis, i.e., given an observable φ in BV , and the projection ε we need to compute rigorously the coefficients {v 0 , . . . , v n } such that To do so, we will use rigorous integration through interval arithmetics, as explained in the book [35].
Once an observable is projected on the Ulam basis, many operations involved in the computation of the diffusion coefficient become componentwise operations on vectors; we explain this point in more details.
The first operation is the integral with respect to Lebesgue measure of an observable projected on the Ulam basis. This is given by the following formula: Suppose now we have computed an approximation h ε of the invariant density with respect to the partition, i.e., 1 0 h ε dx = 1. In the following we will denote its coefficients on the Ulam basis by {w 0 , . . . w n }. Note that the i-th component, w i , is the measure of I i with respect to the measure h ε dm.
The second operation we are interested in is the pointwise product of functions and the relation of the projection ε with this operation. We claim that: We will prove this by expressing the components of ε (φ · h ε ) as a function of the components {w 0 , . . . , w n } of h ε and the components {v 0 , . . . , v n } of ε φ. We claim that First of all recalling that χ 2 I i = χ I i and that χ I i · χ I j = 0 for i = j we have: On the right hand side, since h ε is constant on each I i and equal to w i , we have: These identities simplify the computations when dealing with the Ulam basis. It is worth noting that these identities imply that: Moreover, it is worth observing that, if P ε is the Ulam approximation and φ is an observable:

Item (1) in Algorithm 3.1
In this step, we find l * such that item (1) of Algorithm 3.1 is satisfied. In particular we want to find l * such that As explained in Remark 3.3, instead of verifying item (1) to be smaller than τ 2 , we verify that it is smaller than τ 256 . This will give us more room in verifying item (2) so that the sum of the errors from both items is smaller than τ . Since the system satisfies (A2), there exist 0 < ρ * < 1, and C * > 0 such that for any g ∈ BV 0 , and any k ∈ N, We want to find a 0 < ρ * < 1 and a C * > 0 so that (4.2) is satisfied. Once these two numbers are computed, we can easily find l * (see (3.5)) so that item (1) is satisfied. To compute ρ * and C * we follow [16] whose main idea is to build a system of iterated inequalities governed by a positive matrix M such that: where means component-wise inequalities, e.g. for vectors − → By using Lemma 2.2 and Appendix, we get that, if ||P n ε | BV 0 || 1 ≤ α 2 , the following inequalities are satisfied: (4.4) Using the inequalities above we have that: Following the ideas of [16] we have that where ρ * is the dominant eigenvalue of M and (a, b) is the corresponding left eigenvector. Thus, our main task now is to identify all the entries of the above matrix. The first one is M, a bound on the L 1 norm of the iterates of P and P ε ; by definition, we have that ||P n || ≤ 1 and ||P ε || 1 ≤ 1, therefore M = 1. The two constants α 2 and n 1 in M are two constants that give us an estimate of the speed at which P ε contracts the space BV 0 . Let P ε be the discretized Ulam operator and fix α 2 < 1; we want to find and n 1 ≥ 0 such that ∀v ∈ BV 0 with α 2 < 1. We follow the idea of [15] and use the computer to estimate n 1 with a rigorous computation; we refer to their paper for the algorithm used to certify n 1 and the corresponding numerical estimates and methods. Consequently, (4.3) is satisfied with n 1 = 28 , α ≤ Remark 4.1 Using equation (4.5) and supposing l * = k · n 1 we see that, for any ψ in BV 0 :

Item (2) of Algorithm 3.1
From now on l * is fixed and it is equal to 112. So far, we executed the first loop of the Algorithm 3.1; i.e., Remark 4.2 Note in our computation above we have obtained l * such that

Item (3) of Algorithm 3.1
In this step, we have to find ε * , a mesh size of the Ulam discretization, such that To bound this term we need a rigorous approximation of the T -invariant density h, in the L 1 -norm; we follow the ideas (and refer to the algorithm) of [15]. Set: (4.8) The following table contains, for different mesh sizes ε, error bounds for the terms in equation (4.7); in particular a bound on κ defined in (4.8):

Remark 4.3
The code implementing rigorous computation of diffusion coefficients for piecewise uniformly expanding maps is avalaible at the research section of the following personal page: http://www.im.ufrj.br/nisoli/

A Non Rigorous Verification
We also perform a non-rigorous experiment to compute σ 2 in the above example. Let F ζ be the set of floating point numbers in [0, 1] with ζ binary digits. Note that the system has high entropy, so we have to be careful in our computation and choose ζ big. Due to high expansion of the system, in few iterations the ergodic average along the simulated orbit may have little in common with the orbit of the real system. So, we have to do computations with a really high number of digits (ζ = 1024 binary digits).
Let {x 0 , . . . , x n−1 } be n random floating points in F l ; fix k and for each i = 0, . . . , n − 1 let Let μ be an approximation of the average of φ with respect to the invariant measure, obtained by integrating the observable using the approximation of the invariant density: The estimatorμ gives a non-rigorous estimate for the average of the observable with respect to the invariant measure, while the estimatorσ 2 is an estimator for the diffusion coefficient.
The table below shows the outcome of the experiment with n = 20, 000. In Fig. 2, a histogram plot of the distribution of A k (x i ) for k = 10, k = 200, n = 20, 000. In red we have the normal distribution with average μ and variance σ 2 ε * ,l * / √ k. The output of this non-rigourous experiment is in line with the output from our rigorous computation in Sect. 4.5. P n 1 = P n ε = 1.