1.1. A brief review of a nonlocal vector calculus

To make some sense of how the model (1.5) arises and why we refer to it as a nonlocal analogue of the model (1.1), we need to first introduce some elements of a nonlocal vector calculus. Note that in this article we also use elements of the fractional calculus. However, because there are many excellent references about the fractional calculus and its applications (see e.g. Baleanu, Diethelm, Scalas and Trujillo Reference Baleanu, Diethelm, Scalas and Trujillo2016, Mainardi Reference Mainardi1997, Meerschaert and Sikorskii Reference Meerschaert and Sikorskii2012), we do not discuss them here.

The classical vector calculus provides a set of tools that are in ubiquitous use for the modelling, analysis, discretization and numerical analysis of PDE models, an obvious example being the use of Green’s first identity to transform the strong formulation of a PDE into a weak formulation. The foundations of that calculus are the familiar divergence, gradient and curl differential operators, upon which an edifice is built that includes vector identities (e.g. div curl ${\boldsymbol {v}}$ = 0), integral theorems (e.g. Gauss’s theorem) and much, much more.

A nonlocal vector calculus has been developed (Gunzburger and Lehoucq Reference Gunzburger and Lehoucq2010, Du et al. Reference Du, Gunzburger, Lehoucq and Zhou2012a , Du, Gunzburger, Lehoucq and Zhou Reference Du, Gunzburger, Lehoucq and Zhou2013a , Alali, Gunzburger and Liu Reference Alali, Gunzburger and Liu2015, Mengesha and Du Reference Mengesha and Du2016, Du Reference Du2019, D’Elia et al. Reference D’Elia, Flores, Li, Radu and Yu2019a ) to deal with nonlocal models such as (1.5) in much the same way as the classical vector calculus is used to deal with PDE models such as (1.1). Here, mostly following Du et al. (Reference Du, Gunzburger, Lehoucq and Zhou2013a ), we provide a brief introduction to the nonlocal vector calculus, including notions that are used in the rest of the article. We remark that certain elements of the nonlocal vector calculus have previously appeared in Coulhon and Grigoryan (Reference Coulhon and Grigoryan1998), Zhou and Schölkopf (Reference Zhou and Schölkopf2005), Gilboa and Osher (Reference Gilboa and Osher2008), Lézoray, Ta and Elmoataz (Reference Lézoray, Ta and Elmoataz2010), Lou, Zhang, Osher and Bertozzi (Reference Lou, Zhang, Osher and Bertozzi2010) and Jiang, Lim, Yao and Ye (Reference Jiang, Lim, Yao and Ye2011), for example. However, the discussions in those papers, compared to that in Du et al. (Reference Du, Gunzburger, Lehoucq and Zhou2013a ), are limited in scope and in application and provide only a partial development of a nonlocal vector calculus that mimics the classical vector calculus.

The foundation of the nonlocal vector calculus is integral operators that mimic the three differential operators upon which the classical vector calculus is built.

Given the vector-valued functions ${{\boldsymbol \nu }_{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})} \colon {({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})\times ({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}\to {{\mathbb R}^{d}}$ and ${\boldsymbol \alpha }_{\delta }{({\boldsymbol {x}},{\boldsymbol {y}})}\colon {({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})\times ({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}\to {{\mathbb R}^{d}}$ , the action of the nonlocal divergence operator ${{\mathcal D}_{\delta }}\colon {\Omega \cup \Omega _{{\mathcal I}_{\delta }}} \to {\mathbb R}$ on ${\boldsymbol \nu }{({\boldsymbol {x}},{\boldsymbol {y}})}$ is defined as

(1.13) $$ \begin{align} ({{\mathcal D}_{\delta}}{\boldsymbol\nu})({\boldsymbol{x}}) := \int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} ({\boldsymbol\nu}({\boldsymbol{y}},{\boldsymbol{x}})+{\boldsymbol\nu}{({\boldsymbol{x}},{\boldsymbol{y}})})\cdot{\boldsymbol\alpha}_{\delta}{({\boldsymbol{x}},{\boldsymbol{y}})} {{\mathcal X}_{{B_{\delta}({\boldsymbol{x}})}}({\boldsymbol{y}})}{{\,{\mathrm{d}}} {\boldsymbol{y}}} \end{align} $$

for all ${\boldsymbol {x}}\in {\Omega \cup \Omega _{{\mathcal I}_{\delta }}}$ , where ${\boldsymbol \alpha }_{\delta }{({\boldsymbol {x}},{\boldsymbol {y}})}$ denotes an antisymmetric function, that is, ${\boldsymbol \alpha }_{\delta }({\boldsymbol {y}},{\boldsymbol {x}})=-{\boldsymbol \alpha }_{\delta }{({\boldsymbol {x}},{\boldsymbol {y}})}$ for all ${\boldsymbol {x}},{\boldsymbol {y}}\in {\Omega \cup \Omega _{{\mathcal I}_{\delta }}}$ , and ${\mathcal X}_{({\,\cdot \,})}({\,\cdot \,})$ denotes the indicator function. Note that unlike its differential counterpart $\nabla \cdot $ , the nonlocal operator ${{\mathcal D}_{\delta }}$ is not uniquely defined, that is, not only do we have an unspecified parameter $\delta>0$ but also we have said nothing about ${\boldsymbol \alpha }_{\delta }{({\boldsymbol {x}},{\boldsymbol {y}})}$ other than it is an antisymmetric function. The choices one makes for $\delta $ and especially ${\boldsymbol \alpha }_{\delta }{({\boldsymbol {x}},{\boldsymbol {y}})}$ can result in operators having very different properties. Thus $\delta $ and ${\boldsymbol \alpha }_{\delta }{({\boldsymbol {x}},{\boldsymbol {y}})}$ are modelling choices dictated by the specific application one considers. We have already mentioned, in (1.4), the wide choices of $\delta $ that are of interest; several choices for ${\boldsymbol \alpha }_{\delta }{({\boldsymbol {x}},{\boldsymbol {y}})}$ are considered in Section 2.2.1.

Simple manipulations show that, under suitable regularity assumptions that we do not dwell on here,

(1.14) $$ \begin{align} \int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} v {{\mathcal D}_{\delta}} {\boldsymbol\nu} {\,{\mathrm{d}}} {\boldsymbol{x}} = \int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}}\int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} {\boldsymbol\nu} {\,\cdot\,} {{\mathcal D}_{\delta}^{\ast}} v {{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}}, \end{align} $$

where

(1.15) $$ \begin{align} ({{\mathcal D}_{\delta}^{\ast}} u){({\boldsymbol{x}},{\boldsymbol{y}})} := - (u({\boldsymbol{y}}) -u({\boldsymbol{x}})) {\boldsymbol\alpha}_{\delta}{({\boldsymbol{x}},{\boldsymbol{y}})} {{\mathcal X}_{{B_{\delta}({\boldsymbol{x}})}}({\boldsymbol{y}})} \end{align} $$

for all ${{\boldsymbol {x}},{\boldsymbol {y}}}\in {\Omega \cup \Omega _{{\mathcal I}_{\delta }}}$ . From (1.14), it is natural to refer to operator ${{\mathcal D}_{\delta }^{\ast }}$ as being the adjoint operator corresponding to ${{\mathcal D}_{\delta }}$ . With ${{\mathcal D}_{\delta }^{\ast }}$ being the adjoint of ${{\mathcal D}_{\delta }}$ , one may formally refer to $-{{\mathcal D}_{\delta }^{\ast }}$ as being a nonlocal gradient operator.

Let ${{\boldsymbol \Theta }_{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}\colon ({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})\times \ ({\Omega \cup \Omega _{{\mathcal I}_{\delta }}}) \to {\mathbb R}^{{d}\times {d}}$ denote a tensor-valued function that is symmetric in the function sense, i.e. ${{\boldsymbol \Theta }_{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}={{\boldsymbol \Theta }_{\delta }} ({\boldsymbol {y}}, {\boldsymbol {x}})$ , and symmetric and non-negative definite in the matrix sense. Then, from the definitions of the nonlocal divergence and gradient operators, we have that

(1.16) $$ \begin{align} &-{{\mathcal D}_{\delta}} ({{\boldsymbol\Theta}_{\delta}}{{\mathcal D}_{\delta}^{\ast}} u)\\ &\quad = 2\int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} (u({\boldsymbol{y}})-u({\boldsymbol{x}})) {\boldsymbol\alpha}_{\delta}{({\boldsymbol{x}},{\boldsymbol{y}})}{\,\cdot\,} ( {{\boldsymbol\Theta}_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})}{\boldsymbol\alpha}_{\delta}{({\boldsymbol{x}},{\boldsymbol{y}})} ){{\mathcal X}_{{B_{\delta}({\boldsymbol{x}})}}({\boldsymbol{y}})} {{\,{\mathrm{d}}} {\boldsymbol{y}}}. \notag \end{align} $$

For all ${{\boldsymbol {x}},{\boldsymbol {y}}}\in {\Omega \cup \Omega _{{\mathcal I}_{\delta }}}$ , we define the kernel ${\gamma _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ as

(1.17) $$ \begin{align} {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} := {\boldsymbol\alpha}_{\delta}{({\boldsymbol{x}},{\boldsymbol{y}})}{\,\cdot\,} ({{\boldsymbol\Theta}_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})}{\boldsymbol\alpha}_{\delta}{({\boldsymbol{x}},{\boldsymbol{y}})}){{\mathcal X}_{{B_{\delta}({\boldsymbol{x}})}}({\boldsymbol{y}})}. \end{align} $$

Note that ${\gamma _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ defined in this way satisfies the symmetry condition (1.7) and the support condition (1.8), where the former follows from the fact that ${\mathcal X}_{{B_{\delta }({\boldsymbol {x}})}}({\boldsymbol {y}})$ is itself a symmetric function, that is, if ${\boldsymbol {y}}\in {B_{\delta }({\boldsymbol {x}})}$ then necessarily ${\boldsymbol {x}}\in B_{\delta }({\boldsymbol {y}})$ .

We then have that (1.16) can be expressed as

(1.18) $$ \begin{align} {{\mathcal D}_{\delta}} ({{\boldsymbol\Theta}_{\delta}}{{\mathcal D}_{\delta}^{\ast}} u)&=-2\int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} (u({\boldsymbol{y}})-u({\boldsymbol{x}})) {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} {{\,{\mathrm{d}}} {\boldsymbol{y}}}\nonumber\\ & = -2\int_{{B_{\delta}({\boldsymbol{x}})}} (u({\boldsymbol{y}})-u({\boldsymbol{x}})) {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} {{\,{\mathrm{d}}} {\boldsymbol{y}}}, \end{align} $$

where ${\gamma _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ is given by (1.17) and ${{\mathcal D}_{\delta }}$ and ${{\mathcal D}_{\delta }^{\ast }}$ are given by (1.13) and (1.15), respectively.

Of course, we recognize that

(1.19) $$ \begin{align} {{\mathcal D}_{\delta}} ({{\boldsymbol\Theta}_{\delta}}{{\mathcal D}_{\delta}^{\ast}} u)=-{{\mathcal L}_{\delta}} u, \end{align} $$

where ${{\mathcal L}_{\delta }} u$ is defined in (1.6) or (1.11). The fact that the operator ${{\mathcal L}_{\delta }}$ can be written as a composition of nonlocal divergence and gradient operators justifies referring to (1.5) as a (variable coefficient) nonlocal Poisson problem, that is, ${{\mathcal L}_{\delta }} u$ may be viewed as a nonlocal analogue of $\nabla {\,\cdot \,} ( {{\boldsymbol {D}}} \nabla u)$ for suitable choices of ${\boldsymbol \alpha }_{\delta }{({\boldsymbol {x}},{\boldsymbol {y}})}$ , and, if ${{\boldsymbol \Theta }_{\delta }}$ is the identity tensor, ${{\mathcal L}_{\delta }} u$ is indeed a nonlocal analogue of the $\Delta u$ .

The final bit of the nonlocal vector calculus that we need is a (generalized) nonlocal Green’s first identity. In the context of this section, that identity is given by (Du et al. Reference Du, Gunzburger, Lehoucq and Zhou2013a )

(1.20) $$ \begin{align} \int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} v{{\mathcal D}_{\delta}} ({{\boldsymbol\Theta}_{\delta}}{{\mathcal D}^{\ast}} u) {{\,{\mathrm{d}}} {\boldsymbol{y}}} = \int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}}\int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} {{\mathcal D}_{\delta}^{\ast}} v {\,\cdot\,} ({{\boldsymbol\Theta}_{\delta}}{{\mathcal D}_{\delta}^{\ast}})u {{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}}. \end{align} $$

We have

(1.21) $$ \begin{align} {{\mathcal D}_{\delta}^{\ast}} v {\,\cdot\,} ({{\boldsymbol\Theta}_{\delta}}{{\mathcal D}_{\delta}^{\ast}})u = (v({\boldsymbol{y}})-v({\boldsymbol{x}}))(u({\boldsymbol{y}})-u({\boldsymbol{x}})) {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})}, \end{align} $$

so substituting (1.18), (1.10) and (1.21) into (1.20) results in the equivalent form of the nonlocal Green’s first identity given by

(1.22) $$ \begin{align} & \int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}}\int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} (v({\boldsymbol{y}})-v({\boldsymbol{x}}))(u({\boldsymbol{y}})-u({\boldsymbol{x}})) {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})}{{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}} \nonumber\\ &\quad =-\int_{{\Omega}} v({\boldsymbol{x}}) {{\mathcal L}_{\delta}} u ({\boldsymbol{x}}){{\,{\mathrm{d}}} {\boldsymbol{x}}} + \int_{{\Omega_{{\mathcal I}_{\delta}}}} v({\boldsymbol{x}}) {{\mathcal N}_{\delta}} u ({\boldsymbol{x}}){{\,{\mathrm{d}}} {\boldsymbol{x}}}. \end{align} $$

We recognize (1.22) as a nonlocal analogue of the (generalized) local Green’s first identity

$$ \begin{align*} \int_{\Omega} \nabla v{\,\cdot\,}({{\boldsymbol{D}}}\nabla u) {{\,{\mathrm{d}}} {\boldsymbol{x}}} = -\int_{\Omega} v\nabla{\,\cdot\,}({{\boldsymbol{D}}}\nabla u) {{\,{\mathrm{d}}} {\boldsymbol{x}}} + \int_{\partial{\Omega}} v {\boldsymbol{n}} {\,\cdot\,}({{\boldsymbol{D}}}\nabla u) {{\,{\mathrm{d}}} {\boldsymbol{x}}}. \end{align*} $$

Derivations in a more rigorous functional analytic setting of nonlocal integral identities similar to (1.22) and the connection to their local analogues can be found in Mengesha and Du (Reference Mengesha and Du2016) and Du (Reference Du2019).

Thus, in addition to ${{\mathcal D}_{\delta }}$ and ${{\mathcal D}_{\delta }^{\ast }}$ , we have the nonlocal gradient operator ${\mathcal G}_{\delta }$ defined by, for a scalar-valued function $\nu {({\boldsymbol {x}},{\boldsymbol {y}})}$ ,

(1.23) $$ \begin{align} ({\mathcal G}_{\delta}\nu)({\boldsymbol{x}}) := \int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} (\nu({\boldsymbol{y}},{\boldsymbol{x}})+\nu{({\boldsymbol{x}},{\boldsymbol{y}})}){\boldsymbol\alpha}_{\delta}{({\boldsymbol{x}},{\boldsymbol{y}})} {{\mathcal X}_{{B_{\delta}({\boldsymbol{x}})}}({\boldsymbol{y}})} {{\,{\mathrm{d}}} {\boldsymbol{y}}} \end{align} $$

for all ${\boldsymbol {x}}\in {\Omega \cup \Omega _{{\mathcal I}_{\delta }}}$ and its adjoint operator ${\mathcal G}^{\ast }_{\delta }$ given by, for a vector-valued function ${\boldsymbol {u}}({\boldsymbol {x}})$ ,

(1.24) $$ \begin{align} ({\mathcal G}^{\ast}_{\delta} {\boldsymbol{u}}){({\boldsymbol{x}},{\boldsymbol{y}})} = - ({\boldsymbol{u}}({\boldsymbol{y}})-{\boldsymbol{u}}({\boldsymbol{x}})) {\,\cdot\,} {\boldsymbol\alpha}_{\delta}{({\boldsymbol{x}},{\boldsymbol{y}})}{{\mathcal X}_{{B_{\delta}({\boldsymbol{x}})}}({\boldsymbol{y}})} \end{align} $$

for all ${{\boldsymbol {x}},{\boldsymbol {y}}}\in {\Omega \cup \Omega _{{\mathcal I}_{\delta }}}$ . Similarly, we have the nonlocal curl operator ${\mathcal C}_{\delta }$ defined by, for a vector-valued function ${\boldsymbol \nu }{({\boldsymbol {x}},{\boldsymbol {y}})}$ ,

(1.25) $$ \begin{align} ({\mathcal C}_{\delta}{\boldsymbol\nu})({\boldsymbol{x}}) := \int_{{{\mathbb R}^{d}}} ({\boldsymbol\nu}({\boldsymbol{y}},{\boldsymbol{x}})+{\boldsymbol\nu}{({\boldsymbol{x}},{\boldsymbol{y}})})\times{\boldsymbol\alpha}_{\delta}{({\boldsymbol{x}},{\boldsymbol{y}})}{{\mathcal X}_{{B_{\delta}({\boldsymbol{x}})}}({\boldsymbol{y}})} {{\,{\mathrm{d}}} {\boldsymbol{y}}} \end{align} $$

for all ${\boldsymbol {x}}\in {\Omega \cup \Omega _{{\mathcal I}_{\delta }}}$ and its adjoint operator ${\mathcal C}^{\ast }_{\delta }$ defined by, for a vector-valued function ${\boldsymbol {u}}({\boldsymbol {x}})$ ,

(1.26) $$ \begin{align} {({\mathcal C}^{\ast}_{\delta}{\boldsymbol{u}}){({\boldsymbol{x}},{\boldsymbol{y}})} = - {\boldsymbol\alpha}_{\delta}{({\boldsymbol{x}},{\boldsymbol{y}})}\times({{\boldsymbol{u}}}({\boldsymbol{y}})-{{\boldsymbol{u}}}({\boldsymbol{x}})){{\mathcal X}_{{B_{\delta}({\boldsymbol{x}})}}({\boldsymbol{y}})} } \end{align} $$

for all ${{\boldsymbol {x}},{\boldsymbol {y}}}\in {\Omega \cup \Omega _{{\mathcal I}_{\delta }}}$ . Thus, in the nonlocal vector calculus we have pairs of divergence operators ${{\mathcal D}_{\delta }}$ and $-{\mathcal G}^{\ast }_{\delta }$ , gradient operators ${\mathcal G}_{\delta }$ and $-{\mathcal D}^{\ast }_{\delta }$ , and curl operators ${\mathcal C}_{\delta }$ and ${\mathcal C}^{\ast }_{\delta }$ , with ${{\mathcal D}_{\delta }}$ , ${\mathcal G}_{\delta }$ , ${\mathcal C}_{\delta }$ acting on functions of two points, i.e. ${\boldsymbol \nu }{({\boldsymbol {x}},{\boldsymbol {y}})}$ and $\nu {({\boldsymbol {x}},{\boldsymbol {y}})}$ , and ${{\mathcal D}_{\delta }^{\ast }}$ , ${\mathcal G}^{\ast }_{\delta }$ , ${\mathcal C}^{\ast }_{\delta }$ acting on functions of one point, i.e. ${\boldsymbol {u}}({\boldsymbol {x}})$ and $u({\boldsymbol {x}})$ .

With these operators in hand, we have the nonlocal vector identities

(1.27) $$ \begin{align} \begin{cases} {{\mathcal D}_{\delta}}({\mathcal C}^{\ast}_{\delta} {\boldsymbol{u}}) = 0, \qquad {\mathcal C}_{\delta}({\mathcal G}^{\ast}_{\delta} u) = 0,\\ {{\mathcal D}_{\delta}}({{\mathcal D}_{\delta}^{\ast}} {\boldsymbol{u}}) = {\mathcal G}_{\delta}({\mathcal G}^{\ast}_{\delta} {\boldsymbol{u}}) + {\mathcal C}_{\delta}({\mathcal C}^{\ast}_{\delta} {\boldsymbol{u}}), \end{cases} \end{align} $$

which are nonlocal analogues of $\nabla {\,\cdot \,}(\nabla \times {\boldsymbol {u}})=0$ , $\nabla \times (\nabla u)=0$ and $\Delta {\boldsymbol {u}} = \nabla {\,\cdot \,} (\nabla {\boldsymbol {u}}) = \nabla ( \nabla {\,\cdot \,} {\boldsymbol {u}}) - \nabla \times (\nabla \times {\boldsymbol {u}})$ , respectively. In particular, we can view $-{{\mathcal D}_{\delta }}({{\mathcal D}_{\delta }^{\ast }} {\boldsymbol {u}}) = -{\mathcal G}_{\delta }({\mathcal G}^{\ast }_{\delta } {\boldsymbol {u}}) - {\mathcal C}_{\delta } ({\mathcal C}^{\ast }_{\delta } {\boldsymbol {u}})$ as a nonlocal vector Laplacian.

In addition, the operators introduced in this remark are used in D’Elia et al. (Reference D’Elia, Flores, Li, Radu and Yu2019a ) to define the nonlocal Hodge–Helmholtz decomposition of vector-valued functions that depend on two points ${\boldsymbol {x}}$ and ${\boldsymbol {y}}$ . Note that the operators ${{\mathcal D}_{\delta }}$ , ${{\mathcal D}_{\delta }^{\ast }}$ , ${\mathcal G}_{\delta }$ and ${\mathcal G}^{\ast }_{\delta }$ , as operators between functions of two points ${\boldsymbol {x}}$ and ${\boldsymbol {y}}$ and functions of one point ${\boldsymbol {x}}$ or ${\boldsymbol {y}}$ , do not recover their classical counterparts in the local limit $\delta \to 0$ . Thus, in parallel, additional examples of nonlocal operators acting on vector-valued functions of one point only are introduced in Du et al. (Reference Du, Gunzburger, Lehoucq and Zhou2013a ) and further analysed in Mengesha and Du (Reference Mengesha and Du2016), Du (Reference Du2019), Du, Tao, Tian and Yang (Reference Du, Tao, Tian and Yang2019c ) and Lee and Du (Reference Lee and Du2020). In fact, a nonlocal Helmholtz decomposition of a vector-valued function ${\boldsymbol {u}}({\boldsymbol {x}})$ can be established which recovers the classical one in the local limit; see Lee and Du (Reference Lee and Du2020). These operators can also be used to construct a nonlocal analogue of the Stokes equation, which was used to provide a foundation for the analysis of smoothed particle hydrodynamics (Du and Tian Reference Du and Tian2019, Lee and Du Reference Lee and Du2020). Some of these developments can be found in Section 20.

2.1. Function spaces, weak formulations and well-posedness

2.1.1. Function spaces and norms

We define an ‘energy’ space

(2.1) $$ \begin{align} {V({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})} = \{ v({\boldsymbol{x}})\in L^2({\Omega\cup\Omega_{{\mathcal I}_{\delta}}}) \colon \| v \|_{{V({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})}} <\infty\}, \end{align} $$

where $\| \cdot \|_{{V({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}} $ denotes the ‘energy’ norm defined as

(2.2) $$ \begin{align} & \| v \|_{{V({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})}}\\ &\quad = \biggl(\int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}}\int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} (v({\boldsymbol{y}}) -v({\boldsymbol{x}})) ^2 {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} {{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}} + \|v\|^2_{L^2({\Omega\cup\Omega_{{\mathcal I}_{\delta}}}) }\biggr)^{1/2}. \nonumber \end{align} $$

The space ${V({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ is by definition a subspace of $L^2({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})$ and can be shown, for suitably chosen kernels, to be a Hilbert space that comes equipped with the inner product

(2.3) $$ \begin{align} &(u,v)_{{V({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})}} \nonumber\\ &\quad = \int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}}\int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} (v({\boldsymbol{y}})-v({\boldsymbol{x}})) (u({\boldsymbol{y}})-u({\boldsymbol{x}})) {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} {{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}} \nonumber\\ &\quad\quad\, + (u,v)_{L^2({\Omega\cup\Omega_{{\mathcal I}_{\delta}}}) }\quad\text{for all } u,v\in{V({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})}. \end{align} $$

We also have the constrained energy space

(2.4) $$ \begin{align} {V_c({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})} = \{ v({\boldsymbol{x}})\in {V({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})} \colon v({\boldsymbol{x}}) = 0 \ \text{for } {\boldsymbol{x}}\in{\Omega_{{\mathcal I}_{\delta}}} \}. \end{align} $$

Under suitable conditions on the kernels, we can prove a nonlocal Poincaré inequality (Du et al. Reference Du, Gunzburger, Lehoucq and Zhou2012a , Mengesha and Du Reference Mengesha and Du2013, Du Reference Du2019). Thus, we have that the seminorm

(2.5) $$ \begin{align} | v |_{{V_c({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})}} = \biggl( \int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}}\int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} (v({\boldsymbol{y}})-v({\boldsymbol{x}})) ^2 {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} {{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}} \biggr)^{1/2} \end{align} $$

is a norm on ${V_c({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ equivalent to the norm (2.2) with equivalence constants that are independent of $\delta $ as $\delta \to 0$ . Correspondingly, ${V_c({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ is a Hilbert space equipped with the inner product

(2.6) $$ \begin{align} (u,v)_{{V_c({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})}} = \int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}}\int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} (v({\boldsymbol{y}})-v({\boldsymbol{x}})) (u({\boldsymbol{y}})-u({\boldsymbol{x}})) {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} {{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}} \end{align} $$

for all $u,v\in {V_c({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ . For more detailed discussions of the conditions on the kernels and rigorous proofs of the properties on these spaces, we refer to Mengesha and Du (Reference Mengesha and Du2015, Reference Mengesha and Du2016) and Du (Reference Du2019).

To allow inhomogeneous data, we define the space ${V_t({\Omega _{{\mathcal I}_{\delta }}})}$ ,

(2.7) $$ \begin{align} {{V_t({\Omega_{{\mathcal I}_{\delta}}})} = \{ v = w|_{{\Omega_{{\mathcal I}_{\delta}}}} \; \text{for some}\; w\in{V({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})} \},} \end{align} $$

which involves restrictions to the domain ${\Omega _{{\mathcal I}_{\delta }}}$ having finite volume in ${{\mathbb R}^{d}}$ . A norm on ${V_t({\Omega _{{\mathcal I}_{\delta }}})}$ can be defined by

(2.8) $$ \begin{align} \| v \|_{{V_t({\Omega_{{\mathcal I}_{\delta}}})}} = \inf\{ \|w\|_{{V({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})}}\mid w\in{V({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})}, w|_{{\Omega_{{\mathcal I}_{\delta}}}}=v\}. \end{align} $$

One can view ${V_t({\Omega _{{\mathcal I}_{\delta }}})}$ , induced by ${V({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ , as a nonlocal analogue of the trace space $H^{1/2}(\partial {\Omega })$ induced by $H^1({\Omega })$ .

As mentioned in Remark 1.2, specific function spaces are defined by making specific choices for ${\boldsymbol \alpha }_{\delta }{({\boldsymbol {x}},{\boldsymbol {y}})}$ . For some choices, the function spaces so defined can be shown to be equivalent to standard function spaces such as Sobolev spaces. See Section 2.2.1 for some examples.

2.1.2. Weak formulations and well-posedness

A weak formulation of problem (1.5) with ${\mathcal V} u=u$ is defined as follows:

Given ${\gamma _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ defined in (1.17) and given $g({\boldsymbol {x}})\in {V_t({\Omega _{{\mathcal I}_{\delta }}})}$ and $f({\boldsymbol {x}})\in {V_{d}({\Omega })}$ , find $u({\boldsymbol {x}})\in {V({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ such that

(2.9) $$ \begin{align} \begin{cases} {{{\mathcal A}_{\delta}}}(u,v) = {{\langle{f},{v}\rangle}} &\text{for all } v\in{V_c({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})},\\ u({\boldsymbol{x}}) = g({\boldsymbol{x}}) &\text{for all } {\boldsymbol{x}}\in{\Omega_{{\mathcal I}_{\delta}}}, \end{cases} \end{align} $$

where, for all $u({\boldsymbol {x}}),v({\boldsymbol {x}})\in {V({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ , we have the bilinear form

(2.10) $$ \begin{align} {{{\mathcal A}_{\delta}}}(u,v) =\int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}}\int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} (v({\boldsymbol{y}})-v({\boldsymbol{x}})) (u({\boldsymbol{y}})-u({\boldsymbol{x}})) {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} {{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}} \end{align} $$

and, for all $f({\boldsymbol {x}})\in {V_{d}({\Omega })}$ , the dual space of ${V_c({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ defined by extending the $L^2$ inner product as the duality pairing and with the induced norm $\|\cdot \|_{{V_{d}({\Omega })}}$ , we have the linear functional

(2.11) $$ \begin{align} {{\langle{f},{v}\rangle}}=\int_{{\Omega}} f({\boldsymbol{x}})v({\boldsymbol{x}}) {{\,{\mathrm{d}}} {\boldsymbol{x}}}. \end{align} $$

That (2.9)–(2.11) are indeed a weak formulation corresponding to (1.5) is an immediate consequence of setting $v({\boldsymbol {x}})=0$ on ${\Omega _{{\mathcal I}_{\delta }}}$ in the nonlocal Green’s first identity (1.22) so that the term involving ${{\mathcal N}_{\delta }} u$ vanishes. Equation (2.9) is a nonlocal analogue of the weak formulation of problem (1.1) with ${\mathcal B} u=u$ given by

$$ \begin{align*} {\int_{\Omega} \nabla v{\,\cdot\,}({{\boldsymbol{D}}} \nabla u){{\,{\mathrm{d}}} {\boldsymbol{x}}} = {{\langle{f},{v}\rangle}},} \end{align*} $$

for all $v\in H^1_0({\Omega })$ along with $u=g$ on $\partial {\Omega }$ .

Under suitable conditions on the kernel ${\gamma _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ , it can be shown that the bilinear form ${{{\mathcal A}_{\delta }}}(u,v)$ is symmetric and is continuous and coercive with respect to ${V_c({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ so that, by the Lax–Milgram theorem, problem (2.9) is well-posed. In particular, we have the a priori estimate with respect to the energy norm for the solution $u({\boldsymbol {x}})$ of (2.9) given by

(2.12) $$ \begin{align} \| u \|_{V({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})} \le {C_{\delta}} ( \|f\|_{V_{d}({\Omega})} + \| g\|_{V_t({\Omega_{{\mathcal I}_{\delta}}})}) , \end{align} $$

where the constant ${C_{\delta }}$ does not depend on $\delta $ . Moreover, for suitably chosen kernels, one can consider its local ( $\delta \to 0$ ) and fractional ( $\delta \to \infty $ ) limits, as discussed in Part 2.

Remark 2.1 (a nonlocal Neumann model)

We can of course also define a weak formulation of (1.5) for ${\mathcal V} u = {\mathcal N}_{{\delta }} u$ , that is, for a nonlocal Neumann problem. In this case a weak formulation is given by: find $u\in {V({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}\setminus {\mathbb R} $ such that

(2.13) $$ \begin{align} \begin{cases} \displaystyle {{{\mathcal A}_{\delta}}}(u,v) = {{\langle{f},{v}\rangle}} + \int_{{\Omega_{{\mathcal I}_{\delta}}}} g({\boldsymbol{x}})v({\boldsymbol{x}}){{\,{\mathrm{d}}} {\boldsymbol{x}}} {\quad\text{{for all} }} v\in{V({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})}\setminus{\mathbb R}, \\ \displaystyle \int_{{\Omega}} f({\boldsymbol{x}}) {{\,{\mathrm{d}}} {\boldsymbol{x}}} + \int_{{\Omega_{{\mathcal I}_{\delta}}}} g({\boldsymbol{x}}){{\,{\mathrm{d}}} {\boldsymbol{x}}} =0, \end{cases} \end{align} $$

where $f\in V_d(\Omega )$ and for the Neumann problem $g({\boldsymbol {x}})$ belongs to the dual space of ${V({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}|_{\Omega _{{\mathcal I}_{\delta }}}$ . In (2.13), the quotient space ${V({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}\setminus {\mathbb R}$ is used to ensure uniqueness of solutions. The second equation in (2.13) is a compatibility condition on the data f and g that is necessary to ensure existence.

Equation (2.13) is a nonlocal analogue of the weak formulation corresponding to (1.1) with ${\mathcal B} u=({{\boldsymbol {D}}} \nabla u)\cdot {\boldsymbol {n}}$ given by

$$ \begin{align*} \int_{\Omega} \nabla v{\,\cdot\,}({{\boldsymbol{D}}} \nabla u){{\,{\mathrm{d}}} {\boldsymbol{x}}} = \int_{\Omega} fv{{\,{\mathrm{d}}} {\boldsymbol{x}}} + \int_{\partial{\Omega}}gv{{\,{\mathrm{d}}} {\boldsymbol{x}}} \end{align*} $$

for all $v\in H^1({\Omega })\setminus {\mathbb R}$ and a pair of prescribed data, also denoted by f and g, that satisfy the compatibility condition

$$ \begin{align*} {\int_{\Omega} f{{\,{\mathrm{d}}} {\boldsymbol{x}}} + \int_{\partial{\Omega}}g{{\,{\mathrm{d}}} {\boldsymbol{x}}}=0.} \end{align*} $$

We note that other ways of defining a nonlocal Neumann problem can also be used. For additional discussions about nonlocal Neumann problems, see e.g. Cortazar, Elgueta, Rossi and Wolanski (Reference Cortazar, Elgueta, Rossi and Wolanski2008), Du (Reference Du2019), Mengesha and Du (Reference Mengesha and Du2015, Reference Mengesha and Du2016), Tao, Tian and Du (Reference Tian and Du2017) and D’Elia, Tian and Yu (Reference D’Elia, Tian and Yu2019d ). One can also find a fractional version of (1.5) with ${\mathcal V} u = {\mathcal N}_{{\delta }} u$ in Deng, Li, Tian and Zhan (Reference Deng, Li, Tian and Zhang2018) for a nonlocal Neumann problem of the fractional Laplacian in the bounded domain. In addition, for specialized power-like kernels, one can also relate this to problems associated with problems related to the regional Laplacian defined by (3.6).□

For the sake of economy of exposition, we will focus on (2.9) in which ${\mathcal V} u = u$ and, at times, we will even focus on the homogeneous version of (2.9) for which $g({\boldsymbol {x}})=0$ .

2.1.3. Relation to an energy minimization principle

The weak formulation (2.9) can also be derived as the Euler–Lagrange equation for an energy minimization principle. Consider the energy functional

(2.14) $$ \begin{align} \begin{aligned} {\mathcal E}_{\delta}(u;f,{\gamma_{\delta}})& =\displaystyle \frac12\int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}}\int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} (u({\boldsymbol{y}})-u({\boldsymbol{x}})) ^2 {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} {{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}} \\ &\quad\, - \int_{{\Omega}} f({\boldsymbol{x}})u({\boldsymbol{x}}) {{\,{\mathrm{d}}} {\boldsymbol{x}}} \end{aligned} \end{align} $$

and the following minimization principle:

(2.15) $$ \begin{align} &\text{Given } {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} \text{ defined in } (1.17) \text{ and given } g({\boldsymbol{x}})\in{V_t({\Omega_{{\mathcal I}_{\delta}}})}\nonumber \\ &\text{and } f({\boldsymbol{x}})\in{V_{d}({\Omega})}, \text{ find } u({\boldsymbol{x}})\in{V({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})} \text{ such that} \nonumber\\[4pt] &\qquad\qquad\qquad\qquad{\mathcal E}_{\delta}(u;f,{\gamma_{\delta}})\ \text{is minimized}\nonumber \\[4pt] &\text{subject to } u({\boldsymbol{x}}) = g({\boldsymbol{x}}) \text{ for } {\boldsymbol{x}}\in{\Omega_{{\mathcal I}_{\delta}}}. \end{align} $$

Using standard techniques from the calculus of variations, one easily sees that (2.9) is indeed the Euler–Lagrange equation corresponding to the energy minimization principle (2.15) under suitable conditions on the kernel ${\gamma _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ and data f and g.

2.2. Kernel choices and the corresponding energy spaces

We assume that ${\gamma _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ given by (1.17) can be written in the form

(2.16) $$ \begin{align} {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})}= {\phi_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} {\theta_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} {\mathcal X}_{{B_{\delta}({\boldsymbol{x}})}}({\boldsymbol{y}}) , \end{align} $$

where ${\theta _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ and ${\phi _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ denote non-negative, symmetric, scalar-valued functions. This form is very general with respect to both operators and constitutive functions. In fact, by setting

(2.17) $$ \begin{align} {\boldsymbol\alpha}_{\delta}{({\boldsymbol{x}},{\boldsymbol{y}})} = \displaystyle \frac{{{\boldsymbol{y}}-{\boldsymbol{x}}}}{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|} \sqrt{{\phi_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})}} \end{align} $$

and

(2.18) $$ \begin{align} {\theta_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} :=\displaystyle \frac{{{\boldsymbol{y}}-{\boldsymbol{x}}}}{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|} {\,\cdot\,} \biggl( {{\boldsymbol\Theta}_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} \displaystyle \frac{{{\boldsymbol{y}}-{\boldsymbol{x}}}}{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|}\biggr) \end{align} $$

in (1.17), we arrive at (2.16). Thus the only assumption made is that ${\boldsymbol \alpha }_{\delta }{({\boldsymbol {x}},{\boldsymbol {y}})}$ is directed along the vector ${{\boldsymbol {y}}-{\boldsymbol {x}}}$ . Note that ${\boldsymbol \alpha }_{\delta }{({\boldsymbol {x}},{\boldsymbol {y}})}$ given by (2.17) is indeed an antisymmetric function, whereas ${\theta _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ given by (2.18) is indeed a non-negative, symmetric function if ${{\boldsymbol \Theta }_{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ is a symmetric in the function sense and symmetric and non-negative definite in the matrix sense.

We do not combine ${\theta _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ and ${\phi _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ into a single function because we want to keep separate the operator definition and constitutive roles played by ${\gamma _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ ; note that ${\phi _{\delta }}$ is only related to ${\boldsymbol \alpha }_{\delta }$ whereas ${\theta _{\delta }}$ is only related to ${{\boldsymbol \Theta }_{\delta }}$ . Thus ${\phi _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ now takes over the role of ${\boldsymbol \alpha }_{\delta }{({\boldsymbol {x}},{\boldsymbol {y}})}$ in this regard, and ${\theta _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ can be viewed as a constitutive function. Note also that ${\boldsymbol \alpha }_{\delta }{({\boldsymbol {x}},{\boldsymbol {y}})}$ as given by (2.17) can be used to define the ‘first-order’ nonlocal operators ${{\mathcal D}_{\delta }}$ and ${{\mathcal D}_{\delta }^{\ast }}$ , whereas ${\gamma _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ is used to define the ‘second-order’ nonlocal operator ${{\mathcal L}_{\delta }}$ and the ‘Neumann’ nonlocal operator ${{\mathcal N}_{\delta }}$ . We will refer to ${\gamma _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ as the kernel, ${\phi _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ as the kernel function and ${\theta _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ as a constitutive function.

Remark 2.2 (scalar constitutive tensors)

The case of ${{\boldsymbol \Theta }_{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ being a scalar tensor (i.e. if ${{\boldsymbol \Theta }_{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})} = {\theta _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})} {{\boldsymbol {I}}}$ , where ${{\boldsymbol {I}}}$ denotes the identity tensor and ${\theta _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ denotes a scalar-valued function) is by far the most common case considered, just as it is for a scalar diffusivity tensor ${{\boldsymbol {D}}} = a({\boldsymbol {x}}){{\boldsymbol {I}}}$ in the local PDE case. In this case, (2.18) becomes a tautology. Also, in this case, (1.17) simplifies to

$$ \begin{align*} {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})}&= {\theta_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} {\boldsymbol\alpha}_{\delta}{({\boldsymbol{x}},{\boldsymbol{y}})}\cdot{\boldsymbol\alpha}_{\delta}{({\boldsymbol{x}},{\boldsymbol{y}})}{\mathcal X}_{{B_{\delta}({\boldsymbol{x}})}}({\boldsymbol{y}}) \\ &= {\theta_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} |{\boldsymbol\alpha}_{\delta}{({\boldsymbol{x}},{\boldsymbol{y}})}|^2{\mathcal X}_{{B_{\delta}({\boldsymbol{x}})}}({\boldsymbol{y}}). \end{align*} $$

2.2.1. A list of kernel functions in common use

Because ${\theta _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ is a constitutive function which is not specific even within a single application, we focus on choices for the function ${\phi _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ that determines generic properties of a model. However, we again note that because we have assumed that ${{\boldsymbol \Theta }_{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ is a symmetric, non-negative definite tensor, it follows that ${\theta _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ is a non-negative function for all ${\boldsymbol {x}},{\boldsymbol {y}}\in {\Omega \cup \Omega _{{\mathcal I}_{\delta }}}$ .

We also assume, as is usually the case, that for all ${\boldsymbol {x}},{\boldsymbol {y}}\in {\Omega \cup \Omega _{{\mathcal I}_{\delta }}}$ , the constitutive function ${\theta _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ is bounded from above and below by positive constants whose values do not depend on $\delta $ . In fact, often ${\theta _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ does not depend on $\delta $ . Thus we focus on a list of kernel functions ${\phi _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ . In Table 2.1 we provide the energy spaces corresponding to the kernel functions we list.

Table 2.1

Kernel functions ${\phi _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ and the corresponding energy spaces ${V({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ .

Translation-invariant, integrable kernel functions

The kernel function is translation-invariant, that is, ${\phi _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}={\phi _{\delta }}({{\boldsymbol {y}}-{\boldsymbol {x}}})$ , and satisfies for some positive constant $C>0$ ,

(2.19) $$ \begin{align} {C\leq \int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} {\phi_{\delta}} ({{\boldsymbol{y}}-{\boldsymbol{x}}}) {{\,{\mathrm{d}}} {\boldsymbol{y}}} < \infty \quad\text{for all } {\boldsymbol{x}}\in{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}.} \end{align} $$

In the design of numerical methods, one may need to differentiate between cases for which

• • ${\phi _{\delta }} {({\boldsymbol {z}})}$ is a smooth, bounded function of ${{\boldsymbol {z}}}$ ,

• • ${\phi _{\delta }} {({\boldsymbol {z}})}$ is a bounded function of ${{\boldsymbol {z}}}$ but is not smooth, e.g. a piecewise constant function,

• • ${\phi _{\delta }} {({\boldsymbol {z}})}$ is a singular but integrable function, e.g. ${\phi _{\delta }}({\boldsymbol {z}})\propto 1/|{\boldsymbol {z}}|^r$ with $r<{d}$ .

A special case of translation-invariant kernel functions is that of radial kernel functions for which ${\phi _{\delta }}({\boldsymbol {z}})={\phi _{\delta }}(|{\boldsymbol {z}}|)$ .

‘Critical’ kernel functions

The kernel function satisfies

(2.20) $$ \begin{align} {\phi_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} \propto \displaystyle \frac1{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|^{d}}. \end{align} $$

We refer to this kernel as critical because it ‘just misses being integrable’, that is, $1/{|{{\boldsymbol {y}}-{\boldsymbol {x}}}|^{{d}-{\epsilon }}}$ is integrable for any ${\epsilon }>0$ and for any ${d}$ , but is not integrable for any ${d}$ if ${\epsilon }=0$ . For the other kernels in this list, the corresponding energy space is precisely known; see Table 2.1. For critical kernels, that is not the case; what is known is that energy space is a Hilbert space and is a strict subspace of $L^2({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})$ .

‘Peridynamic’ kernel functions

The kernel function satisfies

(2.21) $$ \begin{align} {\phi_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} \propto \displaystyle \frac1{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|}. \end{align} $$

We refer to this kernel as peridynamic because it has the same singular behaviour as one of the commonly used kernels in peridynamic models for solid mechanics. Note that for ${d}=1$ this kernel function is a member of the critical class, whereas for ${d}=2$ it is integrable and for ${d}=3$ it is square-integrable.

Fractional kernel functions

The kernel function satisfies

(2.22) $$ \begin{align} {\phi_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} \propto \displaystyle \frac{1}{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|^{{d}+2s}}. \end{align} $$

We refer to this kernel as fractional because the singularity at ${\boldsymbol {y}}={\boldsymbol {x}}$ is that of the kernel for the fractional Laplace operator. This kernel plays a central role in Section 3.

Remark 2.4. We note that many square-integrable kernel functions of the more general form ${\phi _{\delta }}({\boldsymbol {x}},{\boldsymbol {y}})$ also have ${V({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}=L^2({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})$ as the energy space. However, we do not consider this class of kernels in the later sections.□

Remark 2.5 (solution smoothness in nonlocal diffusion)

For translation-invariant integrable kernel functions such as (2.19), the solution may not be smoother than the data. For example, if $g({\boldsymbol {x}})=0$ and $f({\boldsymbol {x}})\in L^2({\Omega })$ , then generically, the solution $u({\boldsymbol {x}})=({-}{{\mathcal L}_{\delta }})^{-1}f\in L^2({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})$ . Because $L^2({\Omega })$ includes, for example, functions with jump discontinuities possibly over some co-dimension $1$ interfaces, this means that if the data $f({\boldsymbol {x}})$ has a jump discontinuity, so may the solution $u({\boldsymbol {x}})$ . For the fractional kernel function, we could have $2s$ derivatives of smoothing, that is, if the data $f({\boldsymbol {x}})\in H^{-s}({\Omega })$ , then $u({\boldsymbol {x}})\in H^{s}({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})$ , where $H^{s}({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})$ denotes a fractional Sobolev space and $H^{-s}({\Omega })$ its dual space. Note that for $s<1/2$ , $H^{s}({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})$ also contains functions having jump discontinuities.□

Remark 2.6 (scaling constants and limit behaviour)

Often kernel functions take the form

$$ \begin{align*} {\phi_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} = C {\varphi}{({\boldsymbol{x}},{\boldsymbol{y}})}, \end{align*} $$

where C denotes a scaling constant that, depending on the specific instance, may depend on $\delta $ and also on other parameters appearing in the definition of ${\phi _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ . One approach towards defining the scaling constant C is to choose it so that in the limit $\delta \to 0$ , a nonlocal model and its solution reduce to a local PDE model and its solution, respectively. Another approach is to choose C so that in the limit $\delta \to \infty $ , the nonlocal model and its solution reduce to a model that is posed on ${{\mathbb R}^{d}}$ . Below, we will have occasion to follow both approaches. In any case, the scaling constant is a modelling choice so that, for example, if neither type of limiting behaviour is of interest, then the scaling constant should be chosen based on the physics being modelled.□

3.1. Integral fractional Laplacian models for diffusion on bounded domains

Letting ${\Omega }\subset {\mathbb R}^{d}$ denote a bounded Lipschitz domain, we define the integral fractional Laplacian $-(-\Delta )^{s}$ to be the restriction of the full-space operator to functions satisfying a volume constraint on $\Omega _{{\mathcal I}_{\infty }}:= {\mathbb R}^{d}\setminus \Omega $ , i.e. on the interaction domain corresponding to ${\Omega }$ . The fractional Poisson problem is then given by

(3.5) $$ \begin{align} \begin{cases} (-\Delta)^{s}u = f &\text{for all }{\boldsymbol{x}}\in{\Omega}, \\ u=g&\text{for all }{\boldsymbol{x}}\in \Omega_{{\mathcal I}_{\infty}} , \end{cases} \end{align} $$

where we have a given source term $f({\boldsymbol {x}})$ and Dirichlet volume constraint data $g({\boldsymbol {x}})$ defined for ${\boldsymbol {x}}\in {\Omega }$ and ${\boldsymbol {x}}\in {\Omega _{{\mathcal I}_{\infty }}}$ , respectively. Problem (3.5) is a nonlocal analogue of the local Poisson problem for the PDE Laplacian $\Delta $ . It is known (see e.g. Biccari and Hernández-Santamaría Reference Biccari and Hernández-Santamaría2018) that as $s\rightarrow 1^-$ , the solution of integral fractional diffusion model (3.5) strongly converges to the solution of the local diffusion problem in $H^{1-\epsilon }(\Omega )$ .

A related operator on ${\Omega }$ is the regional fractional Laplacian (Bogdan, Burdzy and Chen Reference Bogdan, Burdzy and Chen2003, Chen and Kim Reference Chen and Kim2002)

(3.6) $$ \begin{align} - (-\Delta) _{\text{regional}}^{s} u ({{\boldsymbol{x}}}) = {C_{{d},s}}~ \int_{{\Omega}} \displaystyle \frac{u( {{\boldsymbol{y}}})-u( {{\boldsymbol{x}}})}{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|^{{d}+2s}}{{\,{\mathrm{d}}} {\boldsymbol{y}}}, \end{align} $$

which is used in one of several approaches for generalizing the PDE Poisson problem with a homogeneous Neumann boundary condition to the fractional Laplacian case (Dipierro, Ros-Oton and Valdinoci Reference Dipierro, Ros-Oton and Valdinoci2017). Note that in (3.5) the operator has not changed, that is, the operator defined in (3.1) is used; what has changed is that domain of the operator is changed from ${{\mathbb R}^{d}}$ to the bounded domain ${\Omega }$ . On the other hand, in using the operator in (3.6), not only is the domain changed in the same manner, but the operator itself has changed, that is, (3.1) involves an integral over ${{\mathbb R}^{d}}$ whereas (3.6) involves an integral over ${\Omega }$ .

To discuss the variational form of the fractional Poisson problem (3.5), we use the standard fractional Sobolev space $H^{s} ({{\mathbb R}^{d}}) $ defined via the Fourier transform as

$$ \begin{align*} H^s({{\mathbb R}^{d}}) = \biggl\{ u\in L^2({{\mathbb R}^{d}}) \colon \int_{{{\mathbb R}^{d}}} (1+ |{\boldsymbol\xi}|^{2s}) |{\mathcal F} u({\boldsymbol\xi})|^2 {\,{\mathrm{d}}} {\boldsymbol\xi} <\infty \biggr\}. \end{align*} $$

If ${\Omega }\subset {{\mathbb R}^{d}}$ is a bounded domain, we define the Sobolev space $H^{s}({\Omega })$ as (McLean Reference McLean2000)

$$ \begin{align*} H^{s} ({\Omega}):=\{u\in L^{2} ({\Omega}) \colon {\|u\|}_{H^{s} ({\Omega})} < \infty\} \end{align*} $$

equipped with the norm

$$ \begin{align*} {\|u\|}_{H^{s} ({\Omega})}^{2}= |u|^{2}_{H^{s} ({\Omega})} + {\|u\|}_{L^{2} ({\Omega})}^{2}, \end{align*} $$

where we have the seminorm

$$ \begin{align*} |u|^{2}_{H^{s} ({\Omega})}=\int_{{\Omega}} \int_{{\Omega}} \displaystyle \frac{ (u( {{\boldsymbol{y}}})-u( {{\boldsymbol{x}}})) ^{2}}{ |{{\boldsymbol{y}}-{\boldsymbol{x}}}|^{{d}+2s}} {{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}}. \end{align*} $$

Moreover, when imposing a homogeneous Dirichlet volume constraints, e.g. $g=0$ in (3.5), we use the space

$$ \begin{align*} {H}^{s}_c({\Omega}):=\{u\in H^{s}({{\mathbb R}^{d}}) \colon u=0 \ \text{for all } {\boldsymbol{x}}\in \Omega_{{\mathcal I}_{\infty}}\} \end{align*} $$

that is equipped with the norm

$$ \begin{align*} {\|u\|}_{{H}^{s}_c ({\Omega})}^{2} = {\|u\|}_{H^{s} ({{\mathbb R}^{d}}) }^{2}={\|u\|}_{L^{2} ({\Omega})}^{2} + |u|^{2}_{H^{s} ({{\mathbb R}^{d}}) }. \end{align*} $$

For $s>1/2$ , ${H}^{s}_c({\Omega })$ coincides with the space $H_{0}^{s}({\Omega })$ that is the closure of $C_{0}^{\infty }({\Omega })$ with respect to the $H^{s}({\Omega })$ -norm, whereas for $s<1/2$ , ${H}^{s}_c({\Omega })$ is identical to $H^{s}({\Omega })$ . In the critical case $s=1/2$ , . See e.g. McLean (Reference McLean2000, Chapter 3) for a detailed discussion.

A variational form of (3.5) is derived starting from the integration-by-parts formula (e.g. the nonlocal Green’s first identity (1.22))

$$ \begin{align*} &\displaystyle \frac12\int_{\mathbb{R}^d}\int_{\mathbb{R}^d}(u({{\boldsymbol{y}}})-u({{\boldsymbol{x}}})) (v({{\boldsymbol{y}}})-v({{\boldsymbol{x}}})) {{\gamma}_{s}}{({\boldsymbol{x}},{\boldsymbol{y}})} {{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}} \end{align*} $$

$$ \begin{align*} &\quad = \int_{{\Omega}} ((-\Delta)^{s}u({{\boldsymbol{x}}})) v({{\boldsymbol{x}}}) {{\,{\mathrm{d}}} {\boldsymbol{x}}} + \int_{{\Omega_{{\mathcal I}_{\infty}}}} ({\mathcal N}_{s}u({{\boldsymbol{x}}})) v({{\boldsymbol{x}}}) {{\,{\mathrm{d}}} {\boldsymbol{x}}} , \end{align*} $$

where, as in (1.12), we have the nonlocal operator given by

$$ \begin{align*} \mathcal{N}_{s}u({{\boldsymbol{x}}})= \int_{\mathbb{R}^d} (u({{\boldsymbol{x}}})-u({{\boldsymbol{y}}})) {{\gamma}_{s}}{({\boldsymbol{x}},{\boldsymbol{y}})} {{{\,{\mathrm{d}}} {\boldsymbol{y}}}} \quad\text{for all } {\boldsymbol{x}}\in \Omega_{{\mathcal I}_{\infty}}. \end{align*} $$

An alternative definition of $\mathcal {N}_s$ has been presented in Dipierro et al. (Reference Dipierro, Ros-Oton and Valdinoci2017), corresponding to a nonlocal operator $\mathcal {N}_{\delta }$ that is defined on $\Omega _{\delta }$ as an integral over the domain $\Omega $ only, instead of ${\Omega \cup \Omega _{{\mathcal I}_{\delta }}}$ as in the case of (1.10). Substituting the volume constraint in (3.5) and $v({\boldsymbol {x}})=0$ on ${\Omega _{{\mathcal I}_{\infty }}}$ , we then have the following weak formulation:

(3.7) $$ \begin{align} &\text{Find } u({\boldsymbol{x}})\in {H}^{s}(\Omega) \text{ such that } u({\boldsymbol{x}})=g({\boldsymbol{x}}) \text{ for all } {\boldsymbol{x}}\in \Omega_{{\mathcal I}_{\infty}}\text{ and }\nonumber\\ &{\mathcal A}_s(u,v)={{\langle{f},{v}\rangle}} + (g,v)_{\Omega,\Omega_{{\mathcal I}_{\infty}}} \quad\text{ for all } v\in {H}^{s}_c({\Omega}), \end{align} $$

where

(3.8) $$ \begin{align} \begin{aligned} {\mathcal A}_s(u,v) &=\underbrace{\displaystyle \frac{1}{2} \int_{\Omega} \int_{\Omega}\; (u({\boldsymbol{y}}) - u({\boldsymbol{x}})) (v({\boldsymbol{y}}) -v({\boldsymbol{x}})) {\gamma}_{s} ({\boldsymbol{x}},{\boldsymbol{y}}) {\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}}_{{{\mathcal A}_{\Omega,\Omega}} (u,v)}\\ &\quad\, + \underbrace{\int_{{\Omega}}u({{\boldsymbol{x}}})v({{\boldsymbol{x}}})\int_{\Omega_{{\mathcal I}_{\infty}}} {\gamma}_{s} ({\boldsymbol{x}},{\boldsymbol{y}}) {{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}}}_{{\mathcal A}_{{\Omega,\Omega_{{\mathcal I}_{\infty}}}}(u,v)} \end{aligned} \end{align} $$

and

(3.9) $$ \begin{align} (g,v)_{\Omega,\Omega_{{\mathcal I}_{\infty}}} = \int_{\Omega} v({\boldsymbol{x}}) \int_{\Omega_{{\mathcal I}_{\infty}}} g({\boldsymbol{y}}){{\gamma}_{s}}{({\boldsymbol{x}},{\boldsymbol{y}})} {{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}}. \end{align} $$

The bilinear form ${\mathcal A}_s({\,\cdot \,},{\,\cdot \,})$ is ${H}^{s}_c({\Omega })$ -coercive and continuous so that as long as the right-hand side of (3.7) is bounded, the well-posedness of that problem follows from the Lax–Milgram theorem. If $g=0$ , ${\mathcal A}_s({\,\cdot \,},{\,\cdot \,})$ is identical to the bilinear form ${\mathcal A}_{\delta }({\,\cdot \,},{\,\cdot \,})$ defined in (2.10) with $\delta =\infty $ . On the other hand, for $g\neq 0$ , ${\mathcal A}_s({\,\cdot \,},{\,\cdot \,})$ can be seen as the $\delta =\infty $ limit of the bilinear form corresponding to some inhomogeneously constrained nonlocal diffusion problems discussed in Du (Reference Du2019), instead of the bilinear form ${\mathcal A}_{\delta }({\,\cdot \,},{\,\cdot \,})$ given by (2.10). Note that the bilinear form ${\mathcal A}_{{{\Omega },{\Omega }}}({\,\cdot \,},{\,\cdot \,})$ can be seen to correspond to the regional fractional Laplacian.

Note that

$$ \begin{align*} {\mathcal A}_{{\Omega},{\Omega_{{\mathcal I}_{\infty}}}}(u,v) = {C_{{d},s}} \int_{{\Omega}} u({{\boldsymbol{x}}}) v({{\boldsymbol{x}}})\bigg( \int_{{\Omega_{{\mathcal I}_{\infty}}}} \displaystyle \frac{1}{{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|}^{{d}+2s}}{{\,{\mathrm{d}}} {\boldsymbol{y}}}\bigg) {{\,{\mathrm{d}}} {\boldsymbol{x}}}. \end{align*} $$

The identity

$$ \begin{align*} { \displaystyle \frac{1}{{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|}^{{d}+2s}} = - \displaystyle \frac{1}{2s} \nabla_{{{\boldsymbol{y}}}}\cdot \displaystyle \frac{{{\boldsymbol{y}}-{\boldsymbol{x}}}}{{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|}^{{d}+2s}}} \end{align*} $$

and Gauss’s theorem results in

$$ \begin{align*} {\mathcal A}_{{\Omega},{\Omega_{{\mathcal I}_{\infty}}}}(u,v) = \displaystyle \frac{{C_{{d},s}}}{2s} \int_{{\Omega}} u({{\boldsymbol{x}}}) v({{\boldsymbol{x}}})\bigg( \int_{\partial{\Omega}} \displaystyle \frac{({{\boldsymbol{y}}-{\boldsymbol{x}}})\cdot{\boldsymbol{n}}_{\boldsymbol{y}}}{{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|}^{{d}+2s}}{{\,{\mathrm{d}}} {\boldsymbol{y}}}\bigg) {{\,{\mathrm{d}}} {\boldsymbol{x}}}, \end{align*} $$

where ${\boldsymbol {n}}_{\boldsymbol {y}}$ denotes the unit outer normal to $\partial {\Omega }$ at ${{\boldsymbol {y}}}$ . We then have that

(3.10) $$ \begin{align} {\mathcal A}_s(u,v) = {\mathcal A}_{{{\Omega},{\Omega}}}(u,v) + \displaystyle \frac{{C_{{d},s}}}{2s} \int_{{\Omega}} u({{\boldsymbol{x}}}) v({{\boldsymbol{x}}})\bigg( \int_{\partial{\Omega}} \displaystyle \frac{({{\boldsymbol{y}}-{\boldsymbol{x}}})\cdot{\boldsymbol{n}}_{\boldsymbol{y}}}{{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|}^{{d}+2s}}{{\,{\mathrm{d}}} {\boldsymbol{y}}}\bigg) {{\,{\mathrm{d}}} {\boldsymbol{x}}} \end{align} $$

so that an integral over the unbounded domain ${\Omega _{{\mathcal I}_{\infty }}}$ can be avoided in any computation involving ${\mathcal A}_s(u,v)$ . Of course, for a homogeneous volume constraint $g({\boldsymbol {x}})=0$ on ${\Omega _{{\mathcal I}_{\infty }}}$ , the right-hand side in (3.7) also only involves integrals over ${\Omega }$ . For $g({\boldsymbol {x}})\ne 0$ , one does have to evaluate, for ${\boldsymbol {x}}\in {\Omega }$ , the data integral $\int _{\Omega _{{\mathcal I}_{\infty }}} g({\boldsymbol {y}}){{\gamma }_{s}}{({\boldsymbol {x}},{\boldsymbol {y}})} {\,{\mathrm {d}}} {\boldsymbol {y}}$ that, in principle, can be approximated using a quadrature rule.

Above we have incorporated the inhomogeneous volume constraint directly into the weak form. An alternative is to enforce that condition via Lagrange multipliers (Acosta, Borthagaray and Heuer Reference Acosta, Borthagaray and Heuer2019).

3.1.1. Truncated interaction domains

An alternative to the approach based on the bilinear form (3.10) is to truncate the integration domain ${\Omega _{{\mathcal I}_{\infty }}}$ as is considered in D’Elia and Gunzburger (Reference D’Elia and Gunzburger2013). For $\delta> {\text {diam}} ({\Omega })$ , we pose the truncated problem

(3.11) $$ \begin{align} \begin{cases} (-\Delta)^{s}_{\delta} u = f &\text{for all }{\boldsymbol{x}}\in{\Omega},\\ u=g&\text{for all }{\boldsymbol{x}}\in{{\Omega_{{\mathcal I}_{\delta}}}} , \end{cases} \end{align} $$

where

$$ \begin{align*} (-\Delta) ^{s}_{\delta} u ({{\boldsymbol{x}}}) := \int_{{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}} ( u( {{\boldsymbol{x}}})-u( {{\boldsymbol{y}}})) {{\gamma}_{s}}{{({\boldsymbol{x}},{\boldsymbol{y}})}} {{\,{\mathrm{d}}} {\boldsymbol{y}}}. \end{align*} $$

D’Elia and Gunzburger (Reference D’Elia and Gunzburger2013) and Burkovska and Gunzburger (Reference Burkovska and Gunzburger2019b ) have shown that if $u_{\infty }$ and $u_{\delta }$ denote the solutions of (3.5) and (3.11), respectively, then

(3.12) $$ \begin{align} \| u_{\delta} - u_{\infty}\|_{H^s({\Omega})} \le \displaystyle \frac{C}{\delta^{2s}} \quad\text{as } \delta\to\infty, \end{align} $$

where C depends on norms of f and g.

Remark 3.1. Clearly, (3.12) implies that the solution of (3.11) can be viewed as an approximation of the solution of (3.5). However, (3.11) and its solution are useful in their own right. In (3.5), the horizon $\delta $ is infinite, but, in practice, although $\delta $ may be large, it is not likely to be infinite. In (3.11), $\delta $ is assumed large but finite so that it may provide a better model for practical applications that feature large horizons. In this case one can reverse roles and interpret the solution of (3.5) to be an approximation of the solution of (3.11), with (3.12) now telling us something about the error incurred by replacing the perhaps more useful model (3.11) with the much more studied model (3.5).□

3.2. Spectral fractional Laplacian models for diffusion on bounded domains

An alternative definition of a fractional-order Laplacian on a bounded domain ${\Omega }$ makes use of spectral information about the PDE Laplacian. Here we focus on the case of homogeneous constraints.

Let $0<\lambda _{0}\leq \lambda _{1}\leq \cdots $ and ${\varphi }_{0},{\varphi }_{1},\ldots $ denote the eigenvalues and the corresponding eigenfunctions of the PDE Laplacian $-\Delta $ with a homogeneous Dirichlet boundary condition, that is, we have

(3.13) $$ \begin{align} \begin{cases} -\Delta{\varphi}_{m}({{\boldsymbol{x}}}) =\lambda_{m}{\varphi}_{m}({{\boldsymbol{x}}})& \text{for all }{{\boldsymbol{x}}}\in{\Omega} ,\\ {\varphi}_{m} ({{\boldsymbol{x}}}) =0 & \text{for all } {{\boldsymbol{x}}}\in\partial{\Omega}, \end{cases} \end{align} $$

where the orthogonal eigenfunctions are normalized so that ${\|{\varphi }_{m}\|}_{L^{2}({\Omega })}=1$ . Then the eigenfunctions $\{{\varphi }_{m}\}_{m=0}^{\infty }$ form a complete orthonormal basis for $L^{2}({\Omega })$ . As a result, any function $u\in L^{2}({\Omega })$ can be expanded as

$$ \begin{align*} u=\sum_{m=0}^{\infty}u_{m}{\varphi}_{m} \quad\text{with } u_{m}= ({u},{{\varphi}_{m}}), \end{align*} $$

where $(\cdot ,\cdot )$ denotes the inner product of $L^2({\Omega })$ . We then have

$$ \begin{align*} (-\Delta) u({{\boldsymbol{x}}})= \sum_{m=0}^{\infty}u_{m} \lambda_m {\varphi}_{m} ({{\boldsymbol{x}}}) , \end{align*} $$

and the spectral fractional Laplacian of order $s\in (0,1)$ with homogeneous boundary condition is given by

$$ \begin{align*} ({-\Delta}|_{{\Omega},0}) ^{s}u ({{\boldsymbol{x}}}) = \sum_{m=0}^{\infty} u_{m}\lambda_{m}^{s} {\varphi}_{m}({{\boldsymbol{x}}}). \end{align*} $$

As $s\rightarrow 0$ , the identity operator is recovered, whereas the integer order Laplacian is recovered as $s\rightarrow 1$ .

For a given source term f, the spectral fractional Laplacian Poisson problem for $u({\boldsymbol {x}})$ , ${\boldsymbol {x}}\in {\Omega }$ , is then given by

(3.14) $$ \begin{align} ({-\Delta}|_{{\Omega},0}) ^{s}u ({{\boldsymbol{x}}}) = f({\boldsymbol{x}}) \quad\text{for all }{\boldsymbol{x}}\in{\Omega}. \end{align} $$

Using the heat kernel $p_{{\Omega },0}({{\boldsymbol {x}}},{{\boldsymbol {y}}},t)$ , the spectral fractional Laplacian can be rewritten (see e.g. Abatangelo and Dupaigne Reference Abatangelo and Dupaigne2017) in the form of the nonlocal operator

$$ \begin{align*} ({-\Delta}|_{{\Omega},0}) ^{s}u({{\boldsymbol{x}}}) = \int_{{\Omega}} \gamma_{{\Omega},0;s}({{\boldsymbol{x}}},{{\boldsymbol{y}}}) (u({{\boldsymbol{y}}})-u({{\boldsymbol{x}}})) {{\,{\mathrm{d}}} {\boldsymbol{y}}} + \kappa_{{\Omega},0;s}({{\boldsymbol{x}}})u({{\boldsymbol{x}}}), \end{align*} $$

where

$$ \begin{align*} \gamma_{{\Omega},0;s}{({\boldsymbol{x}},{\boldsymbol{y}})} &= \displaystyle \frac{s}{\Gamma(1-s)} \int_{0}^{\infty} p_{{\Omega},0}(t,{{\boldsymbol{x}}},{{\boldsymbol{y}}}) \displaystyle \frac{1}{t^{1+s}}{{\,{\mathrm{d}}} t}, \end{align*} $$

$$ \begin{align*} \kappa_{{\Omega},0;s}({{\boldsymbol{x}}}) &= \displaystyle \frac{s}{\Gamma(1-s)} \int_{0}^{\infty} \biggl(1-\int_{{\Omega}}p_{{\Omega},0}(t,{{\boldsymbol{x}}},{{\boldsymbol{y}}}) {\,{\mathrm{d}}} {\boldsymbol{y}}\biggr) \displaystyle \frac{1}{t^{1+s}} {\,{\mathrm{d}}} t. \end{align*} $$

If ${\Omega }={{\mathbb R}^{d}}$ , the integral and the spectral definitions of the fractional Laplacian coincide, but they are different for bounded domains (Servadei and Valdinoci Reference Servadei and Valdinoci2014). For example, in the particular case of ${\Omega }={{\mathbb R}^{d}}_{+}$ , i.e. the half-plane, the kernel of the spectral fractional Laplacian is given by

$$ \begin{align*} \gamma_{{{\mathbb R}^{d}}_{+},0;s}{({\boldsymbol{x}},{\boldsymbol{y}})} = \gamma_{s}{({\boldsymbol{x}},{\boldsymbol{y}})}-\gamma_{s}({{\boldsymbol{x}}},-{{\boldsymbol{x}}}) \neq \gamma_{s}{({\boldsymbol{x}},{\boldsymbol{y}})} \quad \text{for all } {{\boldsymbol{x}}},{{\boldsymbol{y}}}\in{{\mathbb R}^{d}}_{+}. \end{align*} $$

Note that the spectral definition can allow for the treatment of non-homogeneous boundary conditions of various types; see e.g. Antil, Pfefferer and Rogovs (Reference Antil, Pfefferer and Rogovs2018) and Cusimano, del Teso, Gerardo-Giorda and Pagnini (Reference Cusimano, del Teso, Gerardo-Giorda and Pagnini2018).

3.3. Additional considerations about fractional Laplacian models

3.3.1. Other integral representations of the fractional Laplacian

The Dunford–Taylor integral is a powerful tool in the numerical analysis of the fractional diffusion problem (Bonito and Pasciak Reference Bonito and Pasciak2015, Bonito and Pasciak Reference Bonito and Pasciak2017, Bonito, Lei and Pasciak Reference Bonito, Lei and Pasciak2019). In particular, the solution u of the Poisson problem (3.14) involving the spectral fractional Laplacian can be given by the Dunford–Taylor representation

$$ \begin{align*} u= ({-\Delta}|_{{\Omega},0})^{{-s}} f = \displaystyle \frac{\sin (s\pi) }{\pi} \int_{0}^{\infty} {\mu^{-s} (\mu-\Delta) ^{-1}f {\,{\mathrm{d}}} \mu}. \end{align*} $$

A combination of using a sinc quadrature rule and a finite element discretization of the reaction–diffusion-type term in the integrand was explored by Bonito and Pasciak (Reference Bonito and Pasciak2015, Reference Bonito and Pasciak2017).

Using similar techniques, the bilinear form ${\mathcal A}_s(\cdot ,\cdot )$ associated with the integral fractional Laplacian can be rewritten as (Bonito et al. Reference Bonito, Lei and Pasciak2019)

(3.15) $$ \begin{align} {\mathcal A}_s(u,v) = \frac{2 \sin (\pi s)}{\pi} \int_0^{\infty} {\mu^{1-2s}} \int_{{{\mathbb R}^{d}}} ((-\Delta)(I-\mu^2 \Delta)^{-1} u )({\boldsymbol{x}})\, v({\boldsymbol{x}}) {{\,{\mathrm{d}}} {\boldsymbol{x}}} {\,{\mathrm{d}}}\mu , \end{align} $$

for any $u,v \in H_c^s(\Omega )$ .

3.3.2. Extension representations

An extension representation of the fractional Laplacian was introduced in Molčanov and Ostrovskiĭ (Reference Molčanov and Ostrovskiĭ1969) and popularized by Caffarelli and Silvestre (Reference Caffarelli and Silvestre2007), who showed that the fractional Poisson problem $(3.5)$ posed on ${{\mathbb R}^{d}}$ can be recast as a Neumann-to-Dirichlet mapping over the extended domain ${{\mathbb R}^{d}}\times [0,\infty )$ , that is, we have

(3.16) $$ \begin{align} \begin{cases} -\nabla{\,\cdot\,} z^{\beta} \nabla U ({{\boldsymbol{x}}},z ) = 0 & \text{for all }({{\boldsymbol{x}}},z)\in{{\mathbb R}^{d}}\times[0,\infty) , \\[7pt] \displaystyle \frac{\partial U}{\partial \nu^{\beta}} ({{\boldsymbol{x}}}) = 2^{1-2s}\displaystyle \frac{\Gamma(1-s)}{\Gamma(s)} f({{\boldsymbol{x}}})& \text{for all }{{\boldsymbol{x}}}\in {{\mathbb R}^{d}}, \end{cases} \end{align} $$

where $\beta = 1-2s$ and

$$ \begin{align*} \displaystyle \frac{\partial U}{\partial \nu^{\beta}}({{\boldsymbol{x}}})= -\lim_{{z}\rightarrow 0^{+}}z^{\beta}\displaystyle \frac{\partial U}{\partial z} ({{\boldsymbol{x}}},{z}) , \end{align*} $$

with the solution to (3.5) recovered by taking the trace of U on ${{\mathbb R}^{d}}$ .

The spectral fractional Laplacian can be recovered, by restricting the extension domain from ${{\mathbb R}^{d}}\times [0,\infty )$ to ${\Omega }\times [0,\infty )$ and imposing a homogeneous Dirichlet boundary condition on the lateral surface (Stinga and Torrea Reference Stinga and Torrea2010).

The apparent advantage of the extension problem is that the nonlocal problem is replaced with a classical, integer-order local problem. This comes at the price of having to deal with a singular weight function and an additional spatial dimension.

3.3.3. Regularity of solutions

In the classical integer-order case, smoother domains and smoother right-hand sides result in smoother solutions of the Poisson problem; see e.g. Taylor (Reference Taylor1996). A lifting property of this type does not hold for the fractional-order Poisson problems (3.5) and (3.14). The typical solution behaviour of (3.5) close to the boundary can be characterized (see Ros-Oton and Serra Reference Ros-Oton and Serra2014) as

(3.17) $$ \begin{align} u({{\boldsymbol{x}}})\approx \operatorname{dist}({\boldsymbol{x}},\partial{\Omega})^{s}. \end{align} $$

The Sobolev regularity of the solution of (3.5) was studied in Acosta, Bersetche and Borthagaray (Reference Acosta, Bersetche and Borthagaray2017) as a special case of a result in Grubb (Reference Grubb2015a ). Suppose that $\partial {\Omega }\in C^{\infty }$ and $f\in H^{r} ({\Omega })$ for $r\geq -s$ and let $u\in {H}^{s}_c ({\Omega })$ denote the solution of the fractional Poisson problem (3.7). Then we obtain the regularity estimate

$$ \begin{align*} u\in \begin{cases} H^{2s+r} ({\Omega}) & \text{if } 0<s+r<1/2 , \\ H^{s+1/2-\varepsilon} ({\Omega})\ \text{for all}\ \varepsilon>0 & \text{if } 1/2\leq s+r . \end{cases} \end{align*} $$

This result is extended to the case of non-homogeneous volume condition in Acosta et al. (Reference Acosta, Borthagaray and Heuer2019). Additional regularity results with respect to s and, for the truncated fractional Laplacian, with respect to $\delta $ , are derived in Burkovska and Gunzburger (Reference Burkovska and Gunzburger2019b ).

The solution of the spectral fractional Poisson problem (3.14) displays different properties. Its behaviour close to the boundary is given by (see Caffarelli and Stinga Reference Caffarelli and Stinga2016)

$$ \begin{align*} u({\boldsymbol{x}})\approx \operatorname{dist}({\boldsymbol{x}},\partial{\Omega})^{\min\{2s,1\}}, \quad s\neq 1/2. \end{align*} $$

Grubb (Reference Grubb2016) showed that if $f\in H_c^{r} ({\Omega })$ , $r\geq -s$ , then the weak solution of the spectral fractional Poisson problem (3.14) satisfies $u\in H_c^{r+2s} ({\Omega })$ .

More detailed regularity results for the spectral fractional Laplacian can be found in Grubb (Reference Grubb2016).

3.3.4. Analytic solutions in one and two dimensions

Closed-form solutions of the fractional Poisson problem (3.5) posed on the unit ball are available; see Dyda, Kuznetsov and Kwaśnicki (Reference Dyda, Kuznetsov and Kwaśnicki2017b ) for detailed derivations. These solutions are useful to have in hand, e.g. for verifying numerical results.

In ${d}=1$ dimension, for the source term

$$ \begin{align*} f_{k,0}^{1D}= 2^{2s}\Gamma (1+s) ^{2}\binom{s+k-1/2}{s}\binom{s+k}{s} P_{k}^{(s,-1/2)} (2x^{2}-1) \end{align*} $$

with $k\geq 0$ , the solution is given by

$$ \begin{align*} u_{k,0}^{1D}= P_{k}^{(s,-1/2)} (2x^{2}-1) (1-x^{2}) ^{s}_{+}, \end{align*} $$

where

$$ \begin{align*} \binom{a}{b}=\displaystyle \frac{\Gamma (a+1) }{\Gamma (b+1) \Gamma (a-b+1) } \end{align*} $$

denotes a generalized binomial coefficient, $P_{k}^{(a,b)}$ denote the Jacobi polynomials, and $a_{+}=\max \{0,a\}$ . Moreover, for

$$ \begin{align*} f_{k,1}^{1D} = 2^{2s}\Gamma (1+s) ^{2}\binom{s+k+1/2}{s}\binom{s+k}{s} x P_{k}^{(s,1/2)} (2x^{2}-1) \end{align*} $$

with $k\geq 0$ , the solution is given by

$$ \begin{align*} u_{k,1}^{1D}= x P_{k}^{(s,1/2)} (2x^{2}-1) (1-x^{2}) ^{s}_{+}. \end{align*} $$

Turning to ${d}=2$ dimensions, for

$$ \begin{align*} f_{k,\ell}^{2D} = 2^{2s}\Gamma (1+s) ^{2}\binom{s+k+\ell}{s}\binom{s+k}{s} r^{\ell}\cos (\ell{\varphi}) P_{k}^{(s,\ell)} (2r^{2}-1) \end{align*} $$

with $\ell ,k\geq 0$ and $(r,{\varphi })$ denoting polar coordinates, the solution is given by

$$ \begin{align*} u_{k,\ell}^{2D} = r^{\ell}\cos (\ell{\varphi}) P_{k}^{(s,\ell)} (2r^{2}-1) (1-r^{2}) ^{s}_{+}. \end{align*} $$

Two example solutions are shown in Figure 3.1. Observe that the solutions contain a term that behaves like $\delta ({\boldsymbol {x}}) ^{s}$ , where $\delta ({\boldsymbol {x}}) $ is the distance from ${\boldsymbol {x}}\in {\Omega }$ to the boundary.

Figure 3.1

Analytic solutions $u_{0,0}^{2D}$ for $s=0.4$ and $u_{1,1}^{2D}$ for $s=0.6$ . The behaviour (3.17) close to the boundary is apparent.

3.4. Inhomogeneous constitutive functions in fractional models

So far in this section we have only considered the fractional Laplacian and its variants, which is analogous to considering local diffusion problems having a diffusivity that is not only constant but equal to one. Of course, in the local case, it is an easy matter to treat inhomogeneous diffusivities, that is, we merely replace $\Delta u$ with $\nabla \cdot ({{\boldsymbol {D}}} \nabla u)$ , where ${{\boldsymbol {D}}}({\boldsymbol {x}})$ generally denotes an inhomogeneous constitutive tensor. Using the nonlocal vector calculus, it is almost as easy to define integral fractional models in which the constitutive properties of the media considered are inhomogeneous even if the nonlocal constitutive function is a tensor. Here we only consider generalizations of the integral fractional Laplacian (3.1); generalizations of the other integral operators discussed in this section follow immediately from that for (3.1).

The inhomogeneous fractional Laplacian is simply defined by choosing the kernel in (3.1) to now be

(3.18) $$ \begin{align} {{\gamma}_{s}}{({\boldsymbol{x}},{\boldsymbol{y}})}= {C_{{d},s}}\displaystyle \frac{{\theta_s}{({\boldsymbol{x}},{\boldsymbol{y}})}}{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|^{{d}+2s}} \quad\text{for all } {\boldsymbol{x}},{\boldsymbol{y}}\in{{\mathbb R}^{d}}, \end{align} $$

where

(3.19) $$ \begin{align} {\theta_s}{({\boldsymbol{x}},{\boldsymbol{y}})} :=\displaystyle \frac{{{\boldsymbol{y}}-{\boldsymbol{x}}}}{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|} {\,\cdot\,} \biggl( {{\boldsymbol\Theta}_{s}}{({\boldsymbol{x}},{\boldsymbol{y}})} \displaystyle \frac{{{\boldsymbol{y}}-{\boldsymbol{x}}}}{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|}\biggr) \end{align} $$

with ${{\boldsymbol \Theta }_{s}}{({\boldsymbol {x}},{\boldsymbol {y}})}$ denoting a constitutive tensor. Thus we have defined the inhomogeneous integral fractional Poisson operator

(3.20) $$ \begin{align} (-\Delta)^{s}_{{{\boldsymbol\Theta}_{s}}} u ({{\boldsymbol{x}}}) := \int_{{{\mathbb R}^{d}}} {C_{{d},s}}{\theta_s}{({\boldsymbol{x}},{\boldsymbol{y}})} \displaystyle \frac{u( {{\boldsymbol{x}}})-u( {{\boldsymbol{y}}})}{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|^{{d}+2s}} {{\,{\mathrm{d}}} {\boldsymbol{y}}} \quad\text{for all } {\boldsymbol{x}}\in{{\mathbb R}^{d}}. \end{align} $$

Note that we could have arrived at (3.20) by defining

$$ \begin{align*} (-\Delta)^{s}_{{{\boldsymbol\Theta}_{s}}} u ({{\boldsymbol{x}}}) := \frac12 {{\mathcal D}_{s}} ({{\boldsymbol\Theta}_{s}} {{\mathcal D}^{\ast}_{s}} u), \end{align*} $$

where ${{\mathcal D}_{s}}$ and ${{\mathcal D}^{\ast }_{s}}$ are defined by (1.13) and (1.15), respectively, with

(3.21) $$ \begin{align} {\boldsymbol\alpha}_{s}{({\boldsymbol{x}},{\boldsymbol{y}})} = \displaystyle \frac{{{\boldsymbol{y}}-{\boldsymbol{x}}}}{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|} \sqrt{\displaystyle \frac{{C_{{d},s}}}{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|^{{d}+2s}}}. \end{align} $$

If ${{\boldsymbol \Theta }_{s}}{({\boldsymbol {x}},{\boldsymbol {y}})}={{\boldsymbol \Theta }_{s}}(|{{\boldsymbol {y}}-{\boldsymbol {x}}}|)$ , i.e. ${{\boldsymbol \Theta }_{s}}$ is a radial function, then (3.1) is a model for a homogeneous medium, that is, in this case, (3.1) is a fractional analogue of the PDE $-\nabla \cdot ({{\boldsymbol {D}}}\nabla u)=f$ in which ${{\boldsymbol {D}}}$ is a constant tensor. If in addition ${{\boldsymbol \Theta }_{s}}$ is a scalar tensor, then (3.1) is a fractional analogue of the PDE $-\kappa \Delta u=f$ in which $\kappa $ is a constant. Thus, to obtain an inhomogeneous integral fractional model, ${{\boldsymbol \Theta }_{s}}{({\boldsymbol {x}},{\boldsymbol {y}})}$ cannot be a radial function. See Remark 2.3.

4.1. Asymptotically compatible schemes

It is known in practice that to obtain consistency between nonlocal models and corresponding local PDE models, the mesh or quadrature point spacing may have to be reduced at a faster pace than the reduction of the horizon parameter (Bobaru et al. Reference Bobaru, Yang, Alves, Silling, Askari and Xu2009, Bobaru and Duangpanya Reference Bobaru and Duangpanya2010, Chen and Gunzburger Reference Chen and Gunzburger2011). Otherwise there could potentially be complications, most notably inconsistent limiting solutions when the horizon parameter is coupled proportionally to the discretization parameter (Tian Reference Tian2017, Tian and Du Reference Tian and Du2013, Tian and Du Reference Tian and Du2014, Tian and Du Reference Tian and Du2020). Asymptotically compatible (AC) schemes, motivated by the findings in Tian and Du (Reference Tian and Du2013) and formally introduced in Tian and Du (Reference Du, Huang and Lehoucq2014), are numerical discretizations of nonlocal models that converge to nonlocal continuum models for a fixed horizon parameter and to the local discrete schemes as the horizon vanishes for both discrete schemes with a fixed numerical resolution and for continuum models with increasing numerical resolution.

Let $h>0$ denote a grid size parameter (or particle spacing) and let $\delta $ denote the horizon parameter or even a more generic model parameter. Instead of h, we could use the dimension N of a spectral basis. The implications of the AC property are illustrated in Figure 4.1. There, $u_{\delta }$ denotes the solution of the continuous nonlocal problem with $\delta>0$ , $u_0$ the solution of the corresponding continuous local problem, $u_{\delta ,h}$ the solution of the discretized nonlocal problem, and $u_{0,h}$ the solution of the discretized local problem. Figure 4.1 is meant to illustrate the basic property of AC schemes, namely that for such schemes it does not matter in which order the limits are taken, or even if $\delta $ and h are related in some way so that they simultaneously tend to zero. Note that figures similar to Figure 4.1 can be drawn for $N\to \infty $ and/or $\delta \to \infty $ .

Figure 4.1

An illustrative diagram for AC schemes.

As seen from detailed studies given in Tian and Du (Reference Tian and Du2013), some popular discretization schemes for nonlocal peridynamics fail to be AC. In particular, if $\delta $ is taken to be proportional to h, then as $h\to 0$ , piecewise constant conforming finite element solutions actually converge to the incorrect limit, similarly to those based on simple Riemann sum quadrature approximations to nonlocal operators. Similar discussions can be found in Bobaru et al. (Reference Bobaru, Yang, Alves, Silling, Askari and Xu2009) and Chen and Gunzburger (Reference Chen and Gunzburger2011) for diffusion models and Seleson, Du and Parks (Reference Seleson, Du and Parks2016) for systems of peridynamic equations.

Asymptotically compatible (AC) schemes, such as conforming Galerkin-type approximations of weak forms (Chen and Gunzburger Reference Chen and Gunzburger2011, Tian and Du Reference Tian and Du2014, Xu, Gunzburger, Burkardt and Du Reference Xu, Gunzburger, Burkardt and Du2016b ), nonconforming discontinuous Galerkin approximations (Du, Ju, Lu and Tian Reference Du, Ju, Lu and Tian2020a , Du, Ju and Lu Reference Du, Ju and Lu2019b , Du and Yin Reference Du and Yin2019) or collocation- or quadrature-based approximations of strong forms (Du and Tian Reference Du and Tian2014, Seleson et al. Reference Seleson, Du and Parks2016, Du et al. Reference Du, Tao, Tian and Yang2019c , Zhang, Gunzburger and Ju Reference Zhang, Gunzburger and Ju2016a , Zhang, Gunzburger and Ju Reference Zhang, Gunzburger and Ju2016b ), offer the potential to solve for approximations of a model of interest with different choices of parameters to gain efficiency and to avoid the pitfall of reaching inconsistent limits.

5.1. Asymptotically compatible conforming finite element methods for nonlocal diffusion

To delve deeper into finite element methods for nonlocal diffusion with a finite range of interactions, we first specialize to a more specific setting for which the local diffusion limit has been rigorously established. Consider the model given by (2.9) with $g=0$ . We also choose the kernel ${\gamma _{\delta }}{({\boldsymbol {x}},{\boldsymbol {y}})}$ as a rescaled translation-invariant kernel, that is,

(5.2) $$ \begin{align} {\gamma_{\delta}}{({\boldsymbol{x}},{\boldsymbol{y}})} =\displaystyle \frac{1}{\delta^{d+2}}\gamma\biggl(\displaystyle \frac{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|}{\delta}\biggr) \end{align} $$

for $\delta \in (0,1]$ and some kernel $\gamma ( {|\cdot |})$ that is assumed to be radial, nonnegative, compactly supported in $B_1({{\mathbf {0}}})$ (the unit ball centred at the origin), and has a bounded second-order moment, that is,

(5.3) $$ \begin{align} \widehat \gamma (|{\boldsymbol{\xi}} |)=|{\boldsymbol{\xi}}|^2\gamma(|{\boldsymbol{\xi}}|)\in L^1_{\text{loc}}({{\mathbb R}^{d}})\quad \text{and}\quad \int_{{B_1({{\mathbf{0}}})} } \widehat\gamma(|{\boldsymbol{\xi}}|) {\,{\mathrm{d}}} {\boldsymbol{\xi}}=d. \end{align} $$

For kernels satisfying these conditions, we refer to the discussions in Section 2.2. For such kernels, the local $\delta \to 0$ limit becomes the homogeneous Dirichlet boundary value problem of the Poisson equation. Let us denote the energy space ${V_{c,\delta }({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}={V_c({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ as defined in (2.4) to highlight the dependence on $\delta $ .

Now, for any fixed $\delta \in (0,1]$ , we introduce conforming finite element spaces $\{{V_{\delta , h}} \} \subset {V_{c,\delta }({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ associated with the triangulation $\tau _h=\{ K \}$ of the domain ${\Omega \cup \Omega _{{\mathcal I}_{\delta }}}$ (or ${\Omega }\cup {\Omega }_{{\mathcal I}_1}$ that contains the domain ${\Omega \cup \Omega _{{\mathcal I}_{\delta }}}$ for any $\delta <1$ ). We set

$$ \begin{align*} V_{\delta, h} := \{ v\in {V_{c,\delta}({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})} \colon v|_K \in P_p(K) \ \text{for all}\ K \in \tau_h \}, \end{align*} $$

where $P_p(K)$ denotes the space of polynomials on $K\in \tau _h$ of degree less than or equal to p. Again, for different $\delta $ , in order to have the finite element functions defined on a common spatial domain, we also assume, as for $ {V_{c,\delta }({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ , that any element in $V_{\delta , h}$ vanishes outside ${\Omega }$ .

As $h\to 0$ , we assume that $\{{V_{\delta , h}} \}$ is dense in $ {V_{c,\delta }({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ , that is, for any $v\in {V_{c,\delta }({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ , there exists a sequence $\{v_h\in {V_{\delta , h}}\}$ such that, for a given $\delta>0$ ,

(5.4) $$ \begin{align} \| v_h-v\|_{{V_{c,\delta}({\Omega\cup\Omega_{{\mathcal I}_{\delta}}})}}\to0\quad\text{as } h\to 0. \end{align} $$

These properties are easily satisfied by standard finite element spaces.

The Galerkin approximation is defined by replacing ${V_{c,\delta }({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ with $V_{\delta , h}$ in (2.9):

(5.5) $$ \begin{align} \text{Find } u_{h,\delta}\in {V_{\delta, h}} \text{ such that } {\mathcal A}_{\delta}(u_{h,\delta}, v_h) = {\langle{f},{v_{h}}\rangle} \text{ for all } v_h\in {V_{\delta, h}}. \end{align} $$

The analysis of Galerkin approximations of parametrized nonlocal variational problems can be formulated within the general framework of AC schemes (Tian and Du Reference Tian and Du2014). Other than the necessary properties of functions and operators that guarantee $u_{\delta }\to u_0$ as $\delta \to 0$ , where $u_0$ denotes the solution of the corresponding PDE, and standard density properties such as (5.4), one of the key ingredients for the asymptotic compatibility of the scheme is the asymptotic density property of the finite element space. Rather than a general definition, we provide a specific instance, adapted to the models studied here, as follows. A family of finite-dimensional spaces $\{ V_{\delta ,h}\subset {V_{c,\delta }({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ , $\delta \in (0,1)$ , $h\in (0, h_0]\}$ is asymptotically dense in $H_0^1(\Omega )$ , if, for all $ v\in H_0^1(\Omega )$ , there exists a sequence $\{v_n\in V_{\delta _n, h_n}\}_{h_n\to 0, \delta _n\to 0}$ as $n\to \infty $ such that

(5.6) $$ \begin{align} \| v-v_n\|_{H^1({\Omega})}\to 0 \quad \text{as } n\to\infty. \end{align} $$

Tian and Du (Reference Du, Huang and Lehoucq2014) have shown rigorously that, for scalar nonlocal diffusion equations such as (2.9), all conforming Galerkin approximations of the nonlocal models containing continuous piecewise linear functions are automatically AC in any space dimension. This means that they can recover the correct local limit as long as both $\delta $ and h are decreasing, even if the nonlocal parameter $\delta $ is reduced at a much faster pace than the mesh spacing h. Even though the analysis of the above AC property is highly technical, an intuitive explanation is that even with h larger than $\delta $ , the nonlocal features that ensure the correct local limit are still encoded in the stiffness matrices thanks to higher-order (than constant) basis functions. In fact, even though, for the class of kernels we are considering, discontinuous piecewise constant finite element spaces are conforming, they do not generally result in AC approximations. This was first noticed numerically by Chen and Gunzburger (Reference Chen and Gunzburger2011).

For general nonlocal systems in any dimension, Tian and Du (Reference Tian and Du2014, Reference Tian and Du2020) have shown that as long as the condition $h=o(\delta )$ is met as $\delta \to 0$ , then we obtain the correct local limit even for discontinuous piecewise constant finite element approximations when they are of the conforming type. Practically speaking, this implies that a mild growth of the bandwidth in the finite element stiffness matrix is needed as the mesh is refined in order to recover the correct local limit for piecewise constant finite element schemes. In fact, Tian and Du (Reference Tian and Du2013, Reference Tian and Du2015) have shown that if a constant bandwidth is kept as the mesh is refined, the approximations may converge to an incorrect local limit.

Naturally, schemes using higher-order basis functions tend to provide higher-order accuracy in the nonlocal setting as well, should the solutions enjoy sufficient regularity. At the moment, the theory on AC schemes does not offer any estimate of the order of convergence with respect to different couplings of h and $\delta $ . Preliminary numerical experiments in Tian and Du (Reference Tian and Du2014) offer some insight about the balance of modelling and discretization errors, but additional theoretical analyses need to be carried out, except for the case of Fourier spectral approximations of nonlocal models with periodic boundary conditions, for which precise error estimates can be found in Du and Yang (Reference Du and Yang2016, Reference Du and Yang2017) and Slevinsky, Montanelli and Du (Reference Slevinsky, Montanelli and Du2018).

5.2. Nonconforming and DG FEMs for nonlocal models with sufficiently singular kernels

The framework of AC schemes is very general. For example, the theory also guided the development of nonconforming discontinuous Galerkin approximations (Tian and Du Reference Tian and Du2015) for nonlocal diffusion and nonlocal peridynamic models with sufficiently singular interaction kernels. Consider the scalar model given by (2.9) with $g= 0$ . If the nonlocal interaction kernel satisfies

(5.7) $$ \begin{align} \int_{|{\boldsymbol{x}}|<\epsilon} |{\boldsymbol{x}}|\gamma_{\delta}(|{\boldsymbol{x}}|) {\,{\mathrm{d}}} {\boldsymbol{x}}=\infty\quad\text{for all } \epsilon\in (0,\epsilon_0], \end{align} $$

then it may be the case that discontinuous finite element solutions do not belong to the associated energy space and alternative formulations have to be developed. For technical reasons, Tian and Du (Reference Tian and Du2015) also assumed that

(5.8) $$ \begin{align} \lim_{\epsilon\to 0} \epsilon^2 \biggl( \int_{|{\boldsymbol{x}}|<\epsilon} |{\boldsymbol{x}}|^2\gamma_{\delta}(|{\boldsymbol{x}}|) {\,{\mathrm{d}}} {\boldsymbol{x}} \biggr)^{-1} = 0. \end{align} $$

Note that for kernels $\gamma _{\delta }(r)$ that behave like $1/r^{d+2s}$ as $r\to 0$ with $s\in (1/2,1)$ , conditions (5.7) and (5.8) are satisfied.

Tian and Du (Reference Tian and Du2015) introduced a nonconforming DG scheme based on the removal of the singularity in the nonlocal interaction kernel in a sufficiently small neighbourhood of the origin parametrized by a cut-off level n. In other words, $\gamma _{\delta }(r)$ in (2.9) is replaced with

(5.9) $$ \begin{align} \gamma_{\delta}^n (r) := \begin{cases} \gamma_{\delta}(r) & \text{if } \gamma_{\delta}(r)\leq n ,\\ n & \text{if } \gamma_{\delta}(r)> n. \end{cases} \end{align} $$

For nonlocal problems with the regularized kernel $\gamma _{\delta }^n (r)$ , discontinuous element spaces, such as conventional discontinuous finite element spaces, can be used as conforming Galerkin approximations, leading to a discrete solution $u_{h,\delta ,n}$ . Naturally, there are other ways to define the cut-off. The essential requirement is that the resulting modified kernel is both radial and integrable (for a given n) and it converges pointwise to the original kernel.

The goal is then to demonstrate that $u_{h,\delta ,n}$ converges to the solution $u_{\delta }$ of the original continuous nonlocal model with a singular kernel as $n\to \infty $ and $h\to 0$ . The convergence theory can be established by generalizing the relevant compactness results given in Bourgain, Brezis and Mironescu (Reference Bourgain, Brezis and Mironescu2001) as $n \to \infty $ for a given $\delta>0$ . Indeed, the following generalization is made in Tian and Du (Reference Tian and Du2015).

Compactness results

Given the kernels $\gamma _{\delta }^n$ and $\gamma _{\delta }$ and the corresponding energy spaces ${V^{n}_{c,\delta }({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ and ${V_{c,\delta }({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ , assume that the energy norms of $\{v_n\in {V^{n}_{c,\delta }({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}\}$ have the uniform bound

$$ \begin{align*} {\text{supp}_n} \int_{\Omega\cup\Omega_{{\mathcal I}_{\delta}}}\int_{\Omega\cup\Omega_{{\mathcal I}_{\delta}}} \gamma_{\delta}^n(|{\boldsymbol{y}}-{\boldsymbol{x}}|) (v_n({\boldsymbol{y}})-v_n({\boldsymbol{x}}))^2 {\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}} \leq C_0. \end{align*} $$

Then $\{v_n\}$ is relatively compact in $L^2({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})$ and any limit function v belongs to ${V_{c,\delta }({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ with

(5.10) $$ \begin{align} \int_{\Omega\cup\Omega_{{\mathcal I}_{\delta}}} \int_{\Omega\cup\Omega_{{\mathcal I}_{\delta}}} \gamma_{\delta}(|{\boldsymbol{y}}-{\boldsymbol{x}}|) (v({\boldsymbol{y}})-v({\boldsymbol{x}}))^2 {\,{\mathrm{d}}} {\boldsymbol{x}} {\,{\mathrm{d}}} {\boldsymbol{y}} \leq C_0 . \end{align} $$

The classical Bourgain–Brezis–Mironescu compactness result established in Bourgain et al. (Reference Bourgain, Brezis and Mironescu2001) can be seen as the local limit of the new compactness result for nonlocal spaces in Tian and Du (Reference Tian and Du2015) with a finite $\delta $ . By applying the framework of asymptotically compatible discretization (Tian and Du Reference Tian and Du2014, Tian and Du Reference Tian and Du2020), reviewed earlier with respect to the mesh parameter h and the cut-off level n, we can obtain the convergence of $u_{h,\delta ,n}$ to $u_{\delta }$ unconditionally as $n\to \infty $ and $h\to 0$ if the underlying finite element space contains the continuous piecewise linear finite element space. Moreover, in this case, if $n\to \infty $ , we expect that for any given $\delta $ and h, $u_{h,\delta ,n} \to u_{h,\delta }$ as $n\to \infty $ , where $u_{h,\delta }$ denotes the conforming finite element approximation of $u_{\delta }$ in the space ${V_{\delta ,h}}\cap {V_{c,\delta }({\Omega \cup \Omega _{{\mathcal I}_{\delta }}})}$ .

If piecewise constant finite elements are used, a conditional convergence theorem has been established in Tian and Du (Reference Tian and Du2015), provided the definition of the cut-off is suitably modified. For instance, consider a kernel of the type

(5.11) $$ \begin{align} \displaystyle \frac{\gamma_{\star}}{|{\boldsymbol{y}}-{\boldsymbol{x}}|^{d+2s}}\leq \gamma_{\delta}(|{\boldsymbol{x}} -{\boldsymbol{y}}|)\leq\displaystyle \frac{\gamma^{\star}}{|{\boldsymbol{y}}-{\boldsymbol{x}}|^{d+2s}} \quad \text{ for } {\boldsymbol{x}},\, {\boldsymbol{y}}\in{\Omega\cup\Omega_{{\mathcal I}_{\delta}}} , \end{align} $$

for some $s\in (1/2,1)$ and positive constants $\gamma _{\star }$ and $\gamma ^{\star }$ . Then one may let

(5.12) $$ \begin{align} \gamma_{\delta}^n(r) := \begin{cases} \gamma_{\delta}(r) & \text{for } r\geq1/n ,\\ \gamma_{\delta}(1/n) & \text{for } 0<r<1/n. \end{cases} \end{align} $$

For the approximate solution $u_{h,\delta ,n}$ defined as the Galerkin approximation to the nonlocal problem with kernel $\gamma _{\delta }^n$ using piecewise constants, we have that $\| u_{h,\delta ,n}-u_{\delta } \|_{L^2}\to 0$ if $h=o(1/n)$ as $n\to \infty $ .

We note that the nonconforming approximations discussed in Tian and Du (Reference Tian and Du2015) in the local limit do not yield a standard nonconforming finite element approximation or a DG approximation to the local problem. One may construct other alternative formulations that can give rise to the conventional nonconforming and DG discretization of local PDEs; see e.g. Du et al. (Reference Du, Ju, Lu and Tian2020a ) for a study based on the DG with penalty formulation.

5.3. Adaptive mesh refinement for nonlocal models

Due to reduced sparsity, for instance, nonlocal models generally incur greater computational costs than their local PDE-based counterparts. Thus, designing effective adaptive methodologies is important, and it is an area worthy of attention. So far the work has been limited to nonlocal models with a finite range of interactions.

Du, Ju, Tian and Zhou (Reference Du, Ju, Tian and Zhou2013b ) provided an a posteriori error analysis of conforming finite element methods for solving linear nonlocal diffusion and peridynamic models. The approach adopted is a residual-type error estimator in the $L^2$ -norm, e.g. of the form $\|{-}\mathcal {L}_{\delta }(u_{\delta ,h})-f\|_{L^2}$ , which remains well-defined and can be easily computed for kernel functions that are both radial and integrable from element-wise contributions without worrying about flux jumps across element boundaries. This is in sharp contrast to the case of second-order elliptic PDEs. The theory of a posteriori error analysis has been rigorously derived for nonlocal volume-constrained problems associated with scalar equations. The reliability and efficiency of the estimators are proved, and relationships between nonlocal and classical local a posteriori error estimates are also studied.

Du, Tian and Zhao (Reference Du, Tian and Zhao2013c ) have also developed a convergent adaptive finite element algorithm for the numerical solution of scalar nonlocal models. For problems involving certain radial but non-integrable kernel functions, the convergence of the adaptive algorithm is rigorously derived with the help of several basic ingredients, such as an upper bound on the estimator, the estimator reduction, and the orthogonality property. How these estimators and methods work in the local limit and for general time-dependent and nonlinear peridynamic models remains to be investigated.

For nonlocal problems having solutions with jump discontinuities, an adaptive finite element method is given in Xu et al. (Reference Xu, Gunzburger and Burkardt2016a ). There, an algorithm is developed that first detects the location of the discontinuity and then refines the grid near the discontinuity. To preserve the $h^2$ accuracy possible with the use of piecewise linear elements even when the exact solution contains jump discontinuities at unknown locations, the elements surrounding the discontinuity should have thickness of ${O}(h^4)$ across the discontinuity. This was already observed for the one-dimensional case in Chen and Gunzburger (Reference Chen and Gunzburger2011) and Xu et al. (Reference Xu, Gunzburger, Burkardt and Du2016b ). In higher dimensions, a naive refinement strategy that results in small, well-shaped elements in the vicinity of a $({d}-1)$ -dimensional surface can thus result in an excessive number of degrees of freedom. The adaptive refinement strategy of (Xu et al. Reference Xu, Gunzburger and Burkardt2016a ) instead results in elongated elements having thickness ${O} (h^4)$ across the discontinuity but length ${O} (h)$ along the discontinuity, as illustrated in Figure 5.1. The presence of elongated elements is not harmful to the error because of the anisotropic behaviour of the solution, that is, it is smooth along the discontinuity but discontinuous across the discontinuity. Whereas robust meshing algorithms for this type of anisotropic refinement and accurate predictions of solution jump discontinuities remain computationally challenging in higher space dimensions, the numerical examples in Xu et al. (Reference Xu, Gunzburger and Burkardt2016a ) illustrate that the adaptive strategy developed there does indeed result in $h^2$ convergence, as is also the case for the one-dimensional numerical results in Chen and Gunzburger (Reference Chen and Gunzburger2011) and Xu et al. (Reference Xu, Gunzburger, Burkardt and Du2016b ).

Figure 5.1

A coarse-grid illustration of a mesh resulting from the adaptive strategy discussed in this section. The elongated yellow elements surround the discontinuity in the solution.

Another useful feature of the adaptive strategy developed in Xu et al. (Reference Xu, Gunzburger and Burkardt2016a ) is that the transition between the elongated elements along discontinuities and well-shaped elements with sides of length ${O} (h)$ away from the discontinuity can be abrupt, that is, there is no need to have a transition zone in which elements gradually grow in size. This is also illustrated in Figure 5.1. This feature is already present in the methods developed in Chen and Gunzburger (Reference Chen and Gunzburger2011) and Xu et al. (Reference Xu, Gunzburger, Burkardt and Du2016b ) for the one-dimensional case. The importance of this feature is that it also helps to keep down the number of degrees of freedom.

The next means of savings in the number of degrees of freedom is to use discontinuous basis functions only in the elongated elements, and continuous basis functions in elements in which the exact solution is smooth. The final means is to switch from the nonlocal model to the corresponding PDE model in all elements that do not interact with the elongated elements that surround the discontinuity, i.e. all elements whose vertices are at distance greater than $\delta $ from the vertices of the elongated elements.

The four degrees of freedom-saving features (i.e. elongated elements, abrupt grid transitions, the use of discontinuous basis function only in the elongated elements, and the switch to a local PDE model away from the discontinuity), result in tremendous savings in the number of degrees of freedom. In fact, the resulting number of degrees of freedom is comparable to that for a local PDE model using a regular grid of elements having sides of ${O} (h)$ . Moreover, the relative savings become greater as the grid size is reduced and as the dimension ${d}$ is increased. See Xu et al. (Reference Xu, Gunzburger and Burkardt2016a ) for several numerical illustrations supporting these conclusions.

With regard to more general coupling strategies of local and nonlocal models for computational effectiveness, relevant works are discussed in a recent review on the subject (Du Reference Du2019). AC schemes have also been studied in the context of multiscale problems modelled by local–nonlocal coupling formulations with a spatially heterogeneous horizon (Du and Tian Reference Du and Tian2018, Tao, Tian and Du Reference Tao, Tian and Du2019), based on a new trace theorem for nonlocal function spaces that provides a stronger version of the classical local counterparts (Tian and Du Reference Tian and Du2017).

5.4. Approximations of fractional models as limit of nonlocal models with a finite range of interactions

One of the key messages we want to get across is to show that many popularly studied fractional PDEs are either specialized nonlocal models or can be treated as limiting cases.

Naturally, for fractional differential operators defined in integral form, by truncating the fractional kernel in the fractional operators to a finite range measured by the horizon parameter $\delta $ , we end up with a nonlocal model parametrized by $\delta $ . Hence we see from Section 3.1.1 that on the continuum level, the fractional models may then be viewed as the infinite horizon limit of nonlocal models with a finite $\delta $ , after a suitable scaling of the kernel.

D’Elia and Gunzburger (Reference D’Elia and Gunzburger2013) and Burkovska and Gunzburger (Reference Burkovska and Gunzburger2019b ) have shown that as $\delta \to \infty $ , the Galerkin approximations of the nonlocal model with the parameter $\delta $ can converge to the solution of the fractional equations, provided h is taken to be sufficiently small. One can actually apply the framework of the AC schemes to show that all conforming Galerkin approximations are AC in the limit $\delta \to 0$ (Tian et al. Reference Tian, Du and Gunzburger2016). That is, the discretization of fractional equations associated with the fractional Laplacian can be viewed as the global (with an infinite nonlocal horizon parameter) limit of nonlocal models with properly scaled fractional-type nonlocal interaction kernels (Tian et al. Reference Tian, Du and Gunzburger2016). Moreover, the convergence does not require the dependence of h on $\delta $ .

Let us give an illustration of the work presented in Tian et al. (Reference Tian, Du and Gunzburger2016). Consider the volume-constrained problems defined for the fractional Laplacian in $\Omega $ with homogeneous Dirichlet condition in ${\Omega _{{\mathcal I}_{\infty }}}={{\mathbb R}^{d}}\setminus {\Omega }$ . By truncations of both the spatial domain and the range of nonlocal interactions, we may end up with a class of parametrized problems

(5.13) $$ \begin{align} \begin{cases} -{\mathcal L}_{\delta} u({\boldsymbol{x}})=-\int_{{{\mathbb R}^{d}}}(u({\boldsymbol{y}})-u({\boldsymbol{x}}))\gamma_{\delta}(|{\boldsymbol{y}}-{\boldsymbol{x}}|) {\,{\mathrm{d}}} {\boldsymbol{y}}=f &\text{for all } {\boldsymbol{x}}\in \Omega ,\nonumber\\ u=0 &\text{for all } {\boldsymbol{x}}\in{{\Omega_{{\mathcal I}_{\delta}}}} \end{cases}\\ \end{align} $$

with

(5.14) $$ \begin{align} \gamma_{\delta}(|{\boldsymbol{y}}-{\boldsymbol{x}}|)= \begin{cases} \displaystyle \frac{{C_{d,s,\delta}}}{ |{\boldsymbol{y}}-{\boldsymbol{x}}|^{d+2s} } & {\boldsymbol{y}}\in B_{\delta}({\boldsymbol{x}}) , \\[9pt] 0 & {\boldsymbol{y}}\in {{\mathbb R}^{d}} \backslash B_{\delta}({\boldsymbol{x}}) . \end{cases} \end{align} $$

Here, $C_{d,s,\delta }={C_{{d},s}}$ is given as in (3.3) for the limiting case of $\delta \to \infty $ so that ${\mathcal L}_{\infty } = (-\Delta )^{s}$ as in (3.1). Note that in the limit of $\delta \to 0$ , we take a $C_{d,s,\delta }$ to be a different scaling factor so that its local limit is ${\mathcal L}_0 =-\Delta $ .

Let h denote the discretization parameter associated with Galerkin approximations of (5.13). For example, h could be the mesh parameter for the finite element discretization or the reciprocal of the number of basis functions in spectral approximations. We can then apply the AC framework to obtain the convergence of the Galerkin approximations of (5.13) to the solution of the fractional equation as $h\to 0$ and $\delta \to \infty $ as long as the finite-dimensional approximations spaces are subspaces of the underlying energy space, that is, as long as we adopt conforming discretizations.

As pointed out in Tian et al. (Reference Tian, Du and Gunzburger2016) and Tian and Du (Reference Tian and Du2020), analyses of the fractional and local limits ( $\delta \to \infty $ and $\delta \to 0$ respectively) of nonlocal models parametrized by a finite $\delta $ and with fractional-type kernels demonstrate that nonlocal models are more general than their fractional and local counterparts and they also serve as a bridge between fractional and local models; see Figure 5.2 and, again, Figure 4.1 and similar diagrams given in Du (Reference Du2019) and Tian and Du (Reference Tian and Du2020). For AC schemes, the bridging roles of nonlocal diffusion with a finite range of nonlocal interactions, presented above at the continuum level, remain valid on the discrete level.

Figure 5.2

Limits of numerical solutions of nonlocal models with finite interaction radius $\delta $ as $\delta \to 0$ (local limit, top) and as $\delta \to \infty $ (fractional limit, bottom).

6.1. Quadrature rules

The fractional Poisson equation (3.5) leads to a dense linear algebraic system

(6.2) $$ \begin{align} {{\boldsymbol{A}}}_{s}\vec{u}=\vec{f}, \end{align} $$

in which the entries in the matrix ${{\boldsymbol {A}}}_{s}=\{\mathcal {A}_{s} (\varphi _{i},\varphi _{j}) \}_{ij}$ involve singular integrals. In order to compute these entries, we decompose the integrals into contributions between pairs of elements $K,\widehat {K}\in {{\mathcal T}_{h}}$ and between pairs consisting of elements $K\in {{\mathcal T}_{h}}$ and external edges $e\in \partial {{\mathcal T}_{h}}$ as follows:

$$ \begin{align*} \mathcal{A}_{ss}(\varphi_{i},\varphi_{i}) = \sum_{K}\sum_{\widehat{K}} \mathcal{A}_{s}^{K\times\widehat{K}}(\varphi_{i},\varphi_{j}) + \sum_{K}\sum_{e} \mathcal{A}_{s}^{K\times e}(\varphi_{i},\varphi_{j}). \end{align*} $$

The individual contributions $\mathcal {A}_{s}^{K\times \widehat {K}}$ and $\mathcal {A}_{s}^{K\times e}$ are given by

$$ \begin{align*} \mathcal{A}_{s}^{K\times\widehat{K}}(\varphi_{i},\varphi_{j}) & = \displaystyle \frac{C(d,s)}{2} \int_{K} {\,{\mathrm{d}}} {\boldsymbol{x}} \int_{\widehat{K}} {\,{\mathrm{d}}} {\boldsymbol{x}} \displaystyle \frac{ (\varphi_{i}({\boldsymbol{x}})-\varphi_{i}({\boldsymbol{y}})) (\varphi_{j}({\boldsymbol{x}})-\varphi_{j}({\boldsymbol{y}})) }{{|{\boldsymbol{x}}-{\boldsymbol{y}}|}^{{d}+2s}} , \\ \mathcal{A}_{s}^{K\times e}(\varphi_{i},\varphi_{j}) &= \displaystyle \frac{C(d,s)}{2s} \int_{K} {\,{\mathrm{d}}} {\boldsymbol{x}} \int_{e} {\,{\mathrm{d}}} {\boldsymbol{x}} \displaystyle \frac{\varphi_{i} ({\boldsymbol{x}}) \varphi_{j} ({\boldsymbol{x}}) ~ \vec{n}_{e}\cdot ({\boldsymbol{x}}-{\boldsymbol{y}}) }{ {|{\boldsymbol{x}}-{\boldsymbol{y}}|}^{{d}+2s}}. \end{align*} $$

Contributions from non-disjoint pairs of elements (see Figure 6.1) are not directly amenable to numerical quadrature, due to their singular nature. Fortunately, these can be treated by adapting techniques used in the boundary element method (BEM) literature to address similar issues arising from singular kernels (Sauter and Schwab Reference Sauter and Schwab2010).

Figure 6.1

Classes of element pairs of configurations that need to be handled: pairs of identical elements (red), element pairs with common edge (yellow), with common vertex (blue) and separated element pairs (green).

However, the fractional Laplacian does pose new difficulties beyond those addressed by the BEM literature, but which can be treated as described by Ainsworth and Glusa (Reference Ainsworth, Glusa and Dick2018b ) and Acosta et al. (Reference Acosta, Bersetche and Borthagaray2017). In particular, Ainsworth and Glusa (Reference Ainsworth and Glusa2017, Reference Ainsworth, Glusa and Dick2018b ) have developed non-uniform order Gauss-type quadrature rules. Alternatively, one could take the approach proposed by Chernov, von Petersdorff and Schwab (Reference Chen and Gunzburger2011). It allows for transforming quadrature rules given on the unit hypercube $[0,1]^{2d}$ to any pair of elements $K\times \widehat {K}$ . Here, the singularity is taken into account through the choice of the weight in the quadrature rules.

Ainsworth and Glusa (Reference Ainsworth and Glusa2017) gave the following result. Denote the quadrature approximation to the bilinear form $\mathcal {A}_{s} (\cdot ,\cdot ) $ by $\mathcal {A}_{s}^{Q} (\cdot ,\cdot ) $ . Then, by using ${O} (\log N_h^{2d})$ quadrature points per element pair $K\times \tilde {K}$ , the consistency error due to quadrature ${|\mathcal {A}_{s}(u,v)-\mathcal {A}_{s}^{Q}(u,v)|}$ is dominated by the discretization error.

A brief discussion about solvers and condition numbers for finite element discretizations of the integral fractional Laplacian is given in Section 12.2.

7.1. Truncation in the extended direction

The truncation approach was pioneered by Nochetto, Otárola and Salgado (Reference Nochetto, Otárola and Salgado2015). A first result shows that truncation of the semi-infinite cylinder ${{\Omega }_z}$ only leads to an exponentially small error.

Let $U\in \mathcal {H}^{1}_{\beta }$ denote the solution of (7.2) and let $U_{z_{\text {trun}}}$ be the extension by zero to ${{\Omega }_z}$ of the solution to (7.2) posed on the truncated domain ${{\Omega }}_{z_{\text {trun}}}$ , $z_{\text {trun}}\geq 1$ , with a homogeneous Dirichlet condition enforced at $z=z_{\text {trun}}$ . Then

$$ \begin{align*} {\|U-U_{z_{\text{trun}}}\|}_{\mathcal{H}^{1}_{\beta}({{\Omega}_z})} \leq C {\mathrm{e}}^{-\sqrt{\lambda_{0}}z_{\text{trun}}/4} {\|f\|}_{H^{-s}(\Omega)}, \end{align*} $$

where $\lambda _{0}$ is the smallest eigenvalue of the integer-order Laplacian as given by (3.13).

Now, let ${{\mathcal T}_{h}}$ denote a quasi-uniform mesh of $\Omega $ with mesh size h and define a graded mesh $\mathcal {T}_{z_{\text {trun}}}$ of the interval $[0,z_{\text {trun}}]$ by

$$ \begin{align*} z_{k}= (k/M) ^{\gamma}z_{\text{trun}}, \quad k=0,\ldots,M, \end{align*} $$

with $M\sim h^{-1}$ and $\gamma>3/(1-\beta )$ . Then define a mesh $\mathcal {T}$ on the extended domain ${{\Omega }}_{z_{\text {trun}}}$ as the tensor product of ${{\mathcal T}_{h}}$ and $\mathcal {T}_{z_{\text {trun}}}$ , and let ${V}_{h}$ denote the space of piecewise linear finite element functions on $\mathcal {T}$ . Let $N=\operatorname {dim}{V}_{h}$ be the overall number of degrees of freedom and choose the truncation to be $z_{\text {trun}}\sim {|\log N|}$ . Then the following convergence result holds.

Let $f\in H_c^{1-s}(\Omega )$ , $U\in \mathcal {H}^{1}_{\beta }({{\Omega }_z})$ be the solution of (7.2), and let $U_{h}\in {V}_{h}$ be the solution of the restriction of (7.2) to ${V}_{h}$ . Moreover, let $u=\operatorname {trace}_{\Omega }U$ denote the solution to the spectral fractional-order Poisson problem (3.14), and let $u_{h}=\operatorname {trace}_{\Omega }U_{h}$ be its discrete approximation. Then

$$ \begin{align*} {\|u-u_{h}\|}_{H_c^{s}(\Omega)} &\leq C {\|U-U_{h}\|}_{\mathcal{H}^{1}_{\beta}({{\Omega}_z})}\\[2pt] &\leq C {|\log h|}^{s} h {\|f\|}_{H_c^{1-s}(\Omega)} \\[2pt] &\leq C {|\log N|}^{s} N^{-1/({d}+1)} {\|f\|}_{H_c^{1-s}(\Omega)}. \end{align*} $$

As observed by Nochetto et al. (Reference Nochetto, Otárola and Salgado2015), this shows that the result is optimal in terms of regularity but sub-optimal in terms of complexity. This shortcoming can be addressed using sparse grids (Banjai et al. Reference Banjai, Melenk, Nochetto, Otárola, Salgado and Schwab2019) or $hp$ -finite elements (Meidner, Pfefferer, Schürholz and Vexler Reference Meidner, Pfefferer, Schürholz and Vexler2018, Banjai et al. Reference Banjai, Melenk, Nochetto, Otárola, Salgado and Schwab2019). Using the latter, exponential convergence can be obtained, including for cases with low regularity of the domain and incompatibility of the forcing term (Banjai et al. Reference Banjai, Melenk, Nochetto, Otárola, Salgado and Schwab2019).

7.2. Spectral method in the extended direction

A different approach was taken by Ainsworth and Glusa (Reference Ainsworth and Glusa2018a ). Instead of truncating the domain, a spectral method is used for the extended direction, based on the observation that the eigenfunctions of the extension problem (7.1) are given by

$$ \begin{align*} \varphi_{m}({\boldsymbol{x}})\psi_{m} ({z}) , \quad m\in\mathbb{N}, \end{align*} $$

where

(7.3) $$ \begin{align} \psi_{m} ({z}) = c_{s} (\lambda_{m}^{1/2}{z}) ^{s} K_{s} (\lambda_{m}^{1/2}{z}) , \end{align} $$

and where $c_{s}= 2^{1-s} / \Gamma (s) $ , $ (\varphi _{m},~\lambda _{m}) $ are the eigenpairs of (3.13) and $K_{\nu }$ are the modified Bessel functions of the second kind.

This motivates the choices

$$ \begin{align*} V_{h}&=\{v_{h}\in H^{1}_{0} (\Omega) \colon {v_{h}}|_{K}\in \mathbb{P}_{k} (K) \ \text{for all}\ K\in\mathcal{T}_{h}\}, \\[2pt] {V}_{h,M}&=\biggl\{V_{h,M}=\sum_{m=0}^{M-1}v_{h,m} ({\boldsymbol{x}}) \widehat{\psi}_{m} ({z}) \colon v_{h,m}\in V_{h} \biggr\}\subset \mathcal{H}_{\beta}^{1} ({{\Omega}_z}) , \end{align*} $$

where

(7.4) $$ \begin{align} \widehat{\psi}_{m} ({z}):= c_{s} (\widehat{\lambda}_{m}^{1/2}{z}) ^{s} K_{s} (\widehat{\lambda}_{m}^{1/2}{z}) \end{align} $$

are defined using suitable approximations $\widehat {\lambda }_{m}$ instead of the true eigenvalues $\lambda _{m}$ .

It turns out that using the asymptotic law for the eigenvalues of the Laplacian and using a coarse finite element approximation gives sufficiently good approximations $\widehat {\lambda }_{m}$ . Moreover, by decimation, only ${|\log h|}^{p}$ for some $p>0$ eigenvalue approximations is required. Overall, this results in the following error estimate.

Let $f\in H_c^{r}(\Omega )$ , $r\geq -s$ , let $U\in \mathcal {H}^{1}_{\beta }({{\Omega }_z})$ be the solution of (7.2) and let $U_{h}\in {V}_{h}$ be the solution of the restriction of (7.2) to ${V}_{h}$ . Moreover, let $u=\operatorname {trace}_{\Omega }U$ be the solution to the spectral fractional-order Poisson problem (3.14) and let $u_{h}=\operatorname {trace}_{\Omega }U_{h}$ be its discrete approximation. Then

$$ \begin{align*} {\|u-u_{h}\|}_{H_c^{s}(\Omega)} &\leq C {\|U-U_{h,M}\|}_{\mathcal{H}^{1}_{\beta}({{\Omega}_z})}\\ & \leq C {|f|}_{H_c^{r}(\Omega)}h^{\min\{k,r+s\}}\sqrt{{|\log h|}} \\ &= C {|f|}_{H_c^{r}(\Omega)}N^{-\min\{k,r+s\}/{d}}{|\log N|}^{p}. \end{align*} $$

9.1. Spectral-Galerkin methods in unbounded domains

For $s\in (0,1)$ , we consider the fractional diffusion problem in ${\Omega }={{\mathbb R}^{d}}$ given by

(9.1) $$ \begin{align} \begin{cases} (-\Delta)^s u ({\boldsymbol{x}}) + {c} u({\boldsymbol{x}}) = f({\boldsymbol{x}}) &\text{for all } {\boldsymbol{x}}\in{{\mathbb R}^{d}} ,\\ u({\boldsymbol{x}})\to 0 & \text{as } |{\boldsymbol{x}}|\rightarrow\infty, \end{cases} \end{align} $$

where ${c}$ denotes a given positive constant and f denotes a given source term. A weak formulation of problem (9.1) is as follows:

(9.2) $$ \begin{align} \begin{aligned} &\text{Given a positive constant } {c} \text{ and } f\in H^{-s}({{\mathbb R}^{d}}), \\ &\text{seek } u\in H^s({{\mathbb R}^{d}}) \text{ such that }\\ &\qquad{{\mathcal A}_s}(u,v) = (f,v ) \quad\text{for all } v\in H^s({{\mathbb R}^{d}}), \end{aligned} \end{align} $$

where, for all $u,v\in H^s({{\mathbb R}^{d}})$ ,

$$ \begin{align*} {\mathcal A}_s(u,v) = \displaystyle \frac{{C_{{d},s}}}{2}\int_{{{\mathbb R}^{d}}}\int_{{{\mathbb R}^{d}}} \displaystyle \frac{ (u({\boldsymbol{y}}) - u({\boldsymbol{x}})) (v({\boldsymbol{y}})-v({\boldsymbol{x}})) }{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|^{{d}+2s}}{{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}} + {c}(u,v) \end{align*} $$

with ${C_{{d},s}}$ is defined in (3.3). The well-posedness of problem (9.2) follows directly from the Lax–Milgram theorem.

Several successful strategies have been developed for determining approximations of solutions of local PDE diffusion problems posed on ${{\mathbb R}^{d}}$ ; see e.g. the survey paper by Shen and Wang (Reference Shen and Wang2009) and the references cited therein. One way is to truncate the domain ${{\mathbb R}^{d}}$ , provided that the solution decays rapidly. However, it is known that solution of the fractional diffusion equation decays slowly with a power law at infinity so that the naive truncation approach results in a relatively large error due to domain truncation.Footnote ³ Another approach is to design an artificial boundary condition imposed on the boundary of a truncated domain that results in the same solution within a bounded region; the use of this approach remains largely an open question in the fractional case. A third approach, which is the approach used here, does not involve domain truncation, but instead approximates the solution in terms of orthogonal functions that potentially can well approximate the solution in the unbounded domain.

For simplicity, we restrict our discussion to the one-dimensional case, i.e. ${d}=1$ ; higher-dimensional cases can be treated in a similar manner by using tensor products of one-dimensional orthogonal polynomials.

9.1.1. Approximation by Hermite polynomials

The orthonormal Hermite polynomials $\{H_n(x)\}$ are defined by the three-term recurrence relation (Abramowitz and Stegun Reference Abramowitz and Stegun1972, Mao and Shen Reference Mao and Shen2017)

$$ \begin{align*} H_0 &= \pi^{-1/4},\qquad H_1(x) = \sqrt{2}\pi^{-1/4}x,\\ H_{n+1} &= x\sqrt{\displaystyle \frac{2}{n+1}}H_n(x)- \sqrt{\displaystyle \frac{n}{n+1}}H_{n-1}(x)\quad\text{for } n=1,2,\ldots. \end{align*} $$

It is well known that those polynomials form an orthonormal basis in the weighted space $L_{{w}}^2({\mathbb R})$ , where the weight function is given by ${w}(x)={\mathrm {e}}^{-x^2}$ . Let $V_N^H$ denote the space spanned by the Hermite polynomials of degree less than or equal to N, and let

$$ \begin{align*} X_N^H = \{ v \colon \sqrt{{w}(x)}u(x) \ \text{for all } u \in V_N^H \}. \end{align*} $$

We define the $L^2_{w}({\mathbb R})$ -projection operator $Q^H_N\colon L^2_{w}({\mathbb R})\to X_N^H$ by

$$ \begin{align*} \int_{{\mathbb R}} u(x) v_N(x){w}(x){\,{\mathrm{d}}} x= \int_{{\mathbb R}} (Q^H_N u)(x) v_N(x){w}(x){\,{\mathrm{d}}} x \quad\text{for all } v_N\in X_N^H. \end{align*} $$

We also define the modified operator $\widehat Q^H_N\colon L^2({\mathbb R}) \rightarrow X_N^H$ by

$$ \begin{align*} \widehat Q^H_N u = {w}^{1/2}Q^H_N (u\,{w}^{-1/2}) \quad\text{for all } u\in L^2({\mathbb R}). \end{align*} $$

Then we have the approximation property

(9.3) $$ \begin{align} \| \widehat Q^H_N u - u \|_{H^s({\mathbb R})} \le {C} N^{(\ell-m)/2} \| u \|_{{\widehat{H}}^m({\mathbb R})} \quad\text{for all } s\in[0,m], \end{align} $$

where ${C}>0$ does not depend on u or N and where the modified Sobolev space ${\widehat {H}}^m({\mathbb R})$ comes equipped with the norm and seminorm

$$ \begin{align*} \| u \|_{{\widehat{H}}^m({\mathbb R})}^2 = \sum_{0\le k\le m} \biggl\| \biggl( {\displaystyle \frac{{\mathrm{d}}}{{\mathrm{d}} x}} + x\biggr)^k u \biggr\|_{L^2({\mathbb R})}^2 , \\ | u |_{{\widehat{H}}^m({\mathbb R})} = \biggl\| \biggl( {\displaystyle \frac{{\mathrm{d}}}{{\mathrm{d}} x}} + x\biggr)^m u \biggr\|_{L^2({\mathbb R})}^2, \end{align*} $$

respectively.

The Hermite–Galerkin discretization of the fractional diffusion problem (9.1) is then as follows:

(9.4) $$ \begin{align} &\text{Given a positive constant } {c} \text{ and } f\in H^{-s}({{\mathbb R}^{d}}),\nonumber \\ &\text{seek } u_N(x) \in X_N^H \text{ such that } \nonumber\\ &{\mathcal A}_s(u_N, v_N) = (I_N^H f, v_N ) \quad\text{for all } v_N\in X_N^H, \end{align} $$

where $I_N^H \colon C({\mathbb R})\rightarrow {X_N^H}$ denotes the Gauss–Hermite interpolation operator with respect to the Gauss–Hermite points (Shen, Tang and Wang Reference Shen, Tang and Wang2011). Then the standard energy estimate and the approximation property (9.3) result in the error estimate

(9.5) $$ \begin{align} \| u - u_N \|_{H^s({\mathbb R})} \le {C} N^{s-m} | u |_{{\widehat{H}}^m({\mathbb R})} + {C} N^{1/6-\ell/2} | f |_{{\widehat{H}}^{\ell}({\mathbb R})} , \end{align} $$

provided that $u\in {\widehat {H}}^m({\mathbb R})$ and $f\in {\widehat {H}}^{\ell }({\mathbb R})$ , where again ${C}>0$ does not depend on u, f or N. This estimate indicates that the numerical solution fails to converge exponentially, because either u and or f will decay algebraically as $x\to \infty $ ; see Mao and Shen (Reference Mao and Shen2017, Section 5).

The Hermite–Galerkin discretization (9.4) is simple to implement because the mass matrix is diagonal, whereas the stiffness matrix can be computed efficiently using the Fourier transform. An implementation of the Hermite–Galerkin method (9.4) using Fourier transforms is provided in Mao and Shen (Reference Mao and Shen2017). See also Mao and Shen (Reference Mao and Shen2017) and Tang, Yuan and Zhou (Reference Tang, Yuan and Zhou2018) for related Hermite collocation methods, whose error estimate also depend on the Fourier transform of the solution u and the source term f.

9.1.2. Approximation by modified Gegenbauer polynomials

In Tang, Wang, Yuan and Zhou (Reference Tang, Wang, Yuan and Zhou2019), the approach taken in Section 9.1.1 is extended to approximations in which the basis is constructed by using Gegenbauer polynomials that are modified using a nonlinear singular mapping. The resulting basis has proved to be better suited to approximating functions with algebraic decay rates (see e.g. Boyd Reference Boyd1987, Boyd Reference Boyd2001) when compared to classical bases of orthogonal polynomials such as Hermite or Laguerre polynomials.

The Gegenbauer polynomials, denoted by ${G}_{n}^{{\lambda }}(t)$ with $ t\in (-1,1)$ and scaling parameter ${\lambda }>-1/2$ , are defined by the three-term recurrence relation

$$ \begin{align*} {G}_{0}^{{\lambda}}(t)&=1,\quad {G}_{n}^{{\lambda}}(t)=2{\lambda} t,\\ n{G}_{n}^{{\lambda}}(t)&=2t(n+{\lambda}-1){G}_{n-1}^{{\lambda}}(t)-(n+2{\lambda}-2){G}_{n-2}^{{\lambda}}(t), \;\; n=2,3,\ldots.\\[-10pt] \end{align*} $$

The Gegenbauer polynomials are orthogonal with respect to the weight function ${w}_{\lambda }(t)=(1-t^2)^{{\lambda }-1/2}$ , and they are closely related to the hypergeometric functions (Gradshteyn and Ryzhik Reference Gradshteyn and Ryzhik2007, page 1000) that are widely used for the analysis of fractional differential equations; see e.g. Ervin, Heuer and Roop (Reference Ervin, Heuer and Roop2018).

We seek approximate solutions in the space

$$ \begin{align*} X_{N}^{G} =\text{span}\{R_{n}^{{\lambda}}(x) \colon n=0,1,\ldots,N\}, \end{align*} $$

where the mapped Gegenbauer functions $R_{n}^{{\lambda }}(x)$ are defined by

$$ \begin{align*} R_{n}^{{\lambda}}(x):= (1+x^2) ^{-\tfrac{{\lambda}+1}{2}}{G}_{n}^{{\lambda}}\biggl(\displaystyle \frac{x}{\sqrt{1+x^2}}\biggr)\quad\text{for all } x\in{{\mathbb R}}. \end{align*} $$

The $L^2({\mathbb R})$ -orthogonal projection operator $Q_{N}^G: L^2({\mathbb R})\to X_N^G$ is defined by

$$ \begin{align*} \int_{{\mathbb R}} u(x) v_N(x) {\,{\mathrm{d}}} x = \int_{{\mathbb R}} Q_{N}^G u(x) v_N(x) {\,{\mathrm{d}}} x \quad\text{for all } v_N\in X_N^G . \end{align*} $$

Next, let $B_m^G({\mathbb R})$ denote the weighted Sobolev space defined by the norm and seminorm

$$ \begin{align*} \| u \|_{B_m^G({\mathbb R})}^2 & = \sum_{0\le k\le m} \biggl\| (1+x^2)^{\tfrac 1 4-\tfrac{{\lambda}+m} 2}\biggl( (1+x^2)^{\tfrac32} \displaystyle \frac{{\mathrm{d}}}{{\mathrm{d}} x}\biggr)^k \Bigl((1+x^2)^{\tfrac{\lambda+1}2}u\Bigr) \biggr\|_{L^2({\mathbb R})}^2 , \\ | u |_{B_m^G({\mathbb R})} & = \biggl\| (1+x^2)^{\tfrac 1 4-\tfrac{{\lambda}+m} 2} \biggl( (1+x^2 )^{\tfrac32} \displaystyle \frac{{\mathrm{d}}}{{\mathrm{d}} x}\biggr)^m \Bigl((1+x^2)^{\tfrac{\lambda+1}2}u\Bigr) \biggr\|_{L^2({\mathbb R})}^2, \end{align*} $$

respectively. Then, for any $u\in H^s({\mathbb R})\cap B_m^G({\mathbb R})$ with integer $ 1\le m\le N+1$ , $ s\in (0,1)$ and ${\lambda }>-1/2$ , we have the approximation property

$$ \begin{align*} \\[-10pt]\| Q_{N}^G u-u \|_{H^s({\mathbb R})}\le {C} N^{s-m} |u|_{B_m^G({\mathbb R})},\\[-10pt] \end{align*} $$

where ${C}>0$ does not depend on N and u.

The mapped Gegenbauer–Galerkin discretization of the fractional diffusion problem (9.1) is as follows:

$$ \begin{align*} &\text{ Given a positive constant } {c} \text{ and } f\in H^{-s}({{\mathbb R}^{d}}),\\[-3pt]&\text{ seek } u_N(x) \in X_N^G \text{ such that} \\[-25pt] \end{align*} $$

$$ \begin{align*} &{\mathcal A}_s(u_N, v_N) = (I_N^G f, v_N ) \quad\text{for all } v_N\in X_N^G, \end{align*} $$

where $I_N^G\colon C({\mathbb R})\rightarrow X_N^G$ denotes the mapped Gegenbauer–Gauss interpolation operator based on the Gegenbauer–Gauss quadrature points (Shen et al. Reference Shen, Tang and Wang2011, Chapter 3). Then, for any $u\in H^s(\mathbb R)\cap {B}_{m}^G(\mathbb R)$ and $f\in {B}_{\ell }^G(\mathbb R) $ with integers $ 1\le m,\ell \le N+1$ and ${\lambda }>-1/2$ , we have the error estimate

(9.6) $$ \begin{align} \| u-u_N\|_{H^s({\mathbb R})}\le {C} N^{s-m} |u|_{{ B}_{m}^G({\mathbb R})}+{C} N^{-\ell} |f|_{{B}_{\ell}^G({\mathbb R})}, \end{align} $$

where ${C}$ denotes a positive constant having value independent of N, u and f. The proof relies on the standard energy estimate as well as the approximation properties of $Q_N^G$ and $I_N^G$ . Compared to the Hermite–Galerkin spectral method of Section 9.1.1, we have that $(-\Delta )^s R_n^{\lambda }(x)$ is not as compact as for the Hermite polynomials, but can be easily approximated by a finite series of hypergeometric functions.

9.1.3. Approximation by modified Chebyshev polynomials

It is noteworthy that $|{\boldsymbol \xi }|^{2s}$ , the symbol of the Fourier transform of the fractional Laplacian, is non-separable so that direct extensions to higher dimensions of methods developed for the one-dimensional setting can be very complicated. To overcome this computational difficulty, Sheng et al. (Reference Sheng, Shen, Tang, Wang and Yuan2019) proposed a reformulation of the weak formulation (9.2) by the Dunford–Taylor formula (see Theorem 4.1 of Bonito, Lei and Pasciak Reference Bonito, Lei and Pasciak2019, and also Section 3.3.1), and constructed a spectral-Galerkin method having Fourier-like bi-orthogonal mapped Chebyshev polynomials as basis functions. This method leads to a diagonalized system that can be efficiently solved via the fast Fourier transform (FFT), and the approach can be easily extended to higher-dimensional problems via tensor products of the basis functions.

For any $u,v\in H^{s}({\mathbb R})$ , by the Dunford–Taylor formula, we have (see (3.15))

$$ \begin{align*} &\displaystyle \frac{{C_{{d},s}}}{2}\int_{{{\mathbb R}^{d}}}\int_{{{\mathbb R}^{d}}} \displaystyle \frac{ (u({\boldsymbol{y}}) - u({\boldsymbol{x}})) (v({\boldsymbol{y}})-v({\boldsymbol{x}})) }{|{{\boldsymbol{y}}-{\boldsymbol{x}}}|^{{d}+2s}} {{\,{\mathrm{d}}} {\boldsymbol{y}} {\,{\mathrm{d}}} {\boldsymbol{x}}} \\ & \quad =\displaystyle \frac{2 \sin (\pi s)}{\pi} \int_0^{\infty} t^{1-2s} \int_{{{\mathbb R}^{d}}} ((-\Delta)( I-t^2 \Delta)^{-1} u ) ({\boldsymbol{x}})\, v({\boldsymbol{x}}) {{\,{\mathrm{d}}} {\boldsymbol{x}}} {\,{\mathrm{d}}} t, \end{align*} $$

where I denotes the identity operator. Letting $w({\boldsymbol {x}}):= ( I-t^2 \Delta )^{-1} u({\boldsymbol {x}})$ , we can rewrite the weak formulation (9.2) as: seek $u\in H^s({{\mathbb R}^{d}})$ such that for all $v\in H^s({{\mathbb R}^{d}})$

(9.7) $$ \begin{align} {\mathcal A}_s(u,v)=\displaystyle \frac{2 \sin (\pi s)}{\pi} \int_0^{\infty} t^{-1-2s} (u-w, v) {\,{\mathrm{d}}} t+ {c}\,(u,v) =(f,v), \end{align} $$

where $w=w({\boldsymbol {x}};u,t)\in H^{1}({{\mathbb R}^{d}})$ is, for any $t>0$ , the solution of

(9.8) $$ \begin{align} t^{2} (\nabla w,\nabla \psi) +(w, \psi) =(u, \psi) \quad\text{for all } \psi \in H^{1}({{\mathbb R}^{d}}). \end{align} $$

The spectral-Galerkin approximation to (9.7)–(9.8) is then given as follows: find $u_N\in {X_N^C}$ such that

(9.9) $$ \begin{align} {\mathcal A}_s(u_N,v_N)&=\displaystyle \frac{2 \sin (\pi s)}{\pi} \int_0^{\infty} t^{-1-2s} (u_N-w_N, v_N) \, {\text d}t+{c}(u_N,v_N) \nonumber \\ &=(I_N f,v_N) \quad\text{for all } v_N\in X_{N}^C, \end{align} $$

where $w_N({\boldsymbol {x}};u_N,t)\in X_{N}^C$ is, for any $t>0$ , the solution of

(9.10) $$ \begin{align} t^{2} (\nabla w_N,\nabla \psi_{{N}}) +(w_N, \psi_{{N}}) =(u_N, \psi_{{N}}) \quad\text{for all } \psi_{{N}} \in X_{N}^C. \end{align} $$

To solve the numerical scheme (9.9)–(9.10), Sheng et al. (Reference Sheng, Shen, Tang, Wang and Yuan2019) used a basis of Fourier-like mapped Chebyshev polynomials. Recall that the classical Chebyshev polynomials are defined by, for $n=0,1,\ldots ,$

$$ \begin{align*} T_n(y)=\cos(n\, {\text { arc}\,\cos}(y))\quad\text{with } y\in (-1,1). \end{align*} $$

Then we define the mapped Chebyshev polynomials by

(9.11) $$ \begin{align} R_{n} (x) := (1+x^2) ^{-{1}/{2}}T_n\biggl(\displaystyle \frac{x}{\sqrt{1+x^2}}\biggr),\quad x\in{\mathbb R}. \end{align} $$

It is well known that the mapped Chebyshev functions are orthonormal in $L^2(\mathbb R)$ but are not orthogonal in $H^1({{\mathbb R}^{d}})$ . Following the spirit of Shen and Wang (Reference Shen and Wang2009), in Sheng et al. (Reference Sheng, Shen, Tang, Wang and Yuan2019) the approximation space is defined by $X_N^C =\text {span}\{R_{n}\colon n=0,1,\ldots ,N\}$ for which a Fourier-like basis that is bi-orthogonal can be constructed, that is, the basis functions are orthogonal with respect to both the $L^2$ and $H^1$ inner products. As a result, (9.10) is diagonalizable, as is (9.9), and hence they can be efficiently solved by an FFT related to Chebyshev polynomials.

The integration in terms of t in (9.9) can be evaluated exactly by using a sinc quadrature scheme (Bonito et al. Reference Bonito, Lei and Pasciak2019) whose computational cost is negligible compared with that of an FFT. Furthermore, this fast algorithm can be easily applied in higher-dimensional settings by using tensor products of the one-dimensional basis functions. The resulting complexity of solving (9.9) is then of ${O} (N(\log N)^{d})$ , where N is the number of degrees of freedom. The relevant error estimate is similar to that in (9.6) with ${\lambda }=0$ .

Remark 9.1. Although the solution to the fractional diffusion problem (9.1) has a simple and closed form in the Fourier domain, all existing spectral methods fail to achieve at an exponential convergence rate. This is because the solution generally decays at an algebraic rate at infinity, even if the source term is smooth and compactly supported. The proposed Hermite and the mapped rational/Chebyshev polynomials can only approximate exponentially decaying functions with spectral accuracy. How one designs a spectral method that accurately captures the algebraic decaying behaviour of the solution to (9.1) remains a largely open problem.

All the estimates stated in this subsection are derived with respect to the energy norm, e.g. $H^s({{\mathbb R}^{d}})$ . There is, of course, interest in obtaining higher approximation rates with respect to the $L^2({{\mathbb R}^{d}})$ -norm. However, the approximation properties of orthogonal polynomials are derived in weighted Sobolev spaces, so the usual duality arguments are not directly applicable here. Thus optimal $L^2({{\mathbb R}^{d}})$ -norm error estimates are still not well understood and warrant further investigation.□

9.2. Spectral-Galerkin methods in bounded domains

Now we turn to fractional diffusion in a bounded domain. Owing to the low regularity of the solution of the underlying problems involving the fractional Laplacian (see e.g. Ervin, Heuer and Roop Reference Ervin, Heuer and Roop2018, Grubb Reference Grubb2015b , Ros-Oton and Serra Reference Ros-Oton and Serra2014), classical spectral methods cannot approximate the solution very well. In this section we focus on the particular fractional Poisson problem (3.5) with $g({\boldsymbol {x}})=0$ , that is, on the problem

(9.12) $$ \begin{align} \begin{cases} (-\Delta)^{s}u = f &\text{for all }{\boldsymbol{x}}\in{\Omega}, \\ u=0&\text{for all }{\boldsymbol{x}}\in{{\Omega_{{\mathcal I}_{\infty}}}}={{\mathbb R}^{d}}\setminus{\Omega} , \end{cases} \end{align} $$

where $f({\boldsymbol {x}})$ is a smooth function and ${\Omega }$ is a bounded domain with a smooth boundary.

Mao, Chen and Shen (Reference Mao, Chen and Shen2016) considered the one-dimensional setting with ${\Omega }=(-1,1)$ and $s\in (1/2,1)$ . In this case, the fractional Laplacian is equivalent to the Riesz-fractional differential derivative that is defined by, for $u\in C_c^{\infty }(-1,1)$ and $x\in (0,1)$

$$ \begin{align*} -{D^{2s}} u(x) := - \displaystyle \frac{1}{2\cos((1-s)\pi)}\displaystyle \frac{{\mathrm{d}}^2}{{\mathrm{d}} x^2}\int_{-1}^1 |x-t|^{1-2s} u(t) {\,{\mathrm{d}}} t =(-\Delta)^s u (x) . \end{align*} $$

Then problem (9.12) is equivalent to the two-point boundary value problem: find the function u satisfying

(9.13) $$ \begin{align} \begin{cases} -{D^{2s}} u(x) =f(x) \quad\text{for all } x\in(-1,1),\\ u(0)=u(1)= 0. \end{cases} \end{align} $$

The discussion in Mao et al. (Reference Mao, Chen and Shen2016) is motivated by the observation that

(9.14) $$ \begin{align} (-\Delta)^{s}\{(1-x^2)^s_+ P_n^{s,s}(x) \} = \displaystyle \frac{{\Gamma}(2s+n+1)}{{\Gamma}(n+1)} P_n^{s,s}(x), \end{align} $$

where $P_n^{a,b}$ with $a, b>-1$ denotes the classical Jacobi polynomials that are orthogonal with respect to the weight function ${w}^{a,b}(x) = (1-x)^a(1+x)^b$ . This implies that if the source term $f(x)$ can be expanded in terms of the Jacobi polynomials $P_n^{s,s}(x)$ , then the solution of the fractional diffusion problem (9.12) is a series in $(1-x^2)^s_+ P_n^{s,s}(x)$ .

By choosing the test space $V_N = \text {span}\{ P_n^{s,s}(x) \colon n=0,\ldots ,N \}$ and the trial space $X_N = \text {span}\{ (1-x^2)_+^s P_n^{s,s}(x) \colon n=0,\ldots ,N \}$ , a Petrov–Galerkin spectral method is given as follows: find $u_N \in X_N$ such that, for all $v_N\in V_N$ ,

$$ \begin{align*} \int_{-1}^1 (-\Delta)^{s} u_N(x) \, v_N(x) (1-x^2)^s {\,{\mathrm{d}}} x = \int_{-1}^1 f(x) \, v_N(x) (1-x^2)^s {\,{\mathrm{d}}} x. \end{align*} $$

Using the relation (9.14) and the approximation property of Jacobi polynomials (Shen et al. Reference Shen, Tang and Wang2011, Theorem 3.35), the following error estimate in the weighted $L^2$ -norm is derived in Mao et al. (Reference Mao, Chen and Shen2016, Theorem 5):

$$ \begin{align*} \|u_N - u \|_{B^0_{-s}} \le {C} N^{-2s-m} \| f \|_{B^m_s} \quad \text{for any integer} ~~m \ge 0, \end{align*} $$

where the norm $\| \cdot \|_{B^m_{\lambda }}$ is defined by

$$ \begin{align*} \| v \|_{B^m_{\lambda}}^2 = \sum_{j=0}^m \int_{-1}^1 | v^{(j)} |^2 (1-x^2)^{j+\lambda} {\,{\mathrm{d}}} x. \end{align*} $$

If m is not an integer, the space $B^m_{\lambda }$ is defined by interpolation (Guo and Wang Reference Guo and Wang2004). The argument also applies to the case when $s\in (0,1/2]$ and ${\Omega }$ is multi-interval (Acosta, Borthagaray, Bruno and Maas Reference Acosta, Borthagaray, Bruno and Maas2018).

The above estimate indicates the exponential convergence under the assumption that the source term f is smooth. The results in Mao et al. (Reference Mao, Chen and Shen2016) extend the studies of spectral methods for solving one-sided fractional equations found in Chen, Shen and Wang (Reference Chen, Shen and Wang2016b ), Zayernouri and Karniadakis (Reference Zayernouri and Karniadakis2014) and Zayernouri, Ainsworth and Karniadakis (Reference Zayernouri, Ainsworth and Karniadakis2015b ); see also Mao and Karniadakis (Reference Mao and Karniadakis2018) and Samiee, Zayernouri and Meerschaert (Reference Samiee, Zayernouri and Meerschaert2019a , Reference Samiee, Zayernouri and Meerschaert2019b ) for general two-sided fractional equations, Zayernouri, Ainsworth and Karniadakis (Reference Zayernouri, Ainsworth and Karniadakis2015a ) for tempered fractional diffusion equations, and Deng, Zhang and Zhao (Reference Deng, Zhang and Zhao2019) for the study of superconvergence points.

A similar analysis may be extended to higher-dimensional cases. Xu and Darve (Reference Xu and Darve2018b ) studied a spectral method for solving the fractional Poisson problem (9.13) with ${\Omega }=B_1({{\mathbf {0}}})\subset \mathbb {R}^{d}$ , ${d}=2,3$ , that is, the unit ball in two or three dimensions. The basic idea is to construct a spectral basis by using Jacobi polynomials and spherical harmonic polynomials (Dyda, Kuznetsov and Kwaśnicki Reference Dyda, Kuznetsov and Kwaśnicki2017a ). In particular, spherical harmonic polynomials of degree $\ell \ge 0$ form a finite-dimensional space, having dimension

$$ \begin{align*} M_{{d},\ell}=\displaystyle \frac{{d}+2\ell-2}{{d}+\ell-2} \binom{{d} +\ell-2}{\ell}. \end{align*} $$

For each $\ell $ we fix a linear basis for this space, denoted by $V_{\ell ,m}$ with $m = 1,\ldots , M_{{d},\ell }$ , that is orthonormal with respect to the surface measure on the unit sphere. Then the functions

$$ \begin{align*} \phi_{\ell,m,n}({\boldsymbol{x}}) = V_{\ell,m}({\boldsymbol{x}})P_n^{s,{{d}}/2+\ell-1}(2|{\boldsymbol{x}}|^2-1)\quad \text{with } \ell\ge0, \ n\ge 0,\ 0\le m\le M_{{d},\ell} \end{align*} $$

form an orthogonal basis in $L^2_{w}({\Omega })$ with respect to the weight function ${w}({{\boldsymbol {x}}})=(1-|{\boldsymbol {x}}|^2)^s$ (Dyda et al. Reference Dyda, Kuznetsov and Kwaśnicki2017a ). Analogous to (9.14), we have

$$ \begin{align*} (-\Delta)^{s} ((1-|{\boldsymbol{x}}|^2)_+^s \phi_{\ell,m,n} ) = \displaystyle \frac{2^{2s}\Gamma(s+n+1)\Gamma({{d}}/{2}+\ell+s+n)}{n!\Gamma({{d}}/{2}+n)} \phi_{\ell,m,n} ({\boldsymbol{x}}). \end{align*} $$

Xu and Darve (Reference Xu and Darve2018b ) have developed and analysed fast and accurate spectral methods based on this property and the approximation property of the Jacobi polynomials.

Unfortunately, the above argument heavily relies on the relation (9.14) so that, due to the singular behaviour of the solution near the boundary, it cannot be directly applied to the fractional diffusion with lower-order terms. Zhang (Reference Zhang2019) studied the regularity of the one-dimensional fractional diffusion–reaction model

(9.15) $$ \begin{align} \begin{cases} (-\Delta)^{s} u(x) + {c} u(x)=f(x) & \text{for all } x\in (-1,1),\\ u(x)= 0 & \text{for all } x\in {\mathbb R}\setminus(-1,1) \end{cases}\quad \end{align} $$

with the constant $s\in (0,1)$ and ${c}>0$ . In particular, Zhang proved that $(1-x^2)^{-s}u$ belongs to $B_{s}^{(2s\wedge m) +2s }$ if $f\in B_{s}^{m}$ . Moreover, the regularity index could be slightly improved, that is,

$$ \begin{align*} (1-x^2)^{-s}u \in B_{s-1}^{(6s-\epsilon)\wedge (1+4s)}\quad \text{provided that } f \in B_{s-1}^{m} \text{ for large } m. \end{align*} $$

The convergence rate of the spectral method certainly deteriorates accordingly, due to a lack of sufficient solution regularity. See also Mao and Shen (Reference Mao and Shen2018) for a spectral element method based on geometric meshes.