Linear and Nonlinear Fractional Elliptic Problems

TL;DR

This paper reviews recent advances in the analysis and numerical methods for linear fractional elliptic problems, focusing on regularity, boundary behavior, and finite element approximations with error estimates and computational challenges.

Contribution

It provides a comprehensive survey of recent analytical and numerical techniques for fractional elliptic problems, highlighting error analysis and computational strategies.

Findings

Error estimates for finite element methods on graded meshes

Analysis of boundary regularity and behavior

Discussion of computational challenges and numerical experiments

Abstract

This paper surveys recent analytical and numerical research on linear problems for the integral fractional Laplacian, fractional obstacle problems, and fractional minimal graphs. The emphasis is on the interplay between regularity, including boundary behavior, and approximability by piecewise linear finite element methods. We discuss several error estimates on graded meshes, and computational challenges associated to implementing and solving efficiently the ensuing integral equations, along with numerical experiments.

Tables

Table 1. Table 1. Computational rates of convergence (with respect to h ℎ h ) for the problem from Example   2.2 in d = 2 𝑑 2 d=2 dimensions. Rates using quasi-uniform meshes are listed in the second row, while rates using graded meshes, with μ = 2 𝜇 2 \mu=2 in ( 2.17 ), are reported in the third row.

Value of $s$	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
Uniform meshes	0.497	0.496	0.498	0.500	0.501	0.505	0.504	0.503	0.532
Graded meshes	1.066	1.040	1.019	1.002	1.066	1.051	0.990	0.985	0.977

Table 2. Table 2. Number of iterations for Gauss-Seidel (GS), conjugate gradient (CG) and preconditioned CG with BPX preconditioner (PCG). Stopping criteria is ‖ 𝐀𝐔 − 𝐅 ‖ 2 ‖ 𝐅 ‖ 2 < 10 − 6 subscript norm 𝐀𝐔 𝐅 2 subscript norm 𝐅 2 superscript 10 6 \frac{\|\mathbf{A}\mathbf{U}-\mathbf{F}\|_{2}}{\|\mathbf{F}\|_{2}}<10^{-6} .

$\overline{J}$	$N$	$s=0.9$			$s=0.5$			$s=0.1$
		GS	CG	PCG	GS	CG	PCG	GS	CG	PCG
8	209	101	18	17	22	14	16	9	21	18
9	413	172	24	20	30	18	18	9	25	21
10	821	314	32	22	42	20	19	9	29	23
11	1357	494	41	23	55	24	20	9	31	24
12	2753	792	55	25	72	30	21	9	36	26
13	4977	1391	73	26	94	34	22	10	37	27
14	9417	2357	95	27	133	39	23	10	40	28

Equations

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

$$\left\{\begin{array}[]{rl}(-\Delta)^{s}u=f&\mbox{ in }\Omega,\\ u=0&\mbox{ in }\Omega^{c}:={\mathbb{R}^{d}}\setminus\Omega.\end{array}\right.$$

$${\tilde{H}}^s(\Omega ):=\{v\in H^s({\mathbb{R}}^d):\mathrm{supp}v\subset \overline{\Omega }\}.$$

$$(v,w{)}_s:=\frac{C_{d,s}}{2}{\iint }_{Q_{\Omega }}\frac{(v(x)-v(y))(w(x)-w(y))}{|x-y{|}^{d+2s}}dxdy,$$

$$(v,w{)}_s=⟨(-\Delta {)}^{s}v,w⟩,\forall w\in {\tilde{H}}^s(\Omega ).$$

$$(u,v{)}_s=⟨f,v⟩\forall v\in {\tilde{H}}^s(\Omega ).$$

$$\Vert u{\Vert }_{{\tilde{H}}^{s+\gamma }(\Omega )}\le C(\Omega ,d,s,\gamma )\Vert f{\Vert }_{H^r(\Omega )}.$$

$$u(x)=\frac{\Gamma (\frac{d}{2})}{2^{2s}\Gamma (\frac{d+2s}{2})\Gamma (1+s)}(1-|x{|}^2{)}_{+}^s,$$

$$u(x)\approx d(x,\partial \Omega {)}^s\phi (x),$$

$$\delta(x,y):=\min\big{\{}\textrm{dist}(x,\partial\Omega),\textrm{dist}(y,\partial\Omega)\big{\}}$$

$$\Vert v{\Vert }_{H_{\kappa }^{\ell }(\Omega )}^2:=\Vert v{\Vert }_{H^k(\Omega )}^2+\sum _{|\beta |=k}{\iint }_{\Omega \times \Omega }\frac{|D^{\beta }v(x)-D^{\beta }v(y){|}^2}{|x-y{|}^{d+2s}}\delta (x,y{)}^{2\kappa }dydx$$

$$H_{\kappa }^{\ell }(\Omega ):=\left\{v\in H^{\ell }(\Omega ):\Vert v{\Vert }_{H_{\kappa }^{\ell }(\Omega )}<\infty \right\}.$$

$$\Vert u{\Vert }_{{\tilde{H}}_{1/2-\epsilon }^{s+1-2\epsilon }(\Omega )}\le \frac{C(\Omega ,d,s)}{\epsilon }\Vert f{\Vert }_{C^{1-s}(\overline{\Omega })}.$$

$$\Vert u{\Vert }_{B_{2,\infty }^{s+t}(\Omega )}\le C(\Omega ,d,s,t)\Vert f{\Vert }_{H^r(\Omega )}.$$

$$\Vert u{\Vert }_{H^{s+\gamma }(\Omega )}\le \frac{C(\Omega ,d,s,\gamma )}{\epsilon }\Vert f{\Vert }_{H^r(\Omega )},$$

$$[v{]}_{s+t,2,\Omega }:=\underset{h\in D\setminus \{0\}}{\sup }\frac{1}{|h{|}^t}|v-v(\cdot +h){|}_{H^s(\Omega )}$$

$$|T_{h}u-u{|}_s^2\le C|h{|}^s|u{|}_s^2,$$

$$|u{|}_{B_{2,\infty }^{3s/2}(\Omega )}\lesssim \Vert f{\Vert }_{B_{2,1}^{-s/2}(\Omega )},$$

$$\Vert u{\Vert }_{{\tilde{H}}^{r+2s}(\Omega )}\lesssim \Vert f{\Vert }_{H^r(\Omega )}.$$

$$\sigma :=\underset{h>0}{\sup }\underset{T\in {\mathcal{T}}_h}{\max }\frac{h_T}{{\rho }_T}<\infty ,$$

$${\mathbb{V}}_h:=\left\{v\in C_0(\Omega ):v=\sum _{i=1}^N{v}_i{\phi }_i\right\},$$

$$(u_h,v_h{)}_s=⟨f,v_h⟩\forall v_h\in {\mathbb{V}}_h.$$

$$\Vert u-u_h{\Vert }_{{\tilde{H}}^s(\Omega )}=\underset{v_h\in {\mathbb{V}}_h}{\min }\Vert u-v_h{\Vert }_{{\tilde{H}}^s(\Omega )}.$$

$$S_T^1:=\bigcup \left\{T^{\prime }\in {\mathcal{T}}_h:T\cap T^{\prime }\ne \varnothing \right\}.$$

$$S_T^2:=\bigcup \left\{T^{\prime }\in {\mathcal{T}}_h:S_T^1\cap T^{\prime }\ne \varnothing \right\},$$

$$|v{|}_{H^s(\Omega )}^2\le \sum _{T\in {\mathcal{T}}_h}\left[{\int }_T{\int }_{S_T^1}\frac{|v(x)-v(y){|}^2}{|x-y{|}^{d+2s}}dydx+\frac{C(d,\sigma )}{s{h}_T^{2s}}\Vert v{\Vert }_{L^2(T)}^2\right]\forall v\in H^s(\Omega ).$$

$${\int }_T{\int }_{S_T^1}\frac{|(v-{\Pi }_{h}v)(x)-(v-{\Pi }_{h}v)(y){|}^2}{|x-y{|}^{d+2s}}dydx\le C{h}_T^{2(\ell -s)}|v{|}_{H^{\ell }(S_T^2)}^2.$$

$$|v-{\Pi }_{h}v{|}_{H^s(\Omega )}\le C(\Omega ,d,s,\sigma ,\ell )h^{\ell -s}|v{|}_{H^{\ell }(\Omega )}.$$

$$h_T\le C(\sigma )\{\begin{array}{ll}{h}^{\mu },& T\cap \partial \Omega \ne \varnothing ,\\ h\text{dist}(T,\partial \Omega {)}^{(\mu -1)/\mu },& T\cap \partial \Omega =\varnothing .\end{array}$$

$$N=\dim {\mathbb{V}}_h\approx \{\begin{array}{ll}{h}^{-2},& \mu <2,\\ h^{-2}|\log h|,& \mu =2,\\ h^{-\mu },& \mu >2.\end{array}$$

Full text

JPB has been supported in part by NSF grant DMS-1411808

WL has been supported in part by NSF grant DMS-1411808 and the Patrick and Marguerite Sung Fellowship in Mathematics.

RHN has been supported in part by NSF grant DMS-1411808

1. Introduction

Diffusion, which is one of the most common physical processes, is the net movement of particles from a region of higher concentration to a region of lower concentration. The assumption that particles respond to Brownian motion leads to classical models of diffusion [36], that have been well studied for a long time. Fick’s first law states that the magnitude of the diffusive flux is proportional to the concentration gradient; by now, it is clear that such a constitutive relation is a questionable model for numerous phenomena [55]. When the associated underlying stochastic process is not given by Brownian motion, the diffusion is regarded as anomalous. In particular, anomalous superdiffusion refers to situations that can be modeled using fractional spatial derivatives or fractional spatial differential operators.

Integer-order differentiation operators are local because the derivative of a function at a given point depends only on the values of the function in an infinitesimal neighborhood of it. In contrast, fractional-order derivatives are nonlocal, integro-differential operators. A striking example of such an operator is the fractional Laplacian of order $s\in(0,1)$, which we will denote by $(-\Delta)^{s}$, and is given by

[TABLE]

We refer to [63] for an illustration of how the heat equation involving the fractional Laplacian arises from a simple random walk with jumps. The nonlocal structure of the operator (1.1) is apparent: to evaluate $(-\Delta)^{s}u$ at a spatial point, information involving all spatial points is needed.

This work deals with fractional diffusion. Our main goal is to review finite element methods (FEMs) to approximate solution of elliptic problems involving $(-\Delta)^{s}$ or related operators on bounded domains. We shall not discuss methods for the spectral fractional Laplacian; the surveys [10, 52] offer comparison between such an operator and the fractional Laplacian (1.1) and review other numerical methods. We point out that the fractional Laplacian (1.1) of order $s\in(0,1)$ is the infinitesimal generator of a $2s$-stable Lévy process. In this regard, problems on a bounded domain with homogeneous Dirichlet boundary conditions arise when the process is killed upon exiting the domain.

Throughout this work, we assume that $\Omega$ is a bounded Lipschitz domain. Whenever additional assumptions on $\partial\Omega$ are required, we shall state them explicitly. Even though there is a wide variety of numerical methods for fractional-order problems available in the literature [52], in this work we shall focus on piecewise linear finite element methods as in [10]. We emphasize the interplay between regularity, including boundary behavior, and approximability. In fact, the convergence rates achievable for the fractional elliptic PDEs discussed below, both linear and nonlinear, are limited by the presence of an algebraic boundary layer regardless of the regularity of $\partial\Omega$ and the polynomial degree for shape regular elements.

The paper is organized as follows. Section 2 deals with the homogeneous Dirichlet problem for the fractional Laplacian in $\Omega$; we discuss regularity of solutions and discuss both theoretical and computational aspects of conforming finite element discretizations. We also comment on some recent applications of this approach and on an alternative nonconforming FEM based on a Dunford-Taylor representation of the weak form of $(-\Delta)^{s}$. Afterwards, in Section 3 we address the obstacle problem for the fractional Laplacian. To derive optimal convergence estimates, we focus on weighted Sobolev regularity, where the weight is a power of the distance to the boundary of $\Omega$. These estimates follow from a precise quantification of boundary regularity of solutions and how solutions detach from the obstacle. Finally, Section 4 deals with fractional minimal graphs, which in fact are subgraphs that minimize a suitable nonlocal perimeter. We formulate a variational form for this problem, which is nonlinear and degenerate. We report on approximation properties of a conforming finite element scheme, and show convergence rates with respect to a novel geometric quantity. The paper concludes with a couple of computational explorations of the behavior of fractional minimal graphs for $d=2$.

2. Linear Problems

In this section we consider the homogeneous Dirichlet problem for the fractional Laplacian (1.1). Given $f:\Omega\to\mathbb{R}$, one seeks a function $u$ such that

[TABLE]

2.1. Variational Formulation

The natural variational framework for (2.1) is within the fractional Sobolev space ${\widetilde{H}}^{s}(\Omega)$, that is defined by

[TABLE]

We refer to [13] for definitions and elementary properties of fractional Sobolev spaces. Here we just state that on the space ${\widetilde{H}}^{s}(\Omega)$, because of the fractional Poincaré inequality, the natural inner product is equivalent to

[TABLE]

where $Q_{\Omega}=({\mathbb{R}^{d}}\times{\mathbb{R}^{d}})\setminus(\Omega^{c}\times\Omega^{c})$. The corresponding norm is $\|v\|_{{\widetilde{H}}^{s}(\Omega)}:=(v,v)_{s}^{1/2}$.

The duality pairing between ${\widetilde{H}}^{s}(\Omega)$ and its dual $H^{-s}(\Omega)$ shall be denoted by $\langle\cdot,\cdot\rangle$. In view of (2.2) we see that, whenever $v\in{\widetilde{H}}^{s}(\Omega)$ then $(-\Delta)^{s}v\in H^{-s}(\Omega)$ and that

[TABLE]

Therefore, given $f\in H^{-s}(\Omega)$, the weak formulation of (2.1) reads: find $u\in{\widetilde{H}}^{s}(\Omega)$ such that

[TABLE]

Existence and uniqueness of weak solutions, and stability of the solution map $f\mapsto u$, follow straightforwardly from the Lax-Milgram Theorem.

2.2. Regularity

A priori, it is not clear how smooth weak solutions are. If $f$ is more regular than $H^{-s}(\Omega)$, then $u$ could be expected to be more regular than ${\widetilde{H}}^{s}(\Omega)$. We now review some results regarding regularity of solutions to problem (2.4). Since our main interest is to derive convergence rates for finite element discretizations, we shall focus on Sobolev regularity estimates.

Recently, using Fourier analytical tools, Grubb [47] obtained estimates of solutions in terms of the so-called Hörmander $\mu$-spaces [48], but such estimates can be reinterpreted in terms of standard Sobolev spaces. A drawback of the following result from [47] is that it assumes the domain $\Omega$ to have smooth boundary, which is a too restrictive condition for finite element applications.

Theorem 2.1** (regularity on smooth domains).**

Let $\Omega$ be a domain with $\partial\Omega\in C^{\infty}$, $s\in(0,1)$, $f\in H^{r}(\Omega)$ for some $r\geq-s$, $u$ be the solution of (2.4) and $\gamma=\min\{s+r,1/2-\varepsilon\}$, with $\varepsilon>0$ arbitrarily small. Then, $u\in\widetilde{H}^{s+\gamma}(\Omega)$ and the following regularity estimate holds:

[TABLE]

A rather surprising feature of the previous result is that no matter how smooth the right hand side $f$ is, we cannot guarantee that solutions are any smoother than $\widetilde{H}^{s+1/2-\varepsilon}(\Omega)$. Indeed, because the fractional Laplacian is an operator of order $2s$, it could be expected to have a lift of order $2s$. As the following example [42] shows, such a reasoning is incorrect, and Theorem 2.1 is sharp.

Example 2.2** (limited regularity).**

Consider $\Omega=B(0,1)\subset{\mathbb{R}^{d}}$ and $f\equiv 1$. Then, the solution to (2.1) is given by

[TABLE]

where $t_{+}=\max\{t,0\}$.

As Example 2.2 illustrates, rough boundary behavior causes the reduced Sobolev regularity of solutions. Ros-Oton and Serra [59] studied problem (2.1) using potential theory tools, and were able to obtain a fine characterization of boundary behavior of solutions, that led them to deduce Hölder regularity estimates. In particular, it turns out that the asymptotic expansion

[TABLE]

where $\varphi$ is a smooth function, is generic.

The Hölder estimates from [59] give rise to Sobolev estimates for solutions in terms of Hölder norms of the data. To capture the boundary behavior, reference [2] introduced fractional weighted norms, where the weight is a power of the distance to the boundary, and developeds estimates in such norms. We denote

[TABLE]

and, for $\ell=k+s$, with $k\in\mathbb{N}$ and $s\in(0,1)$, and $\kappa\geq 0$, we define the norm

[TABLE]

and the associated space

[TABLE]

The regularity estimate in the weighted Sobolev scale (2.7) reads as follows.

Theorem 2.3** (weighted Sobolev estimate).**

Let $\Omega$ be a bounded, Lipschitz domain satisfying the exterior ball condition, $s\in(0,1)$, $f\in C^{1-s}(\overline{\Omega})$ and $u$ be the solution of (2.4). Then, for every $\varepsilon>0$ we have $u\in\widetilde{H}^{s+1-2\varepsilon}_{1/2-\varepsilon}(\Omega)$ and

[TABLE]

For simplicity, the theorem above was stated using the weight $\kappa=1/2-\varepsilon$; a more general form of the result can be found also in [10]. We point out that, for its application in finite element analysis, such a choice is optimal. In principle, increasing the exponent $\kappa$ of the weight allows for a higher differentiability in the solution, with no restriction on $\kappa$ above (as long as $f$ is sufficiently smooth). However, when exploiting this weighted regularity by introducing approximations on a family of shape-regular and graded meshes, the order of convergence (with respect to the number of degrees of freedom) is only incremented as long as $\kappa<1/2$.

It is worth pointing out that Theorem 2.3 is valid for Lipschitz domains satisfying an exterior ball condition. Although such a condition on the domain is much less restrictive than the $C^{\infty}$ requirement in Theorem 2.1, for polytopal domains it implies convexity. For that reason, we present here a result of an ongoing work [15], that characterizes regularity of solutions in terms of Besov norms.

Theorem 2.4** (regularity on Lipschitz domains).**

Let $\Omega$ be a bounded Lipschitz domain, $s\in(0,1)$ and $f\in H^{r}(\Omega)$ for some $r\in(-s,0]$. Then, the solution $u$ to (2.1) belongs to the Besov space $B^{s+t}_{2,\infty}(\Omega)$, where $t=\min\{s+r-\varepsilon,1/2\}$, with $\varepsilon>0$ arbitrarily small, with

[TABLE]

Consequently, by an elementary embedding, we deduce

[TABLE]

where $\gamma=\min\{s+r,1/2\}-\varepsilon$ is ‘almost’ as in Theorem 2.1.

We briefly outline the main ideas in the proof of Theorem 2.4, which follows a technique proposed by Savaré [61] for local problems. The point is to use the classical Nirenberg difference quotient method, and thus bound a certain Besov seminorm of the solution $u$.

Let $t\in(0,1)$ and $D$ be a set generating ${\mathbb{R}^{d}}$ and star-shaped with respect to the origin (for example, a cone). Then the functional

[TABLE]

induces the standard seminorm in the Besov space $B^{s+t}_{2,\infty}(\Omega)$. Because $\Omega$ is a Lipschitz domain, it satisfies a uniform cone property; upon a partition of unity argument, this gives (finitely many) suitable sets $D$ where translations can be taken. For a localized translation operator $T_{h}$, it is possible to prove a bound of the form

[TABLE]

which, in view of (2.10), yields $u\in B^{3s/2}_{2,\infty}(\Omega)$. Once this estimate has been obtained, a bootstrap argument leads to (2.8). Moreover, a refined estimate in $B^{3s/2}_{2,\infty}(\Omega)$ reads

[TABLE]

and interpolation with $|u|_{\widetilde{H}^{s}(\Omega)}\lesssim\|f\|_{H^{-s}(\Omega)}$ yields the following result [15].

Theorem 2.5** (lift theorem on Lipschitz domains).**

Let $\Omega$ be a bounded Lipschitz domain, $s\in(0,1)$ and $f\in H^{r}(\Omega)$ for some $r\in(-s,-s/2]$. Then, the solution $u$ to (2.1) belongs to the Sobolev space $\widetilde{H}^{r+2s}(\Omega)$, with

[TABLE]

2.3. Finite element discretization

In this section, we discuss a direct finite element method to approximate (2.4) using piecewise linear continuous functions. We consider a family $\{\mathcal{T}_{h}\}_{h>0}$ of conforming and simplicial meshes of $\Omega$, that we assume to be shape-regular, namely,

[TABLE]

where $h_{T}=\mbox{diam}(T)$ and $\rho_{T}$ is the diameter of the largest ball contained in $T$. As usual, the subindex $h$ denotes the mesh size, $h=\max_{T\in\mathcal{T}_{h}}h_{T}$ and we take elements to be closed sets. We denote by $\mathcal{N}_{h}$ the set of interior vertices of $\mathcal{T}_{h}$, by $N$ the cardinality of $\mathcal{N}_{h}$, and by $\{\varphi_{i}\}_{i=1}^{N}$ the standard piecewise linear Lagrangian basis, with $\varphi_{i}$ associated to the node $\texttt{x}_{i}\in\mathcal{N}_{h}$. Thus, the set of discrete functions is

[TABLE]

and is clearly conforming: $\mathbb{V}_{h}\subset{\widetilde{H}}^{s}(\Omega)$ for all $s\in(0,1)$.

With the notation described above, the discrete counterpart to (2.4) reads: find $u_{h}\in\mathbb{V}_{h}$ such that

[TABLE]

Because $u_{h}$ is the projection of $u$ onto $\mathbb{V}_{h}$ with respect to the ${\widetilde{H}}^{s}(\Omega)$-norm, we have the best approximation property

[TABLE]

Therefore, in order to obtain a priori rates of convergence in the energy norm, it suffices to bound the distance between the discrete spaces and the solution. Although the bilinear form $(\cdot,\cdot)_{s}$ involves integration on $\Omega\times{\mathbb{R}^{d}}$, one can apply an argument based on a fractional Hardy inequality, to prove that the energy norm may be bounded in terms of fractional–order norms on $\Omega$ (see [2]). It follows that bounding errors within $\Omega$ leads to error estimates in the energy norm.

A technical aspect of fractional-order seminorms is that they are not additive with respect to domain decompositions. With the goal of deriving interpolation estimates, we define the star or first ring of an element $T\in\mathcal{T}_{h}$ by

[TABLE]

We also introduce the star of $S^{1}_{T}$ (or second ring of $T$),

[TABLE]

and the star of the node $\texttt{x}_{i}\in\mathcal{N}_{h}$, $S_{i}:=\mbox{supp}(\varphi_{i})$.

Faermann [37] proved the localization estimate

[TABLE]

This inequality shows that to estimate fractional seminorms over $\Omega$, it suffices to compute integrals over the set of patches $\{T\times S^{1}_{T}\}_{T\in\mathcal{T}_{h}}$ plus local zero-order contributions. Bearing this in mind, one can prove the following type of estimates for suitable quasi-interpolation operators (see, for example, [16, 27]).

Proposition 2.6** (interpolation estimates on quasi-uniform meshes).**

Let $T\in{\mathcal{T}_{h}}$, $s\in(0,1)$, $\ell\in(s,2]$, and $\Pi_{h}$ be a suitable quasi-interpolation operator. If $v\in H^{\ell}(S^{2}_{T})$, then

[TABLE]

where $C=C(\Omega,d,s,\sigma,\ell)$. Therefore, for all $v\in H^{\ell}(\Omega)$, it holds

[TABLE]

The statement (2.16) in Proposition 2.6 could have also been obtained by interpolation of standard integer-order interpolation estimates. However, the technique of summing localized fractional-order estimates also works for graded meshes (cf. (2.18) and (2.19) below).

Combining estimate (2.16) with the best approximation property (2.14) and the regularity estimates described in Section 2.2, we can derive convergence rates. Concretely, the estimates from Theorem 2.1 and Theorem 2.4 translate into a priori rates for quasi-uniform meshes. However, optimal application of Theorem 2.3 requires a certain type of mesh grading. In two-dimensional problems ($d=2$), this can be attained by constructing graded meshes in the spirit of [46, Section 8.4]. In addition to shape regularity, we assume that the family $\{\mathcal{T}_{h}\}$ satisfies the following property: there is a number $\mu\geq 1$ such that given a parameter $h$ representing the meshsize at distance $1$ to the boundary $\partial\Omega$ and $T\in\mathcal{T}_{h}$, we have

[TABLE]

The number of degrees of freedom is related to $h$ by means of the parameter $\mu$ because (recall that $d=2$)

[TABLE]

Also, $\mu$ needs to be related to the exponent $\kappa$ used in the Sobolev regularity estimate (cf. Theorem 2.3). It can be shown that the choice $\mu=2$, that corresponds to $\kappa=1/2$, yields optimal convergence rates in terms of the dimension of $\mathbb{V}_{h}$. We also remark that, as discussed in [13], for three-dimensional problems ($d=3$), the grading strategy (2.17) becomes less flexible, and yields lower convergence rates. For optimal mesh grading beyond $\mu=2$ for both $d=2,3$ we need to break the shape regularity assumption and resort to anisotropic finite elements. They in turn are less flexible in dealing with the isotropic fractional norm of $H^{s}(\Omega)$ and its localization [37]. This important topic remains open.

Quasi-interpolation estimates in weighted Sobolev spaces (2.7) can be derived in the same way as in Proposition 2.6. More precisely, the weighted counterparts to (2.15) and (2.16) read

[TABLE]

for all $v\in H^{\ell}_{\kappa}(S^{2}_{T})$ and

[TABLE]

respectively. The constants above are $C=C(\Omega,d,s,\sigma,\ell,\kappa)$.

We collect all the convergence estimates in the energy norm –involving quasi-uniform and graded meshes– in a single statement [2].

Theorem 2.7** (energy error estimates for linear problem).**

Let $u$ denote the solution to (2.4) and denote by $u_{h}\in\mathbb{V}_{h}$ the solution of the discrete problem (2.13), computed over a mesh ${\mathcal{T}_{h}}$ consisting of elements with maximum diameter $h$. If $f\in L^{2}(\Omega)$, we have

[TABLE]

where $\alpha=\min\{s,1/2\}$. Additionally, if $\Omega\subset\mathbb{R}^{2}$, $f\in C^{1-s}(\overline{\Omega})$ and the family $\{\mathcal{T}_{h}\}$ satisfies (2.17) with $\mu=2$, we have

[TABLE]

To illustrate that Theorem 2.7 is sharp, we solve the problem from Example 2.2 on the discrete spaces (2.12) using a family of quasi-uniform meshes and a family of meshes graded according to (2.17). In Table 1, we report computational convergence rates in the energy norm for several values of $s$. We observe good agreement with the rates predicted by Theorem 2.7.

2.4. Computational challenges

Having at hand theoretical estimates for finite element discretizations of (2.4), we still need to address how to compute discrete solutions and, in particular, how to accelerate the assembly and solution of the discrete system that arises.

Matrix assembly. We first comment on key aspects of the finite element implementation for problems in dimension $d=2$. If $\mathbf{U}=(u_{i})_{i=1}^{N}$ and $u_{h}=\sum_{i=1}^{N}u_{i}\varphi_{i}$, it follows from (2.12) and (2.13) that the linear finite element system reads $\mathbf{A}\mathbf{U}=\mathbf{F}$ with stiffness matrix $\mathbf{A}$ and right-hand side vector $\mathbf{F}$ given by

[TABLE]

Computation of the stiffness matrix is not an easy task. There are two numerical difficulties in taking a direct approach. In first place, the bilinear form $(\cdot,\cdot)_{s}$ requires integration on unbounded domains; we point out that –at least for homogeneous problems as the ones considered here– integration over $\Omega\times\Omega^{c}$ can be reduced to a suitable integration over $\Omega\times\partial\Omega$ by using the Divergence Theorem [4]. Secondly, suitable quadrature rules to compute the stiffness matrix entries are required. To handle the singular (non-integrable) kernel $|x|^{-d-2s}$, one could use techniques from the boundary element method [26, 60]; we refer to [1] for details.

Compression. Note that, independently of $s$, finite element spaces (2.12) give rise to full stiffness matrices. Indeed, for any pair of nodes $\texttt{x}_{i},\texttt{x}_{j}$ such that $S_{i}\cap S_{j}=\emptyset$,

[TABLE]

Thus, computation of the stiffness matrix $\mathbf{A}$ involves a large number of far-field interactions, that is, elements $\mathbf{A}_{ij}$ for $\texttt{x}_{i}$ and $\texttt{x}_{j}$ sufficiently far. However, these elements should be significantly smaller than the ones that involve nodes close to one another. In [4, 64] the cluster method from the boundary element literature was applied and, instead of computing and storing all individual elements from $\mathbf{A}$, far field contributions are replaced by suitable low-rank blocks. The resulting data-sparse representation has $\mathcal{O}(N\log^{\alpha}N)$ complexity for some $\alpha\geq 0$. Reference [50] shows that the inverse of $\mathbf{A}$ can be represented using the same block structure as employed to compress the stiffness matrix.

Preconditioning. There are also issues to be addressed regarding the solution of the dense matrix equation $\mathbf{A}\mathbf{U}=\mathbf{F}$. The use of matrix factorization to solve such a system has complexity $\mathcal{O}(N^{3})$. As an alternative, one can use a conjugate gradient method and thereby the number of iterations needed for a fixed tolerance scales like $\sqrt{\kappa(\mathbf{A})}$, where $\kappa(\mathbf{A})$ is the condition number of $\mathbf{A}$ and satisfies [5]

[TABLE]

Therefore, for two-dimensional problems, we deduce $\kappa(\mathbf{A})=\mathcal{O}(h^{-2s})$ for quasi-uniform meshes, while $\kappa(\mathbf{A})=\mathcal{O}(h^{-2}|\log h|^{s})$ for meshes graded according to (2.17) with $\mu=2$. In the latter case, diagonal preconditioning allows us to recover the same condition number as for uniform meshes [5].

Recently, there have been some advances in the development of preconditioners for fractional diffusion. For instance, multigrid preconditioners were mentioned in [4], while operator preconditioners were studied in [45]. We now briefly comment on some features of an additive Schwarz preconditioner of BPX-type [17] (see also [38]). Assume we have a hierarchy of discrete spaces $\mathbb{V}_{0}\subset\ldots\mathbb{V}_{J}=\mathbb{V}$, with mesh size $h_{j}=\gamma^{2j}$, and let $\iota_{j}:\mathbb{V}_{j}\to\mathbb{V}$ be the inclusion operator. The basic ingredients needed to apply the general theory for additive Schwarz preconditioners are:

• •

Stable decomposition: for every $v\in\mathbb{V}$, there exists a decomposition $v=\sum_{j=0}^{J}v_{j}$ with $v_{j}\in\mathbb{V}_{j}$ such that

[TABLE]

A fundamental ingredient to prove this estimate for polyhedral domains is the optimal regularity pickup estimate for Lipschitz domains of Theorem 2.5, that allows us to perform an Aubin-Nitsche duality argument.

• •

Boundedness: for every $v=\sum_{j=0}^{J}v_{j}$ with $v_{j}\in\mathbb{V}_{j}$,

[TABLE]

As usual, boundedness of multilevel decompositions can be proved by estimating how much scales interact (i.e., using a strengthened Cauchy-Schwarz inequality). Nonlocality adds some difficulties to the derivation of such an estimate, because one cannot integrate by parts elementwise. The argument in [17] is based on the Fourier representation of the fractional Laplacian.

The conditions (2.20) and (2.21) imply that the preconditioner $\mathbf{B}:=\sum_{j=0}^{J}h_{j}^{2s-d}\iota_{j}\iota^{\prime}_{j}$ satisfies $\kappa(\mathbf{B}\mathbf{A})\leq\frac{c_{0}}{c_{1}}$ for graded bisection grids [17]. We illustrate this statement in Table 2 for $\Omega=(-1,1)^{2}$, $f=1$ and $s=0.9,0.5,0.1$. We observe a mild increase of iteration counts but rather robust performance with respect to $s$.

2.5. Applications and related problems

Numerical methods for fractional diffusion models have been extensively studied recently. Let us mention some applications on linear elliptic problems of the approach treated in this section:

• •

A posteriori error analyisis and adaptivity: The reduced regularity of solutions of (2.4) and the high computational cost of assembling the stiffness matrix $\mathbf{A}$ motivate the pursue of suitable adaptive finite element methods. A posteriori error estimates of residual type have been proposed and analyzed in [4, 39, 44, 58], and gradient-recovery based estimates in [64]. Adaptivity is however a topic of current research [38, 39, 44].

• •

Eigenvalue problems: The fractional eigenvalue problem arises, for example, in quantum mechanics problems in which the Brownian-like quantum paths are replaced by Lévy-like ones in the Feynman path integral [51]. Reference [18] studies conforming finite element approximations and applies the Babuška-Osborn theory [7], thereby obtaining convergence rates for eigenfunctions (in the energy and in the $L^{2}$-norms) and eigenvalues. Other methods, implemented on one-dimensional problems, include finite differences [35] and matrix methods [43, 65].

• •

Control problems: Finite element methods for linear-quadratic optimal control problems involving the fractional Laplacian (1.1) have been studied recently. In these problems, the control may be located inside [31] or outside the domain [6].

• •

Non-homogeneous Dirichlet conditions: A mixed method for the non-homogeneous Dirichlet problem for the integral fractional Laplacian was proposed in [3]. Such a method is based on weak enforcement of the Dirichlet condition and using a suitable non-local derivative [32] as a Lagrange multiplier. To circumvent approximating the nonlocal derivative, [6] proposed approximations of the non-homogeneous Dirichlet problem by a suitable Robin exterior value problem.

2.6. Nonconforming FEM: Dunford-Taylor approach

We finally report on a finite element approach for (2.1) proposed in [11] and based on the Fourier representation of the ${\widetilde{H}}^{s}(\Omega)$-inner product:

[TABLE]

This expression can be equivalently written as

[TABLE]

To see this, use Parseval’s formula to obtain

[TABLE]

followed by the change of variables $z=t|\xi|$, which converts the repeated integrals in the expression for $(v,w)_{s}$ into separate integrals, one of them being

[TABLE]

Although identity (2.22) is not an integral representation of the operator $(-\Delta)^{s}$, but rather of the bilinear form $(\cdot,\cdot)_{s}$, we regard it as a Dunford-Taylor representation.

To set up this formal calculation in the correct functional framework, given $u\in{\widetilde{H}}^{s}(\Omega)\subset L^{2}({\mathbb{R}}^{d})$ and $t>0$, let $v(u,t)\in H^{2+s}({\mathbb{R}}^{d})$ be the solution to $v-t^{2}\Delta v=-u$ in ${\mathbb{R}^{d}}$, or equivalently $v=-(I-t^{2}\Delta)^{-1}u$. Therefore, $\Delta v=t^{-2}(v+u)$ and

[TABLE]

This representation is the starting point of a three-step numerical method [11].

• •

Sinc quadrature: the change of variables $t=e^{-y/2}$ yields

[TABLE]

Thus, given an integer $N>0$ and a set of points $\{y_{j}\}_{j=-N}^{N}$ with uniform spacing $\approx N^{-1}$, the sinc quadrature $Q_{s}(u,w)$ approximation of $(\cdot,\cdot)_{s}$ is given by

[TABLE]

• •

Domain truncation: We stress that, in spite of $u$ being supported in $\Omega$, $v(u,t)$ is supported in all of ${\mathbb{R}^{d}}$ for all $t$; hence, some truncation is required. The method from [11] considers, for given $M>0$, a family of balls $B^{M}(t)$ that contain $\Omega$ and whose radius depends on $M$ and $t$ and can be computed a priori.

• •

Finite element approximation: Finally, a standard finite element discretization on $B^{M}(t)$ is performed. This requires meshes that fit $\Omega$ and $B^{M}(t)\setminus\Omega$ exactly – a non-trivial task; let us denote the discrete spaces on $\Omega$ and $B^{M}(t)$ by $\mathbb{V}_{h}$ and $\mathbb{V}_{h}^{M}$, respectively. Given $\psi\in L^{2}({\mathbb{R}^{d}})$, $t>0$ and $M>0$, we define $v_{h}^{M}=v_{h}^{M}(\psi,t)\in\mathbb{V}_{h}^{M}$ to be the unique solution of

[TABLE]

The fully discrete bilinear form reads:

[TABLE]

Using a Strang’s type argument to quantify the consistency errors generated by the three steps above, one obtains the a priori estimate [11, Theorem 7.7]

[TABLE]

where $\beta\in(s,3/2)$. Choosing $\beta=s+1/2-\varepsilon$, which is consistent with Theorem 2.1 and Theorem 2.4, $M=\mathcal{O}(|\log h|)$ and $N=\mathcal{O}(|\log h|^{2})$ gives the convergence rate

[TABLE]

This is similar to the rate obtained in Theorem 2.7 for quasi-uniform meshes. To the best of the authors’ knowledge, implementation of this approach over graded meshes, while feasible in theory, has not yet been pursued in practice.

3. Fractional Obstacle Problem

In this section we review finite element methods for the solution of the obstacle problem for the integral fractional Laplacian which, from now on, we shall simply refer to as the fractional obstacle problem.

The fractional obstacle problem appears, for example, in optimal stopping times for jump processes. In particular, it is used in the modeling of the rational price of perpetual American options [28]. More precisely, if $u$ represents the rational price of a perpetual American option where the assets prices are modeled by a Lévy process $X_{t}$, if $\chi$ denotes the payoff function, then $u$ solves a fractional obstacle problem with obstacle $\chi$.

An a posteriori error analysis of approximations of variational inequalities involving integral operators on arbitrary bounded domains was performed in [58]. We also comment on two recent works related to the approach we review here. Reference [20] deals with finite element discretizations to obstacle problems involving finite and infinite-horizon nonlocal operators. The experiments shown therein were performed on one-dimensional problems with uniform meshes, and indicate convergence with order $h^{1/2}$ in the energy norm. A theoretical proof of that convergence order was obtained in [12], where approximations using the approach discussed in Section 2.6 were considered. We also refer to [44] for computational comparisons between adaptive strategies and uniform and graded discretizations in two-dimensional problems.

3.1. Variational formulation

As before, we assume that $\Omega\subset{\mathbb{R}^{d}}$ is an open and bounded domain and, for the sake of applying weighted regularity estimates, we assume that $\Omega$ has a Lipschitz boundary and satisfies the exterior ball condition. Given $s\in(0,1)$ and functions $f:\Omega\to\mathbb{R}$ and $\chi:\overline{\Omega}\to\mathbb{R}$, with $\chi<0$ on $\partial\Omega$, the obstacle problem is a constrained minimization problem on ${\widetilde{H}}^{s}(\Omega)$ associated with a quadratic functional. Defining the admissible convex set

[TABLE]

the solution to the fractional obstacle problem is $u=\textrm{argmin}_{v\in{\mathcal{K}}}\mathcal{J}(v)$ where

[TABLE]

Existence and uniqueness of solutions is standard. Taking first variation of $\mathcal{J}$, we deduce that such a minimizer $u\in{\mathcal{K}}$ solves the variational inequality

[TABLE]

It can be shown [57] that, if $f\in L^{p}(\Omega)$ for $p>d/2s$, then the solution to the obstacle problem is indeed a continuous function, and that it satisfies the complementarity condition

[TABLE]

For our discussion, we assume that $f$ is such that the solution is defined pointwise, and consequently, we define the coincidence (or contact) and non-coincidence sets,

[TABLE]

The complementarity condition (3.2) can be succinctly expressed as $\lambda\geq 0$ in $\Lambda$ and $\lambda=0$ in $N$. The set $\partial\Lambda$, where the solution detaches from the obstacle, is the free boundary.

3.2. Regularity

The following regularity results for solutions to the fractional obstacle problem are instrumental for error analysis. We recall our assumption that the obstacle $\chi$ is a continuous function and strictly negative on $\partial\Omega$:

[TABLE]

Furthermore, we shall assume that $f\geq 0$. Heuristically, these assumptions should guarantee that the behavior of solutions near $\partial\Omega$ is dictated by a linear problem and that the nonlinearity is confined to the interior of the domain. Finally, to derive regularity estimates, we assume that the data satisfy

[TABLE]

where $\epsilon>0$ is sufficiently small, so that $1-2s+\epsilon$ is not an integer. Under these conditions, Caffarelli, Salsa and Silvestre [23] proved that the solution to the problem posed in the whole space (with suitable decay conditions at infinity) is of class $C^{1,s}({\mathbb{R}}^{d})$. It is worth examining the limiting cases $s=1$ and $s=0$. The former corresponds to the classical obstacle problem whose solutions are of class $C^{1,1}({\mathbb{R}}^{d})$. The latter reduces to $\min\{u-\chi,u-f\}=0$ whose solutions are just of class $C^{0,1}({\mathbb{R}}^{d})$. The regularity of [23] is thus a natural intermediate result.

We emphasize that deriving interior regularity estimates for (3.1) from this result, which is valid for problems posed in ${\mathbb{R}^{d}}$, is not as straightforward as for classical problems. Indeed, the nonlocal structure of $(-\Delta)^{s}$ implies that, if $0\leq\eta\leq 1$ is a smooth cut-off function such that $\eta=1$ in $\{\chi>0\}$, then

[TABLE]

To overcome this difficulty, reference [16] proceeds as follows. Given a set $D$ such that $\{\chi>0\}\subset D\subset\Omega$, one can define a cutoff $\eta$ such that $D\subset\{\eta=1\}$ and split the space roughly into a region where $\eta=1$, a region where $\eta=0$ and a transition region. In the first two regions, $(-\Delta)^{s}(\eta u)$ essentially coincides with a convolution operator with kernel $|z|^{-d-2s}$ but regularized at the origin, while the latter region is contained in the non-contact set $N$ and allows one to invoke interior regularity estimates for linear problems involving $(-\Delta)^{s}$. An important outcome is that solutions to fractional obstacle problems are more regular near the free boundary ($C^{1,s}$) than near the domain boundary ($C^{0,s}$). This is critical for approximation.

Alternatively, one may invoke the Caffarelli-Silvestre extension [24] to obtain local regularity estimates [23]. Since the extension problem involves a degenerate elliptic equation with a Muckenhoupt weight of class $A_{2}$ that depends only on the extended variable, one needs to combine fine estimates for degenerate equations with the translation invariance in the $x$-variable of the Caffarelli-Silvestre weight.

Once the interior regularity of solutions is established, one can invoke the Hölder boundary estimates for linear problems [59] and perform an argument similar to the one in [2] to deduce weighted Sobolev regularity estimates [16].

Theorem 3.1** (weighted Sobolev regularity for the obstacle problem).**

Let $\Omega$ be a bounded Lipschitz domain satisfying the exterior ball condition, $s\in(0,1)$, and $\chi\in C^{2,1}(\Omega)$ satisfying (3.3). Moreover, let $0\leq f\in\mathcal{F}_{s}(\overline{\Omega})$ and $u\in{\widetilde{H}}^{s}(\Omega)$ be the solution to (3.1). For every $\varepsilon>0$ we have that $u\in{\widetilde{H}}^{s+1-2\varepsilon}_{1/2-\varepsilon}(\Omega)$ with the estimate

[TABLE]

We have stated the estimate in Theorem 3.1 in weighted spaces because we are interested in the application of that result for finite element schemes over graded meshes. With the same arguments as in [16], it can be shown that the solution to the fractional obstacle problem (3.1) satisfies $u\in{\widetilde{H}}^{s+1/2-\varepsilon}(\Omega)$ and

[TABLE]

A similar result, for the obstacle problem for a class of integro-differential operators, was obtained in [12]. In the case of purely fractional diffusion (i.e., problems without a second-order differential operator), the estimate builds on [47] (cf. Theorem 2.1). Therefore, we point out that using Theorem 2.4, the requirement that $\Omega$ be a $C^{\infty}$ domain in [12, Cases A and B] can be relaxed to $\Omega$ being Lipschitz.

3.3. Finite element discretization

We consider the same finite element setting as in Section 2.3: let $\mathbb{V}_{h}$ be linear Lagrangian finite element spaces as in (2.12) over a family of conforming and simplicial meshes ${\mathcal{T}_{h}}$. An instrumental tool in the analysis we review here is the interpolation operator $\Pi_{h}:L^{1}(\Omega)\to\mathbb{V}_{h}$ introduced in [25] that, besides satisfying (2.15), is positivity preserving: it satisfies $\Pi_{h}v\geq 0$ for all $v\geq 0$. Such a property yields that, for every $v\in{\mathcal{K}}$,

[TABLE]

We therefore define the discrete admissible convex set

[TABLE]

and consider the discrete fractional obstacle problem: find $u_{h}\in{\mathcal{K}}_{h}$ such that

[TABLE]

We illustrate the delicate interplay between regularity and approximability next. We exploit that $u$ is both globally of class $C^{0,s}(\overline{\Omega})$, via graded meshes as in the linear problem, and locally of class $C^{1,s}(\Omega)$. First, we split the error as $\|u-u_{h}\|_{{\widetilde{H}}^{s}(\Omega)}^{2}=(u-u_{h},u-I_{h}u)_{s}+(u-u_{h},I_{h}u-u_{h})_{s}$, use Cauchy-Schwarz inequality and the interpolation estimate (2.19) to deduce

[TABLE]

This is a consequence of Theorem 2.3 and the use of graded meshes with parameter $\mu=2$ as in the linear theory of Section 2.3. For the remaining term we integrate by parts and utilize the discrete variational inequality (3.7) to arrive at

[TABLE]

Invoking the complementarity condition (3.2), we obtain

[TABLE]

We next observe that the integrand does not vanish only for elements $T$ in the vicinity of the free boundary, namely $T$’s for which $u\neq\chi$ and $(-\Delta)^{s}u\neq f$. Exploiting that $u\in C^{1,s}(\Omega)$, we infer that $(-\Delta)^{s}u-f\in C^{0,1-s}(\Omega)$, whence

[TABLE]

This yields the following optimal energy error estimate. We refer to [16] for details.

Theorem 3.2** (error estimate for obstacle problem).**

Let $u$ be the solution to (3.1) and $u_{h}$ be the solution to (3.7), respectively. Assume that $\chi\in C^{2,1}(\Omega)$ satisfies (3.3) and that $f\in\mathcal{F}_{s}(\overline{\Omega})$. If $d=2$, $\Omega$ is a convex polygon, and the meshes satisfy the grading hypothesis (2.17) with $\mu=2$, then we have that

[TABLE]

where $C>0$ depends on $\chi$, $s$, $d$, $\Omega$, $\varrho$ and $\|f\|_{\mathcal{F}_{s}(\overline{\Omega})}$.

We conclude this section with a computational example illustrating the qualitative behavior of solutions. Further experiments can be found in [16].

Example 3.3** (qualitative behavior).**

Consider problem (3.1), posed in the unit ball $B_{1}\subset\mathbb{R}^{2}$, with $f=0$ and the obstacle

[TABLE]

Figure 1 shows solutions for $s\in\{0.1,0.5,0.9\}$ on meshes graded according to (2.17) with $\mu=2$. The coincidence set $\Lambda$, which contains a neighborhood of the singular point $(1/4,0)$ is displayed in color in the bottom view. It can be observed that, while for $s=0.9$ the discrete solution resembles what is expected for the classical obstacle problem, the solution for $s=0.1$ is much flatter in the non-coincidence set $N$. Because $f\geq 0$, the solution $u$ satisfies $u\geq 0$. Therefore, the solution $u$ approaches $\chi_{+}$ in the limit $s\to 0$, while $u$ is expected to touch the obstacle only at the singular point and detach immediately in the diffusion limit $s=1$. The discrete nonlinear system has been solved using a semi-smooth Newton method.

4. Fractional minimal graphs

In this section we discuss the fractional minimal graph problem. The line of study of this nonlinear fractional problem, that can be regarded as a nonlocal version of the classical Plateau problem, and other related problems, began with the seminal works by Imbert [49] and Caffarelli, Roquejoffre and Savin [22].

As a motivation for the notion of fractional minimal sets, we show how the fractional perimeter arises in the study of a nonlocal version of the Ginzburg-Landau energy, extending a well-known result for classical minimal sets [56]. Let $\Omega\subset{\mathbb{R}^{d}}$ be a bounded set with Lipschitz boundary, $\varepsilon>0$ and define the energy

[TABLE]

where $W(t)=\frac{1}{4}(1-t^{2})^{2}$ is a double-well potential. Then, for every sequence $\{u_{\varepsilon}\}$ of minimizers of $\mathcal{J}_{\varepsilon}[u;\Omega]$ with uniformly bounded energies there exists a subsequence $\{u_{\varepsilon_{k}}\}$ such that

[TABLE]

where $E$ is a set with minimal perimeter in $\Omega$. We now consider a different regularization term: given $s\in(0,1/2)$, we set

[TABLE]

where ${Q_{\Omega}=\left({\mathbb{R}^{d}}\times{\mathbb{R}^{d}}\right)\setminus\left({\Omega}^{c}\times{\Omega}^{c}\right)}$ as in (2.2). The first term in the definition of $\mathcal{J}^{s}_{\varepsilon}$ involves the $H^{s}({\mathbb{R}^{d}})$-norm of $u$, except that the interactions over $\Omega^{c}\times\Omega^{c}$ are removed; for a minimization problem in $\Omega$, these are indeed fixed. As proved in [62], for every sequence $\{u_{\varepsilon}\}$ of minimizers of $\mathcal{J}^{s}_{\varepsilon}$ with uniformly bounded energies there exists a subsequence $\{u_{\varepsilon_{k}}\}$ such that

[TABLE]

However, instead of minimizing the perimeter in $\Omega$, here the set $E$ is a $s$-minimal set in $\Omega$, because it minimizes the so-called fractional perimeter $P_{s}(E,\Omega)$ among all measurable sets $F\subset{\mathbb{R}^{d}}$ such that $F\setminus\Omega=E\setminus\Omega$. In [22] this notion of fractional perimeter (also known as nonlocal perimeter) was proposed, and nonlocal minimal set problems were studied. We refer to [19, Chapter 6] and [29] for nice introductory expositions to the topic and applications.

4.1. Formulation of the problem and regularity

Our goal is to compute fractional minimal graphs, that is, to study the nonlocal minimal surface problem under the restriction of the domain being a cylinder. Concretely, from now on we consider $\Omega^{\prime}=\Omega\times\mathbb{R}$ with $\Omega\subset{\mathbb{R}^{d}}$ being a bounded Lipschitz domain. We assume that the exterior datum is the subgraph of some uniformly bounded function $g:{\mathbb{R}^{d}}\setminus\Omega\to\mathbb{R}$,

[TABLE]

The fractional minimal graph problem consists in finding a locally $s$-minimal set $E$ in $\Omega^{\prime}$ such that $E\setminus\Omega^{\prime}=E_{0}$. We refer to [53] for details on why the notion of locally $s$-minimality is the ‘correct’ one. Under the conditions described above, it can be shown that minimal sets need to be subgraphs, that is,

[TABLE]

for some function $u$ (cf. [54, Theorem 4.1.10]). We shall refer to such a set $E$ as a nonlocal minimal graph in $\Omega$.

In order to find nonlocal minimal graphs, we introduce the space

[TABLE]

(we write $\mathbb{V}^{0}$ whenever $g\equiv 0$) and, considering the weight function $F_{s}\colon\mathbb{R}\to\mathbb{R}$,

[TABLE]

we define the energy functional

[TABLE]

In [54, Chapter 4] it is shown that finding nonlocal minimal graphs is equivalent to minimizing the energy $I_{s}$ over $\mathbb{V}^{g}$. Existence of solution $u$ follows from the existence of locally $s$-minimal sets [53], while uniqueness is a consequence of $I_{s}$ being strictly convex. We also point out that for any function $v\colon{\mathbb{R}^{d}}\to\mathbb{R}$, its energy $I_{s}[v]$ is closely related to certain ${W^{2s}_{1}}$-seminorms [14, Lemma 2.5]:

[TABLE]

To give a clearer picture of the nonlocal minimal graph problem, we compare it to its classical counterpart. Given a bounded domain $\Omega\subset{\mathbb{R}^{d}}$ with sufficiently smooth boundary, and $g\colon\partial\Omega\to{\mathbb{R}^{d}}$, the classical Plateau problem consists in finding $u\colon\Omega\to{\mathbb{R}^{d}}$ that minimizes the graph surface area functional

[TABLE]

among those functions $u\in H^{1}(\Omega)$ satisfying $u=g$ on $\partial\Omega$. By taking first variation of $I$, it follows that the minimizer $u$ satisfies

[TABLE]

The left hand side in (4.7) consists of an $H^{1}$-inner product between $u$ and $v$, with a possibly degenerate weight that depends on $u$. For the nonlocal problem, after taking first variation of $I_{s}$ in (4.4), we obtain that $u$ is a minimizer if and only if

[TABLE]

where the bilinear form $a_{u}\colon\mathbb{V}^{g}\times\mathbb{V}^{0}\to\mathbb{R}$ is given by

[TABLE]

and $\widetilde{G}_{s}(\rho)=\int_{0}^{1}(1+\rho^{2}r^{2})^{-(d+1+2s)/2}dr$ and hence it satisfies $\rho\widetilde{G}_{s}(\rho)=G_{s}(\rho)=F^{\prime}_{s}(\rho)$. Similar to (4.7), the left hand side $a_{u}(u,v)$ in (4.8) is a weighted $H^{s+\frac{1}{2}}$-inner product with a possibly degenerate weight depending on $u$.

As for the regularity of nonlocal minimal graphs, the following result is stated in [21, Theorem 1.1] and builds on the arguments in [9, 41].

Theorem 4.1** (interior smoothness of nonlocal minimal graphs).**

Assume $E\subset\mathbb{R}^{d+1}$ is a locally $s$-minimal set in $\Omega^{\prime}=\Omega\times\mathbb{R}$, given by the subgraph of a measurable function $u$ that is bounded in an open set $\Lambda\supset\Omega$. Then, $u\in C^{\infty}(\Omega)$.

Remark 4.2* (stickiness).*

Theorem 4.1 does not address boundary regularity. By using (4.5), it can be easily proved that $u\in W^{2s}_{1}(\Omega)$ but, because $2s<1$, this does not even guarantee that $u$ has a trace on $\partial\Omega$. In fact, nonlocal minimal graphs can develop discontinuities across $\partial\Omega$. Furthermore, nonlocal minimal graphs generically exhibit this sticky behavior [33, 34].

4.2. Finite element discretization

In this section, we review the finite element discretization of the nonlocal minimal graph problem proposed in [14] and discuss its convergence and error estimates.

For simplicity, we assume that $\mbox{supp}(g)\subset\Lambda$ for some bounded set $\Lambda$. As before, we take a family $\{\mathcal{T}_{h}\}_{h>0}$ of conforming, simplicial and shape-regular meshes on $\Lambda$, which we impose to mesh $\Omega$ exactly. To account for non-zero boundary data, we make a slight modification on (2.12) to define the discrete spaces

[TABLE]

and we define

[TABLE]

where $\Pi_{h}^{c}$ denotes the Clément interpolation operator in $\Omega^{c}$.

With the notation introduced above, the discrete problem seeks $u_{h}\in\mathbb{V}^{g}_{h}$ such that

[TABLE]

Existence and uniqueness of solutions to this discrete problem follow directly from (4.5) and the strict convexity of $I_{s}$. To prove the convergence of the finite element scheme, the approach in [14] consists in proving that the discrete energy is consistent and afterwards using a compactness argument.

Theorem 4.3** (convergence for the nonlocal minimal graph problem).**

Let $s\in(0,1/2)$, $\Omega$ be a bounded Lipschitz domain and $g$ be uniformly bounded and satisfying $\mbox{supp}(g)\subset\Lambda$ for some bounded set $\Lambda$. Let $u$ and $u_{h}$ be, respectively, the solutions to (4.8) and (4.10). Then, it holds that

[TABLE]

The theorem above has the important feature of guaranteeing convergence without any regularity assumption on the solution. However, it does not offer convergence rates. We now show estimates for a geometric notion of error that mimics the one analyzed in [40] for the classical Plateau problem (see also [8, 30]). Such a notion of error is given by

[TABLE]

where $Q(\boldsymbol{a})=\sqrt{1+|\boldsymbol{a}|^{2}}$, $\widehat{\nu}(\boldsymbol{a})=\frac{(\boldsymbol{a},-1)}{Q(\boldsymbol{a})}$. Because $\widehat{\nu}(\nabla u)$ is the normal unit vector on the graph of $u$, the quantity $e(u,u_{h})$ is a weighted $L^{2}$-discrepancy between the normal vectors. For the nonlocal minimal graph problem, [14] introduced

[TABLE]

where $G_{s}(\rho)=F_{s}^{\prime}(\rho)$, the constant $\widetilde{C}_{d,s}=\frac{1-2s}{\alpha_{d}}$, $\alpha_{d}$ is the volume of the $d$-dimensional unit ball and $d_{v}$ is the difference quotient of the function $v$,

[TABLE]

In [14], this novel quantity $e_{s}(u,u_{h})$ is shown to be connected with a notion of nonlocal normal vector, and its asymptotic behavior as $s\to 1/2^{-}$ is established.

Theorem 4.4** (asymptotics of $e_{s}$).**

For all $u,v\in H^{1}_{0}(\Lambda)$, we have

[TABLE]

A simple ‘Galerkin orthogonality’ argument allows to derive an error estimate for $e_{s}(u,u_{h})$ (cf. [14, Theorem 5.1]).

Theorem 4.5** (geometric error).**

Under the same hypothesis as in Theorem 4.3, it holds that

[TABLE]

Therefore, to obtain convergence rates with respect to $e_{s}(u,u_{h})$, it suffices to prove interpolation estimates for the nonlocal minimizer. Although minimal graphs are expected to be discontinuous across the boundary, we still expect that $u\in BV(\Lambda)$ in general. Under this circumstance, the error estimate (4.14) leads to

[TABLE]

4.3. Numerical experiments

We conclude by presenting several numerical experiments and discussing the behavior of nonlocal minimal graphs. The first example we compute is on a one-dimensional domain, and is proposed and theoretically studied in [33, Theorem 1.2] as an illustration of stickiness phenomena.

Example 4.6** (stickiness in $1$D).**

Let $\Omega=(-1,1)$ and $g(x)=\textrm{sign}(x)$ for $x\in\Omega^{c}$. Discrete nonlocal minimal graphs for $s\in\{0.1,0.25,0.4\}$ are shown in Figure 2 (left).

Although our method requires discrete functions to be continuous across the boundary of $\Omega$, the presented $1$D picture clearly suggests a stickiness phenomenon. In addition, the plot also indicates that stickiness becomes more observable when $s$ gets closer to [math].

Example 4.7** (stickiness in an annulus).**

Let $\Omega=B_{1}\setminus\overline{B}_{1/2}\subset\mathbb{R}^{2}$, where $B_{r}$ denotes an open ball with radius $r$ centered at the origin, and let $g=\chi_{B_{1/2}}$. The discrete nonlocal minimal graph for $s=0.25$ is plotted in Figure 2 (right). In this example, stickiness is clearly observed on $\partial B_{1/2}$, while the stickiness on $\partial B_{1}$ is less noticeable.

Example 4.8** (effect of $s$).**

Let $\Omega=B_{1}\subset\mathbb{R}^{2}$, $g=\chi_{B_{3/2}\setminus\Omega}$. Figure 3 shows minimizers for several values of $s$. As $s\to 1/2$, the nonlocal minimal graphs get closer to the classical minimal graph, which is trivially constant in $\Omega$. On the other hand, as $s$ decreases we observe a stronger jump across $\partial\Omega$.

In the end, we make brief comments on the computational side for the examples presented above. In all of the experiments for fractional minimal graphs, we use Newton’s method to solve the nonlinear equation (4.10). Although the Jacobian matrix $\mathbf{A}_{u_{k}}$ in the iterative process corresponds to an $H^{s+1/2}-$inner product $a_{u_{k}}(w,v)$ with degenerate weights, our experiments indicate its condition number behaves like

[TABLE]

for $u_{k}$ whose gradient blows up near the boundary $\partial\Omega$, quasi-uniform meshes and dimensions $d=1,2$. This behavior is the same as the one in linear fractional diffusion with order $s+\frac{1}{2}$. However, the degenerate weight does bring more difficulties in preconditioning.

Another thing worth to be pointed out is the use of Dirichlet condition for our discrete space $\mathbb{V}_{h}^{g}$ requires the discrete function $u_{h}$ to be continuous across $\partial\Omega$. Due to the stickiness phenomenon in Remark 4.2, this may not be true for the solution $u$ of the minimal graph problem. Fortunately, this does not preclude the convergence in ‘trace blind’ fractional Sobolev spaces $W^{2s}_{1}(\Omega)$, and we are still able to capture discontinuities across the boundary in practice. While permitting discontinuities would be desirable, it has conflicts with using Newton’s method in solving (4.10) because the bilinear form $a_{u}(w,v)$ in (4.9) may not be well-defined. The question of how to solve the nonlinear equation (4.10) faster when allowing discontinuous across $\partial\Omega$ is still under investigation.