Exponentially Convergent Trapezoidal Rules to Approximate Fractional Powers of Operators

In this paper we are interested in the approximation of fractional powers of self-adjoint positive operators. Starting from the integral representation of the operators, we apply the trapezoidal rule combined with a double-exponential transform of the integrand function. In this work we show how to improve the existing error estimates for the scalar case and also extend the analysis to operators. We report some numerical experiments to show the reliability of the estimates obtained.


Introduction
In this work we are interested in the numerical approximation of L −α , α ∈ (0, 1).
Here L is a self-adjoint positive operator acting in an Hilbert space H in which the eigenfunctions of L form an orthonormal basis of H, so that L −α can be written through the spectral decomposition of L. In other words, for a given g ∈ H, we have (1) L −α g = +∞ j=1 µ −α j g, ϕ j ϕ j where µ j and ϕ j are the eigenvalues and the eigenfunctions of L, respectively, and •, • denotes the H-inner product.Throughout the paper we also assume σ(L) ⊆ [1, +∞), where σ(L) denotes the spectrum of L.
Applications of (1) include the numerical solution of fractional equations involving the anomalous diffusion, in which L is related to the Laplacian operator, and this is the main reason for which in recent years a lot of attention has been placed on the efficient approximation of fractional powers.Among the approaches recently introduced we quote here the methods based on the best uniform rational approximations of functions closely related to λ −α that have been studied in [6,7,8,9].Another class of methods relies on quadrature rules arising from the Dunford-Taylor integral representation of λ −α [1,2,3,4,5,17,18].Very recently, time stepping methods for a parabolic reformulation of fractional diffusion equations, proposed in [19], have been interpreted by Hofreither in [10] as rational approximations of λ −α .
In this work, starting from the integral representation (2) L −α = 2 sin(απ) π +∞ 0 t 2α−1 (I + t 2 L) −1 dt, α ∈ (0, 1), where I is the identity operator in H, we consider the trapezoidal rule applied to the single and the double-exponential transform of the integrand function.The former approach has been extensively studied in [5], where the authors also provide reliable error estimates.The rate of convergence has been shown to be of type where n is closely related to the number of nodes.Our aim here is just to review some theoretical aspects in order to refine the choice of the parameters that allow faster convergence, even if, still of type (3).As for the double-exponential transform, widely investigated in [12,13,14,15,16] for general scalar functions, in this work we show how to improve the existing error estimates for the function λ −α .We also extend the analysis to operators, showing that it is possible to reach a convergence rate of type exp −c n ln n .
While theoretically disadvantageous with respect to the single-exponential approach, we show that the double-exponential approach is actually faster at least for α ∈ [1/2, 1).
The paper is organized as follows.In Section 2 we make a short background concerning the trapezoidal rule with particular attention to functions that decay exponentially at infinity.In Section 3 we review the existing convergence analysis for the trapezoidal rule combined with a single-exponential transform and we refine the choice of the parameters that allow a faster convergence.Section 4 is devoted to the trapezoidal rule combined with a double-exponential transform.Here the convergence analysis is derived for the approximation of the scalar function λ −α and is then extended to the case of the operator L −α .Some concluding remarks are finally reported in Section 5.

A general convergence result for the trapezoidal rule
Given a generic continuous function f : R → R, in this section we make a short background concerning the trapezoidal approximation where h is a suitable positive value.Given M and N positive integers, we denote the truncated trapezoidal rule by For the error we have where The quantities E D and E T := E TL + E TR are referred to as the discretization error and the truncation error, respectively.
Definition 1. [11, Definition 2.12] For d > 0, let D d be the infinite strip domain of width 2d given by Let B(D d ) be the set of functions analytic in D d that satisfy and The next theorem gives an estimate for the discretization error of the trapezoidal rule when applied to functions in B(D d ). Then, Proof.By (5) we immediately have Using Theorem 1 we obtain (6).
The above result states that for functions that decay exponentially for x → ±∞ it may be possible to have exponential convergence after a proper selection of h.When working with the more general situation and, through the change of variable t = ψ(x), transform (7) to A suitable choice of the mapping ψ may allow to work with a function g ψ that fulfills the hypothesis of Theorem 2 so that I(g) can be evaluated with an error that decays exponentially.Since the aim of the paper is the computation of L −α with σ(L) ⊆ [1, +∞), for λ ≥ 1 we consider now the integral representation (2) and a change of variable t = ψ(x), ψ : (−∞, +∞) → (0, +∞), let (10) Let moreover be the truncated trapezoidal rule for the computation of λ −α , that is, for the computation of We denote the error by and for operator argument ) H→H .The remainder of the paper is devoted to the analysis of two special choices for ψ, the single-exponential (SE) and the double-exponential (DE) transforms.

The Single-Exponential transform
The SE transform is defined by (13) ψ SE (x) = exp(x), so that from ( 9) and ( 10) we get Since the poles of this function are given by we have that g λ,ψSE is analytic in that is, the strip domain with d = π/2, independently of α and λ.Now, in order to prove that g λ,ψSE belongs to B(D ψSE ) (see Definition 1), following the analysis given in [5] we first note that for η ∈ R, |η| < π/2 and λ ≥ 1, Therefore, ( 14) This implies that and also that Therefore, by (1) we can conclude that g λ,ψSE belongs to B(D ψSE ).By ( 4), for the discretization error we have , for h ≤ 2πd we obtain After choosing h, the contribute of the three exponentials can be equalized by taking M and N such that where ⌈•⌉ is the ceiling function.Denoting by n = M + N + 1 the total number of inversions we have that and therefore by (15) we obtain By (12), since L is assumed to be self-adjoint and σ(L) ⊆ [1, +∞), we have that and then, finally, taking d = π/2 we obtain The analysis just presented is almost identical to the one given in [5].Nevertheless in that paper the authors define d = π/4, while here we have shown that one can take d = π/2.This choice produces a remarkable speedup, as shown in Figure 1, where for different values of α we have considered the error in the computation of L −α for the artificial operator 10 16 ].

Double-Exponential transformation
The DE transform we use here is given by ( 18) We consider in (8) the change of variable The function involved in this case is and we employ the trapezoidal rule to compute The parameter τ needs to be selected in some way and the analysis is provided in Section 4.5.Its introduction is motivated by the fact that, when moving from λ to L, the method (the choice of M , N and h) and the error estimates have to be derived by working uniformly in the interval [1, +∞) containing σ(L).As in the SE case, the function g λ,ψDE (x) exhibits a fast decay for x → ±∞, but the definition of the strip of analyticity is now much more difficult to handle since everything now depends on λ and τ .In addition, we can bound the denominator using the results given in [13, p. 388], that is, From the above relations we find where Let x * be such that π sinh x * cos η = ln(τ /λ); we have

4.2.
Error estimate for the scalar case.The bound (20) implies that In addition, assuming d = d(λ, τ ) < π/2, it can be observed that Using (1), for the discretization error we have .
The remaining task is to estimate the truncation error.Using (20) we obtain Similarly, The above results are summarized as follows.
Proposition 1.Using the double-exponential transform, for the quadrature error it holds where ξ(d) is defined by (22).

Defining (26)
h = ln 4dn µ as in [13, Theorem 2.14], we first observe that (see ( 23)) Setting M = N = n, the choice of h as in (26) leads to a truncation error that decays faster than the discretization one, because for an arbitrary constant c (see ( 24)-( 25)) As consequence the idea is to assume the discretization error as estimator for the global quadrature error, that is, using (23) and ( 27), 28) is very similar to the one given in [13, Theorem 2.14], that reads where ν = max (α, 1 − α) and M = n, N = n − χ (or viceversa depending on α), where χ > 0 is defined in order to equalize the contribute of the truncation errors [13,Theorem 2.11].The important difference is given by the factor λ −α that replaces τ −α , and this is crucial to correctly handle the case of λ → +∞.In this situation the error of the trapezoidal rule goes 0 because g λ,ψDE (x) → 0 as λ → +∞ (see (19)).Anyway, as we shall see, d → 0 as λ → +∞, so that the exponential term itself is not able to reproduce this situation.An example is given in Figure 2 in which we consider λ = 10 12 and τ = 100.

4.3.
The poles of the integrand function.All the analysis presented so far is based on the assumption that the integrand function is analytic in the strip D d , for a certain d = d(λ, τ ).Therefore we have to study the poles of this function, that is, we have to study the equation We have By solving the above equation for each k, we obtain the complete set of poles.Assuming to work with the principal value of the logarithm and taking k = 0, we obtain the poles closest to the real axis x 0 and its conjugate x 0 , where In order to apply the bound on the strip we have to define The introduction of the factor r is necessary to avoid ξ(d) → +∞ as Im x 0 → π/2, which verifies for λ → τ (see ( 22)).

Asymptotic behaviors.
Setting we have and therefore we can write (31) as Assuming λ ≫ τ , that is, s ≫ 1, and using (33) we obtain Therefore, for λ ≫ τ , (34) Assume now λ = 1 and τ ≫ 1.By (31) we have we have that finally leads to (35) 4.4.The minimax problem.Let us define the function representing the (λ, τ )-dependent factor of the error estimate given by (28), that is, where K α is defined by (29).Since our aim is to work with a self-adjoint operator with spectrum contained in [1, +∞) the problem consists in defining properly the parameter τ.This can be done by solving As for the true error, experimentally one observes that τ must be taken much greater than 1, independently of α.Therefore, from now on the analysis will be based on the assumption τ ≫ 1. Regarding the function ϕ(λ, τ ), by taking d = d(λ, τ ) as in (32) and n sufficiently large, again, one experimentally observes that with respect to λ the function initially decreases, reaches a local minimum (for λ = τ in which d = rπ/2), then a local maximum (much greater than τ ), and finally goes to 0 for λ → +∞ (see Figure 3).In this view, denoting by λ the local maximum, for n sufficiently large the problem (36) reduces to the solution of (37) ϕ(1, τ ) = ϕ(λ, τ ).

4.4.1.
Evaluating the local maximum.Since 0 < d ≤ rπ/2, 0 < r < 1, we have where C is a constant depending on r.Therefore by (22), so that we neglect the contribute of this function in what follows.
Since the maximum is seen to be much larger than τ , we consider the approximation (34).Therefore we have to solve that, after some manipulation leads to We find the equation and since we finally have to solve For large n we have so that the solution of (39) can be approximated by For any given τ ≥ 1, it can be observed experimentally that λ * is a very good approximation of the local maximum (see Figure 3).We also remark that the assumption on n stated in (26), that leads to the error estimate (28), is automatically fulfilled for λ = λ * , at least for α not too small.Indeed, using (32) and (34) we first observe that (41) d(λ * , τ ) ≈ rπ Then by (38), using µ ≤ 1/2 and assuming for instance 0.9 < r < 1, after some simple computation we find µe 4d(λ * , τ ) ≈ µe 4πr Therefore the condition (26) holds true for n ≥ 1/(9α).
By defining (42) and (41) we have and hence, after some computation By (42) we have Joining the above approximations we finally obtain we have again ξ(d(1, τ )) ≈ 2 and therefore Using (38) we find ) 4.5.Approximating the optimal value for τ .We need to solve (37).Using the approximations (44) and (43) we impose Solving the above equation we find represents an approximate solution of (37).
In Figure 3 we plot the function ϕ(λ, τ * ) for λ ∈ [1, 10 20 ], in an example in which n = 40, α = 1/2, and τ * ∼ = 84.4defined by (45).Moreover we show the results of the approximations (40), ( 43) and (44), for τ = τ * .Clearly the ideal situation would be to have τ * such that ϕ(1, τ * ) = ϕ(λ, τ * ), but notwithstanding Remembering that In Figure 4 we show the behavior of the method for the computation of L −α , with L defined in (17), together with the estimate (47).For comparison, in the same pictures we also plot the error of the SE approach.As mentioned in the introduction, the DE approach appears to be faster for 1/2 ≤ α < 1.

Conclusions
In this work we have analyzed the behavior of the trapezoidal rule for the computation of L −α , in connection with the single and the double-exponential transformations.All the analysis has been based on the assumption of L unbounded, so that the results can be applied even to discrete operators, with spectrum arbitrarily large, without the need to know its amplitude, that is, the largest eigenvalue.In particular we have revised the analysis for the single-exponential transform and we have introduced new error estimates for the scalar and the operator case for the double-exponential transform.The sharp estimate obtained for the scalar case has been fundamental for the proper selection of the parameter τ that is necessary to obtain good results also for the operator case.

Figure 1 .
Figure 1.Error for the trapezoidal rule applied with the singleexponential transform with d = π/4 and d = π/2, and error estimate given by (16).

Figure 4 .
Figure 4. Error for the trapezoidal rule applied with the doubleexponential transform (error DE), with the single-exponential transform (error SE) and error estimate given by (47).