Finite element algorithms for nonlocal minimal graphs

We discuss computational and qualitative aspects of the fractional Plateau and the prescribed fractional mean curvature problems on bounded domains subject to exterior data being a subgraph. We recast these problems in terms of energy minimization, and we discretize the latter with piecewise linear finite elements. For the computation of the discrete solutions, we propose and study a gradient flow and a Newton scheme, and we quantify the effect of Dirichlet data truncation. We also present a wide variety of numerical experiments that illustrate qualitative and quantitative features of fractional minimal graphs and the associated discrete problems.

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 ∈ (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 (−∆) s denotes the fractional Laplacian of order s ∈ (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 ∈ (0, 1/2) converge (again, after rescaling) to motion by fractional mean curvature. Fractional minimal sets also arise in the Γ-limit of nonlocal Ginzburg-Landau energies [28]. We now make the definition of fractional perimeter precise. Let s ∈ (0, 1/2) and two sets A, B ⊂ R d , d ≥ 1. Then, the fractional perimeter of order s of A in B is Given some set A 0 ⊂ R d \ 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 = Ω×R, where Ω is a bounded and sufficiently smooth domain, with exterior data being a subgraph, for some continuous function g : R d \ Ω → R. We briefly remark some key features of this problem: • A technical difficulty arises immediately: all sets A that coincide with A 0 in R d+1 \B have infinite s-perimeter in B. To remedy this issue, one needs to introduce the notion of locally minimal sets [24]. • 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 Ω c . Therefore, an equivalent formulation to the Plateau problem for nonlocal minimal graphs consists in finding a function u : R d → R, with the constraint u = g in Ω c , such that it minimizes the strictly convex energy 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 ∂Ω. As shown by Dipierro, Savin and Valdinoci [17], stickiness is indeed the typical behavior of nonlocal minimal graphs in case Ω ⊂ R. When Ω ⊂ 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 ∈ C ∞ (Ω); 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 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 ∈ (0, 1/2) and g ∈ L ∞ (Ω c ) be given. We consider the space The space V g can be understood as that of functions in W 2s 1 (Ω) with 'boundary value' g. The seminorm in V g does not take into account interactions over Ω c ×Ω 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 (Ω c ). In particular, g may not decay at infinity. In case g is the zero function, the space V g coincides with the standard zero-extension Sobolev space W 2s 1 (Ω); for consistency of notation, we denote such a space by V 0 .
For convenience, we introduce the following notation: given a function u ∈ V g , the form a u : It is worth noticing that G s (ρ) → 0 as |ρ| → ∞. Thus, the weight in (2.1) degenerates whenever the difference quotient |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 ∈ V g such that In light of the previous considerations, equation (2.3) can be regarded as a fractional diffusion problem of order s + 1/2 in 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 ∈ g + H 1 0 (Ω) such that The integral above can be interpreted as a weighted H 1 -form, where the weight depends on u and degenerates as |∇u| → ∞.
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 |u(x)−u(y)| |x−y| → ∞. In this sense, it is not surprising that the fractional-order problems converge to the local ones as s → 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: supp(g) ⊂ Λ, for some bounded set Λ.
The discussion of approximations in case of unboundedly supported data is postponed to Section 5.
3.1. Discrete setting. We consider a family {T h } h>0 of conforming and simplicial triangulations of Λ, and we assume that all triangulations in {T h } h>0 mesh Ω exactly. Moreover, we assume {T h } h>0 to be shape-regular, namely: where h T = diam(T ) and ρ T is the diameter of the largest ball contained in the element T ∈ T h . The vertices of T h will be denoted by {x i }, and the star or patch of {x i } is defined as where ϕ i is the nodal basis function corresponding to the node x i . To impose the condition u = g in Ω c at the discrete level, we introduce the exterior interpolation operator . Thus, Π c h g(x i ) coincides with the standard Clément interpolation of g on x i for all nodes x i such that S i ⊂ R d \ Ω. On the other hand, for nodes x i ∈ ∂Ω, Π c h only averages over the elements in S i that lie in Ω c .
We consider discrete spaces consisting of piecewise linear functions over T h , To account for the exterior data, we define the discrete counterpart of V g , With the same convention as before, we denote by V 0 h the corresponding space in case g ≡ 0. Therefore, the discrete weak formulation reads: find 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 ∈ L ∞ (Ω 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 V g h . 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 I h : V g → V g h that combines the exterior Clément interpolation (3.1) with an interior interpolation operator. More precisely, we set h involves averaging over element stars contained in Ω. Because the minimizer u is smooth in the interior of Ω, but we have no control on its boundary behavior other than the global bound u ∈ W 2s 1 (Ω), 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 ∈ (0, 1/2), Ω be a bounded domain, g ∈ C(Ω c ), and u be the solution to (2.3). Then, the interpolation operator (3.3) satisfies Once we have proved the convergence of I h u to u, energy consistency follows immediately. Since the energy dominates the W 2s 1 (Ω)-norm [7, Lemma 2.5], we can prove convergence in W 2r 1 (Ω) for all r ∈ [0, s) by arguing by compactness. Theorem 3.3 (convergence). Assume s ∈ (0, 1/2), Ω ⊂ R d is a C 1,1 domain, and g ∈ C c (R d ). Let u be the minimizer of I s on V g and u h be the minimizer of I s on V g h . Then, it holds that lim h→0 u − u h W 2r 1 (Ω) = 0, ∀r ∈ [0, s).
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 Ω. 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 ∂Ω. Thus, it is convenient to a priori adapt meshes to better capture the jump of u near ∂Ω.
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 µ ≥ 1. Then, we consider meshes T h such that every element T ∈ T h satisfies These meshes, typically with µ = 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 Ω, because we need to impose the exterior condition u = g on Ω c . As shown in [4,9], the construction In our applications, because Theorem 3.3 gives no theoretical convergence rates, we are not restricted to the choice µ = 2 in two dimensions: a higher µ 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 µ increases, cf. Subsection 7.2.

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 (·, ·), 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 (·, ·) appearing in (3.2). We point out that, following the implementation techniques from [1,2], if we are given

is linear in v h and the latter function can be written in the form
We split Q Ω = (Ω × Ω) ∪ (Ω × Ω c ) ∪ (Ω c × Ω) and, because 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 N Ω elements mesh Ω, while the remaining N Λ − N Ω mesh Λ \ Ω, 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 × 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 × Λ c are of the form where the weight function ω is defined as Since the only restriction on the set Λ is that supp(g) ⊂ Λ, without loss of generality we assume that Λ = B R is a d-dimensional ball with radius R. In such a case, the integral over Λ c can be transformed using polar coordinates into: where ρ 0 (e, x) is the distance from x to ∂B R in the direction of e, which is given by the formula The integral over (ρ 0 (e, x), ∞) can be transformed to an integral over (0, 1) by means of the change of variable ρ = ρ 0 (e, x) ρ −1/(2s) , and then approximated by Gaussian quadrature. Combining this approach with suitable quadrature over ∂B 1 and T l , we numerically compute the integral over T l × Λ c for a given u h .

Gradient Flow. Although we can compute
h , a gradient flow is a feasible approach to solve for the unique minimizer u h .
Given α ∈ [0, 1), and with the convention that H 0 = L 2 , we first consider a time-continuous H α -gradient flow for u h (t), namely u j (t)ϕ i , local existence and uniqueness of solutions in time for (4.2) follow from the fact that a u h (u h , ϕ i ) is Lipschitz with respect to u j for any ϕ 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 τ > 0 and iteration h makes (4.3) amenable for its computational solution. The following proposition proves the stability of the semiimplicit scheme. Its proof mimics the one of classical mean curvature flow [18].
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 G s satisfies G s (r) = 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 | G s (|r 0 |) above and use that G s is even, decreasing on [0, ∞) and non-negative, to obtain This proves (4.5).
Combining this with (4.4) finishes the proof.
Upon writing w k h := u k+1 h − u k , the semi-implicit scheme (4.3) becomes (4.6), which is the crucial step of Algorithm 4.1 to solve (3.2). Equation (4.6) boils down to Algorithm 4.1 Semi-implicit gradient flow 1: Select an arbitrary initial u 0 h ∈ V g h , let k = 0, and set w 0 h H α (Ω) = Inf. Choose a time step τ > 0 and a small number ε > 0.
is just a mass matrix, while if α > 0, M is the stiffness matrix for the linear fractional diffusion problem of order α, given by The matrix K k = K k ij is the stiffness matrix for a weighted linear fractional diffusion of order s + 1 2 , whose elements K k ij := a u k h (ϕ i , ϕ j ) are given by and can be computed as described in Subsection 4.1. The right hand side vector is Because of Proposition 4.1 (stability of H α -gradient flow), the loop in Algorithm 4.1 terminates in finite steps. Moreover, using the continuity of a u k h (·, ·) in [H 1 2 +s (Ω)] 2 , which is uniform in u k h , together with an inverse estimate and 0 ≤ α ≤ 1 2 + s gives where the hidden constant depends on the mesh shape-regularity and h min is the minimum element size. Therefore, the last iterate u k h of Algorithm 4.1 satisfies the residual estimate

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.

Algorithm 4.2 Damped Newton Algorithm
1: Select an arbitrary initial u 0 h ∈ V g h and let k = 0. Choose a small number ε > 0.

4:
Determine the minimum n ∈ N such that u k,n h := u k h + 2 −n w k h satisfies

6: end while
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 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 ∈ V g and v, w ∈ V 0 , its discrete counterpart  Proof. Since I s [u h ] is strictly convex, the convergence of u k h 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 4.2 is to solve the equation (4.7). Due to the linearity of with respect to v h and w k h , we just need to solve a linear system 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, is the stiffness matrix for a weighted linear fractional diffusion of order s + 1 2 . Since the only difference with the semi-implicit gradient flow algorithm is the weight, the elements in 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 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 Ω H to be a bounded open domain containing Ω and such that d(x, Ω) H for all x ∈ ∂Ω H , and choose a cutoff function η H ∈ C ∞ (Ω c ) satisfying We replace g by g H := gη H , and consider problem (2.3) using g H as Dirichlet condition. Let u H ∈ 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 ∈ [0, s), u H h → u H in W 2r 1 (Ω) as h → 0. Therefore we only need to show that, in turn, the minimizers of the truncated problems satisfy u H → u as H → ∞ 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 Proof. We prove only the first estimate, as the second one follows in the same fashion.
From definition (1.2), it follows immediately that F s (0) = F s (0) = 0, and thus F s (ρ) ≤ Cρ 2 if ρ 1. Combining this with the fact that |v( for a.e. x ∈ Ω, y ∈ Ω 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. Proof. Since u H is the minimizer over V g H and T H u ∈ V g H , we deduce Conversely, using that u is the minimzer over V g and Eu H ∈ V g , we obtain and thus conclude the proof.
The energy I s is closely related to the W 2s 1 (Ω)-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 (Ω) for all r ∈ [0, s). 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 (Ω). It follows that, up to a subsequence, u H converges in L 1 (Ω) to a limit u. Also, because u H = g in Ω H , we can extend u by g on Ω c , and have u H → u a.e in R d . We then can invoke Fatou's lemma and Corollary 5.2 to deduce that Because u ∈ V g , we deduce that u = u whence u H → u in L 1 (Ω) as H → 0. By interpolation, we conclude that convergence in W 2r 1 (Ω) holds for all r ∈ [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 ∈ Ω, a certain nonlocal mean curvature at x, u(x) is equal to a given function f (x). For a set E ⊂ R d+1 and x ∈ ∂E, such nonlocal mean curvature operator is defined as [12] H s [E]( x) := P.V.
In turn, for 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 → 1 2 − , it is necessary to normalize the operator H s accordingly. Let α d denote the volume of the d-dimensional unit ball, and consider the prescribed nonlocal mean curvature problem (6.1) The scaling factor 1−2s dα d yields [7, Lemma 5.8] where H[E] denotes the classical mean curvature operator. Therefore, in the limit s → 1 2 − , 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 ∈ 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 ∈ V g is the minimizer of K s [·; f ] if and only if it satisfies for every v ∈ 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 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 ε 0 > 0 such that for every measurable set A ⊂ Ω, On the other hand, Hölder's inequality gives This is 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)).

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 Ω, and we show that graphs with prescribed fractional mean curvature may be discontinuous in Ω.
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 Ω = (−1, 1) ⊂ R and g(x) = M sign(x), where M > 0. Reference [15] proves that, for every s ∈ (0, 1/2), stickiness (i.e. the solution being discontinuous at ∂Ω) 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. Subsection 3.3) with parameter µ = 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 Ω = (−1, 1) into |Ω| 1/µ h subintervals and the smallest ones have size h µ . 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 Ω 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 Ω, 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∈Ω 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. 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 , µ = 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.

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 κ(K), where κ(K) is the condition number of the stiffness matrix K. For linear problems of order s involving the fractional Laplacian (−∆) s , the condition number of K satisfies [3] Reference [3] also shows that diagonal preconditioning yields κ(K) = O N 2s/d , where N is the dimension of the finite element space.
Using the Matlab function condest, we estimate the condition number of the Jacobian matrix in the last Newton iteration in the example from Subsection 7.1 with M = 1, with and without diagonal preconditioning. Figure 2 summarizes that, according to (7.1), is expected to be equal to 9/14 ≈ 0.643.  (2.3)). While the conditioning for graded meshes is significantly poorer, when µ = 2 diagonal preconditioning recovers condition numbers comparable to the ones on quasi-uniform meshes.
Let N = dimV g h 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 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 µ = 2, µ = 3), the behavior is less clear. When using diagonal preconditioning for µ = 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 → 0, H → ∞. Here, we study numerically the effect of data truncation by running experiments on a simple two-dimensional problem. Consider Ω = B 1 ⊂ R 2 and g ≡ 1; then, the nonlocal minimal graph u is a constant function. For H > 0, we set Ω H = B H+1 . and compute nonlocal minimal graphs on Ω with Dirichlet data g H = χ Ω H , which is a truncation of g ≡ 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 Ω is strictly less than 1. For s = 0.25, we plot the L 1 (Ω) and L ∞ (Ω) 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 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 Ω = (−1, 1), s = 0.02, g(x) = 1 for a ≤ |x| ≤ 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 Ω and has stickiness on the boundary. In addition, the figure confirms that lim x→1 − u a (x) = ∞, 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 x→1 − u a (x) = −∞ since g(x) = 0 near the boundary of Ω. This fact implies that u(x) cannot be convex near x = 1 for 1 < a < 2. Furthermore, as a → 1 + one expects that u a (x) → u 1 (x) and thus that u a be convex in the interior of Ω = (−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 Ω. In fact, we see from the right panel in Figure 4 that the nonlocal minimal graph u in Ω 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 Ω ⊂ R 2 be the unit ball, s = 0.25, and g(x) = 1 for 129 128 ≤ |x| ≤ 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 ∂Ω. An argument analogous to the one we discussed in the previous paragraph also explains this behavior in a twodimensional 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 Ω ⊂ R 2 , if nonlocal minimal graphs are continuous at some point x ∈ ∂Ω 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.
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 µ = 2 −9 . The flattening of the curves as x → 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    Table 2 Example of Subsection 7.5: experimental slopes ∂yu h (x, 0) for x = 2 −k and k = 6, . . . , 9. As x → 0 + , these slopes become smaller; this geometric rigidity is easier to capture for larger s. 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 Ω = B(0, 1) ⊂ 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 K s [u; f ] goes from 0 initially down to −6.6 × 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 Ω. We present the following two examples for d = 1 and d = 2.
In first place, let Ω = (−1, 1) ⊂ R, s = 0.01, g = 0 and consider f (x) = 1.5 sign(x). We use a mesh graded toward x = 0, ±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 = ±1 and x = 0. As a second illustration of interior discontinuities, let Ω = (−1, 1) 2 ⊂ R 2 , s = 0.01, g = 0 and consider f (x, y) = 4 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 Ω 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 Ω = B(0, 1) \ 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 Ω = (−1, 1) 2 for s = 0.01. Near the boundary of the domain Ω, we observe a steep slope in the middle of the edges; however, stickiness is not observed at the convex corners of Ω.
We finally investigate stickiness at the boundary of the L-shaped domain Ω = (−1, 1) 2 \ (0, 1) × (−1, 0) with s = 0.25, g = 0, f = −1. We observe in Figure 11 that  The plot in the top-right panel corresponds to a radial slice (y = 0, 0.25 ≤ x ≤ 1) of the annulus, while the ones in the bottom-left and bottom-right show slices along the diagonal (0 ≤ y = x ≤ 1) and perpendicular to an edge of the square (y = 0.5, 0 ≤ x ≤ 1), respectively. We observe that stickiness is larger near the concave portions of the boundary than near the convex ones, and that it is absent in the corners of the square.
in order to maintain the balance between the integral in Ω ∩ B r (x) with the one in Ω c ∩ B r (x). Actually, if R(x) < 1, as happens at convex corners x ∈ ∂Ω, 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 ∈ (0, 1/2) on bounded domains Ω 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 (Ω). We discuss two converging approaches for computing discrete minimizers: a semiimplicit gradient flow scheme and a damped Newton method. Both of these algorithms Fig. 11. Stickiness on the L-shaped domain Ω = (−1, 1) 2 \ (0, 1) × (−1, 0) with prescribed fractional mean curvature f = −1 in Ω and Dirichlet condition g = 0 in Ω c . The plots in the middle correspond to slices along y = x and y = −x respectively, while the ones in the bottom are slices along x = 0 or y = 0.5 respectively. We see that the largest stickiness takes place at the reentrant corner while there is no stickiness at the convex corners.
require the computation of a matrix related to weighted linear fractional diffusion problems of order s + 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 Subsection 7.1 show the method's capability of accurately estimating the jump of solutions across the boundary. In Subsection 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 Ω, as indicated by our experiments in Subsection 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 R d \Ω, we study the effect of data truncation, and show that the finite element minimizers u H h computed on meshes T h over computational domains Ω H converge to the minimal graphs as h → 0, H → 0 in W 2r 1 (Ω) for r ∈ [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.