 # Approximation of a fractional power of an elliptic operator

Some mathematical models of applied problems lead to the need of solving boundary value problems with a fractional power of an elliptic operator. In a number of works, approximations of such a nonlocal operator are constructed on the basis of an integral representation with a singular integrand. In the present paper, new and more convenient integral representations are proposed for operators with fractional powers. Approximations are based on the classical quadrature formulas. The results of numerical experiments on the accuracy of quadrature formulas are presented. The proposed approximations are used for numerical solving a model two-dimensional boundary value problem for fractional diffusion.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Applied mathematical modeling of nonlocal processes Baleanu (2012); Uchaikin (2013) is often associated with solving boundary value problems that include fractional powers of elliptic operators Pozrikidis (2018). For example, suppose that in a bounded domain on the set of functions , there is defined the operator : . We seek the solution of the problem for the equation with the fractional power elliptic operator for a given .

Various numerical techniques, such as finite difference or finite volume methods, can be used to approximate problems with fractional power elliptic operators. The implementation of such algorithms require to calculate the action of a matrix (operator) function on a vector

, where is a given matrix (operator) and is a given vector. For instance, in order to calculate the solution of the discrete fractional power elliptic problem, we get , where .

There exist various approaches Higham (2008) how to calculate . They can be divided into two classes. In the first class, we construct one or another approximation of the operator function, i.e., . In this case, the solution of the problem is connected, e.g., with the solution of a number of standard problems . The second class of methods is based on the construction of approximate solutions .

For a functional approximation of fractional powers, best uniform rational approximation Stahl (2003) can be applied. The use of such technology for approximate solving multidimensional problems of fractional diffusion is discussed, for example, in Harizanov et al. (2018).

For a fractional power of a self-adjoint positive operator, the following integral representation holds Balakrishnan (1960); Birman and Solomjak (1987):

 A−α=sin(απ)π∫∞0θ−α(A+θI)−1dθ. (1)

To construct an approximation formula, the corresponding quadrature formulas Bonito and Pasciak (2015) are used. The fractional power of the operator is approximated by the sum of standard operators. To improve the accuracy of approximation, various approaches are applied.

The kernel of the representation (1) has a peculiarity if . To eliminate it, it is possible Bonito and Pasciak (2016) to introduce the new variable . Thus, from (1), we arrive at

 A−α=sin(απ)π∫∞−∞e(1−α)t(A+etI)−1dt.

To obtain the integral on a finite interval, in Frommer et al. (2014), there is employed the transformation

 θ=μ1−η1+η,μ>0.

From (1), we have

 A−α=2μ1−αsin(πα)π∫1−1(1−t)−α(1+t)α−1(μ(1−t)I+(1+t)A)−1dt. (2)

To approximate the right-hand side, we apply the Gauss-Jacobi quadrature formula Ralston and Rabinowitz (2001) with the weight .

In the present paper, a more general integral representation (in comparison with (1)) is used for fractional power operators, and possibilities of constructing new approximations are considered. Section 2 describes a problem for an equation with a fractional power of an elliptic operator of second order. Various types of integral representations of the fractional power operator are discussed in Section 3. In Section 4, we construct the corresponding quadrature formulas. Examples of solving a model two-dimensional problem with the fractional power of the Laplace operator are given in Section 5. The results of the work are summarized in Section 6.

## 2 Problem formulation

In a bounded polygonal domain with the Lipschitz continuous boundary , we search the solution of a problem with a fractional power of an elliptic operator. Here, we use the definition of a fractional power of an elliptic operator that relies on the spectral theory Birman and Solomjak (1987); Kwaśnicki (2017). Introduce the elliptic operator as

with coefficients , . The operator is defined on the set of functions that satisfy on the boundary the following conditions:

 v(x)=0,x∈∂Ω. (4)

In the Hilbert space , we define the scalar product and norm in the standard way:

 (v,w)=∫Ωv(x)w(x)dx,∥v∥=(v,v)1/2.

In the spectral problem

 Aφk=λkφk,x∈Ω,
 φk=0,x∈∂Ω,

we have

 0<λ1≤λ2≤...,

and the eigenfunctions

form a basis in . Therefore,

 v=∞∑k=1(v,φk)φk.

Let the operator be defined in the following domain:

 D(A)={v | v(x)∈L2(Ω), ∞∑k=0|(v,φk)|2λk<∞}.

Under these conditions and the operator is self-adjoint and positive definite:

 A=A∗≥νI,ν>0, (5)

where is the identity operator in . In applications, the value of is unknown (the spectral problem must be solved). Therefore, we suppose that in (5). Let us assume for the fractional power of the operator :

 Aαv=∞∑k=0(v,φk)λαkφk.

More general and mathematically complete definition of fractional powers of elliptic operators is given in the work Carracedo et al. (2001); Yagi (2009). The solution satisfies the equation

 Aαv=f (6)

under the restriction .

We consider the simplest case, where the computational domain is a rectangle

 Ω={x | x=(x1,x2), 0

To solve approximately the problem (6), we introduce in the domain a uniform grid

 ¯¯¯ω={x | x=(x1,x2),xn=inhn,in=0,1,...,Nn,Nnhn=ln, n=1,2},

where and is the set of interior nodes, whereas is the set of boundary nodes of the grid. For grid functions such that , we define the Hilbert space , where the scalar product and the norm are specified as follows:

For , the grid operator can be written as

 Au=−1h21a(x1+0.5h1,x2)(u(x1+h1,h2)−u(x))+1h21a(x1−0.5h1,x2)(u(x)−u(x1−h1,h2))−1h22a(x1,x2+0.5h2)(u(x1,x2+h2)−u(x))+1h22a(x1,x2−0.5h2)(u(x)−u(x1,x2−h2))+c(x)u(x),x∈ω.

For the above grid operators (see Samarskii (2001); Samarskii and Nikolaev (1989)), we have

 A=A∗≥δI,δ>0,

where is the grid identity operator. For problems with sufficiently smooth coefficients and the right-hand side, it approximates the differential operator with the truncation error , .

To handle the fractional power of the grid operator

, let us consider the eigenvalue problem

 Aψk=μkψk.

We have

 0<δ=μ1≤μ2≤...≤μK,K=(N1−1)(N2−1),

where eigenfunctions form a basis in . Therefore

 u=K∑k=1(u,ψk)ψk.

For the fractional power of the operator , we have

 Aαu=K∑k=1(u,ψk)μαkψk.

Using the above approximations, from (6), we arrive at the discrete problem

 u=A−αb. (7)

Taking into account the transformation

 A⟶1δA,b⟶1δαb,

we can assume, without loss of generality, that in (7), we have

 A≥I, (8)

i.e., .

## 3 Integral representation of the fractional power operator

For numerical solving the problem (7), (8), it is necessary to have an integral representation of the function for and . The starting point of our constructions is the formula (see the definite integral 3.141.4 in the book Gradshteyn and Ryzhik (2007))

 ∫∞0θμ−1(p+qθν)n+1dθ=q−μ/νpμ/ν−n−1γ(n,ν,μν), (9)

where

 0<μν

and is the gamma function. Using (9), we can construct different integral representations for the fractional power operator. Let us discuss some possibilities in this direction.

To employ (9), we introduce parameters associated with the operator : . In choosing these dependencies, we must bear in mind that it is necessary to calculate for integers . With this in mind, it is natural to put and consider two possibilities:

 p→A,q→I, (10)
 p→I,q→A. (11)

Using (10), from (9), we get the following integral representation of the fractional power operator:

 Aμ/ν−1=νπsin(πμν)∫∞0θμ−1(A+θνI)−1dθ. (12)

In the particular case of , from (12), it follows the Balakrishnan formula (1). The inverse transition also takes place, i.e., when using the new variable , from (1), we arrive at (12). For (12), from (9), we obtain the representation

 A−μ/ν=νπsin(πμν)∫∞0θμ−1(I+θνA)−1dθ. (13)

This representation corresponds to the use of the new variable in (1). Thus, for , the integral representations (12), (13) can be treated as variants of the Balakrishnan formula (1). Obviously, for , we do not have such the direct relation between (1) and (9), (10) or (9), (11).

In the two-parameter formulas (12) and (13), one parameter is free. If we put , then the integrand in the right-hand side of (12) and (13) have no singularities at the ends of the integration interval. Assuming , from (12), we get

 A−α=sin(πα)(1−α)π∫∞0(A+θ1/(1−α)I)−1dθ. (14)

In case of (13), we have and

 A−α=sin(πα)απ∫∞0(I+θ1/αA)−1dθ. (15)

The ability to use standard quadrature formulas for approximating the fractional power operator is provided by passing from (14) and (15) to a finite integral. Let us consider, e.g., the transition to a new integration variable on the interval such that with a numerical parameter . For (14), we get

 A−α=sin(πα)(1−α)π∫10(1−t)σα/(1−α)−1(1+(σ−1)t)((1−t)σ/(1−α)A+t1/(1−α)I)−1dt. (16)

The integrand in (16) has no singularities for . For , we have

 A−α=sin(πα)α(1−α)π∫10(α+(1−2α)t)((1−t)1/αA+t1/(1−α)I)−1dt. (17)

The same change of the integration variable in the representation (15) results in

 A−α=sin(πα)απ∫10(1−t)σ(1−α)/α−1(1+(σ−1)t)((1−t)σ/αI+t1/αA)−1dt. (18)

We use this representation for , where . For , we get

 A−α=sin(πα)α(1−α)π∫10(1−α+(2α−1)t)((1−t)1/(1−α)I+t1/αA)−1dt. (19)

Within the replacement of by , the formulas (17) and (19) coincide. The appropriate choice of allows us to provide the -th derivative of the integrand without any singularities. In particular, for (16) and (18), it is necessary to take and , respectively. This control of derivatives of the integrand is necessary when constructing quadrature formulas.

It seems reasonable to construct integral representations of the fractional power operator, when we select immediately the integration interval from 0 to 1. We can use, e.g., the definite integral 3.197.10 from the book Gradshteyn and Ryzhik (2007):

 ∫10θμ−1(1−θ)−μ1+pθdθ=(1+p)−μπsin(πμ),0<μ<1,p>−1.

Using the association , we arrive at the following integral representation:

 A−α=sin(πα)π∫10θα−1(1−θ)−α(I+θ(A−I))−1dt. (20)

Introducing the new variable , we see that (20) corresponds to the representation (2) with .

The following formula (see, for example, 3.198 in the book Gradshteyn and Ryzhik (2007)) provides more possibilities:

 ∫10θμ−1(1−θ)ν−1(aθ+b(1−θ+c)μ+νdθ=(a+c)−μ(b+c)−νB(μ,ν),
 a≥0,b≥0,c>0,μ>0,ν>0,

where is the beta function. This representation can be used for positive integers . However, for and even , the integrand has a singularity.

We also highlight the formula 3.234.2 in the book Gradshteyn and Ryzhik (2007):

 ∫10(θμ−11+pθ+θ−μp+θ)dθ=p−μπsin(πμ),0<μ<1,p>0.

It is also inconvenient for the direct application due to the singularity of the integrand for . Using the new variable , we arrive at

 ∫10(tσμ−11+ptσ+tσ(1−μ)−1p+tσ)dt=p−μπσsin(πμ).

The integrand has no sigularities if

 σμ−1≥0,σ(1−μ)−1≥0.

For instance, it is sufficient to take , where . With this in mind, we obtain a new integral representation of the fractional power operator:

 A−α=σsin(πα)π∫10(tσα−1(I+tσA)−1+tσ(1−α)−1(A+tσI)−1)dt. (21)

Selecting , we guarantee that the -th derivative of the integrand has no singularities.

The integral representations (17) (or (19)) and (21) are applied to approximate the fractional power operator using one or another quadrature formulas.

## 4 Numerical integration

The accuracy of quadrature formulas is investigated for the following integral representation (see (16)) of the function with :

 x−α=sin(πα)(1−α)π∫10(1−t)σα/(1−α)−1(1+(σ−1)t)((1−t)σ/(1−α)x+t1/(1−α))−1dt,σ=ϰ1−αα, (22)

with the parameter . When focusing on (21), we are interested in the integral representation

 x−α=σsin(πα)π∫10(tσα−1(1+tσx)−1+tσ(1−α)−1(x+tσ)−1)dt,σ=ϰmax(α−1,(1−α)−1). (23)

We apply the standard quadrature formulas, i.e., the rectangle and Simpson’s quadrature rules Davis and Philip (1984). The integration interval is divided into equal parts. When using the rectangle rule, the integrand is calculated at the midpoint of these parts, whereas Simpson’s rule employs their boundary nodes. For (22) and (23) with the rectangle rule, the approximating function is denoted by and , respectively. Here and

are explicitly highlighted as the key approximation parameters. The approximation error is estimated at

and is evaluated as follows:

 εe(x)=|re(x,α;M,ϰ)−x−α|,e=1,2.

For the rectangle rule, the approximation error of the function is shown in Figure 1 for and various values of the parameter . The calculations were performed at . The significant effect of the parameter is observed, when and . The accuracy of the rectangle quadrature formula is limited by the second derivative of the integrand. Therefore, we must take . To improve the approximation accuracy of the function for large values of , we need to take larger values of the parameter .

Numerical data for the maximal error

 ¯¯¯εe=max1≤x≤1020εe(x),e=1,2,

are presented in Table 1 for various numbers of integration nodes. It is easy to see that when the number of partitions increases, the convergence rate becomes close to theoretical estimates.

The approximation accuracy of the function for the representation (23) is shown in Figure 2. The corresponding numerical data are given in Table 2. We observe similar accuracy results for the representations (22) and (23). When focusing on approximations of the fractional power operator, the choice is in favor of using the representation (22), since both of the representations demonstrate comparable accuracy, whereas the computational complexity of the representation (23) is approximately two times higher.

We also present results on the approximation accuracy for Simpson’s quadrature rule. The approximation error of the function for the representation (22) is given in Table 3. The theoretical dependence of the accuracy on the number of partitions is observed for large and the parameter . Figure 3 shows the error for a fairly detailed partitioning with . It is easy to see that the accuracy increases with increasing .

## 5 Numerical solution of problems with the fractional power elliptic operator

Now we apply the quadrature formulas of the rectangle and Simpson’s rule considered above for (22) to the approximation of the integral representation of the fractional power grid elliptic operator (16) for solving the problem (7).

In the results presented below, the calculations were performed on the grid with . The peculiarities of boundary value problems with the fractional power operator were studied on the problem (3), (4), (6) in the unit square () with the right-hand side

 f(x)=sgn(x1−0.5) sgn(x2−0.5),sgn(x):=⎧⎪⎨⎪⎩−1,x<0,0,x=0,1,x>0. (24)

The numerical solution of the problem (7) with the right-hand side (24) is depicted in Figure 4 for different values of . For convenience of the comparison, there is shown the function

 y(x)=1maxu(x)u(x),x∈ω.

For the discontinuous right-hand side (24), we observe the formation of internal boundary layers with decreasing . Figure 4: Normalized solution of the problem (3), (4), (6), (24) on the grid N1=N2=256.

The accuracy of the approximate solution of the problem with the fractional power of the Laplace operator obtained using the rectangle rule is illustrated by the data in Table 4. Here the relative error is determined as follows:

 ε=∥w−u∥2∥u∥2,∥u∥22=∑x∈ωu2(x)h1h2,ε∞=∥w−u∥∞∥u∥∞,∥u∥∞=maxx∈ω|u(x)|,

where is the exact solution of the problem (7), (24) and is the approximate one. Similar data for Simpson’s quadrature formula are shown in Table 5. In these calculations, we employed in the representation (22) for the rectangle rule and for Simpson’s rule. We observe good accuracy of the approximate solution even for a small number of nodes of the quadrature formula.

For problems with the fractional power elliptic operator, we considered separately the dependence of the accuracy of the approximate solution on the smoothness of the exact solution. Above, we presented the results for the case of the discontinuous right-hand side (see (24)). For

 f(x)=x1x2, (25)

we have a boundary layer on a part of the boundary of the computational domain for small . The results on the accuracy of the approximate solution for this case are given in Table 6. Similar data for the smoother exact solution obtained for

 f(x)=x1(1−x1)x2(1−x2), (26)

are presented in Table 7. The numerical results indicate that the accuracy of the solution of the discrete problem (7) is practically independent from the smoothness of the solution of the continuous problem (6).

## 6 Conclusions

1. Different variants of the integral representation of the function are considered for approximating the fractional power operator. There are highlighted the integral representations, where the integrand has no singularities. By introducing new integration variables, integral representations are proposed for the function , where the derivatives of the integrand are bounded.

2. The rectangle and Simpson quadrature formulas were constructed to approximate the function for various values of . The accuracy of these quadrature formulas is numerically investigated depending on the key approximation parameters.

3. The model problem with the fractional power of the Laplace operator is considered in a rectangle. Its approximate solution is obtained on the basis of the approximation of the function applying the rectangle and Simpson’s quadrature rule. Calculations demonstrate high accuracy of the approximate solution if we use about one hundred quadrature nodes.

## Acknowledgements

The publication was financially supported by the Ministry of Education and Science of the Russian Federation (the Agreement number 02.a03.21.0008).

## References

• Balakrishnan (1960) Balakrishnan, A. V., 1960. Fractional powers of closed operators and the semigroups generated by them. Pacific Journal of Mathematics 10 (2), 419–437.
• Baleanu (2012) Baleanu, D., 2012. Fractional Calculus: Models and Numerical Methods. World Scientific, New York.
• Birman and Solomjak (1987) Birman, M. S., Solomjak, M. Z., 1987. Spectral theory of self-adjoint operators in Hilbert space. Kluwer academic publishers.
• Bonito and Pasciak (2015) Bonito, A., Pasciak, J., 2015. Numerical approximation of fractional powers of elliptic operators. Mathematics of Computation 84 (295), 2083–2110.
• Bonito and Pasciak (2016) Bonito, A., Pasciak, J. E., 2016. Numerical approximation of fractional powers of regularly accretive operators. IMA Journal of Numerical Analysis 37 (3), 1245–1273.
• Carracedo et al. (2001) Carracedo, C. M., Alix, M. S., Sanz, M., 2001. The Theory of Fractional Powers of Operators. Elsevier, Amsterdam.
• Davis and Philip (1984) Davis, P. J., Philip, R., 1984. Methods of Numerical Integration, 2nd Edition. Academic Press.
• Frommer et al. (2014) Frommer, A., Güttel, S., Schweitzer, M., 2014. Efficient and stable Arnoldi restarts for matrix functions based on quadrature. SIAM Journal on Matrix Analysis and Applications 35 (2), 661–683.
• Gradshteyn and Ryzhik (2007) Gradshteyn, I. S., Ryzhik, I. M., 2007. Table of Integrals, Series, and Products. Academic Press.
• Harizanov et al. (2018) Harizanov, S., Lazarov, R., Margenov, S., Marinov, P., Vutov, Y., 2018. Optimal solvers for linear systems with fractional powers of sparse spd matrices. Numerical Linear Algebra with Applications 25 (5), e2167.
• Higham (2008) Higham, N. J., 2008. Functions of Matrices: Theory and Computation. SIAM, Philadelphia.
• Kwaśnicki (2017) Kwaśnicki, M., 2017. Ten equivalent definitions of the fractional Laplace operator. Fractional Calculus and Applied Analysis 20 (1), 7–51.
• Pozrikidis (2018) Pozrikidis, C., 2018. The Fractional Laplacian. CRC Press.
• Ralston and Rabinowitz (2001) Ralston, A., Rabinowitz, P., 2001. A First Course in Numerical Analysis. Dover Publications.
• Samarskii (2001) Samarskii, A. A., 2001. The Theory of Difference Schemes. Marcel Dekker, New York.
• Samarskii and Nikolaev (1989) Samarskii, A. A., Nikolaev, E. S., 1989. Numerical methods for grid equations. Vol. I, II. Birkhauser Verlag, Basel.
• Stahl (2003) Stahl, H. R., 2003. Best uniform rational approximation of on [0, 1]. Acta Mathematica 190 (2), 241–306.
• Uchaikin (2013) Uchaikin, V. V., 2013. Fractional Derivatives for Physicists and Engineers. Higher Education Press.
• Yagi (2009) Yagi, A., 2009. Abstract Parabolic Evolution Equations and Their Applications. Springer, Berlin.