 # Approximate representation of the solutions of fractional elliptical BVP through the solution of parabolic IVP

Boundary value problem for a fractional power of an elliptic operator is considered. An integral representation by means of a standard solution problem for parabolic equations is used to solve such problems. Quadrature generalized Gauss-Laguerre formulas are constructed. We examine the effect of key parameters on the accuracy of the approximate solution: the number of nodes of the quadrature and fractional power of the operator. Computational experiments were performed to model two-dimensional problem with a fractional power of an elliptic operator.

## 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

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 vab2019-10838 .

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 semi-infinite 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 pseudo-parabolic 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 three-level schemes of higher order approximation on regular and irregular grids (see, for example, duan2018numerical ; ciegisvab2019-00201 ) 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 Gauss-Laguerre 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 finite-difference 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 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 birman1987spectral . The elliptic operator is introduced as

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

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

In the Hilbert space , we define the scalar product and norm in a 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, (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 :

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

The solution satisfies the equation

 Aαv=f (4)

under the restriction .

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

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

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

 ¯¯¯ω={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 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

 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 Samarskii1989 ; SamarskiiNikolaev1978 ), we have

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

where is the grid identity operator in . 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. (6)

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 (4), we arrive at the discrete problem

 u=A−αb, (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

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

The approximation of is based on the use of one or another quadrature formulas for the right-hand side of (8).

Here, the possibilities of approximating both the fractional power of the operator and approximating the solution of the problem are realized.

1. Rational function approximation :

 A−α≈m∑i=1γi(A+θiI)−1,

where are the weights, and are the nodes of the quadrature formula for (8).

2. An approximate solution to the problem (4):

 u≈m∑i=1γiui,(A+θiI)ui=b,i=1,2,…,m.

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

where is gamma function. In this case, the solution to the problem (7) is

 u=1Γ(α)∫∞0θα−1e−θAbdθ. (10)

Function is a solution of the following Cauchy problem

 dwdt+Aw=0,0
 w(0)=b. (12)

By the representation (10) we have

 u=1Γ(α)∫∞0θα−1w(θ)dθ. (13)

For (1), (2) the problem (11), (12) is the standard Cauchy problem for a second-order 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 non-stationary 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).

1. With quadrature formulas for (9) we come to the approximation with operator exponentials:

 A−α≈m∑i=1γie−θiA.

Here, as before, , are the weights and nodes of the corresponding quadrature formula.

2. For an approximate solution to the problem (7) from (13) we have

 u≈um=m∑i=1γiw(θi), (14)

where is the solution to the problem (11), (12).

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

 αθα−1w(θ)dθ=d(θαw(θ))−θαdwdθdθ=d(θαw(θ))+θαAw.

Integrating by parts allows writing (13) as

 u=1αΓ(α)∫∞0θαv(θ)dθ, (15)

where . From (11), (12) we have the Cauchy problem for :

 dvdt+Av=0,0
 v(0)=Ab. (17)

We arrive at a similar result based on the representation (9). Replacing by , where has a non-negative integer, we get

For the problem solution with a fractional power of the operator (5), it follows from (18) that

 u=1Γ(α+p)∫∞0θα+p−1e−θAApbdθ. (19)

We define function as a solution of the equation (11), supplemented by the initial condition

 w(0)=c,c=Apb. (20)

To solve the problem (5) from (19) we obtain the representation

 u=1Γ(α+p)∫∞0θα+p−1w(θ)dθ. (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 half-line 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

 w(t)=K∑k=1(c,ψk)e−μktψk.

From this perspective, from (21) we obtain

 u=1Γ(α+p)K∑k=1(c,ψk)ψk∫∞0θα+p−1e−μkθdθ. (22)

Introducing a new integration variable from (22) we get

 u=1δα+pΓ(α+p)K∑k=1(c,ψk)ψk∫∞0ξα+p−1e−ξe−δ−1(μk−δ)ξdξ.

In view of this, it is necessary to construct quadrature formulas for the integrals

 S(α+p,κ)=∫∞0ξα+p−1e−ξϱ(ξ;κ)dξ,0<α<1,p≥0, (23)

in which the integrand is

 ϱ(ξ;κ)=e−κξ,κ≥0, (24)

Concerning (22) we have

 κ=κk,κk=μkδ−1,k=1,2,…,K.

When calculating the integrals (23) with the weight function , we focus on the generalized Gauss-Laguerre quadrature formula davis2007methods . In this case

 S(α+p,κ)≈Sm(α+p,κ)=m∑i=1σiϱ(ξi;κ), (25)

where the nodes of the quadrature formula are the zeros of the generalized Laguerre polynomial , and weights are

 σi=Γ(m+α+p)m!ξi(L(α+p−1)m+1(ξi))2,i=1,2,…,m.

The accuracy of the approximate calculation of

is estimated at

by the value of the relative error

 ε=1qmax0≤κ≤105|S(α+p,κ)−Sm(α+p,κ)|,q=max0≤κ≤105S(α+p,κ).

For convenience of the analysis of the calculated data, the dependence on and is highlighted separately.

The accuracy of the generalized Gauss-Laguerre 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 .

## 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

 a(x)=1,c(x)=0,x∈Ω.

For eigenfunctions and eigenvalues (6) we have (see, e.g., SamarskiiNikolaev1978 ):

 ψk(x)=2∏β=1√2lβsin(kβπxβ),x∈ω,μk=2∑β=14h2βsin2kβπ2Nβ,kα=1,2,...,Nα−1,α=1,2.

Direct calculations yield

 μ1=2∑β=14h2βsin2π2Nβ<8(1l21+1l22),λK=2∑β=14h2βcos2π2Nβ<4(1h21+1h22).

The model problem (4) is solved numerically with two right-hand sides: , when

 f1(x)=x21(1−x1)x22(1−x2),f2(x)=1+sgn(x1x2−0.25).

where is the sign function. We consider the problem with both a smooth right-hand side when choosing and a non-smooth one choosing .

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

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

For the discontinuous right-hand side (, Figure 2), we observe the formation of internal and boundary layers with decreasing . Figure 1: Normalized solution of the problem (6) when f(x)=f1(x) on the grid N1=N2=256. Figure 2: Normalized solution of the problem (6) when f(x)=f2(x) on the grid N1=N2=256.

The relative error the approximate solution of the problem with the fractional power of the elliptic operator is determined as follows:

 ε2=∥˜u−u∥∥u∥,ε∞=∥˜u−u∥∞∥u∥∞,∥u∥∞=maxx∈ω|u(x)|,

where is the exact solution of the problem (7) and is the approximate one. The calculations were performed using the generalized Gauss-Laguerre 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 non-smooth right side.

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.

## 6 Conclusions

1. 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 non-stationary problem, associated with the constant positive definiteness operator, is defined for constructing optimal Gaussian quadratures.

2. An approximate solution is based on the use of the generalized Gauss-Laguerre 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.

3. The model problem with the fractional power of the Laplace operator with a smooth and discontinuous right-hand 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 right-hand 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 (5-6) (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 Self-adjoint 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 reaction-diffusion problems, SIAM Journal on Scientific Computing 39 (1) (2017) A214–A228.
• (13) L. Aceto, P. Novati, Rational approximations to fractional powers of self-adjoint 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.
• (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.
• (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. Gerardo-Giorda, 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. Gerardo-Giorda, 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.