Finite element algorithms for nonlocal minimal graphs

1.   Introduction

This paper is the continuation of ^[7], where the authors proposed and analyzed a finite element scheme for the computation of fractional minimal graphs of order $s \in (0, 1/2)$ over bounded domains. That problem can be interpreted as a nonhomogeneous Dirichlet problem involving a nonlocal, nonlinear, degenerate operator of order $s + 1/2$. In this paper, we discuss computational aspects of such a formulation and perform several numerical experiments illustrating interesting phenomena arising in fractional Plateau problems and prescribed nonlocal mean curvature problems.

The notion of fractional perimeter was introduced in the seminal papers by Imbert ^[22] and by Caffarelli, Roquejoffre and Savin ^[12]. These works were motivated by the study of interphases that arise in classical phase field models when very long space correlations are present. On the one hand, ^[22] was motivated by stochastic Ising models with Kač potentials with slow decay at infinity, that give rise (after a suitable rescaling) to problems closely related to fractional reaction-diffusion equations such as

where $(-\Delta)^s$ denotes the fractional Laplacian of order $s\in (0, 1/2)$ and $f$ is a bistable nonlinearity. On the other hand, reference ^[12] showed that certain threshold dynamics-type algorithms, in the spirit of ^[26] but corresponding to the fractional Laplacian of order $s \in (0, 1/2)$ converge (again, after rescaling) to motion by fractional mean curvature. Fractional minimal sets also arise in the $\Gamma$-limit of nonlocal Ginzburg-Landau energies ^[28].

We now make the definition of fractional perimeter precise. Let $s \in (0, 1/2)$ and two sets $A, B \subset {\mathbb{R}^d}$, $d \ge 1$. Then, the fractional perimeter of order $s$ of $A$ in $B$ is

where $Q_B = ({\mathbb{R}^d} \times {\mathbb{R}^d}) \setminus (B^c \times B^c)$ and $B^c = {\mathbb{R}^d} \setminus B$. Given some set $A_0 \subset {\mathbb{R}^d} \setminus B$, a natural problem is how to extend $A_0$ into $B$ while minimizing the $s$-perimeter of the resulting set. This is the fractional Plateau problem, and it is known that, if $B$ is a bounded set, it admits a unique solution. Interestingly, in such a case it may happen that either the minimizing set $A$ is either empty in $B$ or that it completely fills $B$. This is known as a stickiness phenomenon ^[15].

In this work, we analyze finite element methods to compute fractional minimal graphs on bounded domains. Thus, we consider $s$-minimal sets on a cylinder $B = \Omega \times {\mathbb{R}}$, where $\Omega$ is a bounded and sufficiently smooth domain, with exterior data being a subgraph,

for some continuous function $g: {\mathbb{R}^d} \setminus \Omega \to {\mathbb{R}}$. We briefly remark some key features of this problem:

$\bullet$ A technical difficulty arises immediately: all sets $A$ that coincide with $A_0$ in ${\mathbb{R}^{d+1}} \setminus B$ have infinite $s$-perimeter in $B$. To remedy this issue, one needs to introduce the notion of locally minimal sets ^[24].

$\bullet$ There exists a unique locally $s$-minimal set, and it is given by the subgraph of a certain function $u$, cf. ^[14,25]. Thus, one can restrict the minimization problem to the class of subgraphs of functions that coincide with $g$ on $\Omega^c$.

$\bullet$ If the exterior datum $g$ is a bounded function, then one can replace the infinite cylinder $B = \Omega \times {\mathbb{R}}$ by a truncated cylinder $B_M = \Omega \times (-M, M)$ for some $M > 0$ sufficiently large [25,Proposition 2.5].

$\bullet$ Let $A$ be the subgraph of a certain function $v$ that coincides with $g$ on $\Omega^c$. One can rewrite $P_s(A, B_M)$ as

where $I_s$ is the nonlocal energy functional defined in (1.1) below [25,Proposition 4.2.8], [7,Proposition 2.3].

Therefore, an equivalent formulation to the Plateau problem for nonlocal minimal graphs consists in finding a function $u: {\mathbb{R}^d} \to {\mathbb{R}}$, with the constraint $u = g$ in ${\Omega}^c$, such that it minimizes the strictly convex energy

where $F_s$ is defined as

A remarkable difference between nonlocal minimal surface problems and their local counterparts is the emergence of stickiness phenomena ^[15]. In the setting of this paper, this means that the minimizer may be discontinuous across $\partial\Omega$. As shown by Dipierro, Savin and Valdinoci ^[17], stickiness is indeed the typical behavior of nonlocal minimal graphs in case $\Omega \subset {\mathbb{R}}$. When $\Omega \subset {\mathbb{R}}^2$, reference ^[16] proves that, at any boundary points at which stickiness does not happen, the tangent planes of the traces from the interior necessarily coincide with those of the exterior datum. Such a hard geometric constraint is in sharp contrast with the case of classical minimal graphs. In spite of their boundary behavior, fractional minimal graphs are smooth in the interior of the domain. Indeed, with the notation and assumptions from above it holds that $u \in C^\infty (\Omega)$; see [11,Theorem 1.1], and ^[5,19].

Our previous work ^[7] introduced and studied a finite element scheme for the computation of fractional minimal graphs. We proved convergence of the discrete minimizers as the mesh size tends to zero, both in suitable Sobolev norms and with respect to a novel geometric notion of error ^[7]. Stickiness phenomena was apparent in the experiments displayed in ^[7], even though the finite element spaces consisted of continuous, piecewise linear functions. We also refer the reader to ^[8] for further numerical examples and discussion on computational aspects of fractional minimal graph problems.

This article expands on our previous research ^[7,8] by discussing the convergence properties of a semi-implicit gradient flow and a damped Newton algorithm for computing discrete minimizers. Moreover, we prove convergence of our finite element scheme in $W^{2r}_1(\Omega)$, $r\in[0, s)$, for unboundedly supported data and graphs with prescribed non-zero nonlocal mean curvature. We perform various numerical experiments illustrating special –and sometimes counterintuitive– behaviors of solutions, especially near $\partial\Omega$. We hope these experiments will help understand the peculiarities of graphs with prescribed nonlocal mean curvature.

This paper is organized as follows. Section 2 gives the formulation of the minimization problem we aim to solve, and compares it with the classical minimal graph problem. Afterwards, in Section 3 we introduce our finite element method and review theoretical results from ^[7] regarding its convergence. Section 4 discusses computational aspects of the discrete problem, including the evaluation of the nonlocal form that it gives rise to, and the solution of the resulting discrete nonlinear equation via a semi-implicit gradient flow and a damped Newton method. Because the Dirichlet data may have unbounded support, we discuss the effect of data truncation and derive explicit bounds on the error decay with respect to the diameter of the computational domain in Section 5. Section 6 is concerned with the prescribed nonlocal mean curvature problem. Finally, Section 7 presents a number of computational experiments that explore qualitative and quantitative features of nonlocal minimal graphs and functions of prescribed fractional mean curvature, the conditioning of the discrete problems and the effect of exterior data truncation.

2.   Formulation of the problem

We now specify the problem we aim to solve in this paper and pose its variational formulation. Let $s \in (0, 1/2)$ and $g \in L^\infty(\Omega^c)$ be given. We consider the space

equipped with the norm

where

where $Q_{\Omega} = (\mathbb{R}^d\times\mathbb{R}^d) \setminus (\Omega\times\Omega)$. The space $\mathbb{V}^g$ can be understood as that of functions in $W^{2s}_1 (\Omega)$ with 'boundary value' $g$. The seminorm in $\mathbb{V}^g$ does not take into account interactions over $\Omega^c \times \Omega^c$, because these are fixed for the class of functions we consider; therefore, we do not need to assume $g$ to be a function in $W^{2s}_1(\Omega^c)$. In particular, $g$ may not decay at infinity. In case $g$ is the zero function, the space $\mathbb{V}^g$ coincides with the standard zero-extension Sobolev space $\widetilde{W}^{2s}_1(\Omega)$; for consistency of notation, we denote such a space by $\mathbb{V}^0$.

For convenience, we introduce the following notation: given a function $u \in \mathbb{V}^g$, the form $a_u \colon \mathbb{V}^g \times \mathbb{V}^0 \to \mathbb{R}$ is

where

It is worth noticing that ${\widetilde{G}_s}(\rho) \to 0$ as $|\rho| \to \infty$. Thus, the weight in (2.1) degenerates whenever the difference quotient $\frac{|u(x)-u(y)|}{|x-y|}$ blows up.

The weak formulation of the fractional minimal graph problem can be obtained by the taking first variation of $I_s[u]$ in (1.1) in the direction $v$. As described in ^[7], that problem reads: find $u \in \mathbb{V}^g$ such that

In light of the previous considerations, Eq (2.3) can be regarded as a fractional diffusion problem of order $s+1/2$ in $\mathbb{R}^d$ with weights depending on the solution $u$ and fixed nonhomogeneous boundary data $g$.

Remark 2.1 (comparison with local problems). Roughly, in the classical minimal graph problem, given some boundary data $g$, one seeks for a function $u \in g + H^1_0(\Omega)$ such that

The integral above can be interpreted as a weighted $H^1$-form, where the weight depends on $u$ and degenerates as $|\nabla u| \to \infty$.

In a similar way, problem $(2.3)$ involves a weighted $H^{s+1/2}$-form, in which the weight depends on $u$ and degenerates as $\frac{|u(x)-u(y)|}{|x-y|} \to \infty$. In this sense, it is not surprising that the fractional-order problems converge to the local ones as $s \to 1/2$. We refer to [7,Section 5] for a discussion on this matter.

3.   Finite element discretization

In this section we first introduce the finite element spaces and the discrete formulation of problem (2.3). Afterwards, we briefly outline the key ingredients in the convergence analysis for this scheme. For the moment, we shall assume that $g$ has bounded support:

The discussion of approximations in case of unboundedly supported data is postponed to Section 5.

3.1.   Discrete setting

We consider a family $\{{\mathcal{T}_h} \}_{h > 0}$ of conforming and simplicial triangulations of $\Lambda$, and we assume that all triangulations in $\{{\mathcal{T}_h} \}_{h > 0}$ mesh $\Omega$ exactly. Moreover, we assume $\{{\mathcal{T}_h} \}_{h > 0}$ to be shape-regular, namely:

where $h_T = \text{diam}(T)$ and $\rho_T$ is the diameter of the largest ball contained in the element $T\in {\mathcal{T}_h}$. The vertices of ${\mathcal{T}_h}$ will be denoted by $\{{\texttt{x}}_i\}$, and the star or patch of $\{ {\texttt{x}}_i \}$ is defined as

where $\varphi_i$ is the nodal basis function corresponding to the node ${\texttt{x}}_i$.

To impose the condition $u = g$ in $\Omega^c$ at the discrete level, we introduce the exterior interpolation operator

where $\Pi_h^{{\texttt{x}}_i} g$ is the $L^2$-projection of $g\big|_{S_i \cap \Omega^c}$ onto $\mathcal{P}_1(S_i \cap \Omega^c)$. Thus, $\Pi_h^c g ({\texttt{x}}_i)$ coincides with the standard Clément interpolation of $g$ on ${\texttt{x}}_i$ for all nodes ${\texttt{x}}_i$ such that $S_i \subset {\mathbb{R}^d} \setminus \overline{\Omega}$. On the other hand, for nodes ${\texttt{x}}_i \in\partial\Omega$, $\Pi_h^c$ only averages over the elements in $S_i$ that lie in $\Omega^c$.

We consider discrete spaces consisting of piecewise linear functions over ${\mathcal{T}_h}$,

To account for the exterior data, we define the discrete counterpart of $\mathbb{V}^g$,

With the same convention as before, we denote by $\mathbb{V}_h^0$ the corresponding space in case $g \equiv 0$. Therefore, the discrete weak formulation reads: find $u_h \in \mathbb{V}^g_h$ such that

Remark 3.1 (well-posedness of discrete problem). Existence and uniqueness of solutions to the discrete problem $(3.2)$ is an immediate corollary of our assumption $g \in L^{\infty}({\Omega}^c)$. Indeed, from this condition it follows that $u_h$ is a solution of $(3.2)$ if and only if $u_h$ minimizes the strictly convex energy $I_s[u_h]$ over the discrete space $\mathbb{V}_h^g$.

3.2.   Convergence

In ^[7], we have proved that solutions to (3.2) converge to the fractional minimal graph as the maximum element diameter tends to $0$. An important tool in that proof is a quasi-interpolation operator ${\mathcal I_h} \colon \mathbb{V}^g \to \mathbb{V}_h^g$ that combines the exterior Clément interpolation (3.1) with an interior interpolation operator. More precisely, we set

where $\Pi_h^\circ$ involves averaging over element stars contained in $\Omega$. Because the minimizer $u$ is smooth in the interior of $\Omega$, but we have no control on its boundary behavior other than the global bound $u \in W^{2s}_1(\Omega)$, we can only assert convergence of the interpolation operator in a $W^{2s}_1$-type seminorm without rates.

Proposition 3.2 (interpolation error). Let $s \in (0, 1/2)$, $\Omega$ be a bounded domain, $g \in C(\Omega^c)$, and $u$ be the solution to $(2.3)$. Then, the interpolation operator $(3.3)$ satisfies

Once we have proved the convergence of ${\mathcal I_h} u$ to $u$, energy consistency follows immediately. Since the energy dominates the $W^{2s}_1(\Omega)$-norm [7,Lemma 2.5], we can prove convergence in $W^{2r}_1(\Omega)$ for all $r \in [0, s)$ by arguing by compactness.

Theorem 3.3 (convergence). Assume $s \in (0, 1/2)$, $\Omega \subset {\mathbb{R}^d}$ is a bounded Lipschitz domain, and $g \in C_c({\mathbb{R}^d})$. Let $u$ be the minimizer of $I_s$ on $\mathbb{V}^g$ and $u_h$ be the minimizer of $I_s$ on $\mathbb{V}_h^g$. Then, it holds that

We finally point out that [7,Section 5] introduces a geometric notion of error that mimics a weighted $L^2$ discrepancy between the normal vectors to the graph of $u$ and $u_h$. We refer to that paper for further details.

3.3.   Graded meshes

As mentioned in the introduction, fractional minimal surfaces are smooth in the interior of $\Omega$. The main challenge in their approximation arises from their boundary behavior and concretely, from the genericity of stickiness phenomena, i.e. discontinuity of the solution $u$ across $\partial\Omega$. Thus, it is convenient to a priori adapt meshes to better capture the jump of $u$ near $\partial\Omega$.

In our discretizations, we use the following construction ^[21], that gives rise to shape-regular meshes. Let $h > 0$ be a mesh-size parameter and $\mu \ge 1$. Then, we consider meshes ${\mathcal{T}_h}$ such that every element $T \in {\mathcal{T}_h}$ satisfies

These meshes, typically with $\mu = 2$, give rise to optimal convergence rates for homogeneous problems involving the fractional Laplacian in 2d ^[2,6,9,10]. We point out that in our problem the computational domain strictly contains $\Omega$, because we need to impose the exterior condition $u = g$ on $\Omega^c$. As shown in ^[4,9], the construction (3.4) leads to

In our applications, because Theorem 3.3 gives no theoretical convergence rates, we are not restricted to the choice $\mu = 2$ in two dimensions: a higher $\mu$ allows a better resolution of stickiness. However, our numerical experiments indicate that the condition number of the resulting matrix at the last step of the Newton iteration deteriorates as $\mu$ increases, cf. Section 7.2.

4.   Numerical schemes

Having at hand a finite element formulation of the nonlocal minimal graph problem and proven its convergence as the mesh size tends to zero, we now address the issue of how to compute discrete minimizers in 1d and in 2d. In first place, we discuss the computation of matrices associated to either the bilinear form $a_{u_h} (\cdot, \cdot)$, or related computations. We propose two schemes for the solution of the nonlinear discrete problems (3.2): a semi-implicit gradient flow and a damped Newton method. In this section we also discuss the convergence of these two algorithms.

4.1.   Quadrature

We now consider the evaluation of the forms $a_{u_h}(\cdot, \cdot)$ appearing in (3.2). We point out that, following the implementation techniques from ^[1,2], if we are given $u_h \in \mathbb{V}_h^g$ and $v_h \in \mathbb{V}_h^0$, then we can compute $a_{u_h}(u_h, v_h)$. Indeed, since $a_{u_h}(u_h, v_h)$ is linear in $v_h$ and the latter function can be written in the form $v_h(x) = \sum_{{\texttt{x}}_i \in {\mathcal{N}_h^\circ}} v_i \varphi_i(x)$, we only need to evaluate

We split $Q_{\Omega} = (\Omega \times \Omega)\cup(\Omega \times \Omega^c)\cup(\Omega^c \times \Omega)$ and, because $\widetilde{G}_s$ is an even function (cf. (2.2)), we can take advantage that the integrand is symmetric with respect to $x$ and $y$ to obtain

We assume that the elements are sorted in such a way that the first ${\mathcal{N}_{\Omega}}$ elements mesh $\Omega$, while the remaining ${\mathcal{N}_{\Lambda}} - {\mathcal{N}_{\Omega}}$ mesh $\Lambda \setminus \Omega$, that is,

By doing a loop over the elements of the triangulation, the integrals $a_{i, 1}$ and $a_{i, 2}$ can be written as:

For the double integrals on $T_l \times T_m$ appearing in the definitions of $a_{i, 1}$ and $a_{i, 2}$, we apply the same type of transformations described in ^[1,13,27] to convert the integral into an integral over $[0, 1]^{2d}$, in which variables can be separated and the singular part can be computed analytically. The integrals over $T_l \times \Lambda^c$ are of the form

where the weight function $\omega$ is defined as

Since the only restriction on the set $\Lambda$ is that $\text{supp}(g) \subset \Lambda$, without loss of generality we assume that $\Lambda = B_{R}$ is a $d$-dimensional ball with radius $R$. In such a case, the integral over $\Lambda^c$ can be transformed using polar coordinates into:

where $\rho_0(e, x)$ is the distance from $x$ to $\partial B_{R}$ in the direction of $e$, which is given by the formula

The integral over $(\rho_0(e, x), \infty)$ can be transformed to an integral over $(0, 1)$ by means of the change of variable $\rho = \rho_0(e, x) \widetilde{\rho}^{-1/(2s)}$, and then approximated by Gaussian quadrature. Combining this approach with suitable quadrature over $\partial B_{1}$ and $T_l$, we numerically compute the integral over $T_l \times \Lambda^c$ for a given $u_h$.

4.2.   Gradient Flow

Although we can compute $a_{u_h}(u_h, v_h)$ for any given $u_h \in \mathbb{V}_h^g, \ v_h \in \mathbb{V}_h^0$, the nonlinearity of $a_{u_h}(u_h, v_h)$ with respect to $u_h$ still brings difficulties in finding the discrete solution to (3.2). Since $a_{u_h}(u_h, v_h) = \frac{\delta I_s[u_h]}{\delta u_h}(v_h)$ and $u_h$ minimizes the convex functional $I_s[u_h]$ in the space $\mathbb{V}_h^g$, a gradient flow is a feasible approach to solve for the unique minimizer $u_h$.

Given $\alpha \in [0, 1)$, and with the convention that $H^0 = L^2$, we first consider a time-continuous $H^\alpha$-gradient flow for $u_h(t)$, namely

where $u_h(0) = u_h^0 \in \mathbb{V}^g_h$ (and thus $I_s[u_h^0] < \infty$). Writing $u_h(t) = \sum_{x_j \in {\mathcal{N}_h^\circ}} u_j(t) \varphi_i$, local existence and uniqueness of solutions in time for (4.2) follow from the fact that $a_{u_h}(u_h, \varphi_i)$ is Lipschitz with respect to $u_j$ for any $\varphi_i$. Noticing that the gradient flow (4.2) satisfies the energy decay property

global existence and uniqueness of solutions in time can also be proved.

Similarly to the classical mean curvature flow of surfaces ^[18], there are three standard ways to discretize (4.2) in time: fully implicit, semi-implicit and fully explicit. Like in the classical case, the fully implicit scheme requires solving a nonlinear equation at every time step, which is not efficient in practice, while the fully explicit scheme is conditionally stable, and hence requires the choice of very small time steps. We thus focus on a semi-implicit scheme: given the step size $\tau > 0$ and iteration counter $k\ge0$, find $u_h^{k+1}\in\mathbb{V}_h^g$ that solves

The linearity of $a_{u^k_h}(u^{k+1}_h, v_h)$ with respect to $u^{k+1}_h$ makes (4.3) amenable for its computational solution. The following proposition proves the stability of the semi-implicit scheme. Its proof mimics the one of classical mean curvature flow ^[18].

Proposition 4.1 (stability of $H^\alpha$-gradient flow). Assume $u^{k+1}_h, u^k_h \in \mathbb{V}^g_h$ satisfy $(4.3)$. Then,

Proof. Choose $v_h = u^{k+1}_h - u^k_h \in \mathbb{V}^0_h$ in (4.3) to obtain

Next, we claim that for every pair of real numbers $r_0, r_1$, it holds that

We recall that $F_s$ is defined according to (1.2), that ${\widetilde{G}_s}$ satisfies ${\widetilde{G}_s}(r) = \frac{1}{r} G_s(r)$, and that $G_s = F_s'$. Since $F_s$ is a convex and even function, we deduce

We add and subtract $(|r_1| - |r_0|) \, |r_1| \, {\widetilde{G}_s}(|r_0|)$ above and use that ${\widetilde{G}_s}$ is even, decreasing on $[0, \infty)$ and non-negative, to obtain

This proves (4.5). Finally, define $d_k(x, y) : = \frac{u_h^{k}(x) - u_h^{k}(y)}{|x-y|}$ and set $r_0 = d_k$ and $r_1 = d_{k+1}$ in (4.5) to deduce that

Combining this with (4.4) finishes the proof.

Upon writing $w_h^k : = u_h^{k+1} - u_h^k$, the semi-implicit scheme (4.3) becomes (4.6), which is the crucial step of Algorithm 1 to solve (3.2).

Equation (4.6) boils down to solving the linear system $\left(M + \tau K^{k} \right) W^k = F^{k}$. In case $\alpha = 0$, the matrix $M = (M_{ij})$ is just a mass matrix, while if $\alpha > 0$, $M$ is the stiffness matrix for the linear fractional diffusion problem of order $\alpha$, given by

The matrix $K^{k} = \left(K^{k}_{ij}\right)$ is the stiffness matrix for a weighted linear fractional diffusion of order $s + \frac{1}{2}$, whose elements $K^{k}_{ij} : = a_{u_h^k} (\varphi_i, \varphi_j)$ are given by

and can be computed as described in Section 4.1. The right hand side vector is $F^{k} = - K^k U^k$, where $U^k = \left(U^{k}_i \right)$ is the vector $U^k_i = u_h^k(x_i)$, i.e., $f^{k}_i = -a_{u_h^k}(u_h^k, \varphi_i)$.

Because of Proposition 4.1 (stability of $H^\alpha$-gradient flow), the loop in Algorithm 1 terminates in finite steps. Moreover, {using the continuity of $a_{u_h^k}(\cdot, \cdot)$ in $[H^{\frac12 + s}(\Omega)]^2$, which is uniform in $u_h^k$, together with an inverse estimate and $0\le \alpha\le \frac12 + s$ gives

where the hidden constant depends on the mesh shape-regularity and $h_{\text{min}}$ is the minimum element size. Therefore, the last iterate $u_h^k$ of Algorithm 1 satisfies the residual estimate

4.3.   Damped Newton algorithm

Since the semi-implicit gradient flow is a first order method to find the minimizer of the discrete energy, it may converge slowly in practice. Therefore, it is worth having an alternative algorithm to solve (3.2) faster. With that goal in mind, we present in the following a damped Newton scheme, which is a second order method and thus improves the speed of computation.

To compute the first variation of $a_u(u, v)$ in (2.1) with respect to $u$, which is also the second variation of $I_s[u]$, we make use of $r{\widetilde{G}_s}(s) = G_s(r)$ and obtain

The identity $G'_s(a) = (1+a^2)^{-(d+1+2s)/2}$ can be easily determined from (4.1). Even though this first variation is not well-defined for an arbitrary $u \in \mathbb{V}^g$ and $v, w \in \mathbb{V}^0$, its discrete counterpart $\frac{\delta a_{u_h}(u_h, v_h)}{\delta u_h}(w_h)$ is well-defined for all $u_h \in \mathbb{V}_h^g$, $v_h, w_h \in \mathbb{V}_h^0$ because they are Lipschitz. Our damped Newton algorithm for (3.2) is presented in Algorithm 2.

Lemma 4.2 (convergence of Algorithm 2). The iterates $u_h^k$ of Algorithm 2 converge quadratically to the unique solution of $(3.2)$ from any initial condition.

Proof. Since $I_s[u_h]$ is strictly convex, the convergence of $u_h^k$ to the solution of discrete problem (3.2) is guaranteed by the theory of numerical optimization in finite dimensional spaces (see ^[23], for example).

The critical step in Algorithm 2 is to solve the equation (4.7). Due to the linearity of $\frac{\delta a_{u_h^k}(u_h^k, v_h)}{\delta u_h^k}(w_h^k)$ with respect to $v_h$ and $w_h^k$, we just need to solve a linear system $\widetilde{K}^k W^k = F^k$, where the right hand side $F^k = (f^k_i)$ is the same as the one in solving (4.6), namely, $f^k_i = a_{u_h^k}(u_h^k, \varphi_i)$. The matrix $\widetilde{K}^k = (\widetilde{K}^k_{ij})$, given by

is the stiffness matrix for a weighted linear fractional diffusion of order $s + \frac{1}{2}$. Since the only difference with the semi-implicit gradient flow algorithm is the weight, the elements in $\widetilde{K}^k$ can be computed by using the same techniques as for $K^k$.

5.   Unboundedly supported data

Thus far, we have taken for granted that $g$ has bounded support, and that the computational domain covers $\text{supp}(g)$. We point out that most of the theoretical estimates only require $g$ to be locally bounded. Naturally, in case $g$ does not have compact support, one could simply multiply $g$ by a cutoff function and consider discretizations using this truncated exterior condition. Here we quantify the consistency error arising in this approach. More precisely, given $H > 0$, we consider $\Omega_H$ to be a bounded open domain containing $\Omega$ and such that $d(x, \overline\Omega) \simeq H$ for all $x \in \partial \Omega_H$, and choose a cutoff function $\eta_H \in C^\infty(\Omega^c)$ satisfying

We replace $g$ by $g_H : = g \eta_H$, and consider problem (2.3) using $g_H$ as Dirichlet condition. Let $u^H \in \mathbb{V}^{g_H}$ be the solution of such a problem, and $u_h^H$ be the solution of its discrete counterpart over a certain mesh with element size $h$. Because of Theorem 3.3 we know that, for all $r \in[0, s)$,

Therefore we only need to show that, in turn, the minimizers of the truncated problems satisfy $u^H \to u$ as $H \to \infty$ in the same norm. As a first step, we compare the differences in the energy between truncated and extended functions. For that purpose, we define the following truncation and extension operators:

Proposition 5.1 (truncation and extension). The following estimates hold for every $v \in \mathbb{V}^g \cap L^\infty({\mathbb{R}^d})$, and $w \in \mathbb{V}^{g_H}\cap L^\infty({\mathbb{R}^d})$:

Proof. We prove only the first estimate, as the second one follows in the same fashion. Because $v = T_H v$ in $\Omega_{H}$, we have

From definition (1.2), it follows immediately that $F_s(0) = F'_s(0) = 0$, and thus $F_s (\rho) \le C \rho^2$ if $\rho \lesssim 1$. Combining this with the fact that $|v(x) - v(y)| \le 2 \| v \|_{L^\infty({\mathbb{R}^d})}$ and $|v(x) - T_H v(y)| \le 2 \| v \|_{L^\infty({\mathbb{R}^d})}$ for a.e. $x\in \Omega, y \in \Omega^c$, and integrating in polar coordinates, we conclude

This concludes the proof.

The previous result leads immediately to an energy consistency estimate for the truncated problem.

Corollary 5.2 (energy consistency). The minimizers of the original and truncated problem satisfy

Proof. Since $u^H$ is the minimizer over $\mathbb{V}^{g_H}$ and $T_H u\in\mathbb{V}^{g_H}$, we deduce

Conversely, using that $u$ is the minimzer over $\mathbb{V}^{g}$ and $Eu^H\in\mathbb{V}^{g}$, we obtain

and thus conclude the proof.

The energy $I_s$ is closely related to the $W^{2s}_1(\Omega)$-norm, in the sense that one is finite if and only if the other one is finite [7,Lemma 2.5]. Thus, in the same way as in Theorem 3.3 (convergence), energy consistency yields convergence in $W^{2r}_1(\Omega)$ for all $r \in [0, s)$.

Proposition 5.3 (convergence). Let $\Omega\subset{\mathbb{R}^d}$ be a bounded Lipschitz domain and $g\in L^\infty({\mathbb{R}^d})$. Let $u$ and $u_H$ be minimizers of $I_s$ over $\mathbb{V}^{g}$ and $\mathbb{V}^{g_H}$, respectively. Then for all $r \in [0, s)$, it holds that

Proof. The proof proceeds using the same arguments as in [7,Theorem 4.3]. In fact, from Corollary 5.2 we deduce that $\{ I_s[u^H] \}$ is uniformly bounded and therefore $\{u^H\}$ is bounded in $W^{2s}_1(\Omega)$. It follows that, up to a subsequence, $u^H$ converges in $L^1(\Omega)$ to a limit $\widetilde u$. Also, because $u^H = g$ in $\Omega_H$, we can extend $\widetilde u$ by $g$ on $\Omega^c$, and have $u^H \to u$ a.e in ${\mathbb{R}^d}$. We then can invoke Fatou's lemma and Corollary 5.2 to deduce that

Because $\widetilde u \in \mathbb{V}^g$, we deduce that $\widetilde u = u$ whence $u_H \to u$ in $L^1(\Omega)$ as $H\to0$. By interpolation, we conclude that convergence in $W^{2r}_1(\Omega)$ holds for all $r \in [0, s)$.

6.   Prescribed nonlocal mean curvature

In this section, we briefly introduce the problem of computing graphs with prescribed nonlocal mean curvature. More specifically, we address the computation of a function $u$ such that for a.e. $x \in \Omega$, a certain nonlocal mean curvature at $\big(x, u(x)\big)$ is equal to a given function $f(x)$. For a set $E \subset {\mathbb{R}^{d+1}}$ and $\widetilde{x} \in \partial E$, such nonlocal mean curvature operator is defined as ^[12]

In turn, for $\widetilde{x} = (x, u(x))$ on the graph of $u$, this can be written as [25,Chapter 4]

To recover the classical mean curvature in the limit $s \to \frac12^-$, it is necessary to normalize the operator $H_s$ accordingly. Let $\alpha_{d}$ denote the volume of the $d$-dimensional unit ball, and consider the prescribed nonlocal mean curvature problem

The scaling factor $\frac{1-2s}{d \alpha_d}$ yields [7,Lemma 5.8]

where $H[E]$ denotes the classical mean curvature operator. Therefore, in the limit $s \to \frac12^-$, formula (6.1) formally becomes the following Dirichlet problem for graphs of prescribed classical mean curvature:

An alternative formulation of the prescribed nonlocal mean curvature problem for graphs is to find $u \in \mathbb{V}^g$ minimizing the functional

Because $I_s[u]$ is convex and the second term in the right hand side above is linear, it follows that this functional is also convex. Then, by taking the first variation of (6.4), we see that $u \in \mathbb{V}^g$ is the minimizer of $\mathcal{K}_{s}[\cdot; f]$ if and only if it satisfies

for every $v \in \mathbb{V}^0$. Formally, (6.1) coincides with (6.5) because one can multiply (6.1) by a test function $v$, integrate by parts and take advantage of the fact that $G_s$ is an odd function to arrive at (6.5) up to a constant factor.

One intriguing question regarding the energy $\mathcal{K}_{s}[u; f]$ in (6.4) is what conditions on $f$ are needed to guarantee that it is bounded below. In fact, for the variational formulation of the classical mean curvature problem (6.3), Giaquinta ^[20] proves the following necessary and sufficient condition for well posedness: there exists some ${\varepsilon}_0 > 0$ such that for every measurable set $A \subset \Omega$,

where $\mathcal{H}^{d-1}$ denotes the $(d-1)-$dimensional Hausdorff measure. In some sense, this condition ensures that the function $f$ be suitably small.

Although we are not aware of such a characterization for prescribed nonlocal mean curvature problems, a related sufficient condition for $\mathcal{K}_{s}[u; f]$ to have a lower bound can be easily derived. We discuss this next.

Proposition 6.1 (convergence). Let $s \in (0, 1/2)$, $\Omega \subset {\mathbb{R}^d}$ be a bounded Lipschitz domain, $\Vert f \Vert_{L^{d/(2s)}(\Omega)}$ be sufficiently small, and $g \in C_c({\mathbb{R}^d})$. Then, for all $h > 0$ there exist unique minimizers $u\in\mathbb{V}^g$ and $u_h\in\mathbb{V}_h^g$ of $\mathcal{K}_{s}[\cdot; f]$, and they satisfy

Proof. Exploiting [7,Lemma 2.5 and Proposition 2.7], we deduce that

for suitable constants $C_0 = C_0(d, \Omega, s, \|g\|_{L^\infty(\Omega^c)})$ and $C_1$. On the other hand, Hölder's inequality and the Sobolev embedding $W^{2s}_1(\Omega) \subset L^{d/(d-2s)}(\Omega)$ give

whence $\mathcal{K}_{s}[u; f]$ is bounded from below provided $\Vert f \Vert_{L^{d/(2s)}(\Omega)}$ is suitably small,

Thus, the existence of minimizers of $\mathcal{K}_{s}[\cdot; f]$ on either $\mathbb{V}^g$ or $\mathbb{V}_h^g$ follows by standard arguments. To prove the convergence of discrete minimizers $u_h$, it suffices to extend the arguments of Section 3.2 to $f\ne0$ following [7,Theorem 4.2].

Remark 6.2 (consistency with local problems). The requirement that $\Vert f \Vert_{L^{d/(2s)}(\Omega)}$ be sufficiently small and estimate $(6.7)$ are to some extent consistent with $(6.6)$, because it holds that

due to Hölder's inequality and the isoperimetric inequality, and formally the case $2s = 1$ corresponds to the classical prescribed mean curvature problem (cf. $(6.2)$).

7.   Numerical experiments

This section presents a variety of numerical experiments that illustrate some of the main features of fractional minimal graphs discussed in this paper. From a quantitative perspective, we explore stickiness and the effect of truncating the computational domain. Moreover, we report on the conditioning of the matrices arising in the iterative resolution of the nonlinear discrete equations. Our experiments also illustrate that nonlocal minimal graphs may change their concavity inside the domain $\Omega$, and we show that graphs with prescribed fractional mean curvature may be discontinuous in $\Omega$.

In all the experiments displayed in this section we use the damped Newton algorithm from §4.3. We refer to ^[7] for experiments involving the semi-implicit gradient flow algorithm and illustrating its energy-decrease property.

7.1.   Quantitative boundary behavior

We first consider the example studied in [15,Theorem 1.2]. We solve (3.2) for $\Omega = (-1, 1) \subset {\mathbb{R}}$ and $g(x) = M \text{sign}(x)$, where $M > 0$. Reference ^[15] proves that, for every $s \in (0, 1/2)$, stickiness (i.e., the solution being discontinuous at $\partial \Omega$) occurs if $M$ is big enough and, denoting the corresponding solution by $u^M$, that there exists an optimal constant $c_0$ such that

In our experiments, we consider $s = 0.1, 0.25, 0.4$ and use graded meshes (cf. Section 3.3) with parameter $\mu = 2, h = 10^{-3}$ to better resolve the boundary discontinuity. The mesh size $h$ here is taken in such a way that the resulting mesh partitions $\Omega = (-1, 1)$ into $\lfloor \frac{|\Omega|^{1/\mu}}{h}\rfloor$ subintervals and the smallest ones have size $h^\mu$. Moreover, since this is an example in one dimension and the unboundedly supported data $g$ is piecewise constant, we can use quadrature to approximate the integrals over $\Omega^c$ rather than directly truncating $g$. The left panel in Figure 1 shows the computed solutions with $M = 16$.

In all cases we observe that the discrete solutions $u_h$ are monotonically increasing in $\Omega$, so we let $x_1$ be the free node closest to $1$ and use $u^M_h(x_1)$ as an approximation of $\sup_{x \in \Omega} u^M(x)$. The right panel in Figure 1 shows how $u^M_h(x_1)$ varies with respect to $M$ for different values of $s$.

For $s = 0.1$ and $s = 0.25$ the slopes of the curves are slightly larger than the theoretical rate $M^{\frac{1+2s}{2+2s}}$ whenever $M$ is small. However, as $M$ increases, we see a good agreement with theory. Comparing results for $M = 2^{7.5}$ and $M = 2^8$, we observe approximate rates $0.553$ for $s = 0.1$ and $0.602$ for $s = 0.25$, where the expected rates are $6/11 \approx 0.545$ and $3/5 = 0.600$, respectively. However, the situation is different for $s = 0.4$: the plotted curve does not correspond to a flat line, and the last two nodes plotted, with $M = 2^{7.5}$ and $M = 2^8$, show a relative slope of about $0.57$, which is off the expected $9/14 \approx 0.643$.

We believe this issue is due to the mesh size $h$ not being small enough to resolve the boundary behavior. We run the same experiment on a finer mesh, namely with $h = 10^{-4}, \mu = 2$, and report our findings for $s = 0.4$ and compare them with the ones for the coarser mesh on Table 1. The results are closer to the predicted rate.

7.2.   Conditioning

For the solutions of the linear systems arising in our discrete formulations, we use a conjugate gradient method. Therefore, the number of iterations needed for a fixed tolerance scales like $\sqrt{\kappa(K)}$, where $\kappa(K)$ is the condition number of the stiffness matrix $K$. For linear problems of order $s$ involving the fractional Laplacian $(-\Delta)^s$, the condition number of $K$ satisfies ^[3]

Reference ^[3] also shows that diagonal preconditioning yields $\kappa(K) = \mathcal{O}\left(N^{2s/d} \right)$, where $N$ is the dimension of the finite element space.

Using the Matlab function $\texttt{condest}$, we estimate the condition number of the Jacobian matrix in the last Newton iteration in the example from Section 7.1 with $M = 1$, with and without diagonal preconditioning. Figure 2 summarizes our findings.

Let $N = \text{dim} \mathbb{V}_h^g$ be the number of degrees of freedom. For a fixed $s$ and using uniform meshes, we observe that the condition number behaves like $N^{2s+1} \simeq h^{-2s-1}$: this is consistent with the $s$-fractional mean curvature operator being an operator of order $s+1/2$. For graded meshes (with $\mu = 2, \mu = 3$), the behavior is less clear. When using diagonal preconditioning for $\mu = 2$, we observe that the condition number also behaves like $N^{2s+1}$.

7.3.   Truncation of unboundedly supported data

In Section 5, we studied the effect of truncating unboundedly supported data and proved the convergence of the discrete solutions of the truncated problems $u^H_h$ towards $u$ as $h \to 0$, $H\to \infty$.

Here, we study numerically the effect of data truncation by running experiments on a simple two-dimensional problem. Consider $\Omega = B_1 \subset {\mathbb{R}}^2$ and $g \equiv 1$; then, the nonlocal minimal graph $u$ is a constant function. For $H > 0$, we set $\Omega_H = B_{H+1}$. and compute nonlocal minimal graphs on $\Omega$ with Dirichlet data $g^H = \chi_{\Omega_H}$, which is a truncation of $g \equiv 1$. Clearly, if there was no truncation, then $u_h$ should be constantly $1$; the effect of the truncation of $g$ is that the minimum value of $u^H_h$ inside $\Omega$ is strictly less than $1$. For $s = 0.25$, we plot the $L^1(\Omega)$ and $L^{\infty}(\Omega)$ norms of $u_h - u^H_h$ as a function of $H$ in Figure 3. The slope of the curve is close to $-1.5$ for large $H$, which is in agreement with the $\mathcal{O}(H^{-1-2s})$ consistency error for the energy $I_s$ we proved in Corollary 5.2.

7.4.   Change of convexity

This is a peculiar behavior of fractional minimal graphs. We consider $\Omega = (-1, 1)$, $s = 0.02$, $g(x) = 1$ for $a \le |x| \le 2$ and $g(x) = 0$ otherwise, and denote by $u_a$ the solution of (3.2). For $a = 1$, it is apparent from Figure 4 (left panel) that the solution $u_1$ is convex in $\Omega$ and has stickiness on the boundary. In addition, the figure confirms that $\lim\limits_{x \to 1^{-}} u_a'(x) = \infty$, which is asserted in [17,Corollary 1.3]. On the contrary, for $1 < a < 2$, as can be seen from Figure 4 (right panel), [17,Corollary 1.3] implies that $\lim\limits_{x \to 1^{-}} u_a'(x) = -\infty$ since $g(x) = 0$ near the boundary of $\Omega$. This fact implies that $u(x)$ cannot be convex near $x = 1$ for $1 < a < 2$. Furthermore, as $a \to 1^+$ one expects that $u_a(x) \to u_1(x)$ and thus that $u_a$ be convex in the interior of $\Omega = (-1, 1)$ for $a$ close to $1$. Therefore it is natural that for some values of $a > 1$ sufficiently close to $1$, the solution $u_a$ changes the sign of its second derivative inside $\Omega$. In fact, we see from the right panel in Figure 4 that the nonlocal minimal graph $u$ in $\Omega$ continuously changes from a convex curve into a concave one as $a$ varies from $1$ to $1.5$.

This change of convexity is not restricted to one-dimensional problems. Let $\Omega \subset {\mathbb{R}}^2$ be the unit ball, $s = 0.25$, and $g(x) = 1$ for $\frac{129}{128} \le |x| \le 1.5$ and $g(x) = 0$ otherwise. Figure 5 (right panel) shows a radial slice of the discrete minimal graph, which is a convex function near the origin but concave near $\partial\Omega$. An argument analogous to the one we discussed in the previous paragraph also explains this behavior in a two-dimensional experiment.

7.5.   Geometric rigidity

Stickiness is one of the intrinsic and distintive features of nonlocal minimal graphs. It can be delicate especially in dimension more than one. We now analyze a problem studied in ^[16] that illustrates the fact that for $\Omega \subset {\mathbb{R}}^2$, if nonlocal minimal graphs are continuous at some point $x \in \partial\Omega$ then they must also have continuous tangential derivatives at such a point. This geometric rigidity stands in sharp contrast with the case of either fractional-order linear problems and classical minimal graphs.

Specifically, we consider $\Omega = (0, 1) \times (-1, 1)$ and the Dirichlet data

where $a \in [0, 1]$ and $\gamma > 0$ are parameters to be chosen. We construct graded meshes with $\mu = 2$ and smallest mesh size $h^{\mu} = 2^{-7}$; see Section 3.3. Figure 6 (left panel) displays the numerical solution $u_h$ associated with $s = 0.25, \; \gamma = 2$ and $a = 1/8$.

If one defines the function $u_0(y) = \lim_{x \to 0^+} u(x, y)$, then according to [16,Theorem 1.4], one has $u'_0(0) = 0$ for $a > 0$. We run a sequence of experiments to computationally verify this theoretical result. For meshes with $\mu = 2$ and $h^{\mu} = 2^{-7}, 2^{-8}, 2^{-9}$, the slopes of $u_h$ in the $y$-direction at $(x, 0)$ for $x = 2^{-6}, 2^{-7}, 2^{-8}, 2^{-9}$, are recorded in Table 2 below for $s = 0.1, 0.25, 0.4$. Because computing the slope of $u_h$ at $(x, 0)$ would be meaningless when $x$ is smaller than $h^{\mu}$, we write a $\texttt{N/A}$ symbol in those cases. Our experiments show that the slopes decrease as $x$ approaches $0$.

To further illustrate this behavior, in Figure 6 (right panel) we display the computed solutions $u_h(x, y)$ at $x = 2^{-3}, 2^{-6}, 2^{-9}$, for $s = 0.25$ over a mesh with $h^{\mu} = 2^{-9}$. The flattening of the curves as $x \to 0^+$ is apparent.

7.6.   Prescribed nonlocal mean curvature

This section presents experiments involving graphs with nonzero prescribed mean curvature. We run experiments that indicate the need of a compatibility condition such as (6.6), the fact that solutions may develop discontinuities in the interior of the domain, and point to the relation between stickiness and the nonlocal mean curvature of the domain.

7.6.1.   Compatibility

As discussed in Section 6, the prescribed nonlocal mean curvature problem (6.5) may not have solutions for some functions $f$. To verify this, in Figure 7 we consider $\Omega = B(0, 1) \subset {\mathbb{R}}^2$, $s = 0.25$, $g = 0$ and two choices of $f$. For the picture on the right ($f = -10$), the residue does not converge to $0$, and the energy $\mathcal{K}_{s}[u; f]$ goes from $0$ initially down to $-6.6 \times 10^6$ after $16$ Newton iterations.

7.6.2.   Discontinuities

Another interesting phenomenon we observe is that, for a discontinuous $f$, the solution $u$ may also develop discontinuities inside $\Omega$. We present the following two examples for $d = 1$ and $d = 2$.

In first place, let $\Omega = (-1, 1) \subset {\mathbb{R}}$, $s = 0.01$, $g = 0$ and consider $f(x) = 1.5 \, \text{sign}(x)$. We use a mesh graded toward $x = 0, \, \pm 1$ with $N = 2000$ degrees of freedom and plot the numerical solution $u_h$ in Figure 8. The behavior of $u_h$ indicates that the solution $u$ has discontinuities both at $x = \pm 1$ and $x = 0$.

As a second illustration of interior discontinuities, let $\Omega = (-1, 1)^2 \subset {\mathbb{R}}^2$, $s = 0.01$, $g = 0$ and consider $f(x, y) = 4 \, \text{sign}(xy)$. We use a mesh graded toward the axis and boundary with $N = 4145$ degrees of freedom and plot the numerical solution $u_h$ in Figure 9. The behavior of $u_h$ shows that the solution $u$ has discontinuities near the boundary and across the edges inside $\Omega$ where $f$ is discontinuous.

7.6.3.   Effect of boundary curvature

Next, we numerically address the effect of boundary curvature over nonlocal minimal graphs. For this purpose, we present examples of graphs with prescribed nonlocal mean curvature in several two-dimensional domains, in which we fix $g = 0$ and $f = -1$.

Consider the annulus $\Omega = B(0, 1) \setminus B(0, 1/4)$ and $s = 0.25$. The top row in Figure 10 offers a top view of the discrete solution $u_h$ and a radial slice of it. We observe that the discrete solution is about three times stickier in the inner boundary than in the outer one. The middle and bottom row in Figure 10 display different views of the solution in the square $\Omega = (-1, 1)^2$ for $s = 0.01$. Near the boundary of the domain $\Omega$, we observe a steep slope in the middle of the edges; however, stickiness is not observed at the convex corners of $\Omega$.

We finally investigate stickiness at the boundary of the L-shaped domain $\Omega = (-1, 1)^2 \setminus (0, 1)\times (-1, 0)$ with $s = 0.25, g = 0, f = -1$. We observe in Figure 11 that stickiness is most pronounced at the reentrant corner but absent at the convex corners of $\Omega$.

From these examples we conjecture that there is a relation between the amount of stickiness on $\partial\Omega$ and the nonlocal mean curvature of $\partial\Omega$. Heuristically, let us assume that the Euler-Lagrange equation is satisfied at some point $x \in \partial \Omega$:

where we recall that $G_s$ is defined in (4.1). This fact is not necessarily true, because (6.1) guarantees this identity to hold on $\Omega$ only. Above, we assume that the minimizer is continuous in $\overline{\Omega}$, so that we can set $u(x) : = \lim_{\Omega \ni y \to x} u(y)$. Thus, we can define the stickiness at $x \in \partial \Omega$ as

We point out that in these examples, because the minimizer $u$ attains its maximum on $\Omega^c$ and is constant in that region, we have $M_s \ge 0$. Let $r > 0$ be small, and let us assume that the prescribed curvature is $f(x) = 0$, that we can split the principal value integral in the definition of $H_s$ and that the contribution of the integral on ${\mathbb{R}^d} \setminus B_r(x)$ is negligible compared with that on $B_r(x)$. Then, we must have

If the solution is sticky at $x$, namely $M_s > 0$, then we can approximate

Due to the fact that $G_s\left(\frac{M_s}{|x-y|}\right)$ is strictly increasing with respect to $M_s$, we can heuristically argue that stickiness $M_s(x)$ grows with the increase of the ratio

in order to maintain the balance between the integral in $\Omega \cap B_r(x)$ with the one in $\Omega^c \cap B_r(x)$. Actually, if $R(x) < 1$, as happens at convex corners $x\in\partial\Omega$, it might not be possible for these integrals to balance unless $M_s(x) = 0$. This supports the conjecture that the minimizers are not sticky at convex corners.

8.   Concluding remarks

This paper discusses finite element discretizations of the fractional Plateau and the prescribed fractional mean curvature problems of order $s \in (0, 1/2)$ on bounded domains $\Omega$ subject to exterior data being a subgraph. Both of these can be interpreted as energy minimization problems in spaces closely related to $W^{2s}_1(\Omega)$.

We discuss two converging approaches for computing discrete minimizers: a semi-implicit gradient flow scheme and a damped Newton method. Both of these algorithms require the computation of a matrix related to weighted linear fractional diffusion problems of order $s + \frac{1}{2}$. We employ the latter for computations.

A salient feature of nonlocal minimal graphs is their stickiness, namely that they are generically discontinuous across the domain boundary. Because our theoretical results do not require meshes to be quasi-uniform, we resort to graded meshes to better capture this phenomenon. Although the discrete spaces consist of continuous functions, our experiments in Section 7.1 show the method's capability of accurately estimating the jump of solutions across the boundary. In Section 7.5 we illustrate a geometric rigidity result: wherever the nonlocal minimal graphs are continuous in the boundary of the domain, they must also match the slope of the exterior data. Fractional minimal graphs may change their convexity within $\Omega$, as indicated by our experiments in Section 7.4.

The use of graded meshes gives rise to poor conditioning, which in turn affects the performance of iterative solvers. Our experimental findings reveal that using diagonal preconditioning alleviates this issue, particularly when the grading is not too strong. Preconditioning of the resulting linear systems is an open problem.

Because in practice it is not always feasible to exactly impose the Dirichlet condition on ${\mathbb{R}^d} \setminus \Omega$, we study the effect of data truncation, and show that the finite element minimizers $u_h^H$ computed on meshes ${\mathcal{T}_h}$ over computational domains $\Omega_H$ converge to the minimal graphs as $h \to 0$, $H\to 0$ in $W^{2r}_1(\Omega)$ for $r \in [0, s)$. This is confirmed in our numerical experiments.

Our results extend to prescribed minimal curvature problems, in which one needs some assumptions on the given curvature $f$ in order to guarantee the existence of solutions. We present an example of an ill-posed problem due to data incompatibility. Furthermore, our computational results indicate that graphs with discontinuous prescribed mean curvature may be discontinuous in the interior of the domain. We explore the relation between the curvature of the domain and the amount of stickiness, observe that discrete solutions are stickier on concave boundaries than convex ones, and conjecture that they are continuous on convex corners.

Acknowledgments

JPB has been supported in part by NSF grant DMS-1411808 and Fondo Vaz Ferreira grant 2019-068. WL has been supported in part by NSF grant DMS-1411808 and the Patrick and Marguerite Sung Fellowship in Mathematics of the University of Maryland. RHN has been supported in part by NSF grants DMS-1411808 and DMS-1908267.

Conflict of interest

The authors declare no conflict of interest.