1 Introduction
Various nonlocal process models baleanu2012fractional ; uchaikin , which are based on the use of fractional power of an elliptic operator Pozrikidis16 , are currently actively discussed. In the simplest case, the problem for the fractional Laplacian is formulated as follows. For a bounded domain on the set of functions we define the linear operator . We seek the solution of the problem for the equation with the fractional power elliptic operator for a given .
Various computational approaches for approximate solution of the problems with fractional power of an operator are used. As applied to problems of linear algebra, they are discussed in the work higham2008functions . For the problems with the fractional power of the elliptic operator under the consideration, special attention should be paid to the methods whose computational implementation is based on a solution of standard boundary value problems bonito2017 .
For a functional approximation of fractional powers, best uniform rational approximation stahl2003best can be applied. This approach for approximate solving of multidimensional problems of fractional diffusion is implemented in the works harizanov2018optimal ; HarizanovArxive2019 . The fractional power of the operator is approximated by the sum of standard operators. Similar approximations can be obtained using the appropriate integral representation of fractional powers of the operator and various quadrature formulas. Classical integral Balakrishnan representations balakrishnan1960fractional ; birman1987spectral are the basis of the works bonito2015numerical ; AcetoNovat ; Aceto2019 . Other integral representations are used in our paper vab201910838 .
A solution to a problem with a fractional power of an elliptic operator can be obtained from a solution to a problem of a larger dimension. In Caffarelli it is shown that the solution of the fractional Laplacian problem can be obtained as a solution to the elliptic problem on the semiinfinite cylinder domain. Computational algorithms based on such a relationship are discussed, for example, in nochetto2015pde .
For solving fractional power elliptic problems, we have proposed vabishchevich2014numerical a numerical algorithm based on a transition to a pseudoparabolic equation (Cauchy problem method). The computational algorithm is simple for practical use, robust, and applicable to solving a wide class of problems. To increase the accuracy of the approximate solution, two and threelevel schemes of higher order approximation on regular and irregular grids (see, for example, duan2018numerical ; ciegisvab201900201 ) are used.
Boundary value problems for the fractional power of an elliptic operator are often considered stinga2019user in terms of the methods of semigroups. In the works cusimano2018discretizations ; cusimano2018numerical , a singular integral expression is used to construct an approximate solution by solving the standard Cauchy problem for the parabolic equation. Special quadrature formulas that take into account the singularity and finite element approximations in space are used.
In the present paper, we also use the representation of the solution of the boundary value problem for the fractional power of a positive definite elliptic operator by solving IVP for the parabolic equation. The amplitude dependence of the solution of the auxiliary parabolic problem is highlighted. Approximations of the desired solution of fractional elliptical BVP is based on generalized GaussLaguerre formula davis2007methods , which most fully takes into account the singularity of the integral representation.
The paper is organized as follows. Section 2 describes a problem for an equation with a fractional power of an elliptic operator of second order. The problem is considered in a rectangle using a standard finitedifference approximation on a uniform grid. Singular integral expression through the solution of a standard Cauchy problem for the parabolic equation is discussed in Section 3. In Section 4, we construct the corresponding quadrature formulas. Numerical experiments for a model twodimensional 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 birman1987spectral . The elliptic operator is introduced as
(1) 
with coefficients , . The operator is defined on the set of functions that satisfy the following conditions on the boundary :
(2) 
In the Hilbert space , we define the scalar product and norm in a standard way:
In the spectral problem
we have
and the eigenfunctions
form a basis in . Therefore,Let the operator be defined in the following domain:
Under these conditions and the operator is selfadjoint and positive definite:
(3) 
where is the identity operator in . In applications, the value of is unknown (the spectral problem must be solved). Therefore, we suppose that in (3). Let us assume for the fractional power of the operator :
The solution satisfies the equation
(4) 
under the restriction .
We consider the simplest case where the computational domain is a rectangle
To solve the problem approximately (4), we introduce a uniform grid in the domain
where and is the set of interior nodes, whereas is the set of boundary nodes of the grid. For the grid functions such as , we define the Hilbert space , where the scalar product and the norm are specified as follows:
For , the grid operator can be written as
For the above grid operators (see Samarskii1989 ; SamarskiiNikolaev1978 ), we have
(5) 
where is the grid identity operator in . For problems with sufficiently smooth coefficients and the righthand 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
(6) 
We have
where eigenfunctions form a basis in . Therefore
For the fractional power of the operator , we have
Using the above approximations, from (4), we arrive at the discrete problem
(7) 
where, for example, . This problem is considered for (5), moreover, .
3 Representation by solving parabolic IVP
An approximate solution to the boundary value problem for the fractional degree of an elliptic operator (7) often based on the use of one or another integral representation. For example, in the works mentioned above bonito2015numerical ; AcetoNovat the Balakrishnan formula balakrishnan1960fractional ; birman1987spectral is used when
(8) 
The approximation of is based on the use of one or another quadrature formulas for the righthand side of (8).
Here, the possibilities of approximating both the fractional power of the operator and approximating the solution of the problem are realized.

Rational function approximation :
where are the weights, and are the nodes of the quadrature formula for (8).

An approximate solution to the problem (4):
The inverse operator of the fractional power elliptic problem is treated as a sum of inverse operators of classical elliptic operators and the required approximate solution of the problem is taken as the sum of the solutions of standard elliptical BVP.
The integral representation of the fractional power of the elliptic operator, which is fundamentally different from (8) used in works cusimano2018discretizations ; cusimano2018numerical . It is based on the method of semigroups, когда is defined through an integral formulation as follows
(9) 
where is gamma function. In this case, the solution to the problem (7) is
(10) 
Function is a solution of the following Cauchy problem
(11) 
(12) 
By the representation (10) we have
(13) 
For (1), (2) the problem (11), (12) is the standard Cauchy problem for a secondorder parabolic equation. Based on (13), the fractional elliptical BVP solution is represented by the IVP solution for the parabolic equation. Nonlocal stationary problem is associated with solving a nonstationary local problem and therefore, the considered approach is methodologically close to Cauchy problem method, which is proposed in vabishchevich2014numerical .
We note the most important features of the approximation of the fractional power of the operator and the problem solutions (7) based on the integral representation (9).

With quadrature formulas for (9) we come to the approximation with operator exponentials:
Here, as before, , are the weights and nodes of the corresponding quadrature formula.
Thus, the computational algorithm for solving the problem (7) can be built on a numerical solution of the Cauchy problem (11), (12), according to which we have the approximate solution with the selected quadrature formula (see (14)).
The integral representation (13) is singular, that is, the integrand has a singularity for for . For this reason, we can not rely on the accuracy of quadrature formulas, especially when . We extend the class of integral representations to eliminate the singularities of the integrand.
It is possible to transform (13) taking into account (11), (12). Based on the solutions of the Cauchy problem (11), (12) we have
Integrating by parts allows writing (13) as
(15) 
where . From (11), (12) we have the Cauchy problem for :
(16) 
(17) 
We arrive at a similar result based on the representation (9). Replacing by , where has a nonnegative integer, we get
(18) 
For the problem solution with a fractional power of the operator (5), it follows from (18) that
(19) 
We define function as a solution of the equation (11), supplemented by the initial condition
(20) 
To solve the problem (5) from (19) we obtain the representation
(21) 
For , the integral representation (21) through the solution of the problem (11), (20) corresponds to (16)–(18).
4 Numerical integration
We have two sources of errors in the approximate representation of the solution of the boundary value problem for the fractional power of a discrete elliptic operator through the solution of a parabolic problem. The first is related to the approximate calculation of the integral (21), the second is connected with a numerical solution of the Cauchy problem (11), (20). At this stage of our study, we assume that the problem (11), (20) is solved accurately and we will focus on the errors of approximate integration.
When using quadrature formulas (14) for an approximate solution of the problem (7) it is necessary to take into account the possible singularity of the integral representation: integration on the halfline and the features of the integrand. In the works of Cusimano et al. cusimano2018discretizations ; cusimano2018numerical integration (13) is held at a uniform partitioning of the final interval and the weights of the quadrature formula are associated with the function . Taking into account (6) to solve the problem (11), (20) we have the representation
From this perspective, from (21) we obtain
(22) 
Introducing a new integration variable from (22) we get
In view of this, it is necessary to construct quadrature formulas for the integrals
(23) 
in which the integrand is
(24) 
Concerning (22) we have
When calculating the integrals (23) with the weight function , we focus on the generalized GaussLaguerre quadrature formula davis2007methods . In this case
(25) 
where the nodes of the quadrature formula are the zeros of the generalized Laguerre polynomial , and weights are
The accuracy of the approximate calculation of
is estimated at
by the value of the relative errorFor convenience of the analysis of the calculated data, the dependence on and is highlighted separately.
The accuracy of the generalized GaussLaguerre quadrature formula for calculating the integral with a different number of nodes is given in Table 1. The accuracy decreases with decreasing fractional power and monotonically increases with increasing . A fundamental increase in accuracy is achieved by choosing .
0  5.241730e01  2.617133e01  1.032809e01  4.721032e02  3.089846e02  

1  1.827771e02  1.264074e02  7.111007e03  4.169156e03  3.078424e03  
25  2  2.090947e03  1.583269e03  1.016776e03  6.683762e04  5.249533e04 
3  3.846239e04  3.069550e04  2.136824e04  1.511406e04  1.236642e04  
4  9.538633e05  7.894328e05  5.818673e05  4.339988e05  3.660028e05  
0  4.891202e01  2.202112e01  7.321202e02  2.822765e02  1.669060e02  
1  8.628584e03  5.396535e03  2.569714e03  1.276760e03  8.540722e04  
50  2  5.088380e04  3.493823e04  1.907784e04  1.067459e04  7.615334e05 
3  4.911145e05  3.563271e05  2.118067e05  1.280562e05  9.540530e06  
4  6.498686e06  4.901612e06  3.097237e06  1.982550e06  1.526013e06  
0  4.563886e01  1.852328e01  5.183319e02  1.683123e02  8.980301e03  
1  4.049587e03  2.286550e03  9.186559e04  3.853388e04  2.329379e04  
100  2  1.212817e04  7.529116e05  3.476968e05  1.646274e05  1.062789e05 
3  6.001109e06  3.942100e06  1.986278e06  1.018559e06  6.876088e07  
4  4.108546e07  2.809418e07  1.508100e07  8.206209e08  5.731183e08 
5 The numerical solution of fractional elliptical BVP
We will illustrate the possibilities of an approximate solution to the problem (1), (2), (4) based on the solution of the Cauchy problem (11), (20) using the quadrature formula (25). The model problem for the Laplace operator is considered when
For eigenfunctions and eigenvalues (6) we have (see, e.g., SamarskiiNikolaev1978 ):
Direct calculations yield
The model problem (4) is solved numerically with two righthand sides: , when
where is the sign function. We consider the problem with both a smooth righthand side when choosing and a nonsmooth one choosing .
The numerical solution of the problem (7) with the righthand side is depicted in Figure 1 and Figure 2 for different values of . For convenience of the comparison, the following function is shown
For the discontinuous righthand side (, Figure 2), we observe the formation of internal and boundary layers with decreasing .
The relative error the approximate solution of the problem with the fractional power of the elliptic operator is determined as follows:
where is the exact solution of the problem (7) and is the approximate one. The calculations were performed using the generalized GaussLaguerre quadrature formula with . A fairly detailed grid is used with .
The dependence of the accuracy of the approximate solution of the problem (7) on the key computational parameters is traced according to the data in Table 2 and Table 3. In particular, the most important is the dependence on the number of nodes of the quadrature formula . Note that as the parameter increases, the accuracy of the approximate solution does not decrease. Therefore, we choose . Comparison of Table 2 and Table 3 shows a significant drop in accuracy for a problem with a nonsmooth right side.
error  

0  25  1.867333e07  2.234522e07  9.979937e08  3.388432e08  1.713416e08  
1.101566e03  1.066598e03  6.414770e04  3.516336e04  2.428239e04  
50  2.826563e08  2.756374e08  8.744432e09  2.113860e09  8.727410e10  
5.127673e04  4.435298e04  2.247211e04  1.025590e04  6.344683e05  
100  4.257879e09  3.399113e09  7.664039e10  1.315780e10  4.425344e11  
2.358176e04  1.864016e04  7.926752e05  3.057283e05  1.697664e05  
1  25  1.042484e05  4.322339e06  1.048142e06  2.702530e07  1.233633e07  
5.917636e03  3.726284e03  1.777896e03  8.772466e04  5.808148e04  
50  1.604299e06  5.427969e07  9.390374e08  1.725236e08  6.414784e09  
2.725576e03  1.538243e03  6.132594e04  2.543889e04  1.527295e04  
100  2.460128e07  6.779333e08  8.335361e09  1.089103e09  3.301012e10  
1.267140e03  6.489867e04  2.177537e04  7.583209e05  4.092155e05  
2  25  3.445381e05  1.438506e05  3.540165e06  9.253987e07  4.236814e07  
9.668806e03  6.152098e03  2.969534e03  1.475692e03  9.797893e04  
50  5.381860e06  1.832602e06  3.200697e07  5.929701e08  2.213631e08  
4.393399e03  2.526119e03  1.037259e03  4.408139e04  2.672836e04  
100  8.326382e07  2.311031e07  2.870886e08  3.782843e09  1.151544e09  
2.059608e03  1.068331e03  3.650104e04  1.287746e04  6.980999e05 
error  

0  25  1.084357e02  3.938291e03  6.485652e04  1.177557e04  4.442441e05  
4.370840e01  1.905319e01  6.248551e02  2.430478e02  1.442032e02  
50  6.312810e03  1.936705e03  2.296383e04  2.977748e05  9.176919e06  
4.005614e01  1.588500e01  4.413926e02  1.444312e02  7.761052e03  
100  3.552338e03  9.377168e04  8.059058e05  7.455089e06  1.873778e06  
3.725840e01  1.330906e01  3.095324e02  8.540190e03  4.166878e03  
1  25  4.076650e02  1.271564e02  2.054384e03  3.768885e04  1.428314e04  
5.324292e01  2.576112e01  9.243630e02  3.763498e02  2.275878e02  
50  2.485088e02  6.386131e03  7.389697e04  9.694268e05  3.010355e05  
4.969681e01  2.169114e01  6.579678e02  2.259551e02  1.239358e02  
100  1.494239e02  3.170517e03  2.626296e04  2.456552e05  6.228987e06  
4.616624e01  1.812017e01  4.657524e02  1.352842e02  6.689115e03  
2  25  5.766283e02  1.970018e02  3.575837e03  7.088258e04  2.777896e04  
5.589886e01  2.869214e01  1.107642e01  4.739887e02  2.933699e02  
50  3.549610e02  9.985667e03  1.305410e03  1.877383e04  6.099850e05  
5.225493e01  2.425210e01  7.941726e02  2.898929e02  1.636029e02  
100  2.157698e02  5.001624e03  4.685264e04  4.819297e05  1.280421e05  
4.867483e01  2.038174e01  5.642903e02  1.737444e02  8.854464e03 
The dependence of the accuracy of the approximate solutions on the discrete problem (7) of the number of nodes in space is of great interest as well. The calculation results for various values of are presented in Table 4 and Table 5. The dimension of a discrete problem practically does not affect the accuracy of quadrature formula for an approximate solution of fractional elliptical BVP.
error  

32  25  8.409942e08  1.456405e07  7.988626e08  2.925977e08  1.517098e08  
6.600418e04  8.573971e04  5.919345e04  3.250000e04  2.188000e04  
50  5.068313e09  9.355974e09  4.825193e09  1.451950e09  6.457840e10  
1.553199e04  2.178892e04  1.602107e04  8.719835e05  5.690694e05  
100  1.238361e10  2.538166e10  1.438117e10  4.265139e11  1.791286e11  
2.100015e05  3.149206e05  2.531976e05  1.457506e05  9.665245e06  
64  25  1.610774e07  2.073936e07  9.562513e08  3.287397e08  1.669788e08  
1.032142e03  9.701921e04  6.310096e04  3.372346e04  2.391175e04  
50  2.002888e08  2.324099e08  7.969464e09  1.982252e09  8.265304e10  
4.305003e04  4.156202e04  1.998006e04  9.833823e05  6.224204e05  
100  1.943737e09  2.230604e09  6.161669e10  1.139644e10  3.927857e11  
1.420510e04  1.499068e04  7.332841e05  2.840727e05  1.554351e05  
128  25  1.820715e07  2.204179e07  9.898763e08  3.368579e08  1.704815e08  
1.087668e03  1.060296e03  6.395206e04  3.508900e04  2.412719e04  
50  2.676261e08  2.678691e08  8.598566e09  2.088511e09  8.637732e10  
4.754600e04  4.383331e04  2.231665e04  1.006521e04  6.322556e05  
100  3.761240e09  3.195480e09  7.399354e10  1.283428e10  4.332366e11  
2.236805e04  1.701325e04  7.822322e05  2.923203e05  1.683085e05 
error  

32  25  2.006090e03  1.828621e03  4.691900e04  9.911544e05  3.910657e05  
1.187981e01  1.098096e01  5.279046e02  2.143857e02  1.226111e02  
50  3.760510e04  4.506979e04  1.189540e04  2.059487e05  6.933095e06  
5.566576e02  5.283872e02  2.524998e02  1.103353e02  6.396780e03  
100  3.792739e05  5.832591e05  1.907976e05  3.269577e06  1.009050e06  
1.965240e02  2.111303e02  1.054680e02  4.130081e03  2.183357e03  
64  25  6.360886e03  3.165742e03  5.971673e04  1.123797e04  4.279690e05  
3.258161e01  1.792333e01  5.846466e02  2.324461e02  1.372679e02  
50  2.574703e03  1.312433e03  1.966289e04  2.738334e05  8.597362e06  
2.021267e01  1.333749e01  3.999555e02  1.290906e02  7.372346e03  
100  7.619762e04  4.606884e04  5.933321e05  6.302869e06  1.650515e06  
1.036411e01  7.774238e02  2.644230e02  7.584779e03  3.532676e03  
128  25  9.378730e03  3.740118e03  6.383832e04  1.168124e04  4.414578e05  
4.279399e01  1.880510e01  6.172671e02  2.397168e02  1.435631e02  
50  4.994286e03  1.766193e03  2.226535e04  2.933208e05  9.074022e06  
3.718976e01  1.502341e01  4.274832e02  1.427241e02  7.658653e03  
100  2.419279e03  7.963929e04  7.576465e05  7.224747e06  1.831256e06  
2.838417e01  1.268924e01  2.938177e02  8.281267e03  4.000634e03 
6 Conclusions

To solve the problem with a fractional power of a positive definite operator, we use the known integral representation of the solution by solving the Cauchy problem. The possibilities of such representation are expanded due to the transition to the problem with the new initial conditions. The key dependence of the solution on the time of the auxiliary nonstationary problem, associated with the constant positive definiteness operator, is defined for constructing optimal Gaussian quadratures.

An approximate solution is based on the use of the generalized GaussLaguerre quadrature formula. The accuracy of the quadrature formula is studied on the representative parametric integration problem depending on the key parameters: degrees and the number of nodes of the quadrature formula.

The model problem with the fractional power of the Laplace operator with a smooth and discontinuous righthand side is considered in a rectangle. Its approximate solution is obtained based on the integral representation by solving the parabolic Cauchy problems. The results of the numerical experiments on the influence of the smoothness of the righthand side on the accuracy of the approximate solutions using a different number of integration nodes are presented.
References
 (1) D. Baleanu, Fractional Calculus: Models and Numerical Methods, World Scientific, New York, 2012.
 (2) V. V. Uchaikin, Fractional Derivatives for Physicists and Engineers, Higher Education Press, 2013.
 (3) C. Pozrikidis, The Fractional Laplacian, CRC Press, 2018.
 (4) N. J. Higham, Functions of Matrices: Theory and Computation, SIAM, Philadelphia, 2008.
 (5) A. Bonito, J. P. Borthagaray, R. H. Nochetto, E. Otarola, A. J. Salgado, Numerical methods for fractional diffusion, Computing and Visualization in Science 19 (56) (2018) 19–46.
 (6) H. R. Stahl, Best uniform rational approximation of on [0, 1], Acta Mathematica 190 (2) (2003) 241–306.
 (7) S. Harizanov, R. Lazarov, S. Margenov, P. Marinov, Y. Vutov, Optimal solvers for linear systems with fractional powers of sparse spd matrices, Numerical Linear Algebra with Applications 25 (5) (2018) e2167.
 (8) S. Harizanov, R. Lazarov, P. Marinov, S. Margenov, J. Pasciak, Analysis of numerical methods for spectral fractional elliptic equations based on the best uniform rational approximation, arXiv:1905.08155.
 (9) A. V. Balakrishnan, Fractional powers of closed operators and the semigroups generated by them., Pacific Journal of Mathematics 10 (2) (1960) 419–437.
 (10) M. S. Birman, M. Z. Solomjak, Spectral Theory of Selfadjoint Operators in Hilbert Space, Kluwer academic publishers, 1987.
 (11) A. Bonito, J. Pasciak, Numerical approximation of fractional powers of elliptic operators, Mathematics of Computation 84 (295) (2015) 2083–2110.
 (12) L. Aceto, P. Novati, Rational approximation to the fractional Laplacian operator in reactiondiffusion problems, SIAM Journal on Scientific Computing 39 (1) (2017) A214–A228.
 (13) L. Aceto, P. Novati, Rational approximations to fractional powers of selfadjoint positive operators authors, Numerische Mathematik 143 (1) (2019) 1–16.
 (14) P. N. Vabishchevich, Approximation of a fractional power of an elliptic operator, arXiv:1905.10838.

(15)
L. Caffarelli, L. Silvestre, An extension problem related to the fractional Laplacian, Communications in Partial Differential Equations 32 (8) (2007) 1245–1260.
 (16) R. H. Nochetto, E. Otárola, A. J. Salgado, A PDE approach to fractional diffusion in general domains: a priori error analysis, Foundations of Computational Mathematics 15 (3) (2015) 733–791.
 (17) P. N. Vabishchevich, Numerically solving an equation for fractional powers of elliptic operators, Journal of Computational Physics 282 (1) (2015) 289–302.
 (18) B. Duan, R. Lazarov, J. Pasciak, Numerical approximation of fractional powers of elliptic operators, IMA Journal of Numerical Analysis (drz013) (2019) 1–26. doi:10.1093/imanum/drz013.
 (19) R. Ciegis, P. Vabishchevich, High order numerical schemes for solving fractional powers of elliptic operators, Computers & Mathematics with Applications (2019) 1–11. In Press. doi:10.1016/j.camwa.2019.08.012.
 (20) P. R. Stinga, User’s guide to the fractional Laplacian and the method of semigroups, in: Handbook of Fractional Calculus and Applications. Volume 2: Fractional Differential Equations, de Gruyter, 2019, pp. 235–266.
 (21) N. Cusimano, F. del Teso, L. GerardoGiorda, G. Pagnini, Discretizations of the spectral fractional Laplacian on general domains with Dirichlet, Neumann, and Robin boundary conditions, SIAM Journal on Numerical Analysis 56 (3) (2018) 1243–1272.
 (22) N. Cusimano, F. del Teso, L. GerardoGiorda, Numerical approximations for fractional elliptic equations via the method of semigroups, arXiv:1812.01518.
 (23) P. J. Davis, P. Rabinowitz, Methods of Numerical Integration, Dover Publications, 2007.
 (24) A. A. Samarskii, The Theory of Difference Schemes, Marcel Dekker, New York, 2001.
 (25) A. A. Samarskii, E. S. Nikolaev, Numerical Methods for Grid Equations. Vol. I, II, Birkhauser Verlag, Basel, 1989.