Loading [MathJax]/jax/element/mml/optable/BasicLatin.js
Research article

Finite element algorithms for nonlocal minimal graphs

  • Received: 26 May 2021 Accepted: 07 July 2021 Published: 20 July 2021
  • 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.

    Citation: Juan Pablo Borthagaray, Wenbo Li, Ricardo H. Nochetto. Finite element algorithms for nonlocal minimal graphs[J]. Mathematics in Engineering, 2022, 4(2): 1-29. doi: 10.3934/mine.2022016

    Related Papers:

    [1] Ko-Shin Chen, Cyrill Muratov, Xiaodong Yan . Layered solutions for a nonlocal Ginzburg-Landau model with periodic modulation. Mathematics in Engineering, 2023, 5(5): 1-52. doi: 10.3934/mine.2023090
    [2] Daniele Cerroni, Florin Adrian Radu, Paolo Zunino . Numerical solvers for a poromechanic problem with a moving boundary. Mathematics in Engineering, 2019, 1(4): 824-848. doi: 10.3934/mine.2019.4.824
    [3] Qiang Du, Tadele Mengesha, Xiaochuan Tian . Lp compactness criteria with an application to variational convergence of some nonlocal energy functionals. Mathematics in Engineering, 2023, 5(6): 1-31. doi: 10.3934/mine.2023097
    [4] Serena Della Corte, Antonia Diana, Carlo Mantegazza . Global existence and stability for the modified Mullins–Sekerka and surface diffusion flow. Mathematics in Engineering, 2022, 4(6): 1-104. doi: 10.3934/mine.2022054
    [5] Tommaso Tassi, Alberto Zingaro, Luca Dede' . A machine learning approach to enhance the SUPG stabilization method for advection-dominated differential problems. Mathematics in Engineering, 2023, 5(2): 1-26. doi: 10.3934/mine.2023032
    [6] L. Dieci, Fabio V. Difonzo, N. Sukumar . Nonnegative moment coordinates on finite element geometries. Mathematics in Engineering, 2024, 6(1): 81-99. doi: 10.3934/mine.2024004
    [7] Miyuki Koiso . Stable anisotropic capillary hypersurfaces in a wedge. Mathematics in Engineering, 2023, 5(2): 1-22. doi: 10.3934/mine.2023029
    [8] Anoumou Attiogbe, Mouhamed Moustapha Fall, El Hadji Abdoulaye Thiam . Nonlocal diffusion of smooth sets. Mathematics in Engineering, 2022, 4(2): 1-22. doi: 10.3934/mine.2022009
    [9] Thomas J. Radley, Paul Houston, Matthew E. Hubbard . Quadrature-free polytopic discontinuous Galerkin methods for transport problems. Mathematics in Engineering, 2024, 6(1): 192-220. doi: 10.3934/mine.2024009
    [10] Matteo Novaga, Marco Pozzetta . Connected surfaces with boundary minimizing the Willmore energy. Mathematics in Engineering, 2020, 2(3): 527-556. doi: 10.3934/mine.2020024
  • 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.



    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

    tuε+(Δ)suε+f(uε)ε1+2s=0,

    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,BRd, d1. Then, the fractional perimeter of order s of A in B is

    Ps(A;B):=12QB|χA(x)χA(y)||xy|d+2sdydx,

    where QB=(Rd×Rd)(Bc×Bc) and Bc=RdB. Given some set A0RdB, a natural problem is how to extend A0 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,

    A0={(x,xd+1):xd+1<g(x),xRdΩ},

    for some continuous function g:RdΩR. We briefly remark some key features of this problem:

    A technical difficulty arises immediately: all sets A that coincide with A0 in Rd+1B 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.

    If the exterior datum g is a bounded function, then one can replace the infinite cylinder B=Ω×R by a truncated cylinder BM=Ω×(M,M) for some M>0 sufficiently large [25,Proposition 2.5].

    Let A be the subgraph of a certain function v that coincides with g on Ωc. One can rewrite Ps(A,BM) as

    Ps(A,BM)=Is[v]+C(M,d,s,Ω,g),

    where Is 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:RdR, with the constraint u=g in Ωc, such that it minimizes the strictly convex energy

    Is[u]:=QΩFs(u(x)u(y)|xy|)1|xy|d+2s1dxdy, (1.1)

    where Fs is defined as

    Fs(ρ):=ρ0ρr(1+r2)(d+1+2s)/2dr. (1.2)

    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 ΩR2, 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 uC(Ω); 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 W2r1(Ω), r[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 Ω. 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.

    We now specify the problem we aim to solve in this paper and pose its variational formulation. Let s(0,1/2) and gL(Ωc) be given. We consider the space

    Vg:={v:RdR:v|ΩW2s1(Ω), v=g in Ωc},

    equipped with the norm

    vVg:=vL1(Ω)+|v|Vg,

    where

    |v|Vg:=QΩ|v(x)v(y)||xy|d+2sdxdy,

    where QΩ=(Rd×Rd)(Ω×Ω). The space Vg can be understood as that of functions in W2s1(Ω) with 'boundary value' g. The seminorm in Vg 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 W2s1(Ωc). In particular, g may not decay at infinity. In case g is the zero function, the space Vg coincides with the standard zero-extension Sobolev space ˜W2s1(Ω); for consistency of notation, we denote such a space by V0.

    For convenience, we introduce the following notation: given a function uVg, the form au:Vg×V0R is

    au(w,v):=QΩ˜Gs(u(x)u(y)|xy|)(w(x)w(y))(v(x)v(y))|xy|d+1+2sdxdy, (2.1)

    where

    ˜Gs(ρ):=10(1+ρ2r2)(d+1+2s)/2dr. (2.2)

    It is worth noticing that ˜Gs(ρ)0 as |ρ|. Thus, the weight in (2.1) degenerates whenever the difference quotient |u(x)u(y)||xy| blows up.

    The weak formulation of the fractional minimal graph problem can be obtained by the taking first variation of Is[u] in (1.1) in the direction v. As described in [7], that problem reads: find uVg such that

    au(u,v)=0vV0. (2.3)

    In light of the previous considerations, Eq (2.3) can be regarded as a fractional diffusion problem of order s+1/2 in Rd 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 ug+H10(Ω) such that

    Ω11+|u|2uvdx=0,vH10(Ω).

    The integral above can be interpreted as a weighted H1-form, where the weight depends on u and degenerates as |u|.

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

    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.

    We consider a family {Th}h>0 of conforming and simplicial triangulations of Λ, and we assume that all triangulations in {Th}h>0 mesh Ω exactly. Moreover, we assume {Th}h>0 to be shape-regular, namely:

    σ=sup

    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

    S_i : = \text{supp} ( \varphi_i) ,

    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

    \begin{equation} \Pi_h^c g : = \sum\limits_{{\texttt{x}}_i \in {\mathcal{N}_h^c}} (\Pi_h^{{\texttt{x}}_i} g) ({\texttt{x}}_i) \; \varphi_i, \end{equation} (3.1)

    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} ,

    \mathbb{V}_h : = \{ v \in C(\Lambda) \colon v|_T \in \mathcal{P}_1 \; \forall T \in {\mathcal{T}_h} \}.

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

    \mathbb{V}_h^g : = \{ v \in \mathbb{V}_h \colon \ v|_{\Lambda \setminus \Omega} = \Pi_h^c 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

    \begin{equation} a_{u_h}(u_h, v_h) = 0 \quad \text{for all } v_h \in \mathbb{V}^0_h. \end{equation} (3.2)

    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 .

    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

    \begin{equation} {\mathcal I_h} v : = \Pi_h^\circ \left(v \big|_\Omega \right) + \Pi_h^c g , \end{equation} (3.3)

    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

    \iint_{Q_{\Omega}} \frac{\left| ({\mathcal I_h} u - u)(x) - ({\mathcal I_h} u - u)(y) \right|}{|x-y|^{d+2s}} \; dxdy \to 0 \quad \mathit{\text{as}} \; h \to 0.

    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

    \lim\limits_{h \to 0} \| u - u_h \|_{W^{2r}_1(\Omega)} = 0, \quad \forall r \in [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.

    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

    \begin{equation} h_T \approx \left\lbrace \begin{array}{ll} C(\sigma) h^\mu, & \overline T \cap \partial \Omega \neq \emptyset \\ C(\sigma) h {\text{dist}}(T, \partial\Omega)^{(\mu-1)/\mu}, & \overline T \cap \partial \Omega = \emptyset. \end{array} \right. \end{equation} (3.4)

    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

    \dim \mathbb{V}_h^g \approx \left\lbrace \begin{array}{ll} h^{(1-d)\mu}, & \mu > \frac{d}{d-1}, \\ h^{-d}|\log h|, & \mu = \frac{d}{d-1}, \\ h^{-d}, & \mu < \frac{d}{d-1}. \end{array} \right.

    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.

    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.

    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

    \begin{aligned} a_i & : = a_{u_h}(u_h, \varphi_i) \\ & = \iint_{Q_{\Omega}} {\widetilde{G}_s}\left(\frac{u_h(x)-u_h(y)}{|x-y|}\right) \frac{(u_h(x)-u_h(y))(\varphi_i(x)-\varphi_i(y))}{|x-y|^{d+1+2s}}dxdy . \end{aligned}

    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

    \begin{aligned} a_i = & \iint_{\Omega \times \Omega} {\widetilde{G}_s}\left(\frac{u_h(x)-u_h(y)}{|x-y|}\right) \frac{(u_h(x)-u_h(y))(\varphi_i(x)-\varphi_i(y))}{|x-y|^{d+1+2s}}dxdy \\ & + 2 \iint_{\Omega \times \Omega^c} {\widetilde{G}_s}\left(\frac{u_h(x)-g_h(y)}{|x-y|}\right) \frac{(u_h(x)-g_h(y)) \varphi_i(x)}{|x-y|^{d+1+2s}}dxdy = : a_{i,1} + 2a_{i,2} . \end{aligned}

    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,

    \bigcup\limits_{1 \le i \le {\mathcal{N}_{\Omega}}} \overline{T_i} = \overline{\Omega} \quad \bigcup\limits_{{\mathcal{N}_{\Omega}} + 1 \le i \le {\mathcal{N}_{\Lambda}}} \overline{T_i} = \overline{\Lambda} \setminus \Omega.

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

    \begin{aligned} a_{i,1} = & \sum\limits_{l,m = 1}^{{\mathcal{N}_{\Omega}}} \iint_{T_l \times T_m} {\widetilde{G}_s}\left(\frac{u_h(x)-u_h(y)}{|x-y|}\right) \frac{(u_h(x)-u_h(y))(\varphi_i(x)-\varphi_i(y))}{|x-y|^{d+1+2s}}dxdy ,\\ a_{i,2} = & \sum\limits_{l = 1}^{{\mathcal{N}_{\Omega}}}\sum\limits_{m = {\mathcal{N}_{\Omega}} + 1}^{{\mathcal{N}_{\Lambda}}} \iint_{T_l \times T_m} {\widetilde{G}_s}\left(\frac{u_h(x)-{g_h(y)}}{|x-y|}\right) \frac{(u_h(x)-{g_h(y)}) \varphi_i(x)}{|x-y|^{d+1+2s}}dxdy \\ & + \sum\limits_{l = 1}^{{\mathcal{N}_{\Omega}}} \iint_{T_l \times \Lambda^c} {\widetilde{G}_s}\left(\frac{u_h(x)}{|x-y|}\right) \frac{u_h(x)\varphi_i(x)}{|x-y|^{d+1+2s}}dxdy . \end{aligned}

    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

    \int_{T_l} \varphi_i(x) \omega(x) dx ,

    where the weight function \omega is defined as

    \begin{equation} \begin{split} \omega(x) & : = \int_{\Lambda^c} G_s\left(\frac{u_h(x)}{|x-y|}\right) \frac{1}{|x-y|^{d+2s}} dy , \\ G_s (\rho) & : = \int_0^\rho (1+ r^2)^{-(d+1+2s)/2} dr = \rho \, {\widetilde{G}_s} (\rho). \end{split} \end{equation} (4.1)

    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:

    w(x) = \int_{\partial B_{1}} dS(e) \int_{\rho_0(e,x)}^{\infty} G_s\left(\frac{u_h(x)}{\rho}\right) \rho^{-1-2s} d\rho,

    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

    \rho_0(e,x) = \sqrt{R^2 - |x|^2 + (e \cdot x)^2} - e \cdot x .

    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 .

    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

    \begin{equation} \langle \partial_t u_h, v_h \rangle_{H^{\alpha}(\Omega)} = -\frac{\delta I_s}{\delta u_h}(v_h) = -a_{u_h}(u_h, v_h), \qquad \forall v_h \in \mathbb{V}^0_h , \end{equation} (4.2)

    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

    \frac{d}{dt} I_s[u_h] = \frac{\delta I_s[u_h]}{\delta u_h}(\partial_t u_h) = a_{u_h}(u_h, \partial_t u_h) = -\langle \partial_t u_h, \partial_t u_h \rangle_{H^{\alpha}(\Omega)} \le 0 ,

    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

    \begin{equation} \frac1{\tau} \langle {u^{k+1}_h - u^k_h} \ , \ v_h \rangle_{H^{\alpha}(\Omega)} = -a_{u^k_h}(u^{k+1}_h \ , \ v_h), \qquad \forall v_h \in \mathbb{V}^0_h . \end{equation} (4.3)

    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,

    \begin{equation*} \label{E:stab-semi-implicit} I_s[u^{k+1}_h] + \frac{1}{\tau} \Vert u^{k+1}_h - u^k_h \Vert^2_{H^{\alpha}(\Omega)} \le I_s[u^k_h]. \end{equation*}

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

    \begin{equation} \frac{1}{\tau} \Vert u^{k+1}_h - u^k_h\Vert^2_{H^{\alpha}(\Omega)} = -a_{u^k_h}(u^{k+1}_h \, , \, u^{k+1}_h - u^k_h) . \end{equation} (4.4)

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

    \begin{equation} (r_1^2 - r_1 r_0) \ {\widetilde{G}_s}(r_0) \ge F_s(r_1) - F_s(r_0). \end{equation} (4.5)

    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

    \begin{aligned} F_s(r_1) - F_s(r_0) & = F_s(|r_1|) - F_s(|r_0|) \\ &\le F_s(|r_1|) - \left[ F_s(|r_1|) + (|r_0| - |r_1|) \ G_s(|r_1|) \right] \\ & = (|r_1| - |r_0|) \ |r_1| \ {\widetilde{G}_s}(|r_1|). \end{aligned}

    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

    \begin{aligned} F_s(r_1) - F_s(r_0) & \le (|r_1| - |r_0|) \ |r_1| \ {\widetilde{G}_s}(|r_0|) + |r_1| \ (|r_1| - |r_0|) \left( {\widetilde{G}_s}(|r_1|) - {\widetilde{G}_s}(|r_0|) \right) \\ &\le (|r_1| - |r_0|) \ |r_1| \ {\widetilde{G}_s}(|r_0|) \\ & = (r_1^2 - |r_0|\ |r_1|) \ {\widetilde{G}_s}(|r_0|) \\ &\le (r_1^2 - r_0 r_1) \ {\widetilde{G}_s}(r_0) . \end{aligned}

    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

    \begin{aligned} a_{u^k_h}(u^{k+1}_h \ , \ u^{k+1}_h - u^k_h) & = \iint_{Q_{\Omega}} {\widetilde{G}_s}\left( d_k(x,y) \right) \frac{d_{k+1}(x,y) (d_{k+1}(x,y) - d_{k}(x,y)) }{|x-y|^{d-1+2s}} dxdy \\ &\ge \iint_{Q_{\Omega}} \frac{F_s(d_{k+1}(x,y)) - F_s(d_{k}(x,y))}{|x-y|^{d-1+2s}} dxdy \\ & = I_s[u^{k+1}_h] - I_s[u^{k}_h] . \end{aligned}

    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).

    Algorithm 1 Semi-implicit gradient flow
    1: Select an arbitrary initial u_h^0 \in \mathbb{V}^g_h , let k = 0 , and set \| w_h^0 \|_{H^\alpha(\Omega)} = \texttt{Inf} . Choose a time step \tau > 0 and a small number {\varepsilon} > 0 .
    2: while \| w_h^k \|_{H^\alpha(\Omega)} > {\varepsilon} do
    3:  Find w_h^{k+1} \in \mathbb{V}^0_h such that
                 \begin{equation} \langle w_h^{k+1}, v_h \rangle_{H^{\alpha}(\Omega)} + \tau a_{u_h^k}(w_h^{k+1}, v_h) = -a_{u_h^k}(u_h^k, v_h) \, \quad \forall v_h \in \mathbb{V}^0_h. \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; (4.6) \end{equation}
    4:  Set u_h^{k+1} = u_h^k + \tau \; w_h^{k+1} and k = k+1 .
    5: end while

     | Show Table
    DownLoad: CSV

    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

    M_{ij} : = \iint_{Q_{\Omega}} \frac{(\varphi_i(x)-\varphi_i(y))(\varphi_j(x)-\varphi_j(y))}{|x-y|^{d+2\alpha}}dxdy \quad (\alpha > 0).

    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

    K^{k}_{ij} = \iint_{Q_{\Omega}} {\widetilde{G}_s}\left(\frac{u_h^k(x)-u_h^k(y)}{|x-y|}\right) \frac{(\varphi_i(x)-\varphi_i(y))(\varphi_j(x)-\varphi_j(y))}{|x-y|^{d+1+2s}}dxdy ,

    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

    \big| a_{u_h^k}(w_h^{k+1},v_h) \big| \lesssim |w_h^{k+1}|_{H^{\frac12 + s}(\Omega)} |v_h|_{H^{\frac12 + s}(\Omega)} \lesssim h_{\text{min}}^{-1-2s+2\alpha} |w_h^{k+1}|_{H^{\alpha}(\Omega)} |v_h|_{H^{\alpha}(\Omega)},

    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

    \max\limits_{v_h\in\mathbb{V}^0_h} \frac{\big|a_{u_h^k}(u_h^k, v_h)\big|}{\|v_h\|_{H^{\alpha}(\Omega)}} \lesssim {\varepsilon} \Big( 1 + \tau h_{\text{min}}^{-1-2s+2\alpha} \Big) .

    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 2 Damped Newton Algorithm
    1: Select an arbitrary initial u_h^0 \in \mathbb{V}^g_h and let k = 0 . Choose a small number {\varepsilon} > 0 .
    2: while \Vert \{ a(u_h^k, {\varphi}_i) \}_{i=1}^m \Vert_{l^2} > {\varepsilon} do
    3:  Find w_h^k \in \mathbb{V}^0_h such that
             \begin{equation} \frac{\delta a_{u_h}(u_h^k, v_h)}{\delta u_h^k}(w_h^k) = -a_{u_h^k}(u_h^k, v_h), \qquad \forall v_h \in \mathbb{V}^0_h. \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; (4.7) \end{equation}
    4:  Determine the minimum n \in \mathbb{N} such that u_h^{k, n} := u_h^k + 2^{-n} w_h^k satisfies
             \begin{equation*} \Vert \{ a_{u_h^k}(u_h^{k, n}, {\varphi}_i) \}_{i = 1}^m \Vert_{l^2} \leq (1 - 2^{-n-1}) \Vert \{ a_{u_h^k}(u_h^k, {\varphi}_i) \}_{i = 1}^m \Vert_{l^2} \end{equation*}
    5:  Let u_h^{k+1} = u_h^{k, n} and k = k+1 .
    6: end while

     | Show Table
    DownLoad: CSV

    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

    \begin{equation*} \label{E:second-variation-I_s} \frac{\delta a_u(u,v)}{\delta u}(w) = \iint_{Q_{\Omega}} G'_s \left(\frac{u(x)-u(y)}{|x-y|}\right) \frac{(w(x)-w(y))(v(x)-v(y))}{|x-y|^{d+1+2s}}dxdy. \end{equation*}

    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

    \widetilde{K}^k_{ij} = \iint_{Q_{\Omega}} G'_s\left(\frac{u_h^k(x)-u_h^k(y)}{|x-y|}\right) \frac{(\varphi_i(x)-\varphi_i(y))(\varphi_j(x)-\varphi_j(y))}{|x-y|^{d+1+2s}}dxdy ,

    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 .

    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

    \begin{equation*} \label{eq:cutoff} 0\le \eta_H \le 1, \quad \text{supp}(\eta_H)\subset \overline{\Omega}_{H+1} \setminus \Omega, \quad \eta_H(x) = 1 \quad \text{in} \quad \Omega_{H} \setminus \Omega . \end{equation*}

    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) ,

    u_h^H \to u^H \quad \text{in } W^{2r}_1(\Omega) \quad \text{as } h \to 0.

    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:

    \begin{array}{ll} T_H \colon \mathbb{V}^g \to \mathbb{V}^{g_H}, \quad & T_H v = v \eta_H , \\ E_H \colon \mathbb{V}^{g_H} \to \mathbb{V}^g, \quad & E_H w = w + (1 - \eta_H) g . \end{array}

    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}) :

    \begin{aligned} |I_s[v] - I_s[T_H v]| \lesssim H^{-1-2s}, \\ |I_s[w] - I_s[E_H w]| \lesssim H^{-1-2s}. \end{aligned}

    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

    \begin{aligned} |I_s[v] & - I_s[T_H v]| \\ & \le 2 \int_\Omega \int_{\Omega_H^c} \left| F_s\left(\frac{v(x)-v(y)}{|x-y|}\right) - F_s\left(\frac{v(x)- T_Hv(y)}{|x-y|}\right) \right| \frac{1}{|x-y|^{d+2s-1}} \; dydx . \end{aligned}

    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

    \begin{aligned} |I_s[v] - I_s[T_H v]| & \lesssim \|v\|_{L^\infty(\Omega^c)}^2 \int_\Omega \int_{\Omega_H^c} \frac{1}{|x-y|^{d+2s+1}} \;dxdy \lesssim H^{-1-2s}. \end{aligned}

    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

    \left|I_s[u] - I_s[u^H]\right| \lesssim H^{-1-2s}.

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

    I_s[u^H] - I_s[u] \le I_s[T_H u] - I_s[u] \lesssim H^{-1-2s}.

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

    I_s[u] - I_s[u^H] \le I_s [Eu^H] - I_s [u^H] \lesssim H^{-1-2s},

    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

    \lim\limits_{H \to \infty} \|u - u^H \|_{W^{2r}_1(\Omega)} = 0 .

    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

    I_s [\widetilde u ] \le \liminf\limits_{H \to \infty} I_s [ u^H ] \lesssim \liminf\limits_{H \to \infty} I_s [u] + H^{-1-2s} = I_s [u].

    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) .

    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]

    \begin{equation*} \label{E:NMS-def-NonLocalCurv} H_s[E](\widetilde{x}) : = {\text{P.V.}} \int_{{\mathbb{R}^{d+1}}} \frac{\chi_{E^c}(\widetilde{y}) - \chi_{E}(\widetilde{y})}{|\widetilde{x}-\widetilde{y}|^{d+1+2s}} d\widetilde{y}. \end{equation*}

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

    \begin{aligned} H_s[u](x) = {\text{P.V.}} \int_{{\mathbb{R}}^d} G_s\left(\frac{u(x)-u(y)}{|x-y|}\right) \frac{dy}{|x-y|^{d+2s}}. \end{aligned}

    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

    \begin{equation} \left\{\begin{array}{rl} \frac{1-2s}{d \alpha_d}H_s[u](x) = f(x), & x \in \Omega, \\ u(x) = g(x), & x \in {\mathbb{R}^d} \setminus \Omega. \end{array}\right. \end{equation} (6.1)

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

    \begin{equation} \lim\limits_{s \to \frac12^-}\frac{1-2s}{d \alpha_d}H_s[E](x) = H[E](x), \end{equation} (6.2)

    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:

    \begin{equation} \left\{\begin{array}{rl} \frac{1}{d} \; {\text{div}}\Big( \frac{\nabla u (x)}{\big(1+|\nabla u (x)|^2 \big)^{1/2}} \Big) = f(x), & x \in \Omega, \\ u(x) = g(x), & x \in {\partial\Omega}. \end{array}\right. \end{equation} (6.3)

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

    \begin{equation} \mathcal{K}_{s}[u;f] : = I_s[u] - \frac{d \alpha_d}{1-2s}\int_{\Omega} f(x) u(x) dx. \end{equation} (6.4)

    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

    \begin{equation} \begin{aligned} 0 & = a_u(u,v) - \frac{d \alpha_d}{1-2s}\int_{\Omega} f(x) v(x) dx \\ & = \iint_{Q_{\Omega}} G_s\left(\frac{u(x)-u(y)}{|x-y|}\right) \frac{v(x)-v(y)}{|x-y|^{d+2s}}dx dy - \frac{d \alpha_d}{1-2s}\int_{\Omega} f(x) v(x) dx \end{aligned} \end{equation} (6.5)

    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 ,

    \begin{equation} \Big| \int_A f(x)dx \Big| \le \frac{(1 - {\varepsilon}_0)}{d} \; \mathcal{H}^{d-1}(\partial A), \end{equation} (6.6)

    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

    \lim\limits_{h \to 0} \| u - u_h \|_{W^{2r}_1(\Omega)} = 0, \quad \forall r \in [0,s).

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

    I_s[u] + C_0 \ge C_1 \|u\|_{W^{2s}_1(\Omega)} ,

    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

    \int_{\Omega} f(x)u(x) dx \le \Vert u \Vert_{L^{d/(d-2s)}(\Omega)} \|f\|_{L^{d/(2s)}(\Omega)} \le C_2 \|u\|_{W^{2s}_1(\Omega)} \|f\|_{L^{d/(2s)}(\Omega)},

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

    \begin{equation} \mathcal{K}_{s}[u;f] \ge \|u\|_{W^{2s}_1(\Omega)} \left(C_1 - C_2\|f\|_{L^{d/(2s)}(\Omega)}\right) - C_0. \end{equation} (6.7)

    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

    \Big| \int_A f(x)dx \Big| \le \left( \int_A 1 dx \right)^{\frac{d-1}{d}} \left( \int_A |f(x)|^d dx \right)^{\frac{1}{d}} \lesssim \mathcal{H}^{d-1}(\partial A) \Vert f \Vert_{L^d(\Omega)},

    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)).

    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.

    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

    \begin{equation} \sup\limits_{x \in \Omega} u^M(x) < c_0 M^{\frac{1+2s}{2+2s}}, \quad \inf\limits_{x \in \Omega} u^M(x) > -c_0 M^{\frac{1+2s}{2+2s}}. \end{equation} (7.1)

    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 .

    Figure 1.  Stickiness in 1d. In the setting of Section 7.1, the left panel displays the finite element solutions u^M_h for M = 16 computed over graded meshes with parameters \mu = 2, h = 10^{-3} and s\in \{0, 1, 0.25, 0.4\} . The right panel shows the value of u^M_h(x_1) as a function of M for s \in \{0.1, 0.25, 0.4\} , which is expected to behave according to (7.1).

    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.

    Table 1.  Comparison between computational results for the problem described in Section 7.1 over two different meshes for s = 0.4 . Let M_i be the value of M in the i -th row. In this table, by the slope at M_i we refer to \frac{\log(u_h^{M_i}(x_1)) -\log(u_h^{M_{i-1}}(x_1))}{\log(M_i) - \log(M_{i-1})} that, according to (7.1), is expected to be equal to 9/14 \approx 0.643 .
    Example with h = 10^{-3} Example with h = 10^{-4}
    \log_2(M) u_h^M(x_1) Slope u_h^M(x_1) Slope
    6.0 26.1545 \texttt{N/A} 26.7488 \texttt{N/A}
    6.5 32.4687 0.624 33.4057 0.641
    7.0 40.0845 0.608 41.5497 0.629
    7.5 49.1873 0.590 51.4627 0.617
    8.0 59.9410 0.571 63.4528 0.604

     | Show Table
    DownLoad: CSV

    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]

    \kappa(K) = \mathcal{O}\left( N^{2s/d} \left( \frac{h_{max}}{h_{min}} \right)^{d-2s} \right).

    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.

    Figure 2.  Condition numbers of the Jacobian matrix \widetilde{K} at the last step of Algorithm 2 for the problem described in Section 7.1 with s = 0.1 (left), s = 0.4 (right) and meshes with grading parameters \mu \in \{ 1, 2, 3\} . The condition number on quasi-uniform meshes ( \mu = 1 ) scales as N^{2s+1} , in agreement with the s -fractional mean curvature operator being an operator of order s+1/2 (cf. (2.3)). While the conditioning for graded meshes is significantly poorer, when \mu = 2 diagonal preconditioning recovers condition numbers comparable to the ones on quasi-uniform meshes.

    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} .

    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.

    Figure 3.  Effects of truncation in 2d for s = 0.25 : for g^H = \chi_{\Omega_H} , we compute the L^1 and L^\infty discrepancies between u_h\equiv 1 and u^H_h as a function of H . For both norms we observe a discrepancy of order H^{-1-2s} , in agreement with Corollary 5.2.

    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 .

    Figure 4.  Change of convexity: one-dimensional experiment for s = 0.02 with a = 1 (left panel) and a = 1.0001, 1.01, 1.05, \; 1.1, 1.2, 1.5 (right panel). The solutions u_a exhibit a transition from being convex in \Omega for a = 1 to being concave for a = 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.

    Figure 5.  Change of convexity: one-dimensional experiment with s = 0.02 (left panel) and two-dimensional experiment with s = 0.25 (right panel). The piecewise constant boundary data vanish near the boundary of \Omega and at infinity and are equal to 1 on an intermediate annulus.

    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

    g(x, y) = \gamma \left( \chi_{(-1,-a) \times (0,1)}(x,y) - \chi_{(-1,-a) \times (-1,0)}(x,y) \right)

    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 .

    Figure 6.  Plot of u_h in Section 7.5 for \gamma = 2, a = 1/8 and s = 0.25 . Left panel: top view of the solution. Right panel: slices at x = 2^{-3}, 2^{-6} and 2^{-9} . The fractional minimal graph flattens as x \to 0^+ , in agreement with the fact that for such a minimizer being continuous at some point x\in \partial\Omega implies having continuous tangential derivatives at such a point.

    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 .

    Table 2.  Example of Section 7.5: experimental slopes \partial_y u_h(x, 0) for x = 2^{-k} and k = 6, \ldots, 9 . As x\to 0^+ , these slopes become smaller; this geometric rigidity is easier to capture for larger s .
    s = 0.10
    h^{\mu} x = 2^{-9} x = 2^{-8} x = 2^{-7} x = 2^{-6}
    2^{-7} \texttt{N/A} \texttt{N/A} 8.546\times10^{-2} 1.1945\times10^{-1}
    2^{-8} \texttt{N/A} 5.856\times10^{-2} 8.406\times10^{-2} 1.2140\times10^{-1}
    2^{-9} 3.940\times10^{-2} 5.730\times10^{-2} 8.572\times10^{-2} 1.2332\times10^{-1}
     
    s = 0.25
    h^{\mu} x = 2^{-9} x = 2^{-8} x = 2^{-7} x = 2^{-6}
    2^{-7} \texttt{N/A} \texttt{N/A} 3.466\times10^{-2} 5.473\times10^{-2}
    2^{-8} \texttt{N/A} 2.135\times10^{-2} 3.469\times10^{-2} 5.551\times10^{-2}
    2^{-9} 1.289\times10^{-2} 2.126\times10^{-2} 3.543\times10^{-2} 5.640\times10^{-2}
     
    s = 0.40
    h^{\mu} x = 2^{-9} x = 2^{-8} x = 2^{-7} x = 2^{-6}
    2^{-7} \texttt{N/A} \texttt{N/A} 8.605\times10^{-3} 1.509\times10^{-2}
    2^{-8} \texttt{N/A} 4.763\times10^{-3} 8.613\times10^{-3} 1.540\times10^{-2}
    2^{-9} 2.578\times10^{-3} 4.739\times10^{-3} 8.886\times10^{-3} 1.574\times10^{-2}

     | Show Table
    DownLoad: CSV

    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.

    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.

    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.

    Figure 7.  Compatibility of data: plots of u_h for s = 0.25, f = -1 in \Omega (left) and after 16 Newton iterations for f = -10 in \Omega (right). The right hand side f = -10 turns out to be incompatible for the prescribed nonlocal mean curvature problem in \Omega = B(0, 1) .

    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 .

    Figure 8.  Nonlocal minimal graph with prescribed discontinuous nonlocal mean curvature. Left: plot of u_h in [-1.5, 1.5] , right: plot of u_h near origin.

    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.

    Figure 9.  A graph with prescribed discontinuous nonlocal mean curvature in the square \Omega = (-1, 1)^2 . The left panel displays a top view, while the right panel shows a side view along the slice \{ y = 1/2\} . The solution to (6.1) is discontinuous inside \Omega .

    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 .

    Figure 10.  Top and side views of functions with prescribed fractional mean curvature f = -1 in \Omega that vanish in \Omega^c . Here, \Omega is either an annulus (top row) or a square (middle and bottom row). The plot in the top-right panel corresponds to a radial slice ( y = 0, \, \, 0.25 \le x \le 1 ) of the annulus, while the ones in the bottom-left and bottom-right show slices along the diagonal ( 0 \le y = x \le 1 ) and perpendicular to an edge of the square ( y = 0.5, 0 \le x \le 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.

    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 .

    Figure 11.  Stickiness on the L-shaped domain \Omega = (-1, 1)^2 \setminus (0, 1)\times (-1, 0) with prescribed fractional mean curvature f = -1 in \Omega and Dirichlet condition g = 0 in \Omega^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.

    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 :

    H_s[u](x) = {\text{P.V.}} \int_{{\mathbb{R}}^d} G_s\left(\frac{u(x)-u(y)}{|x-y|}\right) \frac{dy}{|x-y|^{d+2s}} = f(x),

    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

    M_s(x) : = \lim\limits_{\Omega^c \ni y \to x} u(y) - u(x).

    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

    \begin{equation*} \label{eq:stickiness-balance} \int_{\Omega \cap B_r(x)} G_s\left(\frac{u(x)-u(y)}{|x-y|}\right) \frac{dy}{|x-y|^{d+2s}} \approx \int_{\Omega^c \cap B_r(x)} G_s\left(\frac{u(y) - u(x)}{|x-y|}\right) \frac{dy}{|x-y|^{d+2s}}. \end{equation*}

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

    \int_{\Omega^c \cap B_r(x)} G_s\left(\frac{u(y) - u(x)}{|x-y|}\right) \frac{dy}{|x-y|^{d+2s}} \approx \int_{\Omega^c \cap B_r(x)} G_s\left(\frac{M_s}{|x-y|}\right) \frac{dy}{|x-y|^{d+2s}}.

    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

    R(x) : = \frac{|\Omega\cap B_r(x)|}{|\Omega^c\cap B_r(x)|}

    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.

    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.

    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.

    The authors declare no conflict of interest.



    [1] G. Acosta, F. M. Bersetche, J. P. Borthagaray, A short FE implementation for a 2d homogeneous Dirichlet problem of a fractional Laplacian, Comput. Math. Appl., 74 (2017), 784–816. doi: 10.1016/j.camwa.2017.05.026
    [2] G. Acosta, J. P. Borthagaray, A fractional Laplace equation: regularity of solutions and finite element approximations, SIAM J. Numer. Anal., 55 (2017), 472–495. doi: 10.1137/15M1033952
    [3] M. Ainsworth, W. McLean, T. Tran, The conditioning of boundary element equations on locally refined meshes and preconditioning by diagonal scaling, SIAM J. Numer. Anal., 36 (1999), 1901–1932. doi: 10.1137/S0036142997330809
    [4] I. Babuška, R. B. Kellogg, J. Pitkäranta, Direct and inverse error estimates for finite elements with mesh refinements, Numer. Math., 33 (1979), 447–471. doi: 10.1007/BF01399326
    [5] B. Barrios, A. Figalli, E. Valdinoci, Bootstrap regularity for integro-differential operators, and its application to nonlocal minimal surfaces, Ann. Sc. Norm. Super. Pisa Cl. Sci., 13 (2014), 609–639.
    [6] J. P. Borthagaray, P. Ciarlet Jr, On the convergence in {H}^1-norm for the fractional Laplacian, SIAM J. Numer. Anal., 57 (2019), 1723–1743. doi: 10.1137/18M1221436
    [7] J. P. Borthagaray, W. Li, R. H. Nochetto, Finite element discretizations for nonlocal minimal graphs: Convergence, Nonlinear Anal., 189 (2019), 111566. doi: 10.1016/j.na.2019.06.025
    [8] J. P. Borthagaray, W. Li, R. H. Nochetto, Linear and nonlinear fractional elliptic problems, In: 75 Years of Mathematics of Computation, Providence, RI: Amer. Math. Soc., 2020, 69–92.
    [9] J. P. Borthagaray, R. H. Nochetto, A. J. Salgado, Weighted Sobolev regularity and rate of approximation of the obstacle problem for the integral fractional Laplacian, Math. Mod. Meth. Appl. Sci., 29 (2019), 2679–2717. doi: 10.1142/S021820251950057X
    [10] J. P. Borthagaray, L. M. Del Pezzo, S. Martínez, Finite element approximation for the fractional eigenvalue problem, J. Sci. Comput., 77 (2018), 308–329. doi: 10.1007/s10915-018-0710-1
    [11] X. Cabré, M. Cozzi, A gradient estimate for nonlocal minimal graphs, Duke Math. J., 168 (2019), 775–848.
    [12] L. Caffarelli, J.-M. Roquejoffre, O. Savin, Nonlocal minimal surfaces, Commun. Pure Appl. Math., 63 (2010), 1111–1144.
    [13] A. Chernov, T. von Petersdorff, Ch. Schwab, Exponential convergence of hp quadrature for integral operators with Gevrey kernels, ESAIM Math. Mod. Num. Anal., 45 (2011), 387–422. doi: 10.1051/m2an/2010061
    [14] S. Dipierro, O. Savin, E. Valdinoci, Graph properties for nonlocal minimal surfaces, Calc. Var., 55 (2016), 86. doi: 10.1007/s00526-016-1020-9
    [15] S. Dipierro, O. Savin, E. Valdinoci, Boundary behavior of nonlocal minimal surfaces, J. Funct. Anal., 272 (2017), 1791–1851. doi: 10.1016/j.jfa.2016.11.016
    [16] S. Dipierro, O. Savin, E. Valdinoci, Boundary properties of fractional objects: flexibility of linear equations and rigidity of minimal graphs, J. Reine Angew. Math., 2020 (2020), 121–164. doi: 10.1515/crelle-2019-0045
    [17] S. Dipierro, O. Savin, E. Valdinoci, Nonlocal minimal graphs in the plane are generically sticky, Commun. Math. Phys., 376 (2020), 2005–2063. doi: 10.1007/s00220-020-03771-8
    [18] G. Dziuk, Numerical schemes for the mean curvature flow of graphs, In: IUTAM symposium on variations of domain and free-boundary problems in solid mechanics, Springer, 1999, 63–70.
    [19] A. Figalli, E. Valdinoci, Regularity and Bernstein-type results for nonlocal minimal surfaces, J. Reine Angew. Math., 2017 (2017), 263–273.
    [20] M. Giaquinta, On the Dirichlet problem for surfaces of prescribed mean curvature, Manuscripta Math., 12 (1974), 73–86. doi: 10.1007/BF01166235
    [21] P. Grisvard, Elliptic problems in nonsmooth domains, Boston, MA: Pitman (Advanced Publishing Program), 1985.
    [22] C. Imbert, Level set approach for fractional mean curvature flows, Interface. Free Bound., 11 (2009), 153–176.
    [23] C. T. Kelley, Iterative methods for optimization, SIAM, 1999.
    [24] L. Lombardini, Approximation of sets of finite fractional perimeter by smooth sets and comparison of local and global s -minimal surfaces, Interface. Free Bound., 20 (2018), 261–296. doi: 10.4171/IFB/402
    [25] L. Lombardini, Minimization problems involving nonlocal functionals: nonlocal minimal surfaces and a free boundary problem, PhD thesis, Universita degli Studi di Milano and Universite de Picardie Jules Verne, 2018.
    [26] B. Merriman, J. K. Bence, S. Osher, Diffusion generated motion by mean curvature, AMS Selected Lectures in Mathematics Series: Computational Crystal Growers Workshop, 1992.
    [27] S. A. Sauter, C. Schwab, Boundary element methods, Berlin: Springer-Verlag, 2011.
    [28] O. Savin, E. Valdinoci, \Gamma-convergence for nonlocal phase transitions, Ann. Inst. H. Poincaré Anal. Non Linéaire, 29 (2012), 479–500. doi: 10.1016/j.anihpc.2012.01.006
  • This article has been cited by:

    1. Serena Dipierro, Ovidiu Savin, Enrico Valdinoci, Boundary continuity of nonlocal minimal surfaces in domains with singularities and a problem posed by Borthagaray, Li, and Nochetto, 2023, 62, 0944-2669, 10.1007/s00526-023-02606-3
    2. Juan Pablo Borthagaray, Wenbo Li, Ricardo H. Nochetto, 2023, Chapter 2, 978-3-031-34088-8, 27, 10.1007/978-3-031-34089-5_2
  • Reader Comments
  • © 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2262) PDF downloads(108) Cited by(2)

Figures and Tables

Figures(11)  /  Tables(2)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog