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 selfadjoint positive operator, the following integral representation holds Balakrishnan (1960); Birman and Solomjak (1987):
(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
To obtain the integral on a finite interval, in Frommer et al. (2014), there is employed the transformation
From (1), we have
(2) 
To approximate the righthand side, we apply the GaussJacobi 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 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 Birman and Solomjak (1987); Kwaśnicki (2017). Introduce the elliptic operator as
(3) 
with coefficients , . The operator is defined on the set of functions that satisfy on the boundary the following conditions:
(4) 
In the Hilbert space , we define the scalar product and norm in the 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:
(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 :
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
(6) 
under the restriction .
We consider the simplest case, where the computational domain is a rectangle
To solve approximately the problem (6), we introduce in the domain a uniform grid
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
For the above grid operators (see Samarskii (2001); Samarskii and Nikolaev (1989)), we have
where is the grid identity operator. 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
We have
where eigenfunctions form a basis in . Therefore
For the fractional power of the operator , we have
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))
(9) 
where
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:
(10) 
(11) 
Using (10), from (9), we get the following integral representation of the fractional power operator:
(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
(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 twoparameter formulas (12) and (13), one parameter is free. If we put , then the integrand in the righthand side of (12) and (13) have no singularities at the ends of the integration interval. Assuming , from (12), we get
(14) 
In case of (13), we have and
(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
(16) 
The integrand in (16) has no singularities for . For , we have
(17) 
The same change of the integration variable in the representation (15) results in
(18) 
We use this representation for , where . For , we get
(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):
Using the association , we arrive at the following integral representation:
(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:
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):
It is also inconvenient for the direct application due to the singularity of the integrand for . Using the new variable , we arrive at
The integrand has no sigularities if
For instance, it is sufficient to take , where . With this in mind, we obtain a new integral representation of the fractional power operator:
(21) 
Selecting , we guarantee that the th derivative of the integrand has no singularities.
4 Numerical integration
The accuracy of quadrature formulas is investigated for the following integral representation (see (16)) of the function with :
(22) 
with the parameter . When focusing on (21), we are interested in the integral representation
(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: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
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.
1  5.374386e03  3.975611e03  1.777798e03  6.057090e04  1.860264e04  

2  9.694686e05  1.546665e04  6.365306e05  3.000009e05  1.090598e05  
50  3  9.583340e05  8.914710e05  6.363818e05  3.000215e05  1.092388e05 
4  1.264541e04  1.189327e04  8.483208e05  4.000540e05  1.457062e05  
5  1.548898e04  1.485671e04  1.060141e04  5.000617e05  1.821371e05  
6  1.841219e04  1.780018e04  1.271753e04  6.000429e05  2.185621e05  

1  1.294685e03  1.985566e03  8.876394e04  3.029813e04  9.380258e05 
2  2.433513e05  3.853785e05  1.591494e05  7.501982e06  2.730852e06  
100  3  2.426684e05  2.234929e05  1.591401e05  7.502111e06  2.731974e06 
4  3.232626e05  2.983899e05  2.121750e05  1.000297e05  3.642973e06  
5  4.029470e05  3.731919e05  2.652025e05  1.250368e05  4.553745e06  
6  4.814931e05  4.478680e05  3.182170e05  1.500422e05  5.464479e06  
1  1.647562e04  9.922321e04  4.435056e04  1.515251e04  4.710344e05  
2  6.092649e06  9.617959e06  3.978839e06  1.875618e06  6.829863e07  
200  3  6.088573e06  5.596207e06  3.978781e06  1.875626e06  6.830565e07 
4  8.129377e06  7.471178e06  5.304967e06  2.500845e06  9.107634e07  
5  1.016443e05  9.345551e06  6.631108e06  3.126054e06  1.138456e06  
6  1.219226e05  1.121914e05  7.957167e06  3.751252e06  1.366146e06 
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.
1  5.347988e03  3.966817e03  1.772738e03  1.500737e05  1.640307e05  

2  8.982157e05  1.536406e04  6.366643e05  6.004424e05  6.565671e05  
50  3  4.946450e05  4.509810e05  3.184883e05  4.509810e05  4.946450e05 
4  6.626068e05  6.021592e05  4.248002e05  6.021592e05  6.626068e05  
5  8.336423e05  7.541989e05  5.314105e05  7.541989e05  8.336423e05  
6  1.008694e04  9.071839e05  6.381866e05  9.071839e05  1.008694e04  

1  1.289544e03  1.983358e03  8.863758e04  3.751449e06  4.099039e06 
2  1.639890e05  3.840905e05  1.591577e05  1.500672e05  1.639890e05  
100  3  1.231272e05  1.125908e05  7.958861e06  1.125908e05  1.231272e05 
4  1.643508e05  1.501738e05  1.061274e05  1.501738e05  1.643508e05  
5  2.057415e05  1.878095e05  1.326849e05  1.878095e05  2.057415e05  
6  2.473359e05  2.255024e05  1.592525e05  2.255024e05  2.473359e05  
1  1.638621e04  9.916965e04  4.431897e04  9.378377e07  1.024652e06  
2  4.098779e06  9.602172e06  3.978891e06  3.751408e06  4.098779e06  
200  3  3.074926e06  2.813809e06  1.989506e06  2.813809e06  3.074926e06 
4  4.101019e06  3.752073e06  2.652733e06  3.752073e06  4.101019e06  
5  5.128125e06  4.690667e06  3.316076e06  4.690667e06  5.128125e06  
6  6.156442e06  5.629614e06  3.979483e06  5.629614e06  6.156442e06 
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 .
1  1.853227e02  1.799633e02  1.273240e02  6.002109e03  2.185848e03  

2  3.949860e04  2.269479e04  5.687915e05  1.750531e05  5.280150e06  
50  3  1.275793e04  2.277837e05  4.469762e06  1.174967e06  3.586255e07 
4  9.400469e05  5.571588e06  5.614120e07  9.194800e08  3.276035e08  
5  2.658768e04  4.032756e06  4.091810e07  5.210468e08  3.395709e08  
6  9.293393e04  7.461337e06  6.803497e07  7.645749e08  1.132277e07  
1  9.253872e03  8.993163e03  6.366198e03  3.001054e03  1.092924e03  
2  2.695390e05  5.624706e05  1.413623e05  4.328911e06  1.290228e06  
100  3  3.871962e06  2.782215e06  5.489516e07  1.401037e07  3.949594e08 
4  1.981503e06  3.273869e07  3.283343e08  5.249440e09  1.580743e09  
5  2.380378e06  2.949317e07  2.549296e08  3.259146e09  4.821441e10  
6  4.575033e06  4.604162e07  4.247676e08  4.799821e09  6.473959e10  
1  3.535332e03  4.491582e03  3.183099e03  1.500527e03  5.464620e04  
2  1.864876e07  1.399930e05  3.523765e06  1.076370e06  3.189672e07  
200  3  3.715126e08  2.692643e07  6.802848e08  1.710196e08  4.636237e09 
4  7.661090e08  2.266696e08  1.986647e09  3.136784e10  8.664345e11  
5  1.420196e07  2.919238e08  1.591996e09  2.037560e10  3.016176e11  
6  2.398536e07  3.930730e08  2.653151e09  3.000763e10  4.047396e11 
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 righthand side
(24) 
The numerical solution of the problem (7) with the righthand side (24) is depicted in Figure 4 for different values of . For convenience of the comparison, there is shown the function
For the discontinuous righthand side (24), we observe the formation of internal boundary layers with decreasing .
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:
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.
error  

50  1.669630e06  3.172818e06  6.976081e06  9.886524e06  6.807410e06  
1.931587e06  3.483513e06  7.301196e06  1.006909e05  6.868003e06  
100  4.234571e07  7.955079e07  1.749686e06  2.493483e06  1.746540e06  
4.899010e07  8.733658e07  1.830969e06  2.538547e06  1.758218e06  
200  1.062490e07  1.990234e07  4.377758e07  6.247408e07  4.394942e07  
1.229199e07  2.184997e07  4.580968e07  6.359638e07  4.422973e07 
error  

50  1.286770e04  2.641118e07  5.182607e08  2.110766e08  2.297871e07  
1.558345e04  3.437617e07  5.442683e08  3.023970e08  2.475985e07  
100  7.460526e08  9.578675e09  3.256882e09  1.074501e09  2.876888e10  
1.013479e07  1.051652e08  3.407962e09  1.094340e09  2.910291e10  
200  2.528173e09  6.061693e10  2.038772e10  6.772998e11  1.906629e11  
2.924801e09  6.656861e10  2.133260e10  6.895408e11  1.920284e11 
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 righthand side (see (24)). For
(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
(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).
error  

50  1.209581e04  3.444934e07  9.686650e08  2.433874e08  7.815106e08  
1.141529e04  3.688685e07  9.631475e08  2.422181e08  8.320527e08  
100  1.030346e07  2.380002e08  6.073046e09  1.493598e09  3.465175e10  
1.108808e07  2.349288e08  6.059424e09  1.493626e09  3.468453e10  
200  7.245613e09  1.539020e09  3.799707e10  9.359799e11  2.202853e11  
6.960763e09  1.518666e09  3.791153e10  9.359548e11  2.203251e11 
error  

50  5.800855e05  4.000618e07  1.056566e07  2.455535e08  4.044510e08  
6.236828e05  3.995393e07  1.061651e07  2.460929e08  3.914376e08  
100  6.407204e08  2.899259e08  6.627865e09  1.541997e09  3.513618e10  
6.549602e08  2.924586e08  6.659433e09  1.545095e09  3.516994e10  
200  9.952338e09  1.877999e09  4.146733e10  9.659714e11  2.227958e11  
1.006940e08  1.894607e09  4.166460e10  9.678749e11  2.229965e11 
6 Conclusions

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.

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.

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 selfadjoint 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.