Rational approximations to fractional powers of self-adjoint positive operators
Lidia Aceto Paolo Novati

]Rational approximations to fractional powers of self-adjoint positive operators
This work was partially supported by GNCS-INdAM, University of Pisa (Grant
PRA_2017_05) and  FRA-University of Trieste. 

The authors are members of the INdAM research group GNCS.
[2010]47A58, 65F60, 65D32


Lidia Aceto
 Università di Pisa

 Dipartimento di Matematica, via F. Buonarroti 1/C - 56127 Pisa

Italy
lidia.aceto@unipi.it

Paolo Novati 
 Università di Trieste

Dipartimento di Matematica e Geoscienze, via Valerio 12/1,  34127 Trieste

Italy
novati@units.it
 
[
    Paolo Novati
    Received 2018; accepted 2018
================================




We investigate the rational approximation of fractional powers of unbounded
positive operators attainable with a specific integral representation of the
operator function. We provide accurate error bounds by exploiting classical
results in approximation theory involving Padé approximants. The
analysis improves some existing results and the numerical experiments proves
its accuracy.








§ INTRODUCTION
 



Let ℋ be a separable Hilbert space with inner product ⟨· ,·⟩ and corresponding norm ‖
x‖ =⟨ x,x⟩ ^1/2. Let ℒ be a
self-adjoint positive operator with spectrum σ (ℒ)⊆
c,+∞ ), c>0. Moreover, assume that ℒ has compact
inverse. This paper deals with the numerical approximation of ℒ
^-α, 0<α <1, that, in this setting, can be defined through
the spectral decomposition, i.e.,

    ℒ^-αu=∑_s=1^∞μ _s^-α⟨ u,φ_s⟩φ_s,

where {φ _s}_s=1^∞ is the orthonormal system of eigenfunctions of ℒ and {μ _s}_s=1^∞ is the corresponding sequence of positive real
eigenvalues (arranged in order of increasing magnitude and counted according
to their multiplicities). Clearly ℒ^-α is a self-adjoint
compact operator on ℋ. Since the function λ ^-α
is continuous in [c,+∞ ), we have that (see e.g. <cit.>)

    ‖ℒ^-α‖ =sup_λ∈σ (
    ℒ)|λ ^-α| =μ _1^-α.


An important and widely studied example comes from certain fractional models
involving the symmetric space fractional derivative (-Δ )^β /2
of order β (1<β≤ 2) <cit.>; in this situation the
fractional power is generally approximated through the approximation of (-Δ )^β /2-1 <cit.>.

A standard approach to approximate ℒ^-α is by means of ℒ_N^-α where ℒ_N is a finite dimensional
self-adjoint positive operator representing a discretization of ℒ.
Clearly, improving the sharpness of the discretization the typical
situation is that λ _min(ℒ_N)→ c and λ _max(ℒ_N)→ +∞ (λ _min(
ℒ_N) and λ _max(ℒ_N) denoting the
smallest and the largest eigenvalues of ℒ_N).

In this framework, in order to compute ℒ_N^-α
it is quite natural to employ rational forms. For instance, in <cit.>
some rational approximations are obtained by considering 
the best uniform rational approximation of λ^1-α and λ^α on the interval [0,1].
Beside, other well established techniques are the ones
based on existing integral representations of the Markov function λ
^-α and then on the use of suitable quadrature rules that finally
lead to rational approximations of the type

    ℒ_N^-α≈ℛ_k-1,k(ℒ_N),  ℛ_k-1,k(λ )=p_k-1(λ )/q_k(λ ),  
    p_k-1∈Π _k-1, q_k∈Π _k,

(see e.g.  <cit.>).  

In this setting, in <cit.> the rational forms arise from the use
of the Gauss-Jacobi rule for computing the integral representation (see <cit.>) 

    ℒ^-α=sin (απ )/(1-α )π∫_0^∞(ρ ^1/(1-α )I+ℒ)^-1dρ ,

after the change of variable

    ρ ^1/(1-α )=τ1-t/1+t,    τ >0.

Working in finite dimension, the asymptotically optimal choice of the
parameter τ, yields an error of type

    𝒪( exp( -4k√(λ _min(ℒ
    _N)/λ _max(ℒ_N))) ) ,

where k is the number of points of the quadrature rule, corresponding to a
ℛ_k-1,k(λ ) rational form. Of course √(λ
_min(ℒ_N)/λ _max(ℒ_N))→ 0improving the quality of the discretization so that (<ref>) becomes
meaningless whenever ℒ_N represents an arbitrarily sharp
discretization of ℒ.

The basic aim of the present work is to overcome this problem by working in
the infinite dimensional setting. Using the fact that the Gauss-Jacobi
quadrature on Markov functions is related to the Padé approximation, we
derive an expression for the truncation error λ ^-α-ℛ
_k-1,k(λ ):=λ ^-α-τ ^-αR_k-1,k(λ
/τ )  (here  R_k-1,k(λ
/τ) denotes the (k-1,k)-Padé approximant of (λ /τ )^-α),
that leads to an alternative definition of the parameter τ
independent of the discretization and, at the same time, ensuring an
asymptotically optimal rate of convergence. In particular, we are able to
show that the quadrature nodes for (<ref>) can be defined so that the
error for the computation of ℒ^-α decays approximatively
like

    ‖ℒ^-α-τ ^-αR_k-1,k( ℒ/τ) ‖≈sin (απ )c^-α( 2ke^1/2/α) ^-4α,

and therefore sublinearly. Qualitatively, a similar behavior can also been
observed by working with rational Krylov methods to approximate the action
of ℒ^-α (see <cit.>). The sublinearity appears when
considering unbounded spectra. Using the analysis for unbounded operators,
we also show how to improve quantitatively (<ref>) whenever we assume
to work with ℒ_N. The key point consists in taking τ in (
<ref>) dependent on k. 

We remark that all the theory here developed can be easily employed to
compute the action of the unbounded operator ℒ^1-α on a
vector f∈ D(ℒ) (D(ℒ) is the domain of ℒ
), that is ℒ^1-αf. This may occurs for instance when
solving equations involving the above mentioned fractional Laplacian. In
this situation, after evaluating g=ℒf, ℒ^1-αf
can be computed using our analysis on ℒ^-αg.
Nevertheless, the poles of the rational forms here derived can also be used
to compute ℒ^-αg by means of a rational Krylov method.


The paper is organized as follows. In Section <ref> we recall the basic
features of the Gauss-Jacobi based rational forms for computing (<ref>).
Section <ref> contains the error analysis and represents the main
contribute of this paper. In Section <ref> we revisit the error
analysis for the case of bounded spectra. Finally, in Section <ref> we
present some numerical experiments that validate the theoretical results.






§ BACKGROUND ON THE GAUSS-JACOBI APPROACH




Starting from the representation (<ref>), in order to approximate the
fractional Laplacian in <cit.> the authors consider the change of
variable (<ref>), that leads to

    ℒ^-α=2sin (απ )τ ^1-α/π∫_-1^1( 1-t) ^-α( 1+t) ^α
    -2( τ1-t/1+tI+ℒ) ^-1dt.

Using the k-point Gauss-Jacobi rule with respect to the weight function ω (t)=( 1-t) ^-α( 1+t) ^α -1 the above
integral is approximated by the rational form

    ℒ^-α≈∑_j=1^kγ _j(η _jI+ℒ
    )^-1:=τ ^-αR_k-1,k( ℒ/τ) ,

where the coefficients γ _j and η _j are given by

    γ _j=2sin (απ )τ ^1-α/πw_j/
    1+ϑ _j,    η _j=τ (1-ϑ _j)/
    1+ϑ _j;

here w_j and ϑ _j are, respectively, the weights and nodes
of the Gauss-Jacobi quadrature rule.

The choice of τ in (<ref>) is crucial for the quality of the
approximation attainable by (<ref>). As already mentioned in the
Introduction, working with bounded operators, it has been shown in <cit.>
that asymptotically, that is for k→ +∞, the optimal choice
is given by

    τ=√(λ _min(ℒ_N)λ _max(
    ℒ_N)).

With this choice and denoting by κ (ℒ_N) the spectral
condition number of ℒ_N and by ·_2 the induced Euclidean norm, we obtain

    ‖ℒ_N^-α-τ^-αR_k-1,k( ℒ_N/τ)
    ‖ _2≤ C( √(κ (ℒ_N))-1/√(κ (ℒ_N))+1) ^2k,

with C independent of k, which is a sharper version of (<ref>). We remark that τ is independent of k. In what follows we shall
follow a different strategy allowing a dependence on k (in any case the
coefficients γ _j and η _j completely change with k) but
at the same time a `mesh-independence', since we work with the unbounded
operator ℒ.






§ ERROR ANALYSIS
 



Working with the ratio λ /τ , where λ∈
c,+∞ ) and τ >0, as shown in <cit.> the k-point
Gauss-Jacobi quadrature given by (<ref>)-(<ref>) is such that
R_k-1,k( λ/τ) corresponds to the (k-1,k)-Padé approximant of (λ /τ )^-α centered at 1.

In this sense, defining

    z=1-λ/τ,

in what follows we focus the attention on the (k-1,k)-Padé
approximation

    (1-z)^-α≈ R_k-1,k(1-z).

Indicating the truncation error by

    E_k-1,k(1-z):=(1-z)^-α-R_k-1,k(1-z),

we have the following result.

For each integer k≥ 1, the exact representation of the truncation error defined  in (<ref>) is given by

    E_k-1,k(1-z)=Γ (k+1-α )Γ (k+1)/Γ (1-α
    )Γ (2k+1)_2F_1(k+1,k+α ;2k+1;z)/_2F_1(-k,k;α
    ;z^-1)( -z) ^k,

in which Γ denotes  the gamma function and _2F_1
the hypergeometric function.



Since <cit.>

    (1-z)^-α=_2F_1(1,α ;1;z),

the expression for the truncation error is obtained following the analysis given in <cit.>. 
 


For z<1, let v=1-2z^-1 and ξ be defined by

    v±( v^2-1) ^1/2=e^±ξ.

Then, for large values of k we have

    E_k-1,k(1-z)=4sin (απ )v-1/e^(2k+1)ξ(
    1+e^-ξ) ^-2α/( 1-e^-ξ) ^2(1-α )( 1+𝒪( 1/k) ) .




Since z=2/(1-v), using <cit.> we have that

    _2F_1(k+1,k+α ;2k+1;z)=4Γ (2k+1)Γ (1/2)/Γ
    (k+α )Γ (k+1-α )k^-1/2/(-z)^k+1
          × e^-(k+1)ξ(1-e^-ξ)^-3/2+α(1+e^-ξ)^-1/2-α( 1+𝒪( 1/k) ) ,


    _2F_1(-k,k;α ;z^-1)=Γ (k+1)Γ (α )/2Γ
    (1/2)Γ (k+α )k^-1/2
          × (1-e^-ξ)^1/2-α(1+e^-ξ)^α -1/2(
    e^kξ+e^± iπα ^-1/2e^-kξ) ( 1+𝒪( 1/k) ).


Plugging these relations in (<ref>) and using the identities Γ
(1/2)=√(π) and Γ (α )Γ (1-α )=π /sin (πα ), we find the result.  



For large values of k, the following representation for the truncation
error holds

    E_k-1,k( λ/τ) =2sin (απ)(
    λ/τ) ^-α[ λ ^1/2-τ
    ^1/2/λ ^1/2+τ ^1/2] ^2k( 1+𝒪(
    1/k) ) .




Using (<ref>), after some algebra we obtain

    2v-1/e^(2k+1)ξ( 1+e^-ξ) ^-2α/( 1-e^-ξ) ^2(1-α )=( v+1/v-1)
    ^-α1/[ v+( v^2-1) ^1/2] ^2k.

Since z=2/(1-v) and z=1-λ/τ we find

    v=λ +τ/λ -τ.

Substituting this expression in (<ref>) and then the result in (<ref>
), we easily obtain the statement. 


By (<ref>), (<ref>) and (<ref>) we have

    ‖ℒ^-α-τ ^-αR_k-1,k( ℒ/τ) ‖≤max_λ≥ c τ
    ^-α| E_k-1,k( λ/τ)
    | .

As consequence, a suitable value for τ can be found by working with (
<ref>). To this purpose, let us consider the function

    f(λ ,τ ):=( λ/τ) ^-α[
    λ ^1/2-τ ^1/2/λ ^1/2+τ ^1/2] ^2k,

which is the τ-dependent factor of (<ref>). We want to solve

    min_τ >0max_λ≥ cτ ^-αf(λ ,τ ).

For any fixed τ >0, f(λ ,τ )→ +∞ for λ→ 0^+, f(λ ,τ )→ 0 for λ→ +∞, f(λ ,τ )=0 for λ =τ (the
minimum) and, by solving ∂ f(λ ,τ )/∂λ=0, we find a maximum at

    λ=( k+√(k^2+1)) ^2/α ^2τ =s_k^24k^2/α ^2τ ,

where

    1<s_k^2=1+1/2k^2+𝒪(1/k^4).

Clearly λ>τ and hence

    max_λ≥ cτ ^-αf(λ ,τ )=max{τ
    ^-αf(c,τ ),τ ^-αf(λ,τ )} .

Setting

    φ _1( τ) :=τ ^-αf(c,τ ),  φ
    _2( τ) :=τ ^-αf(λ,τ ),

by (<ref>) we find

    φ _1( τ)    =   c^-α[ c^1/2-τ^1/2/c^1/2+τ ^1/2] ^2k,  
    φ _2( τ)    =   τ ^-αf( s_k^24k^2/α ^2τ ,τ)  
       =   τ^-α( s_k^24k^2/α ^2)^-α( 2ks_k-α/2ks_k+α) ^2k
       =   τ^-α( 4k^2e^2/α ^2) ^-α(1+𝒪(1/k^2)),

where the last equality follows from (<ref>) and by considering the
Taylor expansion around y=0 after setting s_k^2=1+y. Since φ
_2( τ) is monotone decreasing, whereas φ
_1( τ) is monotone increasing for τ >c, the solution
of (<ref>) is obtained by solving

    φ_1( τ) =φ_2( τ) τ >c.



Let τ ^∗ be the solution of (<ref>). Then

    τ ^∗≈τ _k:=c( α/2ke) ^2exp( 2W( 4k^2e/α ^2) ) ,

where W denotes the Lambert-W function.



Neglecting the factor (1+𝒪(1/k^2)) in (<ref>), equation (
<ref>) implies

    ( c/τ) ^-α[ τ ^1/2-c^1/2/τ ^1/2+c^1/2] ^2k=( 4k^2e^2/α ^2) ^-α.

Setting x=( c/τ) ^1/2<1 and

    a_k=α/2ke,

by (<ref>) we obtain

    x^-α/k( 1-x/1+x) =a_k^α/k.

Using the approximation

    1-x/1+x≈ e^-2x,

we solve

    e^-2x=(a_kx)^α/k.

Therefore

    -2x=α/kln( a_kx)

which implies

    2k/a_kα=1/a_kxln( 1/a_kx) .

Using the Lambert-W function, the solution for such equation is given by

    1/a_kx=exp( W( 2k/a_kα) ) .

Substituting x by ( c/τ) ^1/2 and using (<ref>)
we obtain the expression of τ _k. 


We remark that since for large z

    W( z) =ln z-ln( ln z) +𝒪(1),

we have (see (<ref>))

    τ _k=c4k^2/α ^2[ ln( 4k^2/α
    ^2e) ] ^-2(1+𝒪(1/k^2)).


By (<ref>) we thus obtain

    φ _2( τ _k)    =   τ _k^-αf(λ,τ _k) 
       =   c^-α( 2ke^1/2/α) ^-4α[
    2ln( 2k/α) +1] ^2α(1+𝒪
    (1/k^2)).

The above analysis yields the following result.


Let τ _k be defined according to (<ref>). Taking τ =τ _k in (
<ref>) we have

    ‖ℒ^-α-τ _k^-αR_k-1,k( 
    ℒ/τ _k) ‖   ≤   2sin (απ) c^-α( 
    2ke^1/2/α) ^-4α
          ×[ 2ln( 2k/α) +1] ^2α(1+𝒪(1/k^2)).




The statement immediately follows from (<ref>), (<ref>), and the
analysis just made. 



The factor c^-α in the bound (<ref>) reveals how the problem
becomes increasingly difficult if the spectrum is close to the branch point
of λ^-α. 







§ THE CASE OF BOUNDED OPERATORS




The theory just developed can be easily adapted to the case of bounded
operators ℒ_N with spectrum contained in [c,λ _N],
where λ _N=λ _max(ℒ_N). In this situation we
want to solve

    min_τ >0max_c≤λ≤λ _Nτ ^-αf(λ ,τ ).

Looking at (<ref>) we have λ=λ
(k)→ +∞ as k→ +∞. As a consequence for, λ≤λ _N (k small) the solution of (<ref>) remains the one approximated by (<ref>) and the bound (<ref>) is still valid. On the contrary, for λ>λ
_N (k large), the bound can be improved as follows.

Remembering the features of the function f(λ ,τ ) introduced in (
<ref>), we have that for λ>λ _N the solution
of (<ref>) is obtained by solving

    φ_1( τ) =φ_3( τ) τ >c,

where φ_1( τ) is defined in (<ref>) and

    φ_3( τ) :=τ ^-αf(λ_N,τ )
    =λ_N^-α[ λ_N^1/2-τ ^1/2/λ_N^1/2+τ ^1/2] ^2k.

It can be easily verified that the equation φ _1( τ)
=φ _3( τ) has in fact two solutions, one in the
interval (0,c) and the other in (c,λ _N). Anyway since φ
_3( τ) is monotone decreasing in [0,λ _N) we
have to look for the one in (c,λ _N) as stated in (<ref>).


Let τ̂^∗ be the solution of (<ref>). Then

    τ̂^∗≈τ _k:=( -αλ
    _N^1/2/8kln( λ _N/c) +√((αλ _N^1/2/8kln( λ _N/c))^2+( c λ _N) ^1/2)) ^2.




From (<ref>) we have

    c^-α[ τ ^1/2-c^1/2/τ ^1/2+c^1/2]
    ^2k=λ _N^-α[ λ _N^1/2-τ ^1/2/λ _N^1/2+τ ^1/2] ^2k.

Setting x=( c/τ) ^1/2<1 and y=( τ /λ
_N) ^1/2<1 by (<ref>) we obtain

    ( 1-x/1+x) =( λ _N/c) ^-α/2k( 1-y/1+y) .

Using (<ref>) we solve

    e^-2x=( λ _N/c) ^-α/2ke^-2y.

Therefore

    -2x=-α/2kln( λ _N/c) -2y

which implies

    x-y=α/4kln( λ _N/c) .

Substituting x by ( c/τ) ^1/2 and y by (
τ /λ _N) ^1/2 after some algebra we obtain

    τ +α/4kλ _N^1/2ln( λ _N/c) τ ^1/2-(c λ _N)^1/2=0.

Then, solving this equation and taking the positive solution, we obtain the
expression of τ _k.   


Observe that by (<ref>), for k→ +∞ we have

    ( τ _k/λ _N) ^1/2   =   -α/8kln( λ _N/c) +√(( α/8kln( λ _N/c) ) ^2+( c/λ_N) ^1/2)
       =   -α/8kln(
    λ _N/c) +( c/λ _N) ^1/4+𝒪( 1/k^2).

Finally, using (<ref>) and the above expression we obtain

    φ_3( τ_k)    =   λ_N^-α[ λ_N^1/2-τ_k^1/2/λ_N^1/2+τ_k^1/2] ^2k
       ≤   λ _N^-αexp( -4k( τ_k/λ_N) ^1/2)  
       =   λ_N^-αexp( -4k[ ( c/λ_N) ^1/4-α/8kln( λ_N/c) 
    ] ) ( 1+𝒪( 1/k) )  
       =   λ _N^-αexp( -4k( c/λ _N) ^1/4) exp( α/2ln( λ
    _N/c) ) ( 1+𝒪( 1/k)
    )  
       =   ( c λ _N) ^-α /2exp( -4k( c/λ _N) ^1/4) ( 1+𝒪( 1/k) ) .

The above analysis yields the following result.


Let k be such that for each k≥k we
have λ=λ(k)>λ _N. Then for
each k≥k, taking in (<ref>) τ =τ _k, where τ _k is given in (<ref>), the following bound holds

    ‖ℒ_N^-α-τ _k^-αR_k-1,k(
    ℒ_N/τ _k) ‖   ≤     2sin (απ )( c λ _N) ^-α /2
          ×exp( -4k( c/λ _N) ^1/4) ( 1+𝒪(1/k) ).



It is important to remark that, qualitatively, we have obtained the same
result of <cit.> and reported in (<ref>) following a completely
different approach. Nevertheless the analysis here presented is
quantitatively more accurate since it provides more information about the
constant preceeding the exponential factor.

Observe moreover that the analysis of this section may be particularly
useful when, in practical situation, one is forced to keep the
discretization quite coarse (so that k may be rather small) and
also to keep small the number of quadrature nodes k. In this case,
defining τ _k as in (<ref>) may provide results much better than
the one attainable with the asymptotically optimal choice τ̃=
√(λ _min(ℒ_N)λ _max(ℒ_N)).







§ NUMERICAL EXPERIMENTS




In this section we present the numerical results obtained by considering two simple cases of self-adjoint positive operators.
In particular, in the first example we try to simulate the behavior of an unbounded operator by working with a
diagonal matrix with a wide spectrum. In the second one we consider the
standard central difference discretization of the one dimensional Laplace
operator with Dirichlet boundary conditions.

We remark that in all the experiments the weights and nodes of the Gauss-Jacobi quadrature rule are computed  by
using the Matlab function  implemented in Chebfun by Hale and Townsend <cit.>.

 
  
 We define A=diag(1, 2, …,N) and ℒ_N=A^p so that σ (ℒ_N)⊆ 1,N^p]. Taking N=100 and p=4,
 in Figure <ref>, for α=0.25, 0.5, 0.75  the error (with respect to the Euclidean norm)
 and the error bound (<ref>) are plotted versus  k, the number of points of the used Gauss-Jacobi quadrature rule.


 
 In Figure <ref>, for  α=0.5 we plot the error obtained using  τ _k taken as in (<ref>) and  τ̃ as in (<ref>), changing the amplitude of the
 spectrum, that is, the value of p. In particular, we fix again N=100 and take p=2,3,4.
 
 
 
 The figure clearly shows the improvement attainable with τ _k for k
 small, and moreover the deterioration of the method for very large spectra
 when using τ̃.
 

 
We consider the linear operator ℒu=-u^'', u:[0, b]→ℝ, 
with Dirichlet boundary conditions u(0)=u(b)=0. It is known that ℒ has a point spectrum consisting
entirely of eigenvalues

    μ_s=π ^2s^2/b^2,    s=1,2,3,….

Using the standard central difference scheme on a uniform grid and setting b=1, in this example we work with the operator

    ℒ_N:=(N+1)^2tridiag(-1,2,-1)∈ℝ^N× N.

The eigenvalues are

    λ_j=4(N+1)^2sin ^2( jπ/2(N+1)) ,    
    j=1,2,… ,N,

so that σ (ℒ_N)⊆π^2,4(N+1)^2].  

The aim of this example is to show the improvement that can be obtained by
using the k-dependent parameter τ _k as in (<ref>) with respect
to the asymptotically optimal one τ̃. By choosing N=500, so that λ _N≈ 10^6, and α =0.5, in Figure <ref> the
errors are reported. In Figure <ref> we also plot the values of τ
_k. We remark that for other choice of α the results are
qualitatively identical.





Finally, still working with this example, we show the accuracy of the bound (
<ref>) for α=0.25, 0.5, 0.75.  The results are reported in Figure <ref>.










§ CONCLUSIONS
 



In this paper we have considered rational approximations of fractional powers of unbounded
positive operators obtained by exploiting the connection between Gauss-Jacobi quadrature on
Markov functions and Padé  approximants. Using classical results in approximation theory, we
have provided very sharp a priori estimates of the truncation errors that allow to properly
define the parameter τ. The numerical experiments confirm that such analysis improves some existing results.















 
















AMN Aceto L., Magherini C., Novati P. , On the construction and
properties of m-step methods for FDEs, SIAM J. Sci. Comput., 37, A653–A675 (2015)

AN Aceto L., Novati P., Rational approximation to the fractional
Laplacian operator in reaction-diffusion problems, SIAM J. Sci. Comput., 39,
A214–A228 (2017)

AN2 Aceto L., Novati P., Efficient implementation of rational
approximations to fractional differential operators, J. Sci. Comput., 76, 651–671 (2018)

Ba Baker G.A., Essentials of Padé approximants. Academic
Press, New York (1975)

RB Bhatia R., Matrix Analysis,  vol. 169 of Graduate Texts in
Mathematics.  Springer-Verlag, New York (1997)

Bo Bonito A., Pasciak J.E., Numerical approximation of fractional
powers of elliptic operators, Math. Comp., 84, 2083–2110 (2015)

E Elliot D., Truncation errors in Padé approximations to
certain functions: an alternative approach, Math. Comp.,
21, 398–406 (1967)

Erd  Erdélyi A., Magnus W., Oberhettinger F., Tricomi F.G.,
Higher trascendental functions, Vol. 1. McGraw-Hill, New York (1953) 

FI Fasi M., Iannazzo B., Computing the weighted geometric mean of
two large-scale matrices and its inverse times a vector, SIAM J. Matrix
Anal. Appl., 39, 178–203 (2018)

fgs Frommer A., Güttel S., Schweitzer M., Efficient and stable
Arnoldi restarts for matrix functions based on quadrature, SIAM J. Matrix
Anal. Appl., 35, 661–683 (2014)

HT Hale N., Townsend  A., Fast and accurate computation of
Gauss-Legendre and Gauss-Jacobi quadrature nodes and weights, SIAM J. Sci.
Comput., 35, A652–A672 (2013)

HetV Harizanov S., Lazarov  R., Margenov  S., Marinov P., Vutov Y., Optimal
solvers for linear systems with fractional powers of sparse SPD matrices, Numerical Linear Algebra with Applications,
25, 115–128  (2018)

HetP Harizanov S., Lazarov  R., Marinov P., Margenov  S., Pasciak J., Comparison
analysis on two numerical methods for fractional diffusion problems based on rational approximations of t^γ,
0 ≤ t ≤ 1, arXiv: 1805.00711v1 (2018)

Il05 Ilić M.,  Liu F.,  Turner I.,  Anh V., Numerical
approximation of a fractional-in-space diffusion equation I, Fract. Calc.
Appl. Anal., 8, 323–341 (2005) 

Ka Kwaśnicki M., Ten equivalent definitions of the fractional
Laplace operator, Fractional Calculus and Applied Analysis, 20, 7–51 (2017)

Mo Moret I., Novati P., Krylov subspace methods for functions of
fractional differential operators, Math. Comp., (2018) DOI:
10.1090/mcom/3332

No Novati P., Numerical approximation to the fractional derivative
operator, Numer. Math., 127, 539–566 (2014)

Ring Ringrose J. R., Compact non-self-adjoint operators. 
Van Nostrand Reinhold Company, London (1971)



