Quasi-optimal convergence rate for an adaptive method for the integral fractional Laplacian

TL;DR

This paper develops and analyzes an adaptive finite element method with a novel weighted residual error estimator for the integral fractional Laplacian, achieving quasi-optimal convergence rates even in challenging regularity regimes.

Contribution

It introduces a weighted residual a posteriori error estimator for the fractional Laplacian and proves its effectiveness in driving an adaptive algorithm with optimal convergence rates.

Findings

The error estimator is reliable and effective for $0<s<1$.

The adaptive algorithm achieves quasi-optimal convergence rates.

The analysis includes new local inverse estimates for the fractional Laplacian.

Abstract

For the discretization of the integral fractional Laplacian $(-\Delta)^s$, $0 < s < 1$, based on piecewise linear functions, we present and analyze a reliable weighted residual a posteriori error estimator. In order to compensate for a lack of $L^2$-regularity of the residual in the regime $3/4 < s < 1$, this weighted residual error estimator includes as an additional weight a power of the distance from the mesh skeleton. We prove optimal convergence rates for an $h$-adaptive algorithm driven by this error estimator. Key to the analysis of the adaptive algorithm are local inverse estimates for the fractional Laplacian.

Equations

$${\Vert {\tilde{h}}_{\ell }^s(-\Delta {)}^s{v}_{\ell }\Vert }_{L^2(\Omega )}\le C{\Vert v_{\ell }\Vert }_{H^s(\Omega )},$$

$$|v{|}_{H^t(\omega )}^2={\int }_{x\in \omega }{\int }_{y\in \omega }\frac{|v(x)-v(y){|}^2}{{\left|x-y\right|}^{d+2t}}dxdy,\Vert v{\Vert }_{H^t(\omega )}^2=\Vert v{\Vert }_{L^2(\omega )}^2+|v{|}_{H^t(\omega )}^2,$$

$${\tilde{H}}^t(\Omega ):=\{u\in H^t({\mathbb{R}}^d):u\equiv 0\text{on}{\mathbb{R}}^d\backslash \Omega \}$$

$${\Vert v\Vert }_{{\tilde{H}}^t(\Omega )}^2:={\Vert v\Vert }_{H^t(\Omega )}^2+{\Vert v/{\rho }^t\Vert }_{L^2(\Omega )}^2,$$

$$(-\Delta {)}^s:H^t({\mathbb{R}}^d)\to H^{t-2s}({\mathbb{R}}^d),t\ge s,$$

$$(-\Delta {)}^{s}u(x):=C(d,s)\text{P.V.}{\int }_{{\mathbb{R}}^d}\frac{u(x)-u(y)}{{\left|x-y\right|}^{d+2s}}dy\text{with}C(d,s):=-2^{2s}\frac{\Gamma (s+d/2)}{{\pi }^{d/2}\Gamma (-s)},$$

$$\alpha :=1-2s\in (-1,1)$$

$${\Vert U\Vert }_{L_{\alpha }^2(\omega )}^2:={\int }_{\omega }{\mathcal{Y}}^{\alpha }{\left|U(x,\mathcal{Y})\right|}^{2}dxd\mathcal{Y}$$

$$\mathrm{div}({\mathcal{Y}}^{\alpha }\nabla U)$$

$$\mathrm{tr}U$$

$$-d_s\underset{\mathcal{Y}\to 0^{+}}{\lim }{\mathcal{Y}}^{\alpha }{\partial }_{\mathcal{Y}}U(x,\mathcal{Y})=(-\Delta {)}^{s}u,d_s=2^{1-2s}\left|\Gamma (s)\right|/\Gamma (1-s).$$

$$(-\Delta {)}^{s}u$$

$$u$$

$$a(u,v):=\frac{C(d,s)}{2}\int {\int }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}\frac{(u(x)-u(y))(v(x)-v(y))}{{\left|x-y\right|}^{d+2s}}dxdy={\int }_{\Omega }fvdx\text{for all }v\in {\tilde{H}}^s(\Omega ).$$

$$h_{\ell }{|}_T:=h_{\ell }(T):=|T{|}^{1/d}\text{for all }T\in {\mathcal{T}}_{\ell }.$$

$$\displaystyle\Omega_{\ell}(T):=\operatorname*{interior}\left(\bigcup_{T^{\prime}\in\mathcal{T}_{\ell}(T)}\overline{T^{\prime}}\right)\subseteq\Omega,\quad\text{where}\quad\mathcal{T}_{\ell}(T):=\big{\{}T^{\prime}\in\mathcal{T}_{\ell}\,:\,\overline{T}\cap\overline{T^{\prime}}\neq\emptyset\big{\}}$$

$$S^{1,1}({\mathcal{T}}_{\ell }):=\{u\in H^1(\Omega ):u{|}_T\in {\mathcal{P}}^1(T)\text{ for all }T\in {\mathcal{T}}_{\ell }\}\text{and}{S}_0^{1,1}({\mathcal{T}}_{\ell }):=S^{1,1}({\mathcal{T}}_{\ell })\cap H_0^1(\Omega ).$$

$$a(u_{\ell },v_{\ell })={\int }_{\Omega }f{v}_{\ell }dx\text{for all }{v}_{\ell }\in S_0^{1,1}({\mathcal{T}}_{\ell }).$$

$${\omega }_{\ell }(x):=\underset{T\in {\mathcal{T}}_{\ell }}{\inf }\mathrm{dist}(x,\partial T)=\underset{T\in {\mathcal{T}}_{\ell }}{\inf }\underset{y\in \partial T}{\inf }\left|x-y\right|.$$

$$\displaystyle\eta_{\ell}(v_{\ell}):=\left\|\widetilde{h}_{\ell}^{s}\big{(}f-(-\Delta)^{s}v_{\ell}\big{)}\right\|_{L^{2}(\Omega)},\,\,\text{where}\,\,\,\widetilde{h}_{\ell}^{s}:=\begin{cases}h_{\ell}^{s}\quad&\text{for $0<s\leq 1/2$},\\ h_{\ell}^{s-\beta}\omega_{\ell}^{\beta}&\text{for $1/2<s<1$ and $\beta:=s-1/2$.}\\ \end{cases}$$

$$\displaystyle\eta_{\ell}(v_{\ell})=\left(\sum_{T\in\mathcal{T}_{\ell}}\eta_{\ell}(T,v_{\ell})^{2}\right)^{1/2}\!,\quad\text{where}\quad\eta_{\ell}(T,v_{\ell}):=\left\|\widetilde{h}_{\ell}^{s}\big{(}f-(-\Delta)^{s}v_{\ell}\big{)}\right\|_{L^{2}(T)}.$$

$${\Vert u-u_{\ell }\Vert }_{{\tilde{H}}^s(\Omega )}\le C_{\mathrm{rel}}{\eta }_{\ell }.$$

$${\eta }_{\ell }^2\le C_{\mathrm{eff}}^2\left({\Vert u-u_{\ell }\Vert }_{{\tilde{H}}^s(\Omega )}^2+\sum _{T\in {\mathcal{T}}_{\ell }}{h}_{\ell }(T){\Vert u-u_{\ell }\Vert }_{H^{1/2+s}({\Omega }_{\ell }^2(T))}^2\right).$$

$$\theta \sum _{T\in {\mathcal{T}}_{\ell }}{\eta }_{\ell }(T{)}^2\le \sum _{T\in {\mathcal{M}}_{\ell }}{\eta }_{\ell }(T{)}^2.$$

$${\eta }_{\ell +n}\le C_{\mathrm{lin}}{q}_{\mathrm{lin}}^n{\eta }_{\ell }\text{for all }\ell ,n\in {\mathbb{N}}_0.$$

$$c_{\mathrm{opt}}{\mathbb{A}}_t(u)\le \underset{\ell \in {\mathbb{N}}_0}{\sup }(\#{\mathcal{T}}_{\ell }{)}^t{\eta }_{\ell }\le C_{\mathrm{opt}}{\mathbb{A}}_t(u),$$

$${\mathbb{A}}_t(u):=\underset{N\in {\mathbb{N}}_0}{\sup }{N}^t\underset{{\mathcal{T}}_{\mathrm{opt}}\in \{\mathcal{T}\in \mathrm{refine}({\mathcal{T}}_0):\#\mathcal{T}-\#{\mathcal{T}}_0\le N\}}{\min }{\eta }_{\mathrm{opt}},$$

$$\displaystyle\left\|h_{\ell}^{s}\,(-\Delta)^{s}v\right\|_{L^{2}(\Omega)}\leq C_{\rm inv}\,\Big{(}\left\|v\right\|_{\widetilde{H}^{s}(\Omega)}^{2}+\sum_{T\in\mathcal{T}_{\ell}}h_{\ell}(T)^{2s}\left\|v\right\|_{H^{2s}(\Omega_{\ell}(T))}^{2}\Big{)}^{1/2}.$$

$${\Vert {\tilde{h}}_{\ell }^s(-\Delta {)}^s{v}_{\ell }\Vert }_{L^2(\Omega )}\le C_{\mathrm{inv}}{\Vert v_{\ell }\Vert }_{{\tilde{H}}^s(\Omega )}.$$

$$\left|{\left(\sum _{T\in {\mathcal{U}}_{\bullet }}{\eta }_{\circ }(T,v_{\circ }{)}^2\right)}^{1/2}-{\left(\sum _{T\in {\mathcal{U}}_{\bullet }}{\eta }_{\bullet }(T,w_{\bullet }{)}^2\right)}^{1/2}\right|\le C_{\mathrm{stab}}{\Vert v_{\circ }-w_{\bullet }\Vert }_{{\tilde{H}}^s(\Omega )}.$$

$${\left(\sum _{T\in {\mathcal{T}}_{\circ }\backslash {\mathcal{T}}_{\bullet }}{\eta }_{\circ }(T,v_{\bullet }{)}^2\right)}^{1/2}\le q_{\mathrm{red}}{\left(\sum _{T\in {\mathcal{T}}_{\bullet }\backslash {\mathcal{T}}_{\circ }}{\eta }_{\bullet }(T,v_{\bullet }{)}^2\right)}^{1/2}.$$

$${\Vert u_{\circ }-u_{\bullet }\Vert }_{{\tilde{H}}^s(\Omega )}\le C_{\mathrm{rel}}{\left(\sum _{T\in {\mathcal{T}}_{\bullet }\backslash {\mathcal{T}}_{\circ }}{\eta }_{\bullet }(T,u_{\bullet }{)}^2\right)}^{1/2}.$$

Full text

††Acknowledgement. The research of JMM and DP is funded by the Austrian Science Fund (FWF) by the special research program Taming complexity in PDE systems (grant SFB F65). Additionally, DP acknowledges support through the FWF research project Optimal adaptivity for BEM and FEM-BEM coupling (grant P27005).

1 Introduction

Fractional differential operators such as the fractional Laplacian $(-\Delta)^{s}$, $0<s<1$, are an increasingly important modeling tool in, e.g., physics, finance, or image processing. In contrast to classical (integer order) differential operators, these fractional operators are nonlocal, which makes both the analysis and the analysis of numerical methods challenging.

The numerical treatment of such non-local operators is currently a very active research field. A variety of approaches for the different versions of the fractional Laplacian in multi-d, or, more generally, fractional powers of differential operators are available: we mention Galerkin/finite element methods (FEMs) (see, e.g., [NOS15, AB17, ABH18, BMN*+*19] and references therein), techniques that exploit the connection of the fractional Laplacian with semigroup theory, [BP15, BP17], and techniques that rely on the connection with eigenfunction expansions, [SXK17, AG18]; for more details, we refer the reader to the recent surveys [BBN*+*18, LPG*+*18]. The vast majority of the numerical analysis literature focuses on a priori error analyses, and few results on a posteriori error analysis are available in spite of the fact that solutions to fractional differential equations typically have singularities (even for smooth input data), which naturally calls for using locally refined meshes; work on a posteriori error estimation in a Galerkin setting includes [NvPZ10, CNOS15, AG17, BBN*+*18, CNOS15] and on gradient recovery [ZHCK17]. One challenge in devising a posteriori error estimators are poor properties of the residual, namely, it is not necessarily in $L^{2}$. One can overcome this lack of $L^{2}$-regularity by measuring the residual in appropriate $L^{p}$-spaces, [NvPZ10, BBN*+*18] (some restrictions on $s$ apply); an alternative route, which is taken in the present work, is to measure the residual in weighted $L^{2}$-spaces, where the weight is given by a power of the distance from the mesh skeleton. The resulting a posteriori error estimator is shown to be reliable in Theorem 2.3 in the full range $0<s<1$. This a posteriori error estimator is a basic building block of the adaptive Algorithm 2.5 that we analyze. We show in Theorem 2.6 that it yields a sequence of approximations that converge at the optimal algebraic rate (with respect to an appropriate nonlinear approximation class). Such an optimal convergence result is well-known for adaptive FEMs for linear second-order elliptic equations (see, e.g., [BDD04, Ste07, CKNS08, FFP14]) or the classical BEM (see [Gan13, FKMP13, FFK*+*14, FFK*+*15]), and the present work extends these results to the integral fractional Laplacian. Our convergence result is obtained using the abstract, general framework of [CFPP14] for proving such optimal convergence results, which reduces the proof of optimal algebraic convergence to verifying four properties ((A1))—((A4)) of the error estimator (cf. Proposition 3.1). A key ingredient to establish the properties (A1)—(A4) in the case of the BEM are local inverse estimates for the nonlocal boundary integral operators, [Gan13, FKMP13, AFF*+*17]. Also in the present case of the fractional Laplacian underlying our analysis is the inverse estimate

[TABLE]

where $v_{\ell}$ is a piecewise polynomial and $\widetilde{h}_{\ell}$ is the local mesh width function (modified by some additional weight function for $s\geq 1/2$; cf. (2.13)), and is essentially sufficient for proving optimal convergence of the adaptive algorithm. In this work, such an inverse estimate for the nonlocal fractional Laplacian is provided in Theorem 2.7. We highlight that such an inverse estimate proves useful in other applications such as the analysis of multilevel preconditioning for the fractional Laplacian on locally refined meshes.

The present paper is structured as follows: In Section 2, we provide the continuous model problem and its Galerkin discretization by piecewise linears. Moreover, both the classical weighted residual a posteriori error indicators (for $0<s<1/2$) and our modified weighted error estimator (for $1/2\leq s<1$) are presented. The adaptive algorithm, the optimal convergence of the algorithm, as well as the inverse inequality, which plays the key role in our proofs, are stated. In Section 3, the four essential properties (A1)—(A4) of the error estimator from [CFPP14] are recalled and verified with the aid of the inverse inequality. The inverse inequality is then proved in Section 4. Finally, Section 5 provides numerical examples that illustrate the optimal convergence of the proposed adaptive method.

Concerning notation: For bounded, open sets $\omega\subset\mathbb{R}^{d}$ integer order Sobolev spaces $H^{t}(\omega)$, $t\in\mathbb{N}_{0}$, are defined in the usual way. For $t\in(0,1)$, fractional Sobolev spaces are given in terms of the seminorm $|\cdot|_{H^{t}(\omega)}$ and the full norm $\|\cdot\|_{H^{t}(\omega)}$ by

[TABLE]

where we denote the Euclidean distance in $\mathbb{R}^{d}$ by $\left|\;\cdot\;\right|$. Moreover, for bounded Lipschitz domains $\Omega\subset\mathbb{R}^{d}$, we define the spaces

[TABLE]

of $H^{t}$-functions with zero extension, equipped with the norm

[TABLE]

where $\rho(x)$ is the distance of a point $x\in\Omega$ to the boundary $\partial\Omega$.

2 Main Results

2.1 The fractional Laplacian and the Caffarelli-Silvestre extension

There are several different ways to define the fractional Laplacian $(-\Delta)^{s}$. A classical definition on the full space ${\mathbb{R}}^{d}$ is in terms of the Fourier transformation ${\mathcal{F}}$, i.e., $({\mathcal{F}}(-\Delta)^{s}u)(\xi)=|\xi|^{2s}({\mathcal{F}}u)(\xi)$. A consequence of this definition is the mapping property, (see, e.g., [BBN*+*18])

[TABLE]

where the Sobolev spaces $H^{t}(\mathbb{R}^{d})$, $t\in\mathbb{R}$, are defined in terms of the Fourier transformation, [McL00, (3.21)]. Alternative, equivalent definitions of $(-\Delta)^{s}$ exist, e.g., via semi-group or operator theory, [Kwa17]. For our purposes, a convenient representation of the fractional Laplacian is as the principal value integral

[TABLE]

where $\Gamma(\cdot)$ denotes the Gamma function. An important observation made in [CS07] is that the fractional Laplacian can be understood as a Dirichlet-to-Neumann operator of a degenerate elliptic PDE, the so-called extension problem on a half space. In order to describe this degenerated elliptic PDE, we need weighted Sobolev spaces. For

[TABLE]

and measurable subsets $\omega\subset\mathbb{R}^{d}\times\mathbb{R}^{+}$, we define the weighted $L^{2}$-norm

[TABLE]

and denote by $L^{2}_{\alpha}(\omega)$ the space of square-integrable functions with respect to the weight $\mathcal{Y}^{\alpha}$. The Caffarelli-Silvestre extension is conveniently described in terms of the Beppo-Levi space ${\mathcal{B}}^{1}_{\alpha}(\mathbb{R}^{d}\times\mathbb{R}^{+}):=\{U\in{\mathcal{D}}^{\prime}(\mathbb{R}^{d}\times\mathbb{R}^{+})\,|\,\nabla U\in L^{2}_{\alpha}(\mathbb{R}^{d}\times\mathbb{R}^{+})\}$. Elements of ${\mathcal{B}}^{1}_{\alpha}(\mathbb{R}^{d}\times\mathbb{R}^{+})$ are in fact in $L^{2}_{\rm loc}(\mathbb{R}^{d}\times\mathbb{R}^{+})$ and one can give meaning to their trace at $\mathcal{Y}=0$, which is denoted $\operatorname{tr}U$. Recalling $\alpha=1-2s$, one has in fact $\operatorname{tr}U\in H^{s}_{\rm loc}(\mathbb{R}^{d})$ (see, e.g., [MK18]). The extension problem by Caffarelli-Silvestre [CS07] reads as follows: For $u\in\widetilde{H}^{s}(\Omega)$ and $\alpha=1-2s$, let $U\in{\mathcal{B}}_{\alpha}^{1}(\mathbb{R}^{d}\times\mathbb{R}^{+})$ solve

[TABLE]

Then, the fractional Laplacian can be recovered as the Neumann data of the extension problem in the sense of distributions, [CS14, Thm. 3.1]:

[TABLE]

Remark 2.1**.**

*For $u\in H^{s}(\mathbb{R}^{d})$, the extension $U$ given by (2.4) can also be characterized as $U=\operatorname{argmin}\{V\in{\mathcal{B}}^{1}_{\alpha}(\mathbb{R}^{d}\times\mathbb{R}^{+})\,:\,\operatorname{tr}V=u\}$.  *

2.2 The model problem

For a bounded Lipschitz domain $\Omega\subset\mathbb{R}^{d}$, we consider the problem

[TABLE]

where $f$ is a given right-hand side. The fractional Laplacian $(-\Delta)^{s}u$ is defined by formula (2.2) for $x\in\Omega$. The weak formulation of (2.6) reads as follows [Kwa17, Thm. 1.1 (e),(g)]: Find $u\in\widetilde{H}^{s}(\Omega)$ such that

[TABLE]

Existence and uniqueness of $u\in\widetilde{H}^{s}(\Omega)$ follow from the Lax–Milgram lemma.

2.3 Discretization

Henceforth, we assume $\Omega$ to be polyhedral111This is not essential but allows us to work with affine elements.. Let $\mathcal{T}_{\ell}$ be a regular (in the sense of Ciarlet) triangulation of $\Omega$ consisting of open simplices that is $\gamma$-shape regular in the sense of \max_{T\in\mathcal{T}_{\ell}}\big{(}\operatorname*{diam}(T)/|T|^{1/d}\big{)}\leq\gamma<\infty. Here, $\operatorname*{diam}(T):=\sup_{x,y\in T}|x-y|$ denotes the Euclidean diameter of $T$, whereas $|T|$ is the $d$-dimensional Lebesgue volume. To ease notation, we define the piecewise constant mesh size function $h_{\ell}\in L^{\infty}(\Omega)$ by

[TABLE]

Moreover, for all elements $T\in\mathcal{T}_{\ell}$, we define the element patch

[TABLE]

consists of $T$ and all of its neighbors. \Omega_{\ell}^{k}(T):={\rm interior}\big{(}\bigcup\big{\{}T^{\prime}\in\mathcal{T}_{\ell}\,:\,\overline{T^{\prime}}\cap\overline{\Omega_{\ell}^{k-1}(T)}\neq\emptyset\big{\}}\big{)} is the $k$-th order patch of $T$ (and $\Omega_{\ell}^{1}(T):=\Omega_{\ell}(T)$). Later on, the index $\ell\in\mathbb{N}_{0}$ indicates the step of an adaptive algorithm, where the triangulations $\mathcal{T}_{\ell}$ are obtained by local mesh refinement based on newest vertex bisection (NVB). We refer, e.g., to [KPP13] for the NVB algorithm in $d=2$ or to [Ste08] for $d\geq 2$.

For $T\subseteq\Omega$, let $\mathcal{P}^{1}(T)$ denote the space of all affine functions on $T$. We define spaces of $\mathcal{T}_{\ell}$-piecewise affine and globally continuous functions

[TABLE]

For the discretization of (2.7), we consider the Galerkin method based on $S^{1,1}_{0}(\mathcal{T}_{\ell})$. The Lax–Milgram lemma provides existence and uniqueness of $u_{\ell}\in S^{1,1}_{0}(\mathcal{T}_{\ell})$ such that

[TABLE]

Throughout the present work, we will need a weight function that measures the distance from the mesh skeleton: For a mesh $\mathcal{T}_{\ell}$ we introduce

[TABLE]

2.4 A posteriori error estimation

Let $v_{\ell}\in S^{1,1}_{0}(\mathcal{T}_{\ell})$. For $0<s<3/4$, due to the mapping properties of the fractional Laplacian given in (2.1), we have $(-\Delta)^{s}v_{\ell}\in L^{2}(\Omega)$. For $3/4<s<1$, the function $(-\Delta)^{s}v_{\ell}$ is (generically) no longer in $L^{2}(\Omega)$ as it has singularities at the mesh skeleton. Its blow-up can be measured in terms of the distance from the mesh skeleton as has been noticed in [BBN*+*18] and is made precise in the following lemma, which is proved in Section 4 below.

Lemma 2.2**.**

*Let $\omega_{\ell}$ be given by (2.12). Given $0<s<1$, fix $\beta>2s-3/2$, e.g., $\beta:=s-1/2$. For any $v_{\ell}\in S^{1,1}_{0}(\mathcal{T}_{\ell})$, we then have $\omega_{\ell}^{\beta}(-\Delta)^{s}v_{\ell}\in L^{2}(\Omega)$. *

To control the discretization error of (2.11), we study the weighted residual error estimator

[TABLE]

Since the $L^{2}$-norm is local, the error estimator can be written as a sum of local contributions

[TABLE]

If $v_{\ell}=u_{\ell}$ is the solution of (2.11), we abbreviate $\eta_{\ell}:=\eta_{\ell}(u_{\ell})$ as well as $\eta_{\ell}(T):=\eta_{\ell}(T,u_{\ell})$ for all $T\in\mathcal{T}_{\ell}$.

Theorem 2.3**.**

For $0<s<1$, the weighted residual error estimator (2.13) is reliable:

[TABLE]

Moreover, for $0<s\leq 1/2$ with $s\neq 1/4$ and $u\in H^{1}_{0}(\Omega)$, the estimator is efficient in the sense that

[TABLE]

The constants $C_{\rm rel},C_{\rm eff}>0$ depend only on $\Omega$, $d$, $s$, and the $\gamma$-shape regularity of $\mathcal{T}_{\ell}$.

Remark 2.4**.**

*The efficiency result (2.16) is weaker than classical efficiency, where the sum on the right-hand side does not appear. However, as this term scales with the correct powers of the mesh width, there is not a large gap between (2.16) and classical efficiency, as for discrete functions $u$ the right-hand side of (2.16) can be bounded by the energy norm with an inverse estimate.  *

2.5 Adaptive mesh refinement

Based on the local contributions of the weighted residual error estimator (2.13), we consider the following standard approach for adaptive mesh refinement of the type SOLVE – ESTIMATE – MARK – REFINE, where the Dörfler criterion (2.17) from [Dör96] is used to select elements for refinement.

Algorithm 2.5**.**

*Input: Initial triangulation $\mathcal{T}_{0}$, adaptivity parameters $0<\theta\leq 1$, $C_{\rm mark}\geq 1$.

For all $\ell=0,1,2,\dots$, iterate the following steps (i)–(iv):*

1. (i)

Compute the solution $u_{\ell}\in S^{1,1}_{0}(\mathcal{T}_{\ell})$ of (2.11). 2. (ii)

Compute refinement indicators $\eta_{\ell}(T)=\eta_{\ell}(T,u_{\ell})$ from (2.14) for all $T\in\mathcal{T}_{\ell}$. 3. (iii)

Determine a set $\mathcal{M}_{\ell}\subseteq\mathcal{T}_{\ell}$ of, up to the multiplicative factor $C_{\rm mark}$, minimal cardinality such that

[TABLE] 4. (iv)

Generate the coarsest NVB refinement $\mathcal{T}_{\ell+1}:={\rm refine}(\mathcal{T}_{\ell},\mathcal{M}_{\ell})$ of $\mathcal{T}_{\ell}$ such that all marked elements $T\in\mathcal{T}_{\ell}$ have been bisected.

Analogous to the works [Ste07, CKNS08, FFP14, CFPP14] for FEM and BEM, the following Theorem 2.6 states linear convergence (2.18) of Algorithm 2.5 with optimal algebraic convergence rates (2.19).

Theorem 2.6**.**

Let $0<\theta\leq 1$ and $1\leq C_{\rm mark}\leq\infty$. Then, there exist constants $0<q_{\rm lin}<1$ and $C_{\rm lin}>0$ such that the sequence $(u_{\ell})_{\ell\in\mathbb{N}}$ generated by Algorithm 2.5 satisfies

[TABLE]

Moreover, if $0<\theta\ll 1$ is sufficiently small and $1\leq C_{\rm mark}<\infty$, and if the initial triangulation $\mathcal{T}_{0}$ satisfies the admissibility condition [Ste08, Section 4] for $d\geq 3$, then the error estimator converges with the best possible algebraic rate: For each $t>0$, there exist $c_{\rm opt}$, $C_{\rm opt}>0$ such that

[TABLE]

where ${\rm refine}(\mathcal{T}_{0})$ is the (infinite) set of all NVB refinements of the initial triangulation $\mathcal{T}_{0}$ and $\eta_{\rm opt}$ is the weighted residual error estimator corresponding to the triangulation $\mathcal{T}_{\rm opt}$.

2.6 Inverse estimates for the nonlocal operator $\boldsymbol{(-\Delta)^{s}}$

A key piece in the axiomatic approach of [CFPP14], which generalizes ideas and techniques developed in [FKMP13, Gan13, FFK*+*14, FFK*+*15] for the convergence analysis of adaptive algorithms involving nonlocal operators, are inverse inequalities for the pertinent operators. The following Theorem 2.7 establishes this inverse estimate for the fractional Laplacian. We mention that the case $s=1/2$ for the fractional Laplacian is closely related to the hyper-singular integral operator in the BEM [Gan13, FFK*+*15], for which the appropriate inverse estimate is established in [AFF*+*17].

Theorem 2.7**.**

1. (i)

Let $0<s\leq 1/2$ with $s\neq 1/4$. Then, there holds for $v\in\widetilde{H}^{2s}(\Omega)$

[TABLE] 2. (ii)

For $0<s<1$ and all $v_{\ell}\in S^{1,1}_{0}(\mathcal{T}_{\ell})$, there holds

[TABLE]

The constant $C_{\rm inv}>0$ depends only on $\Omega$, $d$, $s$, and the $\gamma$-shape regularity of $\mathcal{T}_{\ell}$.

3 Proof of Theorem 2.3 and Theorem 2.6

3.1 Axioms of adaptivity

The following Proposition 3.1 yields validity of the axioms of adaptivity from [CFPP14], where we recall that the present conforming discretization guarantees that $S^{1,1}_{0}(\mathcal{T}_{\bullet})\subseteq S^{1,1}_{0}(\mathcal{T}_{\circ})$ if the triangulation $\mathcal{T}_{\circ}$ is a refinement of $\mathcal{T}_{\bullet}$ (i.e., $\mathcal{T}_{\circ}$ is obtained from $\mathcal{T}_{\bullet}$ by finitely many steps of newest vertex bisection).

Proposition 3.1**.**

There exist constants $C_{\rm stab}$, $C_{\rm rel}>0$ and $0<q_{\rm red}<1$ depending solely on $\Omega$, $d$, $s$, and the $\gamma$-shape regularity of the initial triangulation $\mathcal{T}_{0}$ such that the following properties ((A1))–((A4)) hold for any refinement $\mathcal{T}_{\bullet}$ of the initial triangulation $\mathcal{T}_{0}$ and any refinement $\mathcal{T}_{\circ}$ of $\mathcal{T}_{\bullet}$:

1. (A1)

**Stability: **For all $v_{\circ}\in S^{1,1}_{0}(\mathcal{T}_{\circ}),w_{\bullet}\in S^{1,1}_{0}(\mathcal{T}_{\bullet})$ and any $\mathcal{U}_{\bullet}\subseteq\mathcal{T}_{\bullet}\cap\mathcal{T}_{\circ}$, there holds

[TABLE] 2. (A2)

**Reduction: **For all $v_{\bullet}\in S^{1,1}(\mathcal{T}_{\bullet})$, there holds

[TABLE] 3. (A3)

**Discrete reliability: **The Galerkin approximations $u_{\bullet}\in S^{1,1}_{0}(\mathcal{T}_{\bullet})$ and $u_{\circ}\in S^{1,1}_{0}(\mathcal{T}_{\circ})$ satisfy that

[TABLE] 4. (A4)

**Quasi-orthogonality: **For all $\ell$, $N>0$ and $\varepsilon>0$, it holds that

[TABLE]

Proof of ((A4)).

The quasi-orthogonality follows directly from the fact that the bilinear form $a(\cdot,\cdot)$ of the fractional Laplacian (2.7) is bilinear, elliptic, and symmetric; see [CFPP14, Section 3.6] for details. ∎

Proof of ((A1)).

Let \omega:=\operatorname*{interior}\big{(}\bigcup_{T\in\mathcal{U}_{\bullet}}\overline{T}\big{)}\subseteq\Omega. With Theorem 2.7, we obtain that

[TABLE]

This proves ((A1)) with $C_{\rm stab}=C_{\rm inv}$. ∎

Proof of ((A2)).

Bisection ensures that $|T^{\prime}|\leq|T|/2$ for all $T\in\mathcal{T}_{\bullet}\backslash\mathcal{T}_{\circ}$ and its children/descendants $T^{\prime}\in\mathcal{T}_{\circ}\backslash\mathcal{T}_{\bullet}$ with $T^{\prime}\subset T$. Note that $\omega:=\bigcup_{T^{\prime}\in\mathcal{T}_{\circ}\backslash\mathcal{T}_{\bullet}}\overline{T^{\prime}}=\bigcup_{T\in\mathcal{T}_{\bullet}\backslash\mathcal{T}_{\circ}}\overline{T}$. For $0<s\leq 1/2$, this proves ((A2)) with $q_{\rm red}=2^{-s/d}$, since

[TABLE]

For $1/2<s<1$, we note that $0<\beta=s-1/2<s$ and $s-\beta>0$. Moreover, we have pointwise $\omega_{\circ}(x)\leq\omega_{\bullet}(x)$. Hence, it follows that

[TABLE]

Arguing as before, we prove ((A2)) with $q_{\rm red}=2^{-(s-\beta)/d}$. This concludes the proof. ∎

The proof of discrete reliability ((A3)) relies on the Scott–Zhang projection [SZ90] for quasi-interpolation in $\widetilde{H}^{s}(\Omega)$. While the original work [SZ90] is concerned with the integer-order Sobolev space $H^{1}(\Omega)$, the approach is generalized in [AFF*+*15] to fractional-order Sobolev spaces $\widetilde{H}^{s}(\Omega)$ for $0\leq s\leq 1$. Since the precise construction will matter, we briefly sketch it: Let $\mathcal{N}_{\bullet}$ be the set of all nodes of $\mathcal{T}_{\bullet}$ and let $\mathcal{N}_{\bullet}^{\rm int}$ be the set of nodes of $\mathcal{T}_{\bullet}$, which lie inside $\Omega$ (and not on the boundary $\Gamma$). For each element $T\in\mathcal{T}_{\bullet}$, let $\{\phi^{\star}_{1,T},\ldots,\phi^{\star}_{d,T}\}\subset\mathcal{P}^{1}(T)$ be the dual basis for the nodal basis of $\mathcal{P}^{1}(T)$ with respect to $(\cdot,\cdot)_{L^{2}(T)}$. these functions may be viewed as elements of $S^{1,0}(\mathcal{T}_{\bullet})$ by zero extension outside $T$. For each $z\in\mathcal{N}_{\bullet}$ choose an arbitrary element $T_{z}\in\mathcal{T}_{\bullet}$ with $z\in\overline{T_{z}}$. (This freedom will be exploited later on.) Let $\phi_{z}\in S^{1,1}(\mathcal{T}_{\bullet})$ denote the hat function associated with $z$. For the element $T_{z}$, let $\phi^{\star}_{z}\in\{\phi^{\star}_{1,T_{z}},\ldots,\phi^{\star}_{d,T_{z}}\}\subset\mathcal{P}^{1}(T_{z})$ be such that $\int_{T_{z}}\phi^{\star}_{z}\phi_{z^{\prime}}\,dx=\delta_{zz^{\prime}}$ for all $z^{\prime}\in\mathcal{N}_{\bullet}$. Then, the Scott–Zhang projector $J_{\bullet}:L^{2}(\Omega)\to S^{1,1}_{0}(\mathcal{T}_{\bullet})$ is defined by

[TABLE]

Since the summation is over the interior nodes $\mathcal{N}_{\bullet}^{\rm int}$, we have $J_{\bullet}v\in S^{1,1}_{0}(\mathcal{T}_{\bullet})$. The following lemma collects some of its properties:

Lemma 3.2**.**

For $0<s<1$, there holds for some constant $C_{\rm SZ}>0$ depending only on $s$ and the $\gamma$-shape regularity of $\mathcal{T}_{\bullet}$ and $\Omega$:

1. (i)

$J_{\bullet}V_{\bullet}=V_{\bullet}$*  as well as  $\left\|J_{\bullet}v\right\|_{\widetilde{H}^{s}(\Omega)}\leq C_{\rm SZ}\,\left\|v\right\|_{\widetilde{H}^{s}(\Omega)}$  for all $V_{\bullet}\in S^{1,1}_{0}(\mathcal{T}_{\bullet})$ and all $v\in\widetilde{H}^{s}(\Omega)$.* 2. (ii)

$h_{\bullet}(T)^{-1}\,\left\|(1-J_{\bullet})v\right\|_{L^{2}(T)}+\left\|\nabla(1-J_{\bullet})v\right\|_{L^{2}(T)}\leq C_{\rm SZ}\,\left\|\nabla v\right\|_{L^{2}(\Omega_{\bullet}(T))}$*  for all $v\in H^{1}_{0}(\Omega)$.* 3. (iii)

$\left\|\widetilde{h}_{\bullet}^{-s}\,(1-J_{\bullet})v\right\|_{L^{2}(\Omega)}\leq C_{\rm SZ}\,\left\|v\right\|_{\widetilde{H}^{s}(\Omega)}$*  for all $v\in\widetilde{H}^{s}(\Omega)$.* 4. (iv)

For each refinement $\mathcal{T}_{\circ}$ of $\mathcal{T}_{\bullet}$, the Scott-Zhang operator $J_{\bullet}$ can be chosen such that additionally

[TABLE]

Proof.

(i)–(ii) are proved in [AFF*+*15, Lemma 4]. With $1/\widetilde{h}_{\bullet}^{s}$ replaced by $h_{\bullet}^{-s}$, (iii) is proved in [FFK*+*15, Lemma 3.2], where the estimate

[TABLE]

is established for any $0<t<1$ as well. Hence, it only remains to prove (iii) for $1/2<s<1$, where $\widetilde{h}_{\bullet}^{-s}=h_{\bullet}^{\beta-s}\omega_{\bullet}^{-\beta}$.

With $\widehat{T}$ being the reference element, let $\omega_{\widehat{T}}(\widehat{x}):=\operatorname*{dist}(\widehat{x},\partial\widehat{T})$ be the corresponding weight function. Let $\Phi_{T}:\widehat{T}\to T$ be an affine parametrization of $T\in\mathcal{T}_{\bullet}$. For each $\widehat{x}\in\widehat{T}$, let $x:=\Phi_{T}(\widehat{x})$. Let $x_{\rm min}\in\partial T$ with $\left|x-x_{\rm min}\right|=\operatorname*{dist}(x,\partial T)$. Then,

[TABLE]

and consequently

[TABLE]

According to [Gri85, Thm. 1.4.4.3], we have $\left\|\omega_{\widehat{T}}^{-\beta}\widehat{w}\right\|_{L^{2}(\widehat{T})}\lesssim\left\|\widehat{w}\right\|_{H^{\beta}(\widehat{T})}$ for all $\widehat{w}\in H^{\beta}(\widehat{T})$, since $0<\beta<1/2$. Let $w:=(1-J_{\bullet})v$ and $\widehat{w}:=w\circ\Phi_{T}$. A scaling argument thus proves that

[TABLE]

Hence, by summation over all elements

[TABLE]

This concludes the proof of (iii).

We prove (iv). To ensure the additional property (3.1), we use the freedom still available in the definition of the Scott-Zhang operator $J_{\bullet}$. To that end, let $\mathcal{N}_{\bullet\cap\circ}$ be the nodes of the elements in $\mathcal{T}_{\bullet}\cap\mathcal{T}_{\circ}$. For each $z\in\mathcal{N}_{\bullet\cap\circ}\subset\mathcal{N}_{\bullet}$, select an element $T_{z}\in\mathcal{T}_{\bullet}\cap\mathcal{T}_{\circ}$ such $z$ is a node of $T_{z}$; for the remaining nodes $z\in\mathcal{N}_{\bullet}\setminus\mathcal{N}_{\bullet\cap\circ}$, the element $T_{z}\in\mathcal{T}_{\bullet}$ merely needs to be such that $z$ is a node of $T_{z}$. This defines the Scott-Zhang operator $J_{\bullet}$. This particular choice ensures

[TABLE]

and (3.1) follows from (3.3). ∎

Proof of ((A3)).

Since $a(\cdot,\cdot)$ is elliptic on $\widetilde{H}^{s}(\Omega)$, we have

[TABLE]

In particular, the Cauchy–Schwarz inequality and Lemma 3.2 (iv) show with $\omega:=\bigcup_{T\in\mathcal{T}_{\bullet}\cap\mathcal{T}_{\circ}}\overline{T}$

[TABLE]

Lemma 3.2 (iii) proves that

[TABLE]

The combination of the last two estimates concludes the proof. ∎

3.2 Proof of Theorem 2.3

For reliability (2.15), we sketch the proof from [CFPP14, Section 3.3]. Let $\mathcal{T}_{\bullet}$ be a given triangulation. Recall that uniform refinement leads to convergence. Given $\varepsilon>0$, we may hence choose a refinement $\mathcal{T}_{\circ}$ of $\mathcal{T}_{\bullet}$ such that $\left\|u-u_{\circ}\right\|_{\widetilde{H}^{s}(\Omega)}\leq\varepsilon$. Therefore, the triangle inequality and ((A3)) prove that

[TABLE]

As $\varepsilon\to 0$, we prove (2.15).

To prove the weak efficiency (2.16), we employ the inverse estimates (2.21)–(2.22) from Theorem 2.7. Since $0<s\leq 1/2$ with $s\neq 1/4$, this yields that

[TABLE]

First, Lemma 3.2 (i) guarantees that

[TABLE]

Next, we observe that $0<2s<1$ with $2s\neq 1/2$. Lemma 3.2 (ii) and an interpolation argument yield that

[TABLE]

Since newest vertex bisection is used to generate $\mathcal{T}_{\bullet}$, only a finite number of shapes of patches can occur. In particular, the hidden constant in the last estimate does not depend on $\mathcal{T}_{\bullet}$, but only on $\mathcal{T}_{0}$. This concludes the proof. ∎

3.3 Proof of Theorem 2.6

We note that reduction from [CFPP14] is a consequence of conformity and ((A1))–((A2)), since

[TABLE]

Therefore, Theorem 2.6 immediately follows from [CFPP14, Theorem 4.1].∎

4 Proof of Theorem 2.7

4.1 Proof of Lemma 2.2

For $x\in T$, we split the fractional Laplacian into a principal value part and a smoother, integrable part

[TABLE]

Using polar coordinates $y=x+r\nu$, $\nu\in S^{d-1}$, where $S^{d-1}$ is the $(d-1)$-dimensional unit sphere, and exploiting that $u_{\ell}|_{T}\in\mathcal{P}^{1}(T)$, we may compute the principal value part

[TABLE]

where the last equality follows from interchanging the integration in $\nu$ and $r$. For the second part in the decomposition of $(-\Delta)^{s}u_{\ell}$ in (4.1) a similar computation using the Lipschitz continuity of $u_{\ell}$ provides

[TABLE]

Since $\operatorname*{dist}(x,\partial T)^{\beta+1-2s}=(\omega_{\ell}|_{T})^{\beta+1-2s}$ is square-integrable in view of $\beta>2s-3/2$, Lemma 2.2 follows. ∎

4.2 Proof of Theorem 2.7

As in [AFF*+*17], which provides inverse estimates for the classical boundary integral operators, the main idea of the proof of Theorem 2.7 is a splitting of the operator into a smoother far-field part and a localized near-field part. The near-field and the far-field are treated separately in the following subsections.

The proofs for both statements (2.21)–(2.22) are fairly similar and only differ in the use of the inverse inequality of Lemma 4.2 below.

Let $v\in\widetilde{H}^{2s}(\Omega)$. We start by a localization and a splitting into near-field and far-field. The $L^{2}$-norm in the error estimator can be written as a sum over all elements

[TABLE]

For any constant $c_{T}\in\mathbb{R}$, we have that $(-\Delta)^{s}c_{T}\equiv 0$, which implies $(-\Delta)^{s}v=(-\Delta)^{s}(v-c_{T})$.

Due to the nonlocality of the fractional Laplacian, we need to split the operator into two contributions, a localized near-field part and a smoother far-field part. To this end, we select, for each element $T\in\mathcal{T}_{\ell}$, a cut-off function $\chi_{T}\in C_{0}^{\infty}(\mathbb{R}^{d})$ with the following properties: i) $\operatorname{supp}\chi_{T}\subset\Omega_{\ell}(T)$, where $\Omega_{\ell}(T)$ is the patch of $T$ defined in (2.9); ii) there is a set $B^{\prime}$ with $T\subset B^{\prime}\subset\Omega_{\ell}(T)$ and $\operatorname*{dist}(T,\partial B^{\prime})\sim h_{\ell}(T)$; iii) $\chi_{T}\equiv 1$ in $B^{\prime}$; iv) $\left\|\chi_{T}\right\|_{W^{1,\infty}(\Omega_{\ell}(T))}\leq Ch_{\ell}(T)^{-1}$ as well as $\operatorname*{dist}(\operatorname*{supp}\chi_{T},\partial\Omega_{\ell}(T))>ch_{\ell}(T)$ with constants $C$, $c>0$ depending only on the $\gamma$-shape regularity of $\mathcal{T}_{\ell}$ and $d$.

Then, for each element $T\in\mathcal{T}_{\ell}$, we have the splitting $(-\Delta)^{s}(v-c_{T})=(-\Delta)^{s}(\chi_{T}(v-c_{T}))+\linebreak(-\Delta)^{s}((1-\chi_{T})(v-c_{T}))$ into the near-field $v_{\rm near}^{T}$ and a far-field $v_{\rm far}^{T}$ given by

[TABLE]

We choose

[TABLE]

If $v=v_{\ell}\in S^{1,1}_{0}(\mathcal{T}_{\ell})$, we write $v_{\rm near,\ell}^{T}$, $v_{\rm near,\ell}^{T}$ for the near-field and far-field to emphasize that the fields are related to discrete functions. For the near-field for $v\in\widetilde{H}^{2s}(\Omega)$ with $0<s\leq 1/2$ and $s\neq 1/4$, Lemma 4.2 below provides the estimate

[TABLE]

and for discrete $v_{\ell}\in S^{1,1}_{0}(\mathcal{T}_{\ell})$ and $0<s<1$

[TABLE]

For any $v\in\widetilde{H}^{s}(\Omega)$, Lemma 4.4 gives

[TABLE]

Combining (4.7), (4.9), we prove the inverse inequality (2.21); the estimate (2.22) is obtained from (4.8) and (4.9). ∎

4.2.1 The near-field

In this subsection, we treat the near-field $v_{\rm near}^{T}:=(-\Delta)^{s}(\chi_{T}(v-c_{T}))$. We start with a Poincaré inequality:

Lemma 4.1**.**

Let $c_{T}$ be given by (4.6) and $v\in\widetilde{H}^{s}(\Omega)$. Then, for a constant $C>0$ depending solely on $\Omega$, $s$, the $\gamma$-shape regularity of $\mathcal{T}$, and the fact that NVB is used, we have

[TABLE]

Proof.

If $\overline{T}\cap\partial\Omega=\emptyset$, then the Poincaré inequality (4.10) is shown in [AFF*+*17, Lemma 5.1] and [AB17]. For elements at the boundary, we have chosen $c_{T}=0$. Due to the boundary condition satisfied by $v$, we have the Poincaré inequality $h_{\ell}(T)^{-2s}\left\|v\right\|_{L^{2}(\Omega_{\ell}(T))}^{2}\lesssim\left|v\right|_{H^{s}(\Omega_{\ell}(T))}^{2}$. To see this, one makes four observations: a) the power $h_{\ell}^{-2s}$ is obtained by scaling; b) $\overline{T}\cap\partial\Omega\neq\emptyset$ implies that at least one face of $\partial\Omega_{\ell}(T)$ lies on $\partial\Omega$; c) only finitely many shapes are possible for $\Omega_{\ell}(T)$ since the meshes are created by NVB; d) the Poincaré inequality for the patch $\Omega_{\ell}(T)$ scaled to size $1$ holds with constant $O(1)$ by the compactness of the embedding of $H^{s}$ in $L^{2}$. ∎

Lemma 4.2**.**

There is a constant $C>0$ depending only on $\Omega$, $d$, $s$, and the $\gamma$-shape regularity of $\mathcal{T}_{\ell}$ such that the following holds:

1. (i)

For $v\in\widetilde{H}^{2s}(\Omega)$ with $0<s\leq 1/2$ and $s\neq 1/4$, let the near-field $v_{\rm near}^{T}$ be given by (4.4). Then,

[TABLE] 2. (ii)

For $v_{\ell}\in S^{1,1}_{0}(\mathcal{T}_{\ell})$ and arbitrary $0<s<1$, the near-field $v_{\rm near,\ell}^{T}$ satisfies that

[TABLE]

Proof.

Before proving this lemma, we point out that we will frequently use in sums the shorthand $T^{\prime}\in\Omega_{\ell}(T)$ to mean that the summation is over all $T^{\prime}\in\mathcal{T}_{\ell}$ satisfying $T^{\prime}\subset\Omega_{\ell}(T)$.

Proof of (i): The mapping properties (2.1) of the fractional Laplacian as a pseudo-differential operator of order $2s$ as well as the additional assumption $v\in\widetilde{H}^{2s}(\Omega)\subset H^{2s}(\mathbb{R}^{d})$ imply

[TABLE]

where the last inequality follows from a Hardy inequality that is valid provided that $s\neq 1/4$, see, e.g., [Gri85, Thm. 1.4.4.4].

For $s=1/2$, we have the local $H^{1}$-norm on the right-hand side of (4.11), and we can directly estimate this by $\left\|v\right\|_{H^{1}(\Omega_{\ell}(T))}^{2}$. We note that this case coincides with the result for the hypersingular integral operator in the boundary element method, [AFF*+*17].

The $H^{2s}$-norm is nonlocal for $s\neq 1/2$, but it has a localized upper-bound, cf. [Fae00, Fae02],

[TABLE]

Step 1: (Estimate of $S_{3}$) As $\operatorname*{supp}\chi_{T}\subset\Omega_{\ell}(T)$, the last term sums up to a $L^{2}$-norm on the patch $\Omega_{\ell}(T)$. Since the size of neighboring elements differ only by a constant multiplicative factor, we can estimate $h_{\ell}(T^{\prime})^{-2s}\lesssim h_{\ell}(T)^{-2s}$ for all $T^{\prime}\in\Omega_{\ell}(T)$. The Poincaré inequality on the patch $\Omega_{\ell}(T)$ given in Lemma 4.1 then gives

[TABLE]

Step 2: (Estimate of $S_{1}$) To estimate $S_{1}$, we use that $\operatorname*{supp}\chi_{T}\subset\Omega_{\ell}(T)$. Therefore, only elements $T^{\prime}\in\Omega_{\ell}(T)$ appear. With a hidden constant depending on the $\gamma$-shape regularity of $\mathcal{T}_{\ell}$, this leads to

[TABLE]

The first term is just the $H^{s}(\Omega_{\ell}(T))$-seminorm so that we may concentrate on the second term. Since $\operatorname*{dist}(\operatorname*{supp}\chi_{T},\partial\Omega_{\ell}(T))>ch_{\ell}(T)$, the denominator in the second integrand can be directly estimated by powers of $h_{\ell}(T)$. For arbitrary $T^{\prime}\in\Omega_{\ell}(T)\backslash T$, this immediately leads to

[TABLE]

After summation over $T^{\prime}$, the desired estimate follows from the Poincaré inequality of Lemma 4.1.

Step 3: (Estimate of $S_{2}$) Since the integrand vanishes if $T^{\prime}\notin\Omega_{\ell}(T)$, we may split the sum as in Step 2 into a smooth part and a part defined on the patch $\Omega_{\ell}(T)$, i.e.,

[TABLE]

Using the Lipschitz continuity of $\chi_{T}$, we obtain

[TABLE]

Since we assume $s<1/2$, a direct calculation reveals

[TABLE]

so that $S_{2,1}$ can be estimate in the required way. To estimate $S_{2,2}$ one uses again the Lipschitz continuity of $\chi_{T}$, the observation $\operatorname*{dist}(\operatorname*{supp}\chi_{T},\partial\Omega_{\ell}(T))>ch_{\ell}(T)$, as well as a Poincaré inequality of Lemma 4.1 so that $S_{2,2}$ can be estimated using the same arguments as in (4.2.1). Putting the estimates for $S_{1}$, $S_{2}$, and $S_{3}$ into (4.12) and summing over all elements using the $\gamma$-shape regularity of $\mathcal{T}_{\ell}$, we obtain

[TABLE]

Proof of (ii) in the case $0<s\leq 1/2$: For $2s\leq 1$ and $v_{\ell}\in S^{1,1}_{0}(\mathcal{T}_{\ell})$ we have the classical inverse estimate (e.g., [SS11, Thm. 4.4.2])

[TABLE]

Hence, for $s\in(0,1/2]\setminus\{1/4\}$ and a $v_{\ell}\in S^{1,1}_{0}(\mathcal{T}_{\ell})$, we may directly combine the result of (i) with the inverse estimate (4.17) to get the desired estimate. To obtain the case $s=1/4$, we note that in the above proof of (i), the assumption $s\neq 1/4$ entered only through the estimate (4.12). However, the bound $\|\chi_{T}(v_{\ell}-c_{T})\|_{\widetilde{H}^{2s}(\Omega)}\lesssim\|\chi_{T}(v_{\ell}-c_{T})\|_{H^{2s}(\Omega)}$ is still valid for $v_{\ell}\in S^{1,1}_{0}(\mathcal{T}_{\ell})$.

Proof of (ii) in the case $1/2<s<1$: Here, we cannot use the mapping properties of the fractional Laplacian since $v_{\ell}$ may not be in $H^{2s}(\Omega)$ for $s\geq 3/4$. However, we can directly estimate the fractional Laplacian using that due to $s>1/2$ the function $r^{-2s+1}$ vanishes at infinity, which does not hold for the case $s<1/2$. We write

[TABLE]

The definition of the fractional Laplace leads to

[TABLE]

We treat the integrals $S_{4}$, $S_{5}$ on the right-hand side separately, starting with $S_{5}$. We split the integral $S_{5}$ further into two parts, a principal value integral containing the singularity and a smoother, integrable part:

[TABLE]

In fact, $S_{5,1}=0$. To see this, we use that, for $x$, $y\in T$, we have $v_{\ell}(x)-v_{\ell}(y)=\nabla v_{\ell}|_{T}\cdot(x-y)$, where $\nabla v_{\ell}|_{T}\in\mathbb{R}^{d}$ is constant. Using polar coordinates, the principal value integral can be computed and is equal to zero as in (4.1). To bound $S_{5,2}$ we note that our assumption $s>1/2$ implies that $r^{-2s+1}$ vanishes at infinity. With the Lipschitz continuity of $v_{\ell}$, we estimate $S_{5,2}$ as in (4.1) as

[TABLE]

where the last two estimates follow from a direct computation of the integral over the distance function and an inverse inequality, where we used that the ratio of neighboring elements is bounded by a constant.

It remains to estimate $S_{4}$ in (4.2.1), which is similar to the integral above. Since $\chi_{T}\equiv 1$ on $T$, we get

[TABLE]

Using the Lipschitz continuity of $\chi_{T}$, polar coordinates, and the Poincaré inequality, we therefore may estimate as in (4.2.1)

[TABLE]

Putting the bounds for $S_{4}$ and $S_{5}$ together and summing over all elements, leads to the claimed inverse estimate. ∎

4.2.2 The far-field

In order to treat the far-field part in the proof of the inverse estimate of Theorem 2.7, we utilize the interpretation of the fractional Laplacian as the Dirichlet-to-Neumann map for the extension problem (2.4). This leads us to the problem of controlling second derivatives of the solution to the extension problem. This is done in the following Lemma 4.3, which is proved with the techniques of tangential difference quotients, typical for elliptic PDEs, see, e.g., [Eva98].

Lemma 4.3**.**

Let $\zeta\in C^{\infty}_{0}(\mathbb{R}^{d+1})$ and $B\subset B^{\prime}\subset\mathbb{R}^{d}\times\mathbb{R}^{+}$ with $\zeta|_{B}\equiv 1$ and $\operatorname*{supp}\zeta\cap\mathbb{R}^{d}\times\mathbb{R}^{+}\subset B^{\prime}$. Let $U\in{\mathcal{B}}^{1}_{\alpha}(\mathbb{R}^{d}\times\mathbb{R}^{+})$ satisfy (2.4a) and $\operatorname{tr}(U\zeta)=0$. Then,

[TABLE]

Proof.

With the unit-vector $e_{x_{j}}$ in the $j$-th coordinate and $\tau>0$, we define the difference quotient

[TABLE]

We use the test-function $V=D_{x_{j}}^{-\tau}(\zeta^{2}D_{x_{j}}^{\tau}U)$ in the weak formulation of (2.4). Noting that $V|_{\partial(\Omega\times\mathbb{R}^{+})}=0$ due to the assumption $\operatorname*{tr}(U\zeta)=0$, we therefore obtain

[TABLE]

Therefore, with Young’s inequality, we have the estimate

[TABLE]

Absorbing the first term and taking the limit $\tau\rightarrow 0$, we obtain the sought inequality for all second derivatives in $x$. ∎

The following Lemma 4.4 provides the required estimate for the far-field.

Lemma 4.4**.**

Let $v\in\widetilde{H}^{s}(\Omega)$ and the far-field $v_{\rm far}^{T}$ be given by (4.5). Let $\beta=0$ for $0<s<1/2$ and $\beta=s-1/2$ for $1/2<s<1$. Then, the estimate

[TABLE]

holds with a constant depending only on $\Omega$, $d$, $s$, and the $\gamma$-shape regularity of $\mathcal{T}_{\ell}$.

Proof.

Since $\beta\geq 0$, we have $\widetilde{h}^{s}\lesssim h^{s}$ so that it suffices to show the estimate (4.21) with $\widetilde{h}^{s}$ replaced with $h^{s}$.

Since the fractional Laplacian is the Dirichlet-to-Neumann map for the extension problem (2.4), we need to control the Neumann data of the extension problem, i.e., the trace of $z(x,\mathcal{Y}):=\mathcal{Y}^{\alpha}\partial_{\mathcal{Y}}U(x,\mathcal{Y})$ at $\mathcal{Y}=0$, where $U$ is the solution of (2.4) with boundary data $u=(1-\chi_{T})v+\chi_{T}c_{T}$. We note that the definition of the cut-off function $\chi_{T}$ and the constant $c_{T}$ (in particular: $c_{T}=0$ if $\overline{T}\cap\partial\Omega\neq\emptyset$) imply $v\in\widetilde{H}^{s}$. Moreover, $(-\Delta)^{s}c_{T}\equiv 0$ implies that $v_{\rm far}^{T}=(-\Delta)^{s}u$.

Given $T\in\mathcal{T}_{\ell}$, the $\gamma$-shape regularity of $\mathcal{T}_{\ell}$ implies the existence of sets $T\subset B_{0}\subset B_{0}^{\prime}\subset\Omega_{\ell}(T)$ with $\operatorname{dist}(B_{0},\partial B_{0}^{\prime})\sim h_{\ell}(T)$ and implied constant depending solely on the $\gamma$-shape regularity. Additionally, we may require that $B_{0}^{\prime}\subset\{\chi_{T}\equiv 1\}$. Setting $B:=B_{0}\times(0,h_{\ell}(T))$, $B^{\prime}:=B_{0}^{\prime}\times(0,2h_{\ell}(T))$, we may select a cut-off function $\zeta\in C^{\infty}_{0}(\mathbb{R}^{d+1})$ with $\zeta\equiv 1$ on $B$, $\operatorname*{supp}\zeta\cap\mathbb{R}^{d}\times\mathbb{R}^{+}\subset B^{\prime}$ and $\|\nabla\zeta\|_{L^{\infty}(\mathbb{R}^{d+1})}\lesssim h_{\ell}(T)^{-1}$, where again the implied constant depends solely on the $\gamma$-shape regularity of $\mathcal{T}_{\ell}$. The multiplicative trace inequality from [MK18, Lemma 3.7] applied with $-\alpha$ provides

[TABLE]

The equation $-\operatorname*{div}(\mathcal{Y}^{\alpha}\nabla U)=0$ implies that $-\partial_{\mathcal{Y}}z=\mathcal{Y}^{\alpha}\Delta_{x}U$, so we obtain with $s=(1-\alpha)/2$

[TABLE]

It remains to control $\Delta_{x}U$ in the weighted $L^{2}$-norm. This is done with Lemma 4.3, and we get

[TABLE]

With a standard a priori estimate, the energy norm on the right-hand side can be estimated by the Dirichlet data $(1-\chi_{T})v+\chi_{T}c_{T}$, and this finally implies

[TABLE]

where the last estimate can be found in the proof of Lemma 4.2. This proves the estimate for the far-field. ∎

5 Numerics

We illustrate our theoretical results of the previous sections with some numerical examples in two dimensions. For more numerical examples about adaptive methods for fractional diffusion, we refer to [AG17].

5.1 Implementational Issues

We implement the crucial steps SOLVE (step (i)) and ESTIMATE (step (ii)) in Algorithm 2.5 in MATLAB R2018a as follows:

• •

SOLVE: To obtain the lowest-order discrete solution $u_{\ell}\in S^{1,1}_{0}(\mathcal{T}_{\ell})$, we use the existing MATLAB-code from [ABB17], where the unbounded domain $\mathbb{R}^{d}$ is replaced by a (large) circle around the computational domain $\Omega$. The integrals of the system matrix are computed with high precision quadrature rules. For domains $\Omega$ with curved boundary $\partial\Omega$, the boundary is approximated by a piecewise linear interpolant, which introduces a variational crime. Nevertheless, this approximation improves as the mesh size near $\partial\Omega$ is reduced.

• •

ESTIMATE: For the error estimator, we need to compute the local contributions

[TABLE]

with $\beta=0$ for $0<s\leq 1/2$ and $\beta=s-1/2$ for $1/2<s<1$. In Section 2, we noted that $(-\Delta)^{s}u_{\ell}(x)$ may blow up as $x$ tends to the mesh skeleton. The integral on the triangle $T$ is transformed to an integral on a square by means of the Duffy transformation. There, the integral is approximated by a quadrature on a tensor product composite Gauss rule on meshes that are refined geometrically towards the edges of the square. Evaluating the residual requires evaluating the fractional Laplacian applied to the discrete solution in each quadrature point. For $x\in T$, this is done with the pointwise evaluation formula from [AG17, Lemma 4]

[TABLE]

Here, only integrals over the boundary of each element appear, which are approximated using the standard 1D-MATLAB quadrature function.

As two contributions $\eta_{\ell}(T_{1})$, $\eta_{\ell}(T_{2})$ are independent for $T_{1}\neq T_{2}$, the computation of the error estimator is easily parallelized, which leads to considerable speed-up in the computations.

Remark 5.1**.**

The computation of the error estimator is fairly expensive due to the need for an accurate quadrature rule for the residual. One way to circumvent this is to use a different error estimator inspired by the near-field and far-field splitting from Section 4, as well as the observation that the near-field behaves like a power of the distance to the skeleton. We illustrate our ideas in the one dimensional case: As in formula (4.1), the principal value integral over one element $T=[x_{0},x_{1}]$ can be computed exactly for $x\in T$ as

[TABLE]

For all other elements $T^{\prime}\in\mathcal{T}_{\ell}\backslash\{T\}$ a similar computation can be made. However, only elements in the patch of $T$ are truly of interest, since the integrand is smooth for all other elements. For the elements in the patch, again, only the terms containing $x_{0},x_{1}$ are relevant, which leads to the splitting

[TABLE]

where $[u_{\ell}^{\prime}]_{T,T^{\prime}}$ denotes the jump of the derivative across the elements $T$ and $T^{\prime}$. Similarly as in Theorem 2.3, the residual can then be bounded by

[TABLE]

*The advantage of this approach is that the far-field, written in terms of formula (• ‣ 5.1), can be cheaply computed using standard Gaussian quadrature and for the near-field only the jumps across the elements need to be computed.  *

We present convergence plots in the energy norm, where we use that — due to the Galerkin orthogonality and symmetry of $a(\cdot,\cdot)$ — the error can be computed as

[TABLE]

5.2 Example 1 – Circle

We start with an example on the unit circle $\Omega=B_{1}(0)$ and choose a constant right-hand side $f=2^{2s}\Gamma(1+s)^{2}$ with exact solution given by (see, e.g., [AG17, BBN*+*18])

[TABLE]

Using polar coordinates, we can easily compute the exact energy as

[TABLE]

The exact solution is smooth inside the unit circle, but non-smooth towards the whole boundary, which is typical for solutions of our model problem (2.6). Therefore, we expect that the adaptive algorithm refines the mesh towards the whole boundary, which is indeed the case as shown in Figure 1.

Figure 2 shows the values of the error estimator (2.13) (brown triangles, blue squares) as well as the error (black diamonds, red stars) in the energy norm for uniform and adaptive mesh refinement by Algorithm 2.5 steered by the estimator (2.13) for $s=0.25$ (with $\theta=0.3$) and $s=0.75$ (with $\theta=0.5$).

As expected, uniform mesh refinement leads to a reduced rate of $N^{-1/4}$ due to the singularity of the solution at the boundary of the circle. However, adaptive refinement restores the optimal algebraic rate of $N^{-1/2}$ both for the error and the estimator as predicted by Theorem 2.6.

In Figure 3, we vary the adaptivity parameter $\theta\in(0,1]$ for fixed $s=0.25$. The dashed lines indicate the values of the error estimator and the solid ones the true errors in the energy norm. We observe optimal algebraic rates for all choices $\theta=0.3,0.5,0.7$ that are smaller than 1. $\theta=1$ leads to uniform refinement and therefore non-optimal convergence. As the error and estimator curves are not distinguishable for $\theta<1$, this experiment suggests, that the condition $\theta\ll 1$ in Theorem 2.6 is not necessary.

5.3 Example 2 – L-shaped domain

As a second example, we consider the L-shaped domain $\Omega=(-1,1)^{2}\backslash[0,1)^{2}$ depicted in Figure 1. We prescribe a constant right-hand side $f\equiv 1$. As the exact solution is unknown, we extrapolate the energy from the computed discrete energy.

Again, the adaptive algorithm refines the meshes towards the whole boundary due to the singularity of the solution at the whole boundary. This is in contrast to the known results for the integer Laplacian ($s=1$), where only refinement towards the reentrant corner is observed.

Figure 4 shows the convergence of the error estimator and the error for the cases $s=0.25$ ($\theta=0.3$) and $s=0.75$ ($\theta=0.5$). We observe optimal convergence rates for the adaptive method and reduced rates for uniform refinement just as in Example 1.

5.4 Example 3 - Circle, discontinuous right-hand side

As a final example, we again consider the unit circle $\Omega=B_{1}(0)$, but choose a discontinuous right-hand side

[TABLE]

For this problem, again an exact solution is known, see, e.g., [AG17]. The energy of this solution can be computed using the Meijer G-function as

[TABLE]

In Figure 5, an adaptively generated mesh and the discrete solution are plotted. In contrast to the previous examples, the adaptive algorithm does not only refine the mesh at the boundary, but also along the discontinuity of $f$ at the line $x=0$. However, the refinement towards the boundary tends to be stronger than towards the singularity of the right-hand side.

In Figure 5, convergence rates for the error estimator and the error are depicted. The empirical results are the same as for the previous two examples.