On Rational Krylov and Reduced Basis Methods for Fractional Diffusion
OOn Rational Krylov and Reduced BasisMethods for Fractional Diffusion
Tobias Danczul ∗ Clemens Hofreither † March 1, 2021
We establish an equivalence between two classes of methods for solving fractionaldiffusion problems, namely, Reduced Basis Methods (RBM) and Rational KrylovMethods (RKM). In particular, we demonstrate that several recently proposedRBMs for fractional diffusion can be interpreted as RKMs. This changed pointof view allows us to give convergence proofs for some methods where none werepreviously available.We also propose a new RKM for fractional diffusion problems with poles chosenusing the best rational approximation of the function x − s in the spectral interval ofthe spatial discretization matrix. We prove convergence rates for this method anddemonstrate numerically that it is competitive with or superior to many methodsfrom the reduced basis, rational Krylov, and direct rational approximation classes.We provide numerical tests for some elliptic fractional diffusion model problems. The area of numerical methods for diffusion problems with a fractional in space diffusionoperator has seen intensive development recently. In the present work, our interest lies inthe spectral definition of the fractional diffusion operator in bounded domains with homoge-neous Dirichlet boundary conditions. Numerical treatment of such problems by extension toa higher-dimensional, but local diffusion problem was proposed and analyzed in [35]. Severaldifferent methods [6, 2] make use of quadrature formulae for Dunford-Taylor or related integralrepresentations of the fractional power of the diffusion operator. A reformulation of the frac-tional problem as a pseudo-parabolic equation and solving it via a time-stepping scheme hasbeen proposed in [42, 43]. Methods based on best uniform rational approximation (BURA) ofcertain functions in the spectral domain were developed in [26, 27].It is remarkable that all the above-mentioned methods can viewed as rational approximationmethods, where the exact fractional power of the involved operator is approximated by a ∗ Institute for Analysis and Scientific Computing (TU Wien), Wiedner Hauptstrasse 8-10, 1040 Wien, Austria.
4s the standard (polynomial) Krylov space and we denote by P k ( L ) the space of polynomialsin L of degree at most k . Some fundamental properties of rational Krylov spaces are given inthe following lemma. Lemma 1.
The rational Krylov space Q k +1 ( L, b ) has the properties1. Q k +1 ( L, b ) = K k +1 ( L, q k ( L ) − b ) = P k ( L ) q k ( L ) − b ,2. if d j = −∞ for some j ∈ { , . . . , k } , then b ∈ Q k +1 ( L, b ) ,3. dim( Q k +1 ( L, b )) = dim( K k +1 ( L, b )) = min { k + 1 , M } , where M is the invariance index (see [21]) of the Krylov space K k +1 ( L, b ) .Proof. See [21, Lemma 4.2] for the case d j = −∞ for one j ∈ { , . . . , k } . If all poles are finite,the claims are validated analogously.The connection to the rational approximation methods sketched in Section 3.1 is easilyestablished. If we let p i ( z ) := (cid:89) j (cid:54) = id j (cid:54) = −∞ ( z − d j ) ∈ P k , i = 0 , . . . , k, and agree on the convention ( L − d j I ) − b := b for d j = −∞ , we see that the vectors introducedin (6) satisfy w j = ( L − d j I ) − b = p j ( L ) q k ( L ) − b and therefore w j ∈ Q k +1 due to the first property. Thus, if the vectors ( w , . . . , w k ) arelinearly independent, it follows from the third property that dim( Q k +1 ) = k + 1 ≤ M and Q k +1 = span { w , . . . , w k } . (8)In other words, if the chosen poles are pairwise distinct, the rational Krylov space is identicalto the space spanned by the solutions of the shifted problems (6).An orthonormal basis for the rational Krylov space Q k +1 is typically computed using the rational Arnoldi method [38, 21]. This algorithm requires as its input L , the right-handside b and the poles ( d j ) kj =1 . It entails solving shifted problems similar to (6) and thenorthonormalizes the resulting vectors, resulting in a matrix W ∈ R n × ( k +1) with orthonormalcolumns which spans Q k +1 . For a given scalar function f defined over Λ , an approximation u k +1 to the vector f ( L ) b within this subspace is then found via Rayleigh-Ritz extraction ,namely u k +1 := W f ( L k +1 ) W T b ∈ Q k +1 ( L, b ) , L k +1 := W T LW ∈ R ( k +1) × ( k +1) . (9)The matrix L k +1 is typically much smaller than L , and thus f ( L k +1 ) can be computed, e.g.,by diagonalization. Güttel [21] proves that this procedure is basis-independent, that is, u k +1 depends only on the space Q k +1 ( L, b ) , not the matrix W itself. Furthermore he points out thatRayleigh-Ritz extraction is equivalent to Galerkin projection in the special case f ( z ) = z − ,i.e., when solving a linear system with the matrix L .The rational Krylov space Q k +1 and its basis representation W depend on L and b , render-ing the above procedure nonlinear, but the same is true for standard Krylov space methods. In5ontrast, the direct rational approximation methods described in Section 3.1 are linear sincethey are given by u r = r ( L ) b with r fixed a priori . Nevertheless, we can also find a rationalrepresentation of this form for the rational Krylov method if we allow r to depend on theinput data, as the following result shows. Theorem 2.
The solution obtained by the rational Krylov method satisfies u k +1 = r ( L ) b , where r = p/q k and p ∈ P k is a polynomial such that r satisfies the interpolation conditions r ( µ j ) = f ( µ j ) , j = 1 , . . . , k + 1 , where the rational Ritz values ( µ j ) k +1 j =1 are the eigenvalues of L k +1 .Proof. See [21, Theorem 4.8] for the case d j = −∞ for one j = 0 , . . . , k . If all poles are finite,the proof follows analogously.Since the denominator q k of r is fixed, p is determined by the polynomial interpolationproblem p ( µ j ) = q k ( µ j ) f ( µ j ) for j = 1 , . . . , k + 1 . If the rational Ritz values µ j are pairwisedistinct, p is uniquely determined by these conditions.Clearly, the quality of the approximation u k +1 ∈ Q k +1 to f ( L ) b depends on the rationalKrylov space and therefore on a suitable choice of the poles ( d j ) kj =0 . However, the followingpowerful result shows that within this space, the approximation is quasi-optimal. Here wewrite W ( L ) for the numerical range of L which, in particular, contains the spectrum of L , and (cid:107)·(cid:107) refers to the Euclidean vector norm. Theorem 3.
Let W be an orthonormal basis of Q k +1 ( L, b ) and L k +1 := W T LW . Let f beanalytic in a neighborhood of W ( L ) and u k +1 = W f ( L k +1 ) W T b . For every set Σ ⊇ W ( L ) there holds (cid:107) f ( L ) b − u k +1 (cid:107) ≤ C (cid:107) b (cid:107) min p ∈P k (cid:107) f − p/q k (cid:107) L ∞ (Σ) with a constant C ≤ . . If L is self-adjoint, the result holds with C = 1 .Proof. See [21, Theorem 4.10] and [32, Proposition 3.2].Note the close relation of this result to Theorem 1: roughly, the error obtained using therational Krylov method with given poles ( d j ) is not much larger than the error obtained usingthe best possible rational approximation method with a rational function r having these samepoles. In this section, we show that several recently proposed schemes which are based on RBMsadmit a representation in the rational Krylov framework. To make matters precise, we considerthe discrete parametric reaction-diffusion equation ( tI + L ) w ( t ) = b (10)6or a prescribed right-hand side b ∈ R n and a parameter t ∈ R +0 := R +0 ∪{∞} that encodes thevariability of the problem. We set w ( ∞ ) := b by convention. The RBM seeks to approximatethe manifold of solutions ( w ( t )) t ∈ R +0 in the low-dimensional space V k +1 := V k +1 ( L, b ) := span { w ( t ) , . . . , w ( t k ) } , (11)where ≤ t < · · · < t k are particular parameters which we refer to as snapshots throughoutthis manuscript. The reduced basis analogon of the last claim in Lemma 1 is provided in [11,Lemma 3.5]: it states that dim( V k +1 ) = k + 1 if b is excited by sufficiently many eigenfunc-tions of L . The reduced basis surrogate w k +1 ( t ) ∈ V k +1 for w ( t ) is computed via Galerkinprojection, w k +1 ( t ) := V ( tI k +1 + L k +1 ) − V T b ∈ V k +1 ( L, b ) , L k +1 := V T LV ∈ R ( k +1) × ( k +1) , (12)where I k +1 ∈ R ( k +1) × ( k +1) denotes the identity matrix and V ∈ R n × ( k +1) a matrix whosecolumns form an orthonormal basis of V k +1 ( L, b ) . After an initial computational investment,the reduced space (11) allows us to evaluate the coefficient vector of w k +1 ( t ) in the basis { w ( t ) , . . . , w ( t k ) } for arbitrary t with complexity only depending on k . Due to (8), weimmediately obtain the following result. Lemma 2.
Let ( t j ) kj =0 ⊂ R +0 be pairwise distinct. Then the reduced space V k +1 ( L, b ) withsnapshots ( t j ) kj =0 and the rational Krylov space Q k +1 ( L, b ) with poles ( − t j ) kj =0 coincide. In the following two subsections, we study two classes of reduced basis methods which havebeen applied to the fractional diffusion problem, namely ones based on interpolation and onquadrature, and establish their connection to rational Krylov methods.
Two different model order reduction strategies have been recently proposed in [12] whichcouple interpolation theory with reduced basis technology. In line with [11], the (forward)fractional operator with positive exponent s ∈ (0 , is reinterpreted as a weighted integralover parametrized reaction-diffusion problems L s b = 2 sin( πs ) π (cid:90) ∞ t − s − ( b − v ( t )) dt, v ( t ) := ( I + t L ) − b . (13)Invoking v ( t ) = t − w ( t − ) and the substitution y = t − , where we rename the substitutedvariable t again, we observe that the integrand can be expressed in terms of the parameterfamily w ( t ) via L s b = sin( πs ) π (cid:90) ∞ t s ( t − b − w ( t )) dt. Based on a selection of snapshots ( t j ) kj =0 ⊂ R +0 , the integrand is approximated using a RBM,yielding u Int k +1 := sin( πs ) π (cid:90) ∞ t s ( t − b − w k +1 ( t )) dt. Our terminology differs from standard RBM notation, where the term snapshot is typically employed torefer to the discrete solution w ( t j ) instead of the parameter t j itself.
7s shown in [11, Theorem 4.3], the surrogate evaluates to u Int k +1 = V L sk +1 V T b (14)where V refers to a matrix of orthonormal basis vectors of V k +1 . In [11] it was proven thatthe scheme approximates L s b at exponential convergence rates. Motivated by these results,the authors of [12] proposed a version of (14) for the backward operator. They confirmedexperimentally that u RB k +1 := V f ( L k +1 ) V T b (15)converges exponentially to L − s b if f ( z ) = z − s , but no rigorous proof was known so far. Thefollowing theorem provides the essential tool to close this gap in the literature and allows usto establish a connection to RKMs. Theorem 4.
Let f be analytic in a neighbourhood of W ( L ) and ( t j ) kj =0 ⊂ R +0 pairwise distinct.Then the reduced basis approximation u RB k +1 with snapshots ( t j ) kj =0 ⊂ R +0 coincides with therational Krylov approximation (9) with poles ( − t j ) kj =0 .Proof. As pointed out by Güttel in [21, Lemma 3.3], the rational Krylov approximation isindependent of the choice of the particular basis. In view of (9) and (15), it thus suffices toverify that the corresponding search spaces Q k +1 ( L, b ) and V k +1 ( L, b ) coincide. This is truedue to Lemma 2.The second method presented in [12], also referred to as dual reduced basis approximation,follows a similar idea but is based on L − . Due to Theorem 2.2, 2.3, and Lemma 2.7 in [12],the negative fractional operator can be expressed as L − s b = 2 sin( πs ) π (cid:90) ∞ t − s − ( L − b − ( L − + t I ) − L − b ) dt. Utilizing y = t − while renaming the substituted variable t again, we obtain L − s b = sin( πs ) π (cid:90) ∞ t s ( t − L − b − ( tL − + I ) − L − b ) dt = sin( πs ) π (cid:90) ∞ t s ( t − L − b − w ( t )) dt. The latter is again approximated utilizing reduced basis technology with prescribed snapshots ( t j ) kj =0 ⊂ R +0 by means of u Dual k +1 := sin( πs ) π (cid:90) ∞ t s ( t − L − b − w k +1 ( t )) dt. (16)In [12, Theorem 3.4] it has been shown that (16) can be computed via u Dual k +1 = L − V L s − ∗ ,k +1 V T b , L ∗ ,k +1 := V T L − V. (17)If t j = ∞ for one j ∈ { , . . . , k } , the surrogate can be interpreted as a post-processed rationalKrylov approximation as follows. 8 heorem 5. Let s ∈ (0 , , f ( z ) = z s − , u Dual k +1 the dual reduced basis approximation (16) with snapshots ( t j ) kj =0 , W an orthonormal basis of Q k +1 ( L − , L − b ) with poles d j = − t − j , L k +1 = W T LW , and u k +1 = W f ( L ∗ ,k +1 ) W T L − b . Assume t j = ∞ for one j ∈ { , . . . , k } ,such that d j = − / ∞ = 0 . Then there holds u Dual k +1 = L − u k +1 . Proof.
We deduce V k +1 ( L, b ) = span { ( t I + L ) − b , . . . , ( t k I + L ) − b } = span { ( t L − + I ) − L − b , . . . , ( t k L − + I ) − L − b } = span { ( L − + t − I ) − L − b , . . . , ( L − + t − k I ) − L − b } , which affirms that the reduced space V k +1 ( L, b ) with snapshots ( t j ) kj =0 coincides with therational Krylov space Q k +1 ( L − , L − b ) with poles d j = − t − j , j = 0 , . . . , k . Let w.l.o.g. t = ∞ , or equivalently, d = 0 . Then, by definition, b ∈ Q k +1 ( L − , L − b ) such that b = W W T b = V V T b . Since u k +1 is independent of the particular basis, we have u k +1 = V L s − ∗ ,k +1 V T L − V V T b = V L s − ∗ ,k +1 L ∗ ,k +1 V T b = V L s − ∗ ,k +1 V T b = L u Dual k +1 . This yields u Dual k +1 = L − u k +1 as claimed. Remark 1.
From the rational Krylov perspective, a more natural approach to approximate L − s b would be to directly extract the surrogate from Q k +1 ( L − , L − b ) using the poles d = 0 , d j = − t − j for j ≥ , and f ( z ) = z s − , or equivalently, Q k +1 ( L − , b ) with d = −∞ , d j = − t − j for j ≥ , and f ( z ) = z s . In this way, the post-processing step, i.e., the finalmultiplication with L − , could be avoided. Based on the well-known Dunford-Taylor integral representation L − s = sin( πs ) π (cid:90) ∞ t − s ( t I + L ) − dt (18)for arbitrary positive definite operators L whose domain is contained in a Hilbert space,Bonito and Pasciak [6] presented an exponentially convergent sinc quadrature approximationfor L − s b . Using the substitution y = ln t , the method can be summarized as L − s b = sin( πs ) π (cid:90) ∞−∞ e (1 − s ) y w ( e y ) dy ≈ k ∗ sin( πs ) π N s (cid:88) j = − M s e (1 − s ) y j w ( e y j ) , (19)where k ∗ > is a parameter controlling the accuracy of the quadrature, y j = jk ∗ , and M s = (cid:24) π (1 − s ) k ∗ (cid:25) , N s = (cid:24) π sk ∗ (cid:25) . (20)As pointed out in [29], the method fits in the class of direct rational approximation techniquespresented in Section 3.1. In every quadrature node a parametric reaction-diffusion problem9f the form (10) must be approximated, which turns out to be the method’s bottleneck. Toalleviate the computational expenses, the authors of [7] propose to add an additional layer ofapproximation in the form of a RBM. Given a collection of snapshots ( t j ) kj =0 , the surrogateis defined by u Sinc k +1 := k ∗ sin( πs ) π N s min (cid:88) j = − M s max e (1 − s ) y j w k +1 ( e y j ) , where < s min ≤ s max < describes an interval for s ∈ [ s min , s max ] in which we wish toapproximate L − s b efficiently. Due to Theorems 2 and 4, we have u Sinc k +1 = k ∗ sin( πs ) π N s min (cid:88) j = − M s max e (1 − s ) y j r j ( L ) b , where r j = p j /q k , p j ∈ P k , interpolates f j ( z ) := ( e y j + z ) − in the eigenvalues of L k +1 = V T LV and q k is defined by (7) with d j = − t j . The authors of [14] pursue a similar approach. After algebraic manipulations of (18), aGauss-Laguerre quadrature is proposed to discretize the integral, which reads L − s b = sin( πs − ) πs − (cid:90) ∞ e − y w ( e − ys − ) dy + sin( πs + ) πs + (cid:90) ∞ e − y e ys + w ( e ys + ) dy ≈ sin( πs − ) πs − M − (cid:88) j =1 τ j, − w ( e − yj, − s − ) + sin( πs + ) πs + M + (cid:88) j =1 τ j, + e yj, + s + w ( e yj, + s + ) . Here, s ± := ± ( s − ) and ( τ j, ± , y j, ± ) M ± j =1 are the weights and nodes defining the quadraturerule, respectively. The choice M + = (cid:24) π sk ∗ (cid:25) M − = (cid:24) π − s ) k ∗ (cid:25) is suggested with parameter k ∗ > as in (20). A RBM strategy is applied to each of the twosums to reduce the computational costs. Based on two different distributions of snapshots ( t − j ) k − j =0 and ( t + j ) k + j =0 , k ± ∈ N , together with their respective reduced basis approximations w − k − +1 and w + k + +1 , the surrogate is defined by u Gauss k +1 := sin( πs − ) πs − M − (cid:88) j =1 τ j, − w − k − +1 ( e − yj, − s − ) + sin( πs + ) πs + M + (cid:88) j =1 τ j, + e yj, + s + w + k + +1 ( e yj, + s + ) , where k := k − + k + and t − = t +0 is assumed for simplicity. Interpreting the RBM in the aboveprocedure as a corresponding RKM, Theorems 2 and 4 yield the representation u Gauss k +1 = sin( πs − ) πs − M − (cid:88) j =1 τ j, − r − j ( L ) b + sin( πs + ) πs + M + (cid:88) j =1 τ j, + e yj, + s + r + j ( L ) b , where r ± j = p ± j /q ± k ± , p ± j ∈ P k , interpolates f ± j ( z ) := ( e ± yj, ± s ± + z ) − in the eigenvalues of thecorresponding projected operator, and q ± k ± is defined as in (7) with d ± j = − t ± j , respectively.10e conclude that each of the two quadrature schemes listed above admits a representationas a matrix-vector product of the form r ( L ) b , where r is a rational function determined by theunderlying RBM. Even though the RBM itself allows the interpretation as RKM, u Sinc k +1 and u Gauss k +1 cannot be extracted from a rational Krylov space via Rayleigh-Ritz extraction. Thelatter can be compensated by applying the RBM directly to the integrand of interest withoutdiscretization of the integral itself. To see this, let V k +1 ( L, b ) refer to an arbitrary reducedspace with basis V and L k +1 = V T LV . We apply (18) to L k +1 to deduce L − sk +1 = sin( πs ) π (cid:90) ∞ t − s ( tI k +1 + L k +1 ) − dt, where I k +1 ∈ R ( k +1) × ( k +1) denotes the identity matrix. Hence, V L − sk +1 V T b = sin( πs ) π (cid:90) ∞ t − s V ( tI k +1 + L k +1 ) − V T b dt = sin( πs ) π (cid:90) ∞ t − s w k +1 ( t ) dt. (21)Again, we make use of the transformation y = ln( t ) and invoke (15) to conclude u RB k +1 = sin( πs ) π (cid:90) ∞−∞ e (1 − s ) y w k +1 ( e y ) dy, (22)if f ( z ) = z − s in (15). Similarly, following the idea in [14, Lemma 3.1], one verifies that forthis particular choice of f u RB k +1 = sin( πs − ) πs − (cid:90) ∞ e − y w k +1 ( e − ys − ) dy + sin( πs + ) πs + (cid:90) ∞ e − y e ys + w k +1 ( e ys + ) dy, (23)which shows that the quadrature discretization can be omitted when using RBMs. Mostnotably, this allows us to spare the choice of the particular quadrature as well as the tuningof its associated parameters. Remark 2.
The presented classification of RBMs in fractional diffusion problems is far fromcomplete. E.g., in [5], the authors propose to apply a RBM to the extension framework [9].The elliptic problem on the artificially extended domain is approximated in a way that makesit amenable to reduced basis technology. It is yet unclear whether this approach allows theinterpretation as RKM and requires further investigation.
It is evident that the performance of all algorithms hinges on a good selection of snapshots(or poles) which determine the underlying matrix V . Weak greedy algorithms are amongthe most popular strategies to provide a good choice for ( t j ) kj =0 , see, e.g., [13]. Provideda computationally efficient error estimator, their aim is to iteratively add those parameterswhich seemingly yield the largest discrepancy to the exact solution. The authors of [7] and[14] advocate the implementation of such an algorithm combined with a residual-based errorestimator to extract the snapshots from the desired parameter domain Ξ ⊂ R . This approachcomes with the benefit of nested spaces, i.e., V k ⊂ V k +1 . A difficulty, however, is the fact thatthe efficient query of s (cid:55)→ u k +1 ≈ L − s b , u k +1 ∈ { u Sinc k +1 , u Gauss k +1 } , requires an s -independentselection of snapshots and is thus either limited to proper subsets s ∈ [ s min , s max ] ⊂ (0 , , ornecessitates Ξ to be unbounded. The latter is difficult to tackle numerically. Motivated byour analysis provided in Section 4, these inconveniences might be overcome if one omits thequadrature discretization as in (21) and chooses Ξ = Λ = [ λ , λ n ] .11 number of algorithms for the adaptive choice of poles in rational Krylov methods havebeen proposed as well [18, 16, 23]. They generally rely on the spectral rational interpolantdescribed in Theorem 2 and, unlike the greedy methods described in the previous paragraph,do not require an error estimator in the spatial domain, which typically makes their imple-mentation more efficient. To the best of our knowledge, the performance of these adaptivepole selection rules for fractional diffusion problems has not been studied.In contrast, the choice for ( t j ) kj =0 proposed in [11, 12] is given in closed form independently of b and s ∈ (0 , . It is based on the so-called Zolotarëv points and will be discussed in Section 4in more detail. Their computation only requires the knowledge of the extremal eigenvalues of L . The resulting spaces are not nested, that is, V k (cid:54)⊂ V k +1 . However, this drawback can beavoided by constructing the hierarchical sequence of sampling points proposed in [17], whichasymptotically yields the same convergence rates as the ones obtained by Zolotarëv.If the goal is to approximate L − s b for one fixed value of s ∈ (0 , , it might be more efficientto choose the low-dimensional space accordingly. Several RKMs have been proposed whichchoose poles in dependence of the fractional order s ; see, e.g., [3]. Theorem 3 makes it clearthat the question of optimal poles directly relates to the best uniform rational approximation(BURA) of f ( z ) = z − s in the spectral interval, which has been comprehensively studiedthroughout the last years in the framework of fractional diffusion [24, 25, 27, 29]. Untilrecently, numerical instabilities while computing the BURA were a major obstacle in theavailability of optimal poles. A remedy for this problem was recently proposed in the formof a novel algorithm for the fast and robust computation of BURAs using only standarddouble-precision arithmetic [30]. The quasi-optimality result Theorem 3 suggests the use of the poles of the best uniform rationalapproximation to f ( z ) = z − s in the spectral interval Λ as the poles of the rational Krylovmethod. Let r ∗ k be the rational function of degree at most k which minimizes the maximumerror, (cid:107) z − s − r ∗ k ( z ) (cid:107) L ∞ (Λ) = min p,q ∈P k (cid:107) z − s − p ( z ) q ( z ) (cid:107) L ∞ (Λ) , and ( d j ) kj =0 its poles. It is a classical result that r ∗ k exists and is unique (see, e.g., [4]). Wecan obtain results on its approximation quality from the work of Stahl [41], who has shownthat the best rational approximation ˜ r ∗ k to z s in [0 , satisfies the error estimate (cid:107) z s − ˜ r ∗ k ( z ) (cid:107) L ∞ [0 , ≤ C s exp( − π √ ks ) with a constant C s > which depends on s . Let λ > and λ n be the smallest and largesteigenvalues of L , respectively, and define, as in [27], r k ( z ) := λ − s ˜ r ∗ k ( λ z − ) . Then (cid:107) z − s − r k ( z ) (cid:107) L ∞ ( λ ,λ n ) = λ − s (cid:13)(cid:13)(cid:13)(cid:13)(cid:18) λ z (cid:19) s − ˜ r ∗ k (cid:18) λ z (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) L ∞ ( λ ,λ n ) ≤ C s λ − s exp( − π √ ks ) . One easily sees that r k satisfies the requisite equioscillation conditions and thus is the bestrational approximation to z − s in [ λ , ∞ ) . By definition, it follows that (cid:107) z − s − r ∗ k ( z ) (cid:107) L ∞ (Λ) ≤ (cid:107) z − s − r k ( z ) (cid:107) L ∞ (Λ) ≤ C s λ − s exp( − π √ ks ) . L − s b using the poles ( d j ) kj =0 of thebest rational approximation, Theorem 3 yields the error estimate (cid:107) L − s b − u k +1 (cid:107) ≤ CC s λ − s exp( − π √ ks ) (cid:107) b (cid:107) . (24)Typically, the smallest eigenvalue (which is closely related to the Poincaré constant of Ω )satisfies λ ≥ but is uniformly bounded with respect to the discretization parameters, andwe can thus ignore the dependence on λ .The above estimate makes use only of information on λ . If λ n (or a good bound for it)is known as well, we can directly use the best rational approximation of z − s on Λ = [ λ , λ n ] whose error is smaller than that of the best approximation on [ λ , ∞ ) and base a rationalKrylov method on its poles. To the best of our knowledge, the analytic behavior of theerror of the best rational approximation to z − s on a finite interval is not known. For thespecial case z (cid:55)→ z − / , the rational function which minimizes the relative maximum errorin a finite interval is explicitly known in terms of elliptic functions [10]. It does not seemthat this construction generalizes to different exponents, however. Error estimates for certain ( k − , k ) -Padé approximations to z − s are given in [2]; very roughly speaking, the authors giveestimates of the order ∼ ( k/s ) − s in the case of unbounded spectrum and ∼ exp( − k/ √ κ ) with the condition number κ = λ n /λ for bounded spectrum. The latter bound becomes pooras κ → ∞ and does not have the root-exponential bound by Stahl cited above as its limitingcase.Since best rational approximations are usually not known explicitly, they have to be approx-imated numerically. The most commonly used algorithm for this task is the rational Remezalgorithm, which is based on the equioscillation property of the best-approximation error, butis highly numerically unstable in its classical formulation. To mitigate this problem, extendedarithmetic precision has often been employed (see, e.g., [44]), which however has the draw-back of high computational effort due to the lack of hardware support for extended precisionarithmetic. Alternate approaches were recently proposed in [31, 19], where new formulationsof the Remez algorithm based on the so-called barycentric rational formula were given, sig-nificantly improving the numerical stability. In particular, the minimax routine in the latestversion of the Chebfun software package [15] is based on [19]. Unfortunately, this routine stilldoes not work well for functions of the type we are interested here. For this purpose, the sec-ond author has recently proposed a novel algorithm for best rational approximation based onbarycentric rational interpolation called BRASIL [30] which can compute the needed best ra-tional approximations rapidly, to very high degrees, and using only standard double-precisionarithmetic.
Each of the algorithms presented in the previous sections is directly related to rational Krylovor, in the broader sense, rational approximation methods. This changed point of view allowsus to use standard techniques from these fields to either provide novel convergence results orilluminate available proofs from a different perspective. We start with the following lemmawhich is instrumental in the analysis of one of the aforementioned reduced basis schemes.
Lemma 3.
Let t ∈ R +0 , w k +1 ( t ) the reduced basis approximation of w ( t ) with snapshots ( t j ) kj =0 ⊂ R +0 , and C as in Theorem 3. Then there holds for every set Σ := [ λ l , λ u ] ⊃ Λ , l > , (cid:107) w k +1 ( t ) − w ( t ) (cid:107) ≤ Ct + λ l (cid:107) Θ (cid:107) L ∞ (Σ) (cid:107) b (cid:107) , Θ( z ) := k (cid:89) j =0 t j (cid:54) = ∞ z − t j z + t j . Proof.
The proof follows the outline of [11, Lemma 5.12]. Due to Lemma 2 and Theorem 3we have (cid:107) w k +1 ( t ) − w ( t ) (cid:107) ≤ C (cid:107) b (cid:107) min p ∈P k (cid:107) ( t + z ) − − p ( z ) /q k ( z ) (cid:107) L ∞ (Σ) , where q k is defined as in (7) with d j = − t j . Assume for now t j = ∞ for some j ∈ { , . . . , k } ;without loss of generality, we choose j = 0 . The right-hand side can be bounded by min p ∈P k (cid:107) ( t + z ) − − p ( z ) /q k ( z ) (cid:107) L ∞ (Σ) ≤ (cid:107) ( t + z ) − − ¯ p ( z ) /q k ( z ) (cid:107) L ∞ (Σ) , where ¯ p ∈ P k − is uniquely defined by ¯ p ( t i ) = ( t + t i ) − q k ( t i ) , i = 1 , . . . , k. (25)Thanks to this interpolation property, we have that ¯ p/q k interpolates ( t + z ) − in t j , j =1 , . . . , k . Moreover, the difference of both functions is a rational function of degree ( k, k + 1) ,such that ( t + z ) − − ¯ p ( z ) q k ( z ) = c ( t ) (cid:81) kj =1 ( z − t j )( t + z ) (cid:81) kj =1 ( z + t j ) for some t -dependent constant c ( t ) ∈ R . Multiplying both sides with ( t + z ) and setting z = − t reveals c ( t ) = k (cid:89) j =1 t − t j t + t j with absolute value smaller than , which is why the claim holds if t = ∞ . Otherwise, wecan choose ¯ p ∈ P k according to ¯ p ( t i ) = ( t + t i ) − q k ( t i ) , i = 0 , . . . , k. Similarly to before, one confirms | ( z + t ) − − ¯ p ( z ) /q k ( z ) | ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:81) kj =0 ( z − t j )( t + z ) (cid:81) kj =0 ( z + t j ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , which proves the claim.Instead of applying Theorem 3 directly to f ( z ) = z − s , the authors of [11, 12] aim for aselection of poles according to a (uniform in the parameter t ) rational approximation of theresolvent function f ( z ) = ( t + z ) − . They propose to choose t = ∞ and t j = λ n dn (cid:18) k − j ) + 12 k K ( δ (cid:48) ) , δ (cid:48) (cid:19) , δ (cid:48) = (cid:112) − δ , δ = λ λ n , (26)14or j = 1 , . . . , k , where dn denotes the Jacobi elliptic function and K the elliptic integralof first kind; see [1, Section 16 & 17]. These snapshots are a scaled version of the so-calledZolotarëv points [46, 20, 36], which are known to minimize the maximal deviation of Θ( z ) in Lemma 3 over the spectral interval of L . As a direct consequence, we obtain exponentialconvergence for the reduced basis approximation (15) when using (26) in the case f ( z ) = z − s ,where no analytical result has been available yet. Theorem 6.
Let t = ∞ , ( t j ) kj =1 the scaled Zolotarëv points from (26) , C as in Theorem 3,and f ( z ) = z − s in (15) . Then there holds for all s ∈ (0 , (cid:107) u − u RB k +1 (cid:107) ≤ Cλ − s e − C ∗ k (cid:107) b (cid:107) , where C ∗ = πK ( µ )4 K ( µ ) , µ = (cid:32) − √ δ √ δ (cid:33) , µ = (cid:112) − µ , δ = λ λ n . Proof.
Due to (18) and (21), we have that u = sin( πs ) π (cid:90) ∞ t − s ( tI + L ) − b dt, u RB k +1 = sin( πs ) π (cid:90) ∞ t − s V ( tI k +1 + L k +1 ) − V T b dt. Lemma 3 yields (cid:107) u − u RB k +1 (cid:107) ≤ sin( πs ) π (cid:90) ∞ t − s (cid:107) ( tI + L ) − b − V ( tI k +1 + L k +1 ) − V T b (cid:107) dt = sin( πs ) π (cid:90) ∞ t − s (cid:107) w ( t ) − w k +1 ( t ) (cid:107) dt ≤ C (cid:107) Θ (cid:107) L ∞ (Λ) (cid:107) b (cid:107) sin( πs ) π (cid:90) ∞ t − s ( t + λ ) − dt = 2 Cλ − s (cid:107) Θ (cid:107) L ∞ (Λ) (cid:107) b (cid:107) , where the last equality follows from the scalar version of (18). With ( t j ) kj =1 as in (26), wemake use of the bound (cid:107) Θ (cid:107) L ∞ (Λ) ≤ e − C ∗ k , a well-known property of the Zolotarëv points [46, 20, 36], to complete the proof. This section is devoted to a numerical comparison of the algorithms discussed above, incorpo-rating efficiency, similarities, and performance with respect to several values of the parameter s . All methods are implemented in the open source finite element library Netgen/NGSolve [39, 40]. We consider the fractional diffusion model problem ( − ∆) s u = 1 on Ω , u = 0 on ∂ Ω , (27)on the unit square Ω := (0 , for s ∈ { . , . , . } . To discretize (27), we use a finiteelement space V h ⊂ H (Ω) constructed over a quasi-uniform triangulation of maximal mesh https://ngsolve.org/ h = 0 . and polynomial order p = 1 . The resulting extremal eigenvalues of L satisfy λ ≈ . and λ n ≈ . . For the sake of presentation, we consider only one algorithmfrom the class of quadrature-based RBMs, namely u Sinc k +1 , and omit the dual reduced basisapproximation. For a detailed investigation of u Gauss k +1 and u Dual k +1 we refer to [14] and [12],respectively. In favour of comparability, we choose t = − d = ∞ in all methods underconsideration. The remaining parameters are specified as follows.• For the quadrature-based RBM u Sinc k +1 , we set s min := 0 . , s max := 0 . , k ∗ := 0 . , M s max and N s min according to (20), and ( t j ) kj =1 ⊂ [ e − M s max k ∗ , e N s min k ∗ ] according to theresidual-based weak greedy algorithm proposed in [7].• For the rational Krylov approximation, we choose f ( z ) = z − s and investigate, in viewof Theorem 4, four different configurations of poles or rather snapshots. – For u Zolo k +1 , we choose the snapshots ( t j ) kj =1 as scaled Zolotarëv points (26). Thisconfiguration corresponds to one of the interpolation-based RBM proposed in [12].The evaluations of the Jacobi elliptic function and the elliptic integral is performedby means of the special function library provided by Scipy . – For u Greedy k +1 , we choose the same snapshots as for u Sinc k +1 . – For u Jac k +1 , we choose the poles ( d j ) kj =1 according to the distribution proposed by[3], which is based on a Gauss-Jacobi quadrature approximation for z − s . – For u Bura k +1 , we choose ( d j ) kj =1 according to the BURA poles of f ( z ) = z − s in [ λ , λ n ] obtained by the BRASIL algorithm [30], which is contained in the baryrat open-source Python package developed by the second author.• For the direct rational approximation method u r =: u Direct k +1 , presented in Section 3.1,we choose r ∈ R k,k as the best uniform rational approximation of f ( z ) = z − s on [ λ , λ n ] obtained by the BRASIL algorithm [30].The L -errors between the expensive discrete solution u in the sense of (4) to (27) and itslow-dimensional surrogates obtained by the six methods listed above are reported in Figures1, 2, and 3 for the values of s = 0 . , . , . , respectively.• In all cases, exponential convergence can be observed. For the RKM with BURA poles,the rate of convergence is significantly better than predicted by (24). One reason for thisis the fact that the error bound does not exploit the knowledge of the largest eigenvalueof L , but is based on a selection of poles on the unbounded interval [ λ , ∞ ) , as discussedin Section 3.4.• For s = 0 . , u Sinc k +1 , u Zolo k +1 , u Greedy k +1 , and u Jac k +1 satisfy exponential convergence of order O ( e − C ∗ k ) , with C ∗ as in Theorem 6. In the case of u Zolo k +1 and u Sinc k +1 , this is in accordancewith Theorem 6 and [7, Lemma 3.3], respectively. For s ∈ { . , . } , better convergencerates can be observed. A possible explanation for this is the particular choice of b .In [12] it has already been observed experimentally that for some configurations of theright-hand side, u Zolo k +1 converges with the predicted convergence rate irrespectively ofthe fractional order. https://docs.scipy.org/doc/scipy/reference/special.html https://github.com/c-f-h/baryrat
16 Those two methods which rely on the BURA, that is, u Bura k +1 and u Direct k +1 , providethe best approximation among all tested methods irrespectively of the fractional order.The observed rate of convergence is between O ( e − . C ∗ k ) and O ( e − . C ∗ k ) . In view ofTheorem 1 and 3, it is not surprising that these methods perform qualitatively similar.What stands out, however, is the observation that the quasi-optimal extraction of u Bura k +1 from Q k +1 ( L, b ) yields slightly better results than u Direct k +1 , which is based on the trueBURA of z − s . This is due to the fact that the former incorporates information aboutthe right-hand side, which allows the RKM to bias the surrogate towards the particularchoice of b . The discrepancy between u Bura k +1 and u Direct k +1 becomes more significant if wechoose b in (1) sufficiently smooth with homogeneous boundary conditions. In this case,the excitations b j := ( b , u j ) of b , with u j as in (5), decay quickly, such that u Direct k +1 ,which assumes a uniform distribution of excitations, requires substantially more linearsolves to reach a prescribed accuracy compared to its rational Krylov competitor.• The approximations u Sinc k +1 and u Greedy k +1 coincide for all values of k and s . The additionalquadrature discretization appears to have no impact on the quality of u Sinc k +1 at all. Apossible reason for this might be the fact that the sinc quadrature in (19) is (close to)exact if we replace w by w k +1 . Indeed, we observe numerically that for any λ ∈ [ λ , λ n ] (cid:90) ∞−∞ e (1 − s ) y r y ( λ ) dy ≈ k ∗ N s min (cid:88) j = − M s max e (1 − s ) y j r y j ( λ ) up to machine precision, where r y is the rational function from Theorem 2 with poles inthe negative snapshots that interpolates the resolvent function f y ( z ) := ( e y + z ) − in therational Ritz values of the underlying RKM. This exactness property of the quadraturecan also be observed for u Gauss k +1 . That is, if we greedily sample the snapshots ( t j ) kj =1 from a sufficiently large interval such that all three approximations u Greedy k +1 , u Sinc k +1 , and u Gauss k +1 are built upon that same search space, we observe numerically for, e.g., s = 0 . ,that u Greedy k +1 = u Sinc k +1 = u Gauss k +1 .• As discussed above, the methods based on the BURA provide the most accurate approxi-mation across all scenarios and are specifically tailored towards the fractional parameter s . If, however, solutions to (27) for several values of s are required, u Sinc , u Zolo ,and u Greedy outperform their competitors in terms of efficiency since they allow directquerying of the solution for arbitrary s after an initial offline computation phase. Acknowledgements
The first author has been funded by the Austrian Science Fund (FWF) through grant numberF 65 and W1245. The second author has been partially supported by the Austrian ScienceFund (FWF) grant P 33956-NBL.
References [1] M. Abramowitz and I. A. Stegun.
Handbook of mathematical functions with formulas,graphs, and mathematical tables , volume 55. National Bureau of Standards Applied Math-ematics Series, 1964. 17 ( e − C ∗ k )0 5 10 15 2010 − − − − − k e rr o r SincZoloGreedyJacBuraDirect
Figure 1: Discrete L -error (cid:107) u − u k +1 (cid:107) M , M as in (3), for s = 0 . , where u is the discrete high-fidelity solution of (27) and u k +1 the solution obtained by the respective method. O ( e − C ∗ k )0 5 10 15 2010 − − − − − k e rr o r SincZoloGreedyJacBuraDirect
Figure 2: Discrete L -error (cid:107) u − u k +1 (cid:107) M , M as in (3), for s = 0 . , where u is the discrete high-fidelity solution of (27) and u k +1 the solution obtained by the respective method.[2] L. Aceto and P. Novati. Rational approximations to fractional powers of self-adjoint positive operators. Numerische Mathematik , 143(1):1–16, 2019. doi: 10.1007/s00211-019-01048-4.[3] L. Aceto, D. Bertaccini, F. Durastante, and P. Novati. Rational Krylov methods forfunctions of matrices with applications to fractional partial differential equations.
Journal ( e − C ∗ k )0 5 10 15 2010 − − − − − k e rr o r SincZoloGreedyJacBuraDirect
Figure 3: Discrete L -error (cid:107) u − u k +1 (cid:107) M , M as in (3), for s = 0 . , where u is the discrete high-fidelity solution of (27) and u k +1 the solution obtained by the respective method. of Computational Physics , 396:470–482, 2019. doi: 10.1016/j.jcp.2019.07.009.[4] N.I. Achieser. Theory of Approximation . Dover books on advanced mathematics. DoverPublications, 1992. ISBN 9780486671291.[5] H. Antil, Y. Chen, and A. C. Narayan. Reduced basis methods for fractional Laplaceequations via extension.
SIAM J. Scientific Computing , 41:A3552–A3575, 2018.[6] A. Bonito and J. E. Pasciak. Numerical approximation of fractional powers of ellip-tic operators.
Mathematics of Computation , 84(295):2083–2110, 2015. doi: 10.1090/s0025-5718-2015-02937-8.[7] A. Bonito, D. Guignard, and A. R. Zhang. Reduced basis approximations of the solutionsto spectral fractional diffusion problems.
Journal of Numerical Mathematics , 28(3):147–160, 2020. doi: 10.1515/jnma-2019-0053.[8] O. Burkovska and M. Gunzburger. Affine approximation of parametrized kernels andmodel order reduction for nonlocal and fractional laplace models.
SIAM Journal onNumerical Analysis , 58(3):1469–1494, 2020. doi: 10.1137/19M124321X.[9] L. Caffarelli and L. Silvestre. An extension problem related to the fractional Laplacian.
Communications in Partial Differential Equations , 32(8):1245–1260, 2007. doi: 10.1080/03605300600987306.[10] T. Chiu, T. Hsieh, C. Huang, and T. Huang. Note on the Zolotarev optimal rationalapproximation for the overlap Dirac operator.
Physical Review D , 66(11), 2002. doi:10.1103/physrevd.66.114502. 1911] T. Danczul and J. Schöberl. A reduced basis method for fractional diffusion operators I,2019. URL https://arxiv.org/abs/1904.05599 .[12] T. Danczul and J. Schöberl. A reduced basis method for fractional diffusion operators II,2020. URL https://arxiv.org/abs/2005.03574 .[13] R. DeVore, G. Petrova, and P. Wojtaszczyk. Greedy algorithms for reduced bases inBanach spaces.
Constructive Approximation , 2013.[14] H. Dinh, H. Antil, Y. Chen, E. Cherkaev, and A. Narayan. Model reduction for fractionalelliptic problems using Kato’s formula. arXiv:1904.09332 [math.NA] , April 2019.[15] T. A Driscoll, N. Hale, and L. N. Trefethen.
Chebfun Guide . Pafnuty Publications, 2014.URL .[16] V. Druskin and V. Simoncini. Adaptive rational Krylov subspaces for large-scale dynam-ical systems.
Systems & Control Letters , 60(8):546–560, 2011. doi: 10.1016/j.sysconle.2011.04.013.[17] V. Druskin, L. Knizhnerman, and M. Zaslavsky. Solution of large scale evolutionary prob-lems using rational Krylov subspaces with optimized shifts.
SIAM Journal on ScientificComputing , 31(5):3760–3780, 2009. doi: 10.1137/080742403.[18] V. Druskin, C. Lieberman, and M. Zaslavsky. On adaptive choice of shifts in rationalKrylov subspace reduction of evolutionary problems.
SIAM Journal on Scientific Com-puting , 32(5):2485–2496, 2010. doi: 10.1137/090774082.[19] S.-I. Filip, Y. Nakatsukasa, L. N. Trefethen, and B. Beckermann. Rational minimax ap-proximation via adaptive barycentric representations.
SIAM Journal on Scientific Com-puting , 40(4):A2427–A2455, 2018. doi: 10.1137/17m1132409.[20] A. A. Gonchar. Zolotarëv problems connected with rational functions.
Mathematics ofthe USSR-Sbornik , 78 (120):640–654, 1969.[21] S. Güttel.
Rational Krylov Methods for Operator Functions . PhD thesis, TechnischeUniversität Bergakademie Freiberg, Germany, 2010. URL http://eprints.ma.man.ac.uk/2586/ . Dissertation available as MIMS Eprint 2017.39.[22] S. Güttel. Rational Krylov approximation of matrix functions: Numerical methods andoptimal pole selection.
GAMM-Mitteilungen , 36(1):8–31, 2013. doi: 10.1002/gamm.201310002.[23] S. Güttel and L. Knizhnerman. A black-box rational Arnoldi variant for Cauchy-Stieltjes matrix functions.
BIT Numerical Mathematics , 53(3):595–616, 2013. doi:10.1007/s10543-013-0420-x.[24] S. Harizanov, R. Lazarov, S. Margenov, P. Marinov, and Y. Vutov. Optimal solvers forlinear systems with fractional powers of sparse SPD matrices.
Numerical Linear Algebrawith Applications , 25(5):e2167, 2018. doi: 10.1002/nla.2167.[25] S. Harizanov, R. Lazarov, P. Marinov, S. Margenov, and J. Pasciak. Comparison analysison two numerical methods for fractional diffusion problems based on rational approxima-tions of t γ , ≤ t ≤ . arXiv e-prints , 2018.2026] S. Harizanov, R. Lazarov, S. Margenov, and P. Marinov. The best uniform rationalapproximation: Applications to solving equations involving fractional powers of ellipticoperators, 2019. URL https://arxiv.org/abs/1910.13865 . arXiv:1910.13865.[27] S. Harizanov, R. Lazarov, S. Margenov, P. Marinov, and J. Pasciak. Analysis of numeri-cal methods for spectral fractional elliptic equations based on the best uniform rationalapproximation. Journal of Computational Physics , 2020. doi: 10.1016/j.jcp.2020.109285.Available online.[28] J. S. Hesthaven, G. Rozza, and B. Stamm.
Certified Reduced Basis Methods forParametrized Partial Differential Equations . Springer, Switzerland, 1 edition, 2015. ISBN978-3-319-22469-5. doi: 10.1007/978-3-319-22470-1.[29] C. Hofreither. A unified view of some numerical methods for fractional diffusion.
Comput-ers & Mathematics with Applications , 80(2):332–350, 2020. doi: 10.1016/j.camwa.2019.07.025.[30] C. Hofreither. An algorithm for best rational approximation based on barycentric rationalinterpolation.
Numerical Algorithms , 2021. doi: 10.1007/s11075-020-01042-0.[31] A. C. Ionit , ă. Lagrange rational interpolation and its applications to approximation oflarge-scale dynamical systems . PhD thesis, Rice University, Houston, TY, 2013.[32] L. Knizhnerman, V. Druskin, and M. Zaslavsky. On optimal convergence rate of therational Krylov subspace reduction for electromagnetic problems in unbounded domains.
SIAM J. Numerical Analysis , 47:953–971, 01 2009. doi: 10.1137/080715159.[33] A. Lischke, G. Pang, M. Gulian, F. Song, C. Glusa, X. Zheng, Z. Mao, W. Cai, M. M.Meerschaert, M. Ainsworth, and G. E. Karniadakis. What is the fractional Laplacian?A comparative review with new results.
Journal of Computational Physics , 404:109009,March 2020. doi: 10.1016/j.jcp.2019.109009. URL https://doi.org/10.1016/j.jcp.2019.109009 .[34] I. Moret and P. Novati. Krylov subspace methods for functions of fractional differentialoperators.
Mathematics of Computation , 88(315):293–312, 2018. doi: 10.1090/mcom/3332.[35] R. H. Nochetto, E. Otárola, and A. J. Salgado. A PDE approach to fractional diffusionin general domains: A priori error analysis.
Foundations of Computational Mathematics ,15(3):733–791, 2015. ISSN 1615-3383. doi: 10.1007/s10208-014-9208-x.[36] I. V. Oseledets. Lower bounds for separable approximations of the Hilbert kernel.
Sbornik:Mathematics , 198(3):425–432, 2007.[37] A. Quarteroni, A. Manzoni, and F. Negri.
Reduced Basis Methods for Partial DifferentialEquations . Springer International Publishing, 2016.[38] A. Ruhe. Rational Krylov sequence methods for eigenvalue computation.
Linear Algebraand its Applications , 58:391–405, 1984. doi: 10.1016/0024-3795(84)90221-0.[39] J. Schöberl. Netgen an advancing front 2d/3d-mesh generator based on abstract rules.
Computing and Visualization in Science , 1:41–52, 1997. ISSN 1432-9360.2140] J. Schöberl. C++11 implementation of finite elements in ngsolve. 2014.[41] H. R. Stahl. Best uniform rational approximation of x α on [0 , . Acta Mathematica , 190(2):241–306, 2003. doi: 10.1007/bf02392691.[42] P. N. Vabishchevich. Numerically solving an equation for fractional powers of ellipticoperators.
Journal of Computational Physics , 282:289–302, 2015. doi: 10.1016/j.jcp.2014.11.022.[43] P. N. Vabishchevich. Numerical solving unsteady space-fractional problems with thesquare root of an elliptic operator.
Mathematical Modelling and Analysis , 21(2):220–238,2016. doi: 10.3846/13926292.2016.1147000.[44] R. S. Varga and A. J. Carpenter. Some numerical results on best uniform ratio-nal approximation of x α on [0 , . Numerical Algorithms , 2(2):171–185, 1992. doi:10.1007/bf02145384.[45] D. R. Witman, M. Gunzburger, and J. Peterson. Reduced-order modeling for nonlocaldiffusion problems.
International Journal for Numerical Methods in Fluids , 83(3):307–327, 2016. doi: 10.1002/fld.4269.[46] E. I. Zolotarëv. Collected works.