Numerical Solution of Fractional Elliptic Problems with Inhomogeneous Boundary Conditions

The numerical solution of fractional-order elliptic problems is investigated in bounded domains. According to real-life situations, we assumed inhomogeneous boundary terms, while the underlying equations contain the full-space fractional Laplacian operator. The basis of the convergence analysis for a lower-order boundary element approximation is the theory for the corresponding continuous problem. In particular, we need continuity results for Riesz potentials and the fractional-order extension of the theory for boundary integral equations with the Laplacian operator. Accordingly, the convergence is stated in fractional-order Sobolev norms. The results were confirmed in a numerical experiment.


Introduction
The numerical solution of fractional-order diffusion problems has a number of applications to simulate real-world phenomena such as groundwater flows [1], the dynamics of some proteins [2], geophysical electromagnetics [3], or population dynamics [4]. In the last few decades-initiating in the work [5]-different approaches were developed to investigate this important problem. In the space-fractional case, the core of any investigation is the study of the underlying fractional elliptic problem.
Physically, this phenomenon (also called anomalous diffusion) can be recognized by the sublinear dependence of the mean-squared displacement of single particles as a function of time. For further explanations, we refer to the review article [6].
The main difficulty is to deal with these problems on bounded domains. While many different definitions of the full-space fractional Laplacian are equivalent [7], for bounded domains, even in the case of homogeneous boundary conditions, substantially different definitions exist [8]. As in many situations in applied mathematics, the functional analytic point of view leads to the most adequate model [9]. Since the discretization of the fractional Laplacian results in a dense matrix in any case, the numerical treatment of the related problems needs special techniques; see [10][11][12]. At the same time, based on the nonlocal property of the fractional Laplacian, the zero-order extension to R d also makes sense. The corresponding numerical approximation-containing integrals on R d -also needs special care [13]. For a review paper on nonlocal models and their numerical approximations with further literature, we refer to [14].
The real challenge, however, is the incorporation of inhomogeneous boundary conditions. Indeed, it seems to be hard to model this nonlocal phenomenon if we only have boundary data. Namely, for the computation of the full-space fractional Laplacian, we would need the values of the underlying function on the whole space R d .
A few attempts were made to resolve this problem in the framework of the conventional PDEs, starting from the work [15], which were summarized in detail in [16]. An alternative approach can be found in [17]. It turns out that the boundary integral form of these problems offers a meaningful alternative.
This idea was initiated in [18] and completed with all possible cases in [19]. This can also constitute the basis of the present study, where boundary integral elements are used.
Accordingly, the objective of this work is to set the theoretical basis of the boundary element method for fractional-order elliptic problems, which can successfully deal with inhomogeneous boundary conditions. Extending the classical approach for the secondorder elliptic case, we state the convergence properties of this approach and present a series of numerical experiments.

Materials and Methods
In practice, the unknown function u has to be computed on a bounded domain Ω ⊂ R d with d = 2, 3, which is assumed to have Lipschitz boundary. The convergence results are given in terms of fractional Sobolev norms. For the full space R d , these are usually defined using the Fourier transform F on the vector space: Here, s is an arbitrary positive index. Using this, one can define the local version with: For the details, we refer to the monograph [20].
In the error analysis, we considered a tessellation: of the boundary such that the surface subdomains {Γ j } N j=1 are disjoint open Lipschitz domains, which were assumed to be a C 1 -diffeomorphic map of a convex polygon for d = 3 or a bounded interval for d = 2.
We also need an estimate concerning piecewise Sobolev spaces.
Lemma 1 ([21], p. 37). For the above tessellation and any v ∈ H s (∂Ω) with s > 0, we have the estimate: For stating the well-posedness of the continuous problem, we also need the Sobolev space: with the corresponding norm. The full-space fractional Laplacian operator [7] can be defined pointwise as: where B(x, r) C = R d \ B(x, r). One can also define its weak form as the function (−∆ R d ) α for which: is satisfied for all v ∈ C ∞ 0 (R d ). The success of the boundary integral methods hinges on the existence of the fundamental solution φ α of (−∆ R d ) α , which is given with: We also used the concept of the trace operator γ, which for a bounded Lipschitz domain Ω is continuous as: for s ∈ ( 1 2 , 3 2 ); see [22], Theorem 3.38. We applied the conventional notation ·, · −β,β for the duality pairing between H −β (∂Ω) and H β (∂Ω) with some positive exponent β.
To streamline the presentation, we use the notation · s,Ω for the Sobolev norm in the space H s (Ω) and similarly with Γ. Furthermore, in the estimates, the relation was used if the left-hand side can be estimated with a domain-dependent constant times the right-hand side.
The main problem in practice is the numerical solution of the elliptic boundary value problem: where f , g are given real functions. At this stage, the differential operator (−∆) α is not yet defined. In any case, it should be linear, such that u can be given as u = u 1 + u 2 , where: and: The problem in (6) is equipped with homogeneous boundary conditions, which has already been studied by many authors, as mentioned in Section 1. Therefore, in this study, we investigate the numerical solution of the problem: corresponding to (7), using boundary element approximation.
The result in [19] ensures the well-posedness of (8) as follows.
The main building block of the proof is the analysis of the surface potential corresponding to φ α on ∂Ω, which is given for any x ∈ R d with: where φ α : R d \ 0 → R is given in (4). In particular, the following result was established in [19].

Error Analysis
The surface potential S α allows us to reformulate the problem in (8). We want to find the unique solution to the equation: To approximate the solution of Equation (12), first, we introduced the trial space for a Galerkin approach. In all cases, the subscript h denotes a parameter corresponding to the discretization.
We made use of the trial space: i.e., S 0 h (Γ) is piecewise constant with respect to the surface mesh. This trial space fulfills the Now, we can apply the Galerkin formulation: which implies the Galerkin orthogonality relation: The starting point of the analysis is to state the quasi-optimality of any Galerkin method as follows.
Lemma 2 (Céa-lemma). For the Galerkin solution u h of Equation (13), we have: Proof. Using the coercivity in (11), the orthogonality relation in (14), and finally, the continuity in (10), we obtain that for any w h ∈ S 0 h , the following estimation is valid: Dividing both sides of the equation by u − u h 1 2 −α,∂Ω , we obtain: as stated in the lemma.
We can verify the stability of the above Galerkin method by applying (11) for u = u h : Dividing both sides by u h 1 2 −α , we obtain stability on the boundary: Now, applying (10), we finally obtain stability for the reconstructed solution as: We need the approximation property of S 0 h on each surface subdomain Γ j separately. For this, we introduced the projection operator Q j h : . An appropriate statement can be found in [21], Theorem 10.14, on p. 237, which was adapted to our notation using 1/2 − α := σ. Lemma 3. For any α ∈ [ 1 2 , 3 2 ] and s ∈ [1/2 − α, 1], the projection operator Q j h has the following approximation property: Note that the condition in Theorem 1 and in Theorem 2 for α is also sufficient in Lemma 3.
We can now prove the following error estimate for the boundary element approximation.

Theorem 4.
Assume that the conditions in Theorems 1 and 2 are satisfied. Then, we have the error estimation: for the boundary element based approximationũ h of the solution of (8).
Proof. We first note thatũ h is obtained via reconstruction by applying the surface potential S α as follows:ũ Applying the same for the analytic solution u of (8), we obtain the following expansion for the error of the reconstructed solution: Using this, then applying (10) for u − u h in Theorem 2, using Lemma 1, and finally, the estimation (17), we obtain: as stated in the theorem.

Numerical Experiments
Let us take the domain Ω ∈ R 2 to be the unit disk. We investigated the model problem: for different values of α ∈ ( 1 2 , 3 4 ] according to Theorem 1, where the analytic solutionũ is given with:ũ where: and: g(x) =ũ(x) |x| = 1.
(24) For the experimental error analysis, we first computed the analytic solutionũ inside the domain given in (22) using f in (23). For this, a quadrature was applied on the unit circle with 10,000 equidistributed points. In this way, even the analytic solution should be approximated for the experimental error analysis.
We also computed the boundary condition g for our model problem. For this, we applied the built-in MATLAB algorithm integral, which is suitable to compute the singular integrals on the boundary. In our experiment, we solved Equation (21)  integral cannot compute the boundary conditions accurately enough, so we omitted these cases.
Note that in a real-life problem, usually, the boundary condition g is given (instead of the function f ), so the method can still be used.
Then, we applied the boundary element method in (13) to compute the approximation u h using the boundary data g in (24). To compute the terms in the corresponding discrete variational formulation (13), we need to compute the entries of the matrix S h given by: To obtain the diagonal elements here, one needs to take special care, as a conventional three-point Gauss quadrature does not deliver sufficient accuracy. Instead, we applied the built-in MATLAB subroutine integral2, which can handle singularities in the twodimensional integrand.
Completing the linear system corresponding to (13), we computed the entries in its right-hand side g h as: which were approximated using a simple midpoint quadrature. Finally, we reconstructedũ h according to (20), which was the most technical and time-consuming part of the computation. This was performed in each point of a grid in the unit disk. This grid size was approximately the same as the one on the boundary. To approximate the integrals in (20) for a fixed x ∈ Γ, we applied a midpoint approximation on the boundary using the same mesh as for the solution of (13).
In the boundary element method, we used different partitions of the unit circle, into N = 8, 16, 32, 64, 128, 256, and 512 uniform parts. In each case, the computational error was calculated in the L 2 (Ω)-norm. After the consecutive refinements, we estimated the corresponding convergence rate, which is called the rate in the tables. Note that in (19), the error estimation is given in H α (Ω)-norm, so we expect the order of convergence to be higher in the L 2 (Ω)-norm. According to Theorem 4, the expected rate of convergence in the H α (Ω)-norm is s − (1/2 − α). Sinceũ| Ω is in H 1 (Ω), this is equal to 1.25 for α = 0.75, 1.2 for α = 0.7, and 1.15 for α = 0.65. In Tables 1-3, the results of the experimental error analysis can be found for different exponents. For relative coarse meshes, in each case, the convergence rate is even faster than expected. At the same time, for a relatively fine mesh, the computation with singular integrals becomes inaccurate, which results in oscillations in u h and deteriorates the accuracy of the reconstructed approximationũ h . At the same time, we could achieve a relatively small error of magnitude 10 −4 before this phenomenon affected the convergence rate.
We can observe in Tables 1 and 2 that the convergence rate decreases already for N = 512 and N = 256, respectively. In these cases, α is closer to 1 2 , so that the singularities at the computation of the boundary condition become sharper. This makes the applied quadratures less accurate. For the coarser meshes, the convergence rate is still better than the expected rate. In the case of α = 0.65, the deterioration of the convergence rates starts earlier at N = 256. For the coarser meshes, the convergence rate is still better than the expected rate 1.15.
In Figure 1, the numerical solution is shown together with the given boundary data, while in Figure   Here, one can already observe some deviation, which is shown in detail in Figure 2, where only the error is presented. At the boundary, indeed, relatively large errors emerge, which were then smoothed out by applying the reconstruction.

214
The idea to convert fractional-order elliptic boundary value problems into boundary 215 integral equations leads to well-posed problems with some extra conditions. Also, as 216 the solution is defined everywhere, this corresponds to the non-local nature of these 217 problems. Their numerical approximations, which were investigated here, offer also 218 efficient methods. To see this, we recall that for conventional elliptic boundary value 219 problems, boundary element methods lead to dimensional reduction but also to full 220 matrices in the resulting linear algebraic problems. In case of fractional order problems, 221 however, boundary element methods offer a clear advantage: we have to compute 222 Figure 4. Pointwise error of the numerical solutionũ h of (21) for α = 0.65 and N = 512 inside the boundary. We can observe that the error is significantly larger than in the case of α = 0.75.

Discussion
The idea to convert fractional-order elliptic boundary value problems into boundary integral equations leads to well-posed problems with some extra conditions. Furthermore, as the solution is defined everywhere, this corresponds to the nonlocal nature of these problems. Their numerical approximations, which were investigated here, offer also efficient methods. To see this, we recall that for conventional elliptic boundary value problems, boundary element methods lead to dimensional reduction, as well as to full matrices in the resulting linear algebraic problems. In the case of fractional-order problems, however, boundary element methods offer a clear advantage: we have to compute with dense matrices anyway, such that using these methods significantly reduces the computational costs. Moreover, as we have pointed out, in this framework, it is possible to model and simulate the presence of inhomogeneous boundary data. An important future research direction is to extend the well-posedness of the original problem in (8) for α > 3 4 . Furthermore, the application of higher-order elements and the performance of a simple approach, the so-called method of fundamental solutions [23], should be investigated.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; nor in the decision to publish the results.
Sample Availability: The MATLAB code for the numerical experiments is available from the authors.