# Asymptotically compatible discretization of multidimensional nonlocal diffusion models and approximation of nonlocal Green’s functions

Asymptotically compatible discretization of multidimensional nonlocal diffusion models and... Abstract Nonlocal diffusion equations and their numerical approximations have attracted much attention in the literature as nonlocal modeling becomes popular in various applications. This paper continues the study of robust discretization schemes for the numerical solution of nonlocal models. In particular, we present quadrature-based finite difference approximations of some linear nonlocal diffusion equations in multidimensions. These approximations are able to preserve various nice properties of the nonlocal continuum models such as the maximum principle and they are shown to be asymptotically compatible in the sense that as the nonlocality vanishes, the numerical solutions can give consistent local limits. The approximation errors are proved to be of optimal order in both nonlocal and asymptotically local settings. The numerical schemes involve a unique design of quadrature weights that reflect the multidimensional nature and require technical estimates on nonconventional divided differences for their numerical analysis. We also study numerical approximations of nonlocal Green’s functions associated with nonlocal models. Unlike their local counterparts, nonlocal Green’s functions might become singular measures that are not well defined pointwise. We demonstrate how to combine a splitting technique with the asymptotically compatible schemes to provide effective numerical approximations of these singular measures. 1. Introduction In this paper, we study numerical approximations of the following multidimensional linear nonlocal diffusion equation: \begin{align} - \mathcal{L}_{\delta}u^{\delta} (\mathbf{x})=f(\mathbf{x}),\quad \mathbf{x}\in\Omega\subset \mathbb{R}^{d}, \end{align} (1.1) with Ω being a given domain in $$\mathbb{R}^{d}$$ (d ≥ 2), f a given right-hand side, uδ the solution to be sought and the nonlocal operator $$\mathcal{L}_{\delta }$$ defined by \begin{align} \mathcal{L}_{\delta} u(\mathbf{x})= 2 \int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\rho_{\delta}(|\mathbf{z}|)\left(u(\mathbf{x}+\mathbf{z})-u(\mathbf{x})\right)\, \mathrm{d}\mathbf{z}. \end{align} (1.2) The operator $$\mathcal{L}_{\delta }$$ is parametrized by a positive horizon parameter δ measuring the range of nonlocal interactions. The specific form of such nonlocal interactions is prescribed by a non-negative kernel function ρδ = ρδ(|z|). The integral in (1.2) is interpreted in the principal value sense whenever needed (Mengesha & Du, 2014). More discussions on the nonlocal model are given in Section 2. A main task of this work is to develop quadrature-based finite difference approximation schemes and associated Green’s functions for (1.1). The approximation schemes are able to preserve the maximum principle at the discrete level and they are shown to be asymptotically compatible in the sense that as the nonlocality vanishes, the numerical solutions can give consistent local limits. They are also shown to be of optimal order in both nonlocal and asymptotically local settings. Previous studies of schemes enjoying such optimal-order error estimates, discrete maximum principles and asymptotical compatibility have been mostly confined to the one-dimensional space (Tian & Du, 2013). While the algorithmic development is in a similar spirit to that for the one-dimensional case, the multidimensional extension made in this work involves both new design elements that reflect the multidimensional nature and more technical convergence analysis relying on new interpolation error estimates of nonconventionally defined weighted divided differences. The motivation for our work is rooted in the need to develop robust approximations of nonlocal models that serve as alternatives to classical partial differential equations (PDEs). Recently, there have been many studies on the application of nonlocal modeling to problems in mechanics, physics and materials, biological and social sciences (Bates & Chmaj, 1999; Applebaum, 2004; Gilboa & Osher, 2008; Bobaru & Duangpanya, 2010; Buades et al., 2010; Lou et al., 2010; Tadmor & Tan, 2014). The integral formulations of spatial interactions in nonlocal models can account for nonlocal effects and allow more singular solutions. An example is the nonlocal peridynamics (PD) theory (Silling, 2000) and its applications in studying cracks and materials failure, as well as other mechanical properties and physical processes (Askari et al., 2008; Silling & Lehoucq, 2010; Silling et al., 2010a; Palatucci et al., 2012). Rigorous mathematics of nonlocal models has proved to be necessary in order to gain fundamental insights and to guide the modeling and simulation efforts (Silling & Askari, 2005; Macek & Silling, 2007; Bobaru et al., 2009; Andreu et al. 2010; Kilic & Madenci, 2010; Zhou & Du, 2010; Chen & Gunzburger, 2011; Du & Zhou, 2011; Du et al., 2012, 2013, 2017b). For nonlocal models like (1.1), as $$\delta \rightarrow 0$$, the nonlocal interactions that define the model can become localized so that the zero-horizon (or local) limit of the nonlocal operator, when valid both physically and mathematically, can be represented by a local differential operator. Given such a scenario the corresponding nonlocal model naturally converges to a conventional differential equation model in the local limit; see for instance, the convergence analysis of linear state-based PD models to the classical Navier equation of linear elasticity (Mengesha & Du, 2014). There have been concerns about whether such consistency can be preserved at the discrete level when the nonlocal models are discretized numerically. Discrete schemes of nonlocal models that preserve the correct limiting behavior are called asymptotically compatible (AC) schemes, a notion developed first in the studies by Tian & Du (2013, 2014). In other words, the numerical approximation given by an AC scheme can reproduce the correct local limiting solution as the horizon δ and the mesh spacing (denoted by h) approach zero, if the convergence to such a limiting local solution is valid on the continuum level. While Tian & Du (2013) presented results concerning the AC property for a one-dimensional scalar model solved by a number of different discretization methods, Tian & Du (2014) managed to develop more general results that are applicable to linear systems in multidimensions based on Galerkin finite element approximations on unstructured meshes. A few subsequent studies have been carried out. For example, Tao et al. (2017) extended the study of AC schemes to nonlocal diffusion equations with Neumann-type volume constraints, Du & Yang (2016) discussed AC schemes based on Fourier spectral methods for problems defined on periodic cells and Du & Yang (2017) proposed an efficient and accurate hybrid algorithm to implement the Fourier spectral methods, while Du et al. (2017b) developed discontinuous-Galerkin-based AC schemes. Moreover, AC schemes were proposed in the study by Chen et al. (2017) for nonlocal time–space models and in the study by Du et al. (2016) for the robust recovery of the nonlocal gradient of the solutions to nonlocal models. For technical reasons, much of the numerical analysis in these works, with the exception of the study by Tian & Du (2014) on conforming Galerkin finite element approximations, has focused on nonlocal operators defined for a one-dimensional spatial variable. In comparison with other AC schemes such as those based on the Galerkin finite element and spectral methods, quadrature-based collocation-type finite difference schemes not only offer optimal-order convergent approximations with simpler implementations but also preserve the discrete maximum principle and asymptotic compatibility; see details in Section 3. Furthermore, these quadrature-based schemes are closely related to other discretizations based on strong forms, such as the original mesh-free discretization developed in the study by Silling (2000). As an application of the discretization scheme, we consider in Section 4 the numerical computation of nonlocal Green’s functions Gδ = Gδ(x, y) for the nonlocal diffusion model (1.1). Nonlocal Green’s functions are important tools for nonlocal continuum models (Weckner et al., 2009; Silling et al., 2010b; Wang et al., 2016, 2017), just like their local versions, i.e., conventional Green’s functions associated with PDEs. Since analytical expressions can be expected only in simple cases, the numerical study of nonlocal Green’s functions has particular significance. We show here that while the local/conventional and nonlocal Green’s functions enjoy many similarities and are intimately connected through the local limiting process, there are also some fundamental differences among them. In particular, for a large class of nonlocal interaction kernels the associated Green’s functions Gδ = Gδ(x, y) on the continuum level for δ > 0 are measure-valued distributions, which cannot be interpreted as conventional functions defined pointwise for x and y. These differences are essential features that, if improperly handled, could present undesirable effects on their numerical approximations. This serves as another reminder of the subtlety involved in the nonlocal modeling. To better deal with these complications we discuss a splitting approach that utilizes the developed AC scheme to get effective and robust approximations of nonlocal Green’s functions. 2. Multidimensional nonlocal diffusion model For the sake of a clear illustration we let u = u(x) denote a scalar deformation field and consider (1.1) subject to a nonlocal constraint of the Dirichlet type, \begin{align}\ u=0 \;\; \textrm{on} \;\Omega_{\mathcal{I}}, \end{align} (2.1) where $$\Omega _{\mathcal{I}}$$ is a boundary layer defined as $$\Omega_{\mathcal{I}}=\{\mathbf{x}\notin \Omega|\;\textrm{dist}(\mathbf{x},\Omega)\le \delta\}.$$ For more discussions on nonlocal constraints defined on a domain with a nonzero volume we refer to Du et al. (2012, 2013). Furthermore, let us consider suitably scaled kernels so that as δ goes to zero, the local limit of the nonlocal operator $$\mathcal{L}_{\delta }$$ is exactly the Laplace operator Δ (which is denoted by $$\mathcal{L}_{0}$$). To this end, it requires that the kernel is a non-negative and nonincreasing function and has a finite second moment in z independent of the horizon parameter δ: $$\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}|\mathbf{z}|^{2}\rho_{\delta}(|\mathbf{z}| ) \, \mathrm{d}\mathbf{z}=d.$$ The limiting local model of (1.1) is then given by \begin{align} - \mathcal{L}_{0}u^{0}=f\;\; \textrm{in } \Omega\quad \mbox{ and }\quad u^{0}=0 \;\;\mbox{ on }\; \partial \Omega, \end{align} (2.2) where, without loss of generality, $$\mathcal{L}_{0}=\Delta$$ for the particular case studied here. We note that a popular choice of ρδ is a rescaled kernel given by \begin{align} \rho_{\delta}(|\mathbf{z}|)=\frac{1}{\delta^{2+d}}\rho\left(\frac{|\mathbf{z}|}{\delta}\right)\!\quad \forall\, \mathbf{z}\in \mathcal{B}_{\delta}(\boldsymbol{0}), \end{align} (2.3) where ρ = ρ(ξ) is a non-negative and nonincreasing function with compact support in [0, 1] and a normalized moment \begin{align} {\int_{0}^{1}}\rho(\xi)\xi^{1+d}\, \mathrm{d}\xi=c_{d} \end{align} (2.4) for a given constant cd, though our discussion here is not restricted to such a rescaled form. 3. Quadrature-based finite difference schemes and asymptotic compatibility AC schemes provide robust numerical approximations of nonlocal models (Tian & Du, 2013, 2014) since the convergence of such schemes is insensitive to the choices of modeling and discretization parameters, in particular with respect to either a sufficiently small horizon δ or a sufficiently refined mesh spacing (denoted by h). 3.1 The discretization scheme We now develop an AC quadrature-based finite difference discretization to the nonlocal diffusion equation (1.1) that extends a similar scheme presented first in the study by Tian & Du (2014) and also studied by Du & Tian (2014) for the one-dimensional case. For simplicity, we let $$\left \{\mathbf{x}_{\mathbf{j}}\in \Omega \cup \Omega _{\mathcal{I}}\right \}$$ be the set of nodes (grid points) of a uniform Cartesian mesh $${\mathcal{T}}_{h}$$ with mesh size h. Here, j denotes a multiindex corresponding to $$\mathbf{x}_\mathbf{j}$$ = hj. We note that the particular ordering of the nodes affects the matrix structure corresponding to the discrete system but does not affect the numerical solution. First, at any node $$\mathbf{x}_\mathbf{i}$$ ∈ Ω, the nonlocal operator can be rewritten as \begin{align} \mathcal{L}_{\delta}u(\mathbf{x}_{\mathbf{i}}) = 2\int_{\mathcal{B}_{\delta}({\boldsymbol{0}})}\dfrac{u(\mathbf{z}+\mathbf{x}_\mathbf{i})-u(\mathbf{x}_{\mathbf{i}})}{W(\mathbf{z})}W(\mathbf{z})\rho_{\delta}(|\mathbf{z} | )\, \mathrm{d}\mathbf{z}, \end{align} (3.1) where W(z) represents a weight function. Then a quadrature-based finite difference scheme of the nonlocal operator (1.2) can be given as \begin{align} \mathcal{L}_{h,\,\delta}u(\mathbf{x}_{\mathbf{i}}) = 2\int_{\mathcal{B}_{\delta}({\boldsymbol{0}})}\mathcal{I}_{h}\left(\dfrac{u(\mathbf{z}+\mathbf{x}_{\mathbf{i}})-u(\mathbf{x}_{\mathbf{i}})}{W(\mathbf{z})}\right)W(\mathbf{z})\rho_{\delta}(|\mathbf{z} | )\, \mathrm{d}\mathbf{z}, \end{align} (3.2) where $$\mathcal{I}_{h}(\cdot )$$ represents the piecewise d-multilinear interpolation operator in z associated with the mesh $${\mathcal{T}}_{h}$$, that is, $$\mathcal{I}_{h}(u)$$ is piecewise linear with respect to each component of the spatial variable. The nonlocal constraint with u = 0 is imposed at nodes in $$\Omega _{\mathcal{I}}$$. For the one-dimensional case, in the study by Tian & Du (2014) the function W(z) = |z| was used for $$z\in \mathbb{R}$$. A key ingredient in its extension to the high-dimensional case is given by \begin{align} W(\mathbf{z})=\dfrac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}, \end{align} (3.3) where the notation |⋅|1 stands for the ℓ1 norm in the d-dimensional vector space while |⋅| denotes the standard Euclidean norm. The auxiliary function W = W(z) plays an important role in ensuring asymptotic compatibility for the multidimensional case. The resulting quadrature-based finite difference scheme for solving the nonlocal diffusion problem can be more conveniently written as the following: find $$\left \{u_{h}^{\delta }(\mathbf{x}_{\mathbf{i}})\right \}$$ such that $$u_{h}^{\delta }(\mathbf{x}_{\mathbf{i}})=0$$ for $$\mathbf{x}_{\mathbf{i}}\in \Omega _{\mathcal{I}}$$ and \begin{align} -\mathcal{L}_{h,\,\delta}u_{h}^{\delta}(\mathbf{x}_{\mathbf{i}}) = \sum_{\mathbf{x}_{\mathbf{j}}\in \mathcal{B}_{\delta}(\mathbf{x }_{\mathbf{i}})}\left(u_{h}^{\delta}(\mathbf{x}_{\mathbf{j}})-u_{h}^{\delta}(\mathbf{x}_{\mathbf{i}})\right)p^{\delta}_{\mathbf{x }_{\mathbf{i}},\mathbf{x}_{\mathbf{j}}} = f(\mathbf{x}_{\mathbf{i}}) \quad \forall\; \mathbf{x}_{\mathbf{i}}\in\Omega, \end{align} (3.4) where the sum over $$\mathbf{x }_\mathbf{j}$$ is for grid points in $$\mathcal{B}_{\delta }(\mathbf{x }_{\mathbf{i}})$$ but not including $$\mathbf{x }_\mathbf{i}$$, $$p^{\delta }_{\mathbf{x}_{\mathbf{i}},\mathbf{x}_{\mathbf{j}}}={2 \beta _{\mathbf{j}\mathbf{i}}}/{W(\mathbf{x}_{\mathbf{j}}-\mathbf{x}_{\mathbf{i}})}$$ with \begin{align} \beta_{\mathbf{j}\mathbf{i}} = \int_{ \mathcal{B}_{\delta}( \boldsymbol{0} ) } \phi_{\mathbf{j}}(\mathbf{z}+\mathbf{x}_{\mathbf{i}}) W(\mathbf{z}) \rho_{\delta}(|\mathbf{z}|) \, \mathrm{d}\mathbf{z}, \end{align} (3.5) and $$\phi_\mathbf{j}$$ is the piecewise multilinear basis function satisfying $$\phi_\mathbf{j} (\mathbf{x}_\mathbf{i})$$ = 0 when $$\mathbf{i}\neq \mathbf{j}$$ and $$\phi_\mathbf{j}(\mathbf{x}_\mathbf{j})$$ = 1. Remark 3.1 The discrete scheme (3.4) is valid for general kernels ρδ = ρδ(z) with bounded second moment in z. This is another benefit of using the weight W = W(z) in (3.3). Indeed, consider a kernel ρδ(|z|) = cδ, α, d|z|−α with α ∈ [d, d + 2) for $$\mathbf{z}\in \mathcal{B}_{\delta }(\boldsymbol{0})$$; then the integral in (3.5) may be unbounded for |i −j| = 1 if we do not use such a weight, or say, if W(z) = 1 is taken instead. Remark 3.2 By the symmetry and translation invariance of the uniform mesh it can be observed that $$\beta _{\mathbf{j}\mathbf{i}}=\beta _{\delta }\left (\mathbf{x}_{\mathbf{i}}-\mathbf{x}_{\mathbf{j}}\right )$$ for some function βδ. In fact, let ϕ0 be the piecewise multilinear basis function located at the origin satisfying ϕ0$$(\mathbf{x}_\mathbf{i})$$ = 0 with $$\mathbf{x}_{\mathbf{i}}\neq{\boldsymbol{0}}$$ and ϕ0(0) = 1; then \begin{align} \beta_{\mathbf{j}\mathbf{i}} = \beta_{\delta}(\mathbf{x}_{\mathbf{i}}-\mathbf{x}_{\mathbf{j}}),\qquad \beta_{\delta}(\mathbf{x})= \int_{ \mathcal{B}_{\delta}( \boldsymbol{0} ) } \phi_{0}(\mathbf{z}+\mathbf{x}) W(\mathbf{z}) \rho_{\delta}(|\mathbf{z}|) \,\mathrm{d}\mathbf{z} . \end{align} (3.6) Moreover, it is easy to see that βδ is symmetric with respect to any axis hyperplane. Also, $$\beta_{\mathbf{j}\mathbf{i}}$$ = 0 if the support of $$\phi_\mathbf{j}$$ does not overlap with the ball Bδ$$(\mathbf{x}_\mathbf{i})$$. In comparison, we give the conventional central finite difference scheme for the local diffusion equation (2.2): \begin{align} - \mathcal{L}_{h,\,0}{u_{h}^{0}}(\mathbf{x}_{\mathbf{i}})=-\sum_{|\mathbf{x}_{\mathbf{j}} - \mathbf{x }_{\mathbf{i}}|=h}\left({u_{h}^{0}}\left(\mathbf{x}_{\mathbf{j}}\right)-{u_{h}^{0}}\left(\mathbf{x}_{\mathbf{i}}\right)\right)\frac{1}{|\mathbf{x }_{\mathbf{j}}-\mathbf{x}_{\mathbf{i}}|^{2}} = f(\mathbf{x}_{\mathbf{i}}) \quad \forall\; \mathbf{x}_{\mathbf{i}}\in\Omega, \end{align} (3.7) that is, it involves only the nearest neighbors along the axes. 3.2 Convergence analysis Due to the symmetry of the kernel we may rewrite the nonlocal operator as $$\mathcal{L}_{\delta} u(\mathbf{x})= \int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\rho_{\delta}(|\mathbf{z}|)\big(u(\mathbf{x}+\mathbf{z})+u(\mathbf{x}-\mathbf{z})-2u(\mathbf{x})\big)\, \mathrm{d}\mathbf{z}$$ and the discrete operator as $$\mathcal{L}_{h,\,\delta} u(\mathbf{x}_{\mathbf{i}})= \sum_{\mathbf{j}\neq\boldsymbol{0}}\frac{u\left(\mathbf{x}_{\mathbf{i}}+\mathbf{z}_{\mathbf{j}}\right)+u\left(\mathbf{x}_{\mathbf{i}}-\mathbf{z}_{\mathbf{j}}\right)-2u(\mathbf{x}_{\mathbf{i}})}{|\mathbf{z}_{\mathbf{j}}|^{2}}|\mathbf{z}_{\mathbf{j}}|_{1}\beta_{\delta}({\mathbf{z}_{\mathbf{j}}}) .$$ We first state a few technical lemmas concerning the nonlocal operators and the nonlocal discrete schemes that are nonlocal analogs of their well-known local versions. First, by the fact that W = W(z), the kernel ρδ, as well as the multilinear basis functions $$\phi_\mathbf{j}$$, is non-negative, we see that βδ, $$\beta_{\mathbf{j}\mathbf{i}}$$ and thus $$p^{\delta }_{\mathbf{x}_{\mathbf{j}}\mathbf{x}_{\mathbf{i}}}$$ are non-negative. This implies the M-matrix property of the coefficient matrix, just like the case corresponding to the central difference approximations of the diffusion operator. Lemma 3.3 (M-matrix). The coefficient matrix, also denoted by $$\mathcal{L}_{h,\delta }$$, corresponding to the linear system (3.4) is an M-matrix. An important stability property can be derived via the discrete maximum principle that follows from the M-matrix property of the coefficient matrix. By constructing a suitable quadratic barrier function similar to the one-dimensional case by Tao et al. (2017), as well as the local counterpart, we can get the following lemma. Lemma 3.4 (Uniform stability). For any δ > 0, $$\mathcal{L}_{h,\delta }^{-1}$$ is bounded. Furthermore, for any 0 < δ < δ0, we have \begin{align} \left\|\mathcal{L}_{h,\delta}^{-1} \right\|_{\infty}\le C(\delta_{0}) \end{align} (3.8) for some constant $$C\!\left(\delta _{0}\right )$$ independent of h and δ as h → 0. Next we consider the uniform consistency derived via the truncation error calculation. Here the uniform consistency means that the truncation error is independent of δ for small δ. We first state the following useful and vital lemma, which is a nonlocal analog of the fact that the central difference approximation to the Laplacian, as given by (3.7), is exact for quadratic polynomials. Lemma 3.5 (Quadratic exactness). For any quadratic polynomial in $$\mathbb{R}^{d}$$ given as u(x) = x ⊗ x : M, where M = (mkj) is a constant matrix, we have \begin{align} \mathcal{L}_{h,\,\delta} u(\mathbf{x}_{\mathbf{i}})=\mathcal{L}_{\delta} u(\mathbf{x}_{\mathbf{i}})=\sum_{k} m_{kk}\quad \forall\, \mathbf{i}. \end{align} (3.9) Proof. The fact that $$\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\mathbf{z}\otimes\mathbf{z}\rho_{\delta}(|\mathbf{z}| ) \,\mathrm{d}\mathbf{z}=I_{d},$$ where Id is the d × d identity matrix, yields $$\mathcal{L}_{\delta} u(\mathbf{x}_{\mathbf{i}})\equiv\sum_{k} m_{kk}.$$ Meanwhile, when u(x) = x ⊗ x : M, it is obvious that $$u\left(\mathbf{x}_{\mathbf{i}}+\mathbf{z}_{\mathbf{j}}\right)+u\left(\mathbf{x}_{\mathbf{i}}-\mathbf{z}_{\mathbf{j}}\right)-2u(\mathbf{x}_{\mathbf{i}})=\mathbf{z}_{\mathbf{j}}\otimes\mathbf{z}_{\mathbf{j}}:M.$$ Thus, $$\mathcal{L}_{h,\,\delta} u(\mathbf{x}_{\mathbf{i}}) = \sum_{\mathbf{j}\neq\boldsymbol{0}}\frac{\mathbf{z}_{\mathbf{j}}\otimes\mathbf{z}_{\mathbf{j}}:M}{\left|\mathbf{z}_{\mathbf{j}}\right|{}^{2}}\left|\mathbf{z}_{\mathbf{j}}\right|{}_{1}\beta_{\delta}\left(\mathbf{z}_{\mathbf{j}}\right)= Q:M .$$ We now need another observation on the coefficient matrix, that is, $$\beta _{\delta }\left (\mathbf{z}_{\mathbf{j}}\right )$$ is an even function; moreover, it depends only on the absolute values of the individual components of $$\mathbf{z}_\mathbf{j}$$. Thus, when doing the summation over j, it is easy to check that Qkl = 0 for $$k\neq l$$. We can further check that Qkk is independent of k, and \begin{eqnarray} \sum_{k} Q_{kk}&=&\int_{\mathcal{B}_{\delta}({\boldsymbol{0}})}\left(\sum_{\mathbf{j}\neq\boldsymbol{0}}\frac{|\mathbf{z}_{\mathbf{j}}|^{2}}{|\mathbf{z}_{\mathbf{j}}|^{2}}|\mathbf{z}_{\mathbf{j}}|_{1}\phi_{\mathbf{j}}(\mathbf{z})\right)\frac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}\rho_{\delta}(|\mathbf{z}|) \, \mathrm{d}\mathbf{z}\nonumber\\ &=&\int_{\mathcal{B}_{\delta}({\boldsymbol{0}})}|\mathbf{z}|^{2}\rho_{\delta}(|\mathbf{z}|) \,\mathrm{d}\mathbf{z}=d,\nonumber \end{eqnarray} where in the last step we have used the property that |z|1 is piecewise linear, and thus $$\sum_{\mathbf{j}\neq\boldsymbol{0}}|\mathbf{z}_{\mathbf{j}}|_{1}\phi_{\mathbf{j}}(\mathbf{z})=\sum_{\mathbf{j}}|\mathbf{z}_{\mathbf{j}}|_{1}\phi_{\mathbf{j}}(\mathbf{z})=|\mathbf{z}|_{1}.$$ Hence, we can derive that Q = Id, which gives $$\mathcal{L}_{\delta} u(\mathbf{x}_{\mathbf{i}})=\mathcal{L}_{h,\,\delta} u(\mathbf{x}_{\mathbf{i}})\equiv\sum_{k} m_{kk}.$$ This completes the proof. Before proceeding to the main theorem of this work we introduce some notation for simplicity of presentation. For any function $$u=u(\mathbf{x}): \mathbb{R}^{d}\mapsto \mathbb{R}$$ and $$\boldsymbol{\alpha }\in \mathbb{N}^{d}$$ let $$\mathbf{x}^{\boldsymbol{\alpha}}=\prod_{i=1}^{d}x_{i}^{\alpha_{i}},\quad \partial^{\boldsymbol{\alpha}}u(\mathbf{x})=\partial_{\alpha_{1}}\cdots\partial_{\alpha_{d}}u(\mathbf{x}),$$ and Dku(x) be the set containing all the kth-order derivatives in the form of ∂αu(x) with |α|1 = k, and $$\left|D^{k}u\right|{}_{\infty}=\max_{|\boldsymbol{\alpha}|_{1}=k}\sup_{\mathbf{x}}\partial^{\boldsymbol{\alpha}}u(\mathbf{x}).$$ We further use D2 to denote the Hessian and define $$H^{\mathbf{x}}(\mathbf{z})=u(\mathbf{x}+\mathbf{z})+u(\mathbf{x}-\mathbf{z})-2u(\mathbf{x}),\quad H^{\mathbf{x}}_{1}(\mathbf{z})=\mathbf{z}\otimes\mathbf{z}:D_{2}u(\mathbf{x})$$ and $$J^{\mathbf{x}}(\mathbf{z})=\frac{H^{\mathbf{x}}(\mathbf{z})}{|\mathbf{z}|^{2}}|\mathbf{z}|_{1},\quad J^{\mathbf{x}}_{1}(\mathbf{z})=\frac{H^{\mathbf{x}}_{1}(\mathbf{z})}{|\mathbf{z}|^{2}}|\mathbf{z}|_{1},$$ and let H2x(z) = Hx(z) − H1x(z) and J2x(z) = Jx(z) − J1x(z). By Taylor’s expansion and the mean value theorem we can show that \begin{align} \left\{ \begin{array}{rcl} \displaystyle H^{\mathbf{x}}_{2}(\mathbf{z})&=& \displaystyle \sum_{|\boldsymbol{\alpha}|_{1}=4}C_{\boldsymbol{\alpha}}^{1}\mathbf{z}^{\boldsymbol{\alpha}}\partial^{\boldsymbol{\alpha}}u(\mathbf{x}_{1}(\mathbf{z})),\\ \displaystyle \frac{\partial}{\partial z_{i}}H^{\mathbf{x}}_{2}(\mathbf{z})&=& \displaystyle\sum_{|\boldsymbol{\alpha}|_{1}=3}C_{\boldsymbol{\alpha}}^{2}\mathbf{z}^{\boldsymbol{\alpha}}\partial^{\boldsymbol{\alpha}+\mathbf{e}_{i}}u(\mathbf{x}_{2}(\mathbf{z})),\\ \displaystyle \frac{\partial^{2}}{\partial z_{i}\partial z_{j}}H^{\mathbf{x}}_{2}(\mathbf{z})&=& \displaystyle \sum_{|\boldsymbol{\alpha}|_{1}=2}C_{\boldsymbol{\alpha}}^{3}\mathbf{z}^{\boldsymbol{\alpha}}\partial^{\boldsymbol{\alpha}+\mathbf{e}_{i}+\mathbf{e}_{j}}u(\mathbf{x}_{3}(\mathbf{z})), \end{array} \right. \end{align} (3.10) where ei and ej are unit vectors with only one nonzero element at the ith and jth positions, respectively, and xk(z) are some points dependent on both x and z. The following is an elementary but technical result to be used later. Lemma 3.6 For J2x(z) defined above, if $$u\in C^{4}(\overline \Omega )$$, for any z and x, we have \begin{align} \big|D^{2}J^{\mathbf{x}}_{2}\big|_{\infty}\le C\big|D^{4}u\big|_{\infty}|\mathbf{z}|_{1}, \end{align} (3.11) where C is generic constant. Proof. To carry out the estimates we first note that $$\frac{\partial |\mathbf{z}|^{2}}{\partial z_{i}}=2z_{i}\quad\mbox{ and }\quad \frac{\partial |\mathbf{z}|^{2}}{\partial z_{i}\partial z_{j}}=2\delta_{ij} .$$ Now, by the definition of J2x(z), we have $$\frac{\partial^{2}}{\partial z_{i}\partial z_{j}}J^{\mathbf{x}}_{2}(\mathbf{z})=\textrm{sign}(z_{i})\frac{\partial}{\partial z_{j}}\frac{H^{\mathbf{x}}_{2}(\mathbf{z})}{|\mathbf{z}|^{2}}+\textrm{sign}(z_{j})\frac{\partial}{\partial z_{i}}\frac{H^{\mathbf{x}}_{2}(\mathbf{z})}{|\mathbf{z}|^{2}}+|\mathbf{z}|_{1}\frac{\partial^{2}}{\partial z_{i}\partial z_{j}}\frac{H^{\mathbf{x}}_{2}(\mathbf{z})}{|\mathbf{z}|^{2}},$$ where sign (⋅) is the sign function. For the first term in the above, we have from (3.10) that $$\left| \textrm{sign}(z_{i}) \frac{\partial}{\partial z_{j}}\frac{H^{\mathrm{x}}_{2}(\mathbf{z})}{|\mathbf{z}|^{2}}\right|\le C\left(\left|\frac{\mathbf{z}^{\boldsymbol{\beta}_{1}}}{|\mathbf{z}|^{2}}\right|+\left|\frac{\mathbf{z}^{\boldsymbol{\beta}_{2}}}{|\mathbf{z}|^{4}}\right|\right)|D^{4}u|_{\infty},$$ where $$\left |\boldsymbol{\beta }_{1}\right |_{1}=3$$ and $$\left |\boldsymbol{\beta }_{2}\right |_{1}=5$$. The second term satisfies a similar estimate. Moreover, by (3.10) again, we get $$\left| \frac{\partial^{2}}{\partial z_{i}\partial z_{j}}\frac{H^{\mathbf{x}}_{2}(\mathbf{z})}{{|\mathbf{z}|^{2}}}\right|\le C\left(\left|\frac{\mathbf{z}^{\boldsymbol{\alpha}_{1}}}{|\mathbf{z}|^{4}}\right|+\left|\frac{\mathbf{z}^{\boldsymbol{\alpha}_{2}}}{|\mathbf{z}|^{8}}\right|\right)|D^{4}u|_{\infty},$$ where |α1|1 = 4 and |α2|1 = 8. For any $$\boldsymbol{\alpha }\in \mathbb{N}^{d}$$ it is obvious that $$\left|\frac{\mathbf{z}^{\boldsymbol{\alpha}}}{|\mathbf{z}|^{|\boldsymbol{\alpha}|_{1}}}\right|\le 1,\quad \left|\frac{\mathbf{z}^{\boldsymbol{\alpha}}}{|\mathbf{z}|^{|\boldsymbol{\alpha}|_{1}-1}}\right|\le |\mathbf{z}|_{1}.$$ Combining the above estimates we can get the desired result (3.11). We now present the truncation error analysis for the quadrature-based difference approximation of the nonlocal model. Lemma 3.7 (Uniform nonlocal truncation error). Assume that $$u\in C^{4}(\overline{\Omega \cup \Omega _{\mathcal{I}}})$$. Then it holds that \begin{align} \max_{1\leq i\leq N_{s}} |\mathcal{L}_{h,\,\delta} u(\mathbf{x}_{\mathbf{i}})-\mathcal{L}_{\delta} u(\mathbf{x}_{\mathbf{i}})|\le C|D^{4}u|_{\infty} h^{2}\!, \end{align} (3.12) where C is a constant independent of δ and h. Proof. Without loss of generality we take xi = 0 to be the origin. Note that $$\mathcal{L}_{\delta}u(\boldsymbol{0}) =\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}J^{\boldsymbol{0}}(\mathbf{z})\frac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}\rho_{\delta}(|\mathbf{z}| ) \, \mathrm{d}\mathbf{z}$$ and $$\mathcal{L}_{h,\delta}u(\boldsymbol{0})=\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\mathcal{I}_{h} \left(J^{\boldsymbol{0}}(\mathbf{z})\right)\frac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}\rho_{\delta}(|\mathbf{z}| ) \, \mathrm{d}\mathbf{z} .$$ Simple subtraction gives us \begin{eqnarray*}\left|\mathcal{L}_{h,\delta}u(\boldsymbol{0})-\mathcal{L}_{\delta}u(\boldsymbol{0})\right|&&=\left|\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\left(\mathcal{I}_{h} \left(J^{\boldsymbol{0}}(\mathbf{z})\right)-J^{\boldsymbol{0}}(\mathbf{z})\right)\frac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}\rho_{\delta}(|\mathbf{z}| ) \, \mathrm{d}\mathbf{z}\right|,\\ && \le\left|\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\left(\mathcal{I}_{h} \left(J_{1}^{\boldsymbol{0}}(\mathbf{z})\right)-J_{1}^{\boldsymbol{0}}(\mathbf{z})\right)\frac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}\rho_{\delta}(|\mathbf{z}| ) \, \mathrm{d}\mathbf{z}\right|\\ &&\quad+\left|\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\left(\mathcal{I}_{h} \left(J_{2}^{\boldsymbol{0}}(\mathbf{z})\right)-J_{2}^{\boldsymbol{0}}(\mathbf{z})\right)\frac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}\rho_{\delta}(|\mathbf{z}| ) \,\mathrm{d}\mathbf{z}\right|\\ && :=E_{1}+E_{2}, \end{eqnarray*} where $$J_{1}^{\boldsymbol{0}}(\mathbf{z})=\frac{\mathbf{z}\otimes\mathbf{z}:D_{2}u(\boldsymbol{0})}{|\mathbf{z}|^{2}}|\mathbf{z}|_{1}\quad\mbox{ and }\quad J_{2}^{\boldsymbol{0}}(\mathbf{z})=J^{\boldsymbol{0}}(\mathbf{z})-J_{1}^{\boldsymbol{0}}(\mathbf{z}).$$ From Lemma 3.5, we know that E1 = 0. For E2, we have $$E_{2}\le Ch^{2}\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\left| D^{2} J_{2}^{\boldsymbol{0}}\right|_{\infty}\frac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}\rho_{\delta}(|\mathbf{z}| ) \,\mathrm{d}\mathbf{z}\,.$$ Thus, by using Lemma 3.6 and the moment condition on the kernel that $$\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}|\mathbf{z}|^{2}\rho_{\delta}(|\mathbf{z}| ) \,\mathrm{d}\mathbf{z}=d,$$ we complete the proof. Combining the stability result in Lemma 3.4 and the consistency result in Lemma 3.7, we immediately get the convergence of the quadrature-based finite difference scheme. Theorem 3.8 (Convergence to nonlocal solution). Assume that the nonlocal exact solution uδ of (1.1) is smooth enough, such as $$u^{\delta }\in C^{4}\left (\overline{\Omega \cup \Omega _{\mathcal{I}}}\right )$$, and uδ, h is the numerical solution obtained by the scheme (3.4). Then for any fixed δ there exists some constant Cδ independent of h as h → 0 such that \begin{align} \left\|{u^{\delta,h}}-I_{h}u^{\delta}\right\|_{\infty}\le Ch^{2}, \end{align} (3.13) where Ihuδ is the interpolation of the exact solution on the grid points. For any δ < δ0 the constant Cδ is uniformly bounded by some generic constant C(δ0). Next we show the asymptotic compatibility of the quadrature-based finite difference scheme to the local limiting model. We assume that the local equation (2.2) is discretized by the classical central finite difference scheme denoted by $$\mathcal{L}_{h,0}$$, whose truncation error is \begin{align} \max_{1\leq i\leq N_{s}} |\mathcal{L}_{h,\,0} u(\mathbf{x}_{\mathbf{i}})-\mathcal{L}_{0} u(\mathbf{x}_{\mathbf{i}})|=\mathcal{O}\big(h^{2}\big), \end{align} (3.14) by the standard numerical PDE analysis, if $$u\in C^{4}(\overline \Omega )$$. We now estimate the modeling error at the discrete level, that is, the error between the discrete nonlocal operator $$\mathcal{L}_{h,\delta }$$ and $$\mathcal{L}_{h,0}$$. Lemma 3.9 (Discrete model error). Assume that $$u\in C^{4}\left (\overline{\Omega \cup \Omega _{\mathcal{I}}}\right )$$. Then it holds that \begin{align} \max_{1\leq i\leq N_{s}} |\mathcal{L}_{h,\,\delta} u(\mathbf{x}_{\mathbf{i}})-\mathcal{L}_{h,0} u(\mathbf{x}_{\mathbf{i}})|=\mathcal{O}\left(\delta^{2}\right)+\mathcal{O}\left(h^{2}\right) \end{align} (3.15) when δ and h are sufficiently small. Proof. The triangle inequality gives \begin{eqnarray*} |\mathcal{L}_{h,\delta} u(\mathbf{x}_{\mathbf{i}})-\mathcal{L}_{h,0} u(\mathbf{x}_{\mathbf{i}})|&\le &|\mathcal{L}_{h,\,\delta} u(\mathbf{x}_{\mathbf{i}})-\mathcal{L}_{\delta} u(\mathbf{x}_{\mathbf{i}})|+ |\mathcal{L}_{0} u(\mathbf{x}_{\mathbf{i}})-\mathcal{L}_{\delta} u(\mathbf{x}_{\mathbf{i}})|\\ &&+ |\mathcal{L}_{0} u(\mathbf{x}_{\mathbf{i}})-\mathcal{L}_{h,\,0} u(\mathbf{x}_{\mathbf{i}})|\\ &:=& R_{1}+R_{2}+R_{3}. \end{eqnarray*} By Lemma 3.7 we have $$R_{1}= \mathcal{O}\big(h^{2}\big).$$ For $$u\in C^{4}(\overline \Omega )$$ the continuum property of nonlocal operators gives $$R_{2}=\mathcal{O}\big(\delta^{2}\big).$$ It is well known that $$R_{3}=\mathcal{O}\big(h^{2}\big).$$ The three equalities above together complete the proof. Combining the consistency in Lemma 3.9 with the stability Lemma 3.4 again, we get convergence (and asymptotic compatibility) to the local limit. Theorem 3.10 (Asymptotic compatibility). Assume that the local exact solution u0 of (2.2) is smooth enough, such as $$u^{0}\in C^{4}\left (\overline{\Omega \cup \Omega _{\mathcal{I}}}\right )$$. Then for any δ < δ0 there is some constant C independent of δ and h as h → 0 such that \begin{align} \left\|{u^{\delta,h}}-I_{h}u^{0}\right\|_{\infty}\le C\left(h^{2}+\delta^{2}\right)\!, \end{align} (3.16) where Ihu0 is the interpolation of the exact solution on the grid points. Remark 3.11 For problems on the same domain but with periodic boundary conditions (PBCs) all similar estimates also hold, which requires only that the exact solution of the local model is a $$C_{\mathrm{per}}^{4}(\overline{\Omega })$$ function. Remark 3.12 We now give a more explicit form of the algebraic system in two dimensions. We let $$u^{\delta ,h}_{j,k}$$ denote the nodal value of the nonlocal difference solution at the mesh point $$\left (x_{j}, y_{k}\right )\in \Omega =(-\pi ,\pi )^{2} \subset \mathbb{R}^{2}$$. Then let r = ⌈δ/h⌉ be the smallest integer larger than or equal to δ/h; we have $$\mathcal{L}_{h,\,\delta} u^{\delta,\,h}_{j,\,k} =\sum_{p=0}^{r}\sum_{q=0}^{r}c_{p,\,q}\left(u^{\delta,\,h}_{j+p,\,k+q}+u^{\delta,\,h}_{j-p,\,k+q}+ u^{\delta,\,h}_{j+p,\,k-q}+u^{\delta,\,h}_{j-p,\,k-q}-4u^{\delta,\,h}_{j,\,k}\right),$$ where c0, 0 = 0 and $$c_{p,\,q}= \frac{ph+qh}{p^{2}h^{2}+q^{2}h^{2}} \iint_{B_{+}(0,\, \delta)} \phi_{p,\,q}(s,t)\rho_{\delta}\left(\sqrt{s^{2}+t^{2}}\right)\frac{s^{2}+t^{2}}{s+t} \,\mathrm{d}s\, \mathrm{d}t\,,$$ with B+(0, δ) denoting the first quadrant of the disc at the origin with radius δ. One may also observe that cp, q = cq, p for any p and q. 4. Green’s functions of a nonlocal operator with integrable kernels By definition, the Green’s function of the nonlocal equation (corresponding to the particular nonlocal operator $$\mathcal{L}_{\delta }$$) solves the nonlocal equation with right-hand side given by the Dirac delta measure δy(x) = δ(x −y) of the variable x for a given y. Discussions of Green’s functions for nonlocal models such as PD have been presented in a number of earlier studies (Weckner et al., 2009; Silling et al., 2010b; Wang et al., 2016, 2017). Here we pay particular attention to the case where the nonlocal interaction kernel is given to be an integrable one. A consequence is that, unlike Green’s functions for local diffusion equations, the nonlocal analog is generically a measure function. Hence, we focus on effective constructions in this case. 4.1 Periodic case For the case of PBCs with Ω = [−π, π]d, we may impose the zero mean compatibility condition for both the right-hand side and the solution over Ω. Thus, we formally have the following modified equation of (1.1): $$-{{\mathcal{L}}_{\delta}} {G^{\delta}}_{\mathbf{y}}(\mathbf{x})={\boldsymbol{\delta}}_{\mathbf{y}}-\frac{1}{|\Omega|},$$ where |Ω| is the area of the domain. Taking the L2-inner product of the above equation with uδ(x) over x yields $$\left(-\mathcal{L}_{\delta} G^{\delta}_{\mathbf{y}},u^{\delta}\right)=u^{\delta}(\mathbf{y })-\frac{1}{|\Omega|}\int_{\Omega} u^{\delta}\, \mathrm{d}\mathbf{x} =u^{\delta}(\mathbf{y}).$$ By the nonlocal Green’s identity (Du et al., 2012, 2013), which remains valid for the periodic case, we have $$\left (-\mathcal{L}_{\delta } G^{\delta }_{\mathbf{y }},u^{\delta }\right )=\left (G^{\delta }_{\mathbf{y}},-\mathcal{L}_{\delta } u^{\delta }\right )$$. We then get $$u(\mathbf{y})=\left(f,G^{\delta}_{\mathbf{y}}\right),$$ which demonstrates how Green’s function $$G^{\delta }_{\mathbf{y}}$$ can be used to represent solutions to the nonlocal equation (1.1). Now, for the case where the kernel is integrable we denote $$C_{\delta}=\int_{\mathcal{B}_{\delta}(\bf{0})} \rho_{\delta}(|\mathbf{z}|) \, \mathrm{d}\mathbf{z}.$$ Then the nonlocal operator $$-\mathcal{L}_{\delta }$$ can be written as $$-\mathcal{L}_{\delta} u = -\rho_{\delta} \ast u + C_{\delta} u\,.$$ Let us perform the splitting $$G^{\delta }_{\mathbf{y}}=S_{1}+g_{0}$$ with S1 and g0 to be determined; we then arrive at $$-{{\mathcal{L}}_{\delta}} S_{1}(\mathbf{x})-{\rho_{\delta}} \ast g_{0} + C_{\delta} g_{0}={\boldsymbol{\delta}}_{\mathbf{y}}-\frac{1}{|\Omega|}.$$ Setting \begin{align} g_{0}=C_{\delta}^{-1}\left(\boldsymbol{\delta}_{\mathbf{y}}-\frac{1}{|\Omega|}\right) \quad \mbox{and}\quad \; g_{1}=C_{\delta}^{-1}\rho_{\delta} \ast g_{0} \end{align} (4.1) leads to \begin{align} -\mathcal{L}_{\delta} S_{1}(\mathbf{x})=\rho_{\delta} \ast g_{0} (\mathbf{x})= C_{\delta} g_{1}\, . \end{align} (4.2) Remark 4.1 Note that while g0 is a singular measure we have generically that g1 ∈ L1 for an integrable ρδ. If additional regularity of ρδ can be assumed then we also get g1 of the same regularity class as ρδ. For example, if ρδ is of the class $$C^{\infty }$$ (as a function in the whole space having compact support in the ball of radius δ) then we also have g1 and S1 in $$C^{\infty }$$. Another example, in the one-dimensional case, is that if ρδ is a bounded variation (BV) function then so are g1 and S1. In turn, if we perform a similar splitting once again as S1 = S2 + g1, we have $$-{{\mathcal{L}}_{\delta}} S_{2}(\mathbf{x})={\rho_{\delta}} \ast g_{1}.$$ Similarly, repeating this splitting technique n times we get recursively that $$g_{n}=C_{\delta}^{-1}\rho_{\delta} \ast g_{n-1},\quad -\mathcal{L}_{\delta} S_{n+1}(\mathbf{x })=\rho_{\delta} \ast g_{n},$$ and generically, the later terms obtained by the recursion are more regular than the preceding terms, similarly to the observation in Remark 4.1. In the end, we get $$G^{\delta}_{\mathbf{y}}=\sum_{n=0}^{\infty}g_{n}$$ with the following recursive relation: $$g_{0}=C_{\delta}^{-1}\left(\boldsymbol{\delta}_{\mathbf{y}}-\frac{1}{|\Omega|}\right)\!,\quad g_{n}=C_{\delta}^{-1}\rho_{\delta} \ast g_{n-1},\; n=1,\;2,\;\dots .$$ Thus, we obtain an explicit series form of $$G^{\delta }_{\mathbf{y}}(\mathbf{x})$$, which is given below. Theorem 4.2 For an integrable kernel ρδ the nonlocal Green’s function $$G^{\delta }_{\mathbf{y}}$$ for equation (1.1) with PBCs can be represented by \begin{align} G^{\delta}_{\mathbf{y}}=\sum_{n=0}^{\infty} C_{\delta}^{-(n+1)} \left(\rho_{\delta}\ast\right)^{n}\ast \left(\boldsymbol{\delta}_{\mathbf{y}}-\frac{1}{|\Omega|}\right). \end{align} (4.4) Moreover, the solution uδ of the nonlocal diffusion equation (1.1) is given by $$u^{\delta}(\mathbf{y})=\left(f,G^{\delta}_{\mathbf{y}}\right).$$ Remark 4.3 Due to the lack of elliptic smoothing of the nonlocal diffusion equation corresponding to the integrable kernels, the nonlocal Green’s function remains a composition of a singular part in the form of a Dirac delta measure and an integrable part. In this case and unlike local Green’s functions for elliptic PDEs, $$G^{\delta }_{\mathbf{y}}$$ should not be viewed as a function defined almost everywhere. Moreover, the regularity of the solution of uδ is in general the same as the right-hand side f, since the regularity $$\big (\,f,G^{\delta }_{\mathbf{y}}\big )$$ is dominated by the first term (g0, f), which is proportional to f. One could also split the solution into two parts, one given by (g0, f) while the other accounts for a more regular part. Remark 4.4 It is worth mentioning that the series expression for the Green’s function is a Neumann series expansion for the equation (I − A)G = f where $$Au = C_{\delta }^{-1}\rho _{\delta }\ast u$$ and $$f = C_{\delta }^{-1}\left (\boldsymbol{\delta }_{\mathbf{y}}-\frac{1}{|\Omega |}\right )$$. 4.2 Other cases of nonlocal constraints The approach of constructing Green’s function through a series expansion is also applicable to other nonlocal constraints. For instance, we consider the Dirichlet-type volume-constrained problems given by (1.1) and (2.1). For any y ∈ Ω let us define g0 = δ(x −y), $$\mathbf{x }\in \Omega \cup \Omega _{\mathcal{I}}$$. It is easy to see that, with gn recursively defined by $$\begin{cases} g^{n}(\mathbf{x}) =C_{\delta}^{-1}\rho_{\delta}\ast g^{n-1}(\mathbf{x}), & \mathbf{x}\in\Omega,\\ g^{n}(\mathbf{x}) =0, & \mathbf{x}\in\Omega_{\mathcal{I}},\end{cases}$$ the nonlocal Green’s function is given by $$G^{\delta}_{\mathbf{y}}(\mathbf{x})=\sum_{n=0}^{\infty} g^{n}(\mathbf{x}).$$ As in the case with PBC, we can also use the above expansion to get a series expression for a general solution, though the regularity pickup in the later terms of the series gets more delicate due to the dependence on matching the boundary data with the right-hand side and remains an interesting topic to be further studied. Remark 4.5 Other types of series expansions, in particular Fourier series expansion, of the nonlocal Green’s functions have been presented in the literature by Weckner et al. (2009), Silling et al. (2010b) and Wang et al. (2016, 2017). Such expansions are valid formally, but for integrable kernels the series do not converge pointwisely. 4.3 Discrete nonlocal Green’s function Following the discussion above on the continuum level we can present a discrete analog with the AC quadrature-based finite difference scheme. Given the uniform mesh $$\left \{\mathbf{x}_{\mathbf{j}}\right \}$$ the discrete Green’s function (matrix) $$G^{h,\delta }_{\mathbf{x}_{\mathbf{i}}}\left (\mathbf{x}_{\mathbf{j}}\right )$$ is given by $$-{{\mathcal{L}}_{h,\,\delta}} {G^{h,\,\delta}}_{{\mathbf{x}}_{\mathbf{i}}}\left(\mathbf{x }_{\mathbf{j}}\right)=\frac{1}{h^{d}} {\delta_{\mathbf{i}\mathbf{j}}}-\frac{1}{|\Omega|},$$ where δij is the Kronecker delta function, that is, δii = 1 and δij = 0 for $$\mathbf{i}\neq \mathbf{j}$$. It follows that the discrete solution to the nonlocal model with a discrete right-hand side $$\left \{f_{\mathbf{j}}\right \}$$ is given by $$u^{h,\,\delta}(\mathbf{x}_{\mathbf{i}})=\sum_{\mathbf{j}\neq\mathbf{i} } f_{\mathbf{j}} G^{h,\delta}_{\mathbf{x}_{\mathbf{i}}}\left(\mathbf{x}_{\mathbf{j}}\right)\!\quad\mbox{ for }\quad \sum_{\mathbf{j}} f_{\mathbf{j}} = 0.$$ The study of discrete nonlocal Green’s functions can be seen as a nonlocal analog of similar works on discrete Green’s functions for local PDEs (Chung & Yau, 2000). Although the nonlocal Green’s function at the continuum level is measure valued for integrable interaction kernels the discrete nonlocal Green’s function given above is well defined at all the grid points for any finite δ and h. Due to the lack of regularity, however, we do not expect strong (or pointwise) convergence of $$G^{h,\delta }_{\mathbf{y}}(\mathbf{x})$$ to $$G^{\delta }_{\mathbf{y}}(\mathbf{x})$$ as h → 0 for a given δ. Instead, we could attempt to get better approximations by the splitting technique, that is, we numerically solve for the regular part while invoking the analytical solution of singular part. Fig. 1. View largeDownload slide The function $$u^{h,\delta }_{1}$$ of the regular part of the nonlocal Green’s function. Fig. 1. View largeDownload slide The function $$u^{h,\delta }_{1}$$ of the regular part of the nonlocal Green’s function. Therefore, we obtain a hybrid approximation of $$G^{\delta }_{\mathbf{x}_{\mathbf{j}}}(\mathbf{x }_{\mathbf{i}})$$ by $$G^{h,\,\delta}_{\mathbf{x}_{\mathbf{j}}}(\mathbf{x}_{\mathbf{i}}) = C_{\delta}^{-1}\left(\boldsymbol{\delta}_{\mathbf{x}_{\mathbf{j}}}(\mathbf{x }_{\mathbf{i}})-\frac{1}{|\Omega|}\right)+u^{h,\,\delta}_{1},$$ where the approximation $$u^{h,\delta }_{1}$$ of the regular part of the nonlocal Green’s function S1 in (4.2) satisfies $$-{{\mathcal{L}}_{h,\,\delta}} {u^{h,\,\delta}}_{1} (\mathbf{x}_{\mathbf{i}})= {\rho_{\delta}} \ast g_{0}(\mathbf{x}_{\mathbf{i}}).$$ Given the higher regularity of S1 we expect to have good convergence of $$u^{h,\delta }_{1}$$ to S1 for any given δ, which is illustrated in Fig. 1 for a kernel that is a positive constant over the horizon and vanishes outside (thus piecewise continuous). Moreover, due to the AC property of the discrete scheme, we also have convergence of $$G^{h,\delta }_{\mathbf{x}_{\mathbf{j}}}(\mathbf{x}_{\mathbf{i}})$$ to the local Green’s function as both h and δ approach zero. Notice that both Theorems 3.8 and 3.10 require that the solutions to be approximated are $$C^{4}\left (\overline{\Omega \cup \Omega _{\mathcal{I}}}\right )$$. Thus, if we subtract enough terms of the Green’s function series expansion from the right-hand side, in principle we should be left with a smooth enough right-hand side that the results of both theorems are applicable. More discussions on this splitting technique are to be presented later. We thus can expect a robust approximation of the local and nonlocal Green’s functions. 4.4 More discussion on the splitting technique We next discuss a potential application of the splitting technique presented above in the numerical solution of nonlocal equations. By avoiding using derivatives, nonlocal models allow more singular solutions, which is a distinct advantage in physical modeling such as the study of material cracks. As an illustration we present a typical example showing how the splitting technique improves the performance of Fourier spectral methods for numerically solving nonlocal equations with discontinuous solutions. For simplicity, we take a one-dimensional nonlocal diffusion model (1.1) over a periodic cell Ω = [0, 2π] with a constant kernel $$\rho _{\delta }(s)=\frac{3}{\delta ^{3}}\chi _{[-\delta ,\delta ]}$$, which is of the rescaled form given in (2.3), and a discontinuous right-hand side $$f(x)=\textrm{sign}(x-\pi), \quad x\in [0,2\pi],$$ where sign(⋅) is the sign function. In this setting, the exact solution has a jump at x = π, e.g., like studied in Du & Yang (2016). One numerical solution is obtained by applying the Fourier spectral method directly. It is not surprising that the well-known Gibbs phenomena occur around x = π when standard Fourier spectral methods are directly used to approximate discontinuous solutions. On the other hand, if we use the splitting technique described above once, by first constructing the discontinuous part of the solution analytically, and then solving a new nonlocal diffusion equation with a smoother right-hand side by Fourier spectral methods, we can eliminate the Gibbs phenomena, as demonstrated in Fig. 2. Fig. 2. View largeDownload slide Numerical comparison. Left: numerical solutions; right: zoom in around $$\pi\!.$$ Fig. 2. View largeDownload slide Numerical comparison. Left: numerical solutions; right: zoom in around $$\pi\!.$$ 5. Conclusion In this work, a quadrature-based finite difference scheme for a nonlocal diffusion model in multidimensions is proposed. Our analysis is given for problems with Dirichlet volumetric constraints. The case of a Neumann-type nonlocal volumetric constraint may involve additional complications, due to the need for ghost points, which require further attention to ensure second-order accuracy. One may adapt the scheme developed here to more general models. In fact, an extension based on the idea presented in a draft version of this work has been given for nonlocal convection diffusion equations in the study by Tian et al. (2018) (where only a partial analysis of the diffusion term has been documented and one should refer to the more complete analysis presented in the current work). In deriving the error analysis our attention is focused on the case where the underlying solutions are smooth. Since nonlocal models are effective in describing physical processes with solution singularities, extending the analysis to cases with minimal regularity will be of interest. One can of course also consider extensions to systems involving vector fields such as the bond-based and state-based PD models. Our study here also illustrated an important difference between the nonlocal Green’s function and their local analog as the former may take on a possibly singular measure form. Nonlocal Green’s functions can be useful in various applications such as the mean exit time computation of jump processes and analysis of mechanical responses to point loads in nonlocal mechanics (Weckner et al., 2009; Silling et al., 2010b; Du et al., 2012; Wang et al., 2016, 2017). The splitting technique presented here offers an effective algorithm for computing singular Green’s functions. Moreover, one may further explore the applications of such techniques to the computation of singular solutions to more general nonlocal models as alluded to in Section 4.4. In comparison with the Galerkin finite element AC scheme developed in the study by Tian & Du (2014) that works for linear systems in multidimensions using unstructured meshes, we see that the quadrature-based finite difference AC schemes given here are developed only for uniform Cartesian meshes with the same mesh size in all directions. However, the schemes given here are simpler to implement as they avoid the evaluation of integrals in $$\mathbb{R}^{2d}$$. Furthermore, they lead to M-matrices and the discrete maximum principle that are in general not the case for the finite element discretization. Hence, an interesting future work is to develop hybrid schemes and particle discretization (Liu et al., 1996; Silling, 2000; Askari et al., 2008; Bessa et al., 2014) that can preserve the nice properties of the quadrature-based finite difference AC schemes while applicable to more general meshes. Funding United States National Science Foundation (DMS-1719699), Air Force Office of Scientific Research MURI center and the Army Research Office MURI Grant (W911NF-15-1-0562). References Andreu , F. , Mazón , J. M. , Rossi , J. D. & Toledo , J. ( 2010 ) Nonlocal Diffusion Problems . Mathematical Surveys and Monographs , vol. 165 . Providence, RI : AMS. Applebaum , D. ( 2004 ) Lévy Processes and Stochastic Calculus . Cambridge Studies in Advanced Mathematics , vol. 93 . Cambridge, UK : Cambridge University Press . Askari , E. , Bobaru , F. , Lehoucq , R. B. , Parks , M. L. , Silling , S. A. & Weckner , O. ( 2008 ) Peridynamics for multiscale materials modeling . J. Phys. Conf. Ser. , 125 , 12 – 78 . Google Scholar CrossRef Search ADS Bates , P. W. & Chmaj , A. ( 1999 ) An integrodifferential model for phase transitions: stationary solutions in higher space dimensions . J. Stat. Phys. , 95 , 1119 – 1139 . Google Scholar CrossRef Search ADS Bessa , M. , Foster , J. , Belytschko , T. & Liu , W. K. ( 2014 ) A meshfree unification: reproducing kernel peridynamics . Comput. Mech. , 53 , 1251 – 1264 . Google Scholar CrossRef Search ADS Bobaru , F. & Duangpanya , M. ( 2010 ) The peridynamic formulation for transient heat conduction . Int. J. Heat Mass Transf. , 53 , 4047 – 4059 . Google Scholar CrossRef Search ADS Bobaru , F. , Yang , M. , Alves , L. F. , Silling , S. A. , Askari , E. & Xu , J. ( 2009 ) Convergence, adaptive refinement, and scaling in 1D peridynamics . Int. J. Numer. Methods Eng. , 77 , 852 – 877 . Google Scholar CrossRef Search ADS Buades , A. , Coll , B. & Morel , J. M. ( 2010 ) Image denoising methods. A new nonlocal principle . SIAM Rev. , 52 , 113 – 147 . Google Scholar CrossRef Search ADS Chen , A. , Du , Q. , Li , C. & Zhou , Z. ( 2017 ) Asymptotically compatible schemes for space–time nonlocal diffusion equations . Chaos Solitons Fractals (in press) . Google Scholar CrossRef Search ADS Chen , X. & Gunzburger , M. ( 2011 ) Continuous and discontinuous finite element methods for a peridynamics model of mechanics . Comput. Methods Appl. Mech. Eng. , 200 , 1237 – 1250 . Google Scholar CrossRef Search ADS Chung , F. & Yau , S-T. ( 2000 ) Discrete Green’s functions . J. Comb. Theory A , 91 , 191 – 214 . Google Scholar CrossRef Search ADS Du , Q. , Gunzburger , M. , Lehoucq , R. B. & Zhou , K. ( 2012 ) Analysis and approximation of nonlocal diffusion problems with volume constraints . SIAM Rev. , 56 , 676 – 696 . Du , Q. , Gunzburger , M. , Lehoucq , R. B. & Zhou , K. ( 2013 ) A nonlocal vector calculus, nonlocal volume-constrained problems, and nonlocal balance laws . Math. Models Methods Appl. Sci. , 23 , 493 – 540 . Google Scholar CrossRef Search ADS Du , Q. , Ju , L. & Lu , J. ( 2017a ) A discontinuous Galerkin method for one-dimensional time-dependent nonlocal diffusion problems , Preprint . Du , Q. , Tao , Y. , Tian , X. & Yang , J. ( 2016 ) Robust a posteriori stress analysis for quadrature collocation approximations of nonlocal models via nonlocal gradients . Comput. Methods Appl. Mech. Eng. , 310 , 605 – 627 . Google Scholar CrossRef Search ADS Du , Q. & Tian , X. ( 2014 ) Asymptotically compatible schemes for peridynamics based on numerical quadratures . ASME 2014 International Mechanical Engineering Congress and Exposition (IMECE2014), vol. 1 , Paper No. IMECE2014–39620, pp . V001T01A058 . Du , Q. & Yang , J. ( 2016 ) Asymptotic compatible Fourier spectral approximations of nonlocal Allen–Cahn equations . SIAM J. Numer. Anal. , 54 , 1899 – 1919 . Google Scholar CrossRef Search ADS Du , Q. & Yang , J. ( 2017 ) Fast and accurate implementation of Fourier spectral approximations of nonlocal diffusion operators and its applications . J. Comput. Phys. , 332 , 118 – 134 . Google Scholar CrossRef Search ADS Du , Q. , Yang , J. & Zhou , Z. ( 2017b ) Analysis of a nonlocal-in-time parabolic equation . Discrete Continuous Dyn. Syst. B. , 22 , 339 – 368 . Google Scholar CrossRef Search ADS Du , Q. & Zhou , K. ( 2011 ) Mathematical analysis for the peridynamic nonlocal continuum theory . Math. Model. Numer. Anal. , 45 , 217 – 234 . Google Scholar CrossRef Search ADS Gilboa , G. & Osher , S. ( 2008 ) Nonlocal operators with applications to image processing . Multiscale Model. Simul. , 7 , 1005 – 1028 . Google Scholar CrossRef Search ADS Kilic , B. & Madenci , E. ( 2010 ) Coupling of peridynamic theory and the finite element method . J. Mech. Mater. Struct. , 5 , 707 – 733 . Google Scholar CrossRef Search ADS Liu , W. K. , Chen , Y. , Jun , S. , Chen , J. , Belytschko , T. , Pan , C. , Uras , R. & Chang , C. ( 1996 ) Overview and applications of the reproducing kernel particle methods . Arch. Comput. Methods Eng. , 3 , 3 – 80 . Google Scholar CrossRef Search ADS Lou , Y. , Zhang , X. , Osher , S. & Bertozzi , A. ( 2010 ) Image recovery via nonlocal operators . J. Sci. Comput. , 42 , 185 – 197 . Google Scholar CrossRef Search ADS Macek , R. & Silling , S. A. ( 2007 ) Peridynamics via finite element analysis . Finite Elem. Anal. Des. , 43 , 1169 – 1178 . Google Scholar CrossRef Search ADS Mengesha , T. & Du , Q. , ( 2014 ) Nonlocal constrained value problems for a linear peridynamic Navier equation . J. Elast. , 116 , 27 – 51 . Google Scholar CrossRef Search ADS Palatucci , G. , Savin , O. & Valdinoci , E. ( 2012 ) Local and global minimizers for a variational energy involving a fractional norm . Ann. Mat. Pura Appl. , 192 , 673 – 718 . Google Scholar CrossRef Search ADS Silling , S. A. ( 2000 ) Reformulation of elasticity theory for discontinuities and long-range forces . J. Mech. Phys. Solids , 48 , 175 – 209 . Google Scholar CrossRef Search ADS Silling , S. A. & Askari , E. ( 2005 ) A meshfree method based on the peridynamic model of solid mechanics . Comput. Struct. , 83 , 1526 – 1535 . Google Scholar CrossRef Search ADS Silling , S. A. & Lehoucq , R. B. ( 2010 ) Peridynamic theory of solid mechanics . Adv. Appl. Mech. , 44 , 73 – 168 . Google Scholar CrossRef Search ADS Silling , S. A. , Weckner , O. , Askari , E. & Bobaru , F. ( 2010a ) Crack nucleation in a peridynamic solid . Int. J. Fract. , 162 , 219 – 227 . Google Scholar CrossRef Search ADS Silling , S. A. , Zimmermann , M. & Abeyaratne , R. ( 2010b ) Deformation of a peridynamic bar . J. Elast. , 73 , 173 – 190 . Google Scholar CrossRef Search ADS Tadmor , E. & Tan , C. ( 2014 ) Critical thresholds in flocking hydrodynamics with non-local alignment . Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. , 372 : 20130401 . Tao , Y. , Tian , X. & Du , Q. ( 2017 ) Nonlocal diffusion and peridynamic models with Neumann type constraints and their numerical approximations . Appl. Math. Comput. , 305 , 282 – 298 . Tian , H. , Ju , L. & Du , Q. ( 2018 ) A conservative nonlocal convection-diffusion model and asymptotically compatible finite difference discretization . Comput. Methods Appl. Mech. Eng. , 320 , 46 – 67 . Google Scholar CrossRef Search ADS Tian , X. & Du , Q. ( 2013 ) Analysis and comparison of different approximations to nonlocal diffusion and linear peridynamic equations . SIAM J. Numer. Anal. , 51 , 3458 – 3482 . Google Scholar CrossRef Search ADS Tian , X. & Du . Q. ( 2014 ) Asymptotically compatible schemes and applications to robust discretization of nonlocal models . SIAM J. Numer. Anal. , 52 , 1641 – 1665 . Google Scholar CrossRef Search ADS Wang , L. , Xu , J. & Wang , J. ( 2016 ) The Green’s functions for peridynamic non-local diffusion . Proc. R. Soc. A , 472 , 20160185 . Wang , L. , Xu , J. & Wang , J. ( 2017 ) Static and dynamic Green’s functions in peridynamics . J. Elast. , 126 , 95 – 125 . Google Scholar CrossRef Search ADS Weckner , O. , Brunk , G. , Epton , M. , Silling , S. & Askari , E. ( 2009 ) Green’s functions in non-local three-dimensional linear elasticity. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. , 53 , 705 – 728 . Zhou , K. & Du , Q. ( 2010 ) Mathematical and numerical analysis of linear peridynamic models with nonlocal boundary conditions . SIAM J. Numer. Anal. , 48 , 1759 – 1780 . Google Scholar CrossRef Search ADS © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) For permissions, please e-mail: journals. permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# Asymptotically compatible discretization of multidimensional nonlocal diffusion models and approximation of nonlocal Green’s functions

, Volume Advance Article – Apr 5, 2018
19 pages

/lp/ou_press/asymptotically-compatible-discretization-of-multidimensional-nonlocal-PnitjMMELr
Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/dry011
Publisher site
See Article on Publisher Site

### Abstract

Abstract Nonlocal diffusion equations and their numerical approximations have attracted much attention in the literature as nonlocal modeling becomes popular in various applications. This paper continues the study of robust discretization schemes for the numerical solution of nonlocal models. In particular, we present quadrature-based finite difference approximations of some linear nonlocal diffusion equations in multidimensions. These approximations are able to preserve various nice properties of the nonlocal continuum models such as the maximum principle and they are shown to be asymptotically compatible in the sense that as the nonlocality vanishes, the numerical solutions can give consistent local limits. The approximation errors are proved to be of optimal order in both nonlocal and asymptotically local settings. The numerical schemes involve a unique design of quadrature weights that reflect the multidimensional nature and require technical estimates on nonconventional divided differences for their numerical analysis. We also study numerical approximations of nonlocal Green’s functions associated with nonlocal models. Unlike their local counterparts, nonlocal Green’s functions might become singular measures that are not well defined pointwise. We demonstrate how to combine a splitting technique with the asymptotically compatible schemes to provide effective numerical approximations of these singular measures. 1. Introduction In this paper, we study numerical approximations of the following multidimensional linear nonlocal diffusion equation: \begin{align} - \mathcal{L}_{\delta}u^{\delta} (\mathbf{x})=f(\mathbf{x}),\quad \mathbf{x}\in\Omega\subset \mathbb{R}^{d}, \end{align} (1.1) with Ω being a given domain in $$\mathbb{R}^{d}$$ (d ≥ 2), f a given right-hand side, uδ the solution to be sought and the nonlocal operator $$\mathcal{L}_{\delta }$$ defined by \begin{align} \mathcal{L}_{\delta} u(\mathbf{x})= 2 \int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\rho_{\delta}(|\mathbf{z}|)\left(u(\mathbf{x}+\mathbf{z})-u(\mathbf{x})\right)\, \mathrm{d}\mathbf{z}. \end{align} (1.2) The operator $$\mathcal{L}_{\delta }$$ is parametrized by a positive horizon parameter δ measuring the range of nonlocal interactions. The specific form of such nonlocal interactions is prescribed by a non-negative kernel function ρδ = ρδ(|z|). The integral in (1.2) is interpreted in the principal value sense whenever needed (Mengesha & Du, 2014). More discussions on the nonlocal model are given in Section 2. A main task of this work is to develop quadrature-based finite difference approximation schemes and associated Green’s functions for (1.1). The approximation schemes are able to preserve the maximum principle at the discrete level and they are shown to be asymptotically compatible in the sense that as the nonlocality vanishes, the numerical solutions can give consistent local limits. They are also shown to be of optimal order in both nonlocal and asymptotically local settings. Previous studies of schemes enjoying such optimal-order error estimates, discrete maximum principles and asymptotical compatibility have been mostly confined to the one-dimensional space (Tian & Du, 2013). While the algorithmic development is in a similar spirit to that for the one-dimensional case, the multidimensional extension made in this work involves both new design elements that reflect the multidimensional nature and more technical convergence analysis relying on new interpolation error estimates of nonconventionally defined weighted divided differences. The motivation for our work is rooted in the need to develop robust approximations of nonlocal models that serve as alternatives to classical partial differential equations (PDEs). Recently, there have been many studies on the application of nonlocal modeling to problems in mechanics, physics and materials, biological and social sciences (Bates & Chmaj, 1999; Applebaum, 2004; Gilboa & Osher, 2008; Bobaru & Duangpanya, 2010; Buades et al., 2010; Lou et al., 2010; Tadmor & Tan, 2014). The integral formulations of spatial interactions in nonlocal models can account for nonlocal effects and allow more singular solutions. An example is the nonlocal peridynamics (PD) theory (Silling, 2000) and its applications in studying cracks and materials failure, as well as other mechanical properties and physical processes (Askari et al., 2008; Silling & Lehoucq, 2010; Silling et al., 2010a; Palatucci et al., 2012). Rigorous mathematics of nonlocal models has proved to be necessary in order to gain fundamental insights and to guide the modeling and simulation efforts (Silling & Askari, 2005; Macek & Silling, 2007; Bobaru et al., 2009; Andreu et al. 2010; Kilic & Madenci, 2010; Zhou & Du, 2010; Chen & Gunzburger, 2011; Du & Zhou, 2011; Du et al., 2012, 2013, 2017b). For nonlocal models like (1.1), as $$\delta \rightarrow 0$$, the nonlocal interactions that define the model can become localized so that the zero-horizon (or local) limit of the nonlocal operator, when valid both physically and mathematically, can be represented by a local differential operator. Given such a scenario the corresponding nonlocal model naturally converges to a conventional differential equation model in the local limit; see for instance, the convergence analysis of linear state-based PD models to the classical Navier equation of linear elasticity (Mengesha & Du, 2014). There have been concerns about whether such consistency can be preserved at the discrete level when the nonlocal models are discretized numerically. Discrete schemes of nonlocal models that preserve the correct limiting behavior are called asymptotically compatible (AC) schemes, a notion developed first in the studies by Tian & Du (2013, 2014). In other words, the numerical approximation given by an AC scheme can reproduce the correct local limiting solution as the horizon δ and the mesh spacing (denoted by h) approach zero, if the convergence to such a limiting local solution is valid on the continuum level. While Tian & Du (2013) presented results concerning the AC property for a one-dimensional scalar model solved by a number of different discretization methods, Tian & Du (2014) managed to develop more general results that are applicable to linear systems in multidimensions based on Galerkin finite element approximations on unstructured meshes. A few subsequent studies have been carried out. For example, Tao et al. (2017) extended the study of AC schemes to nonlocal diffusion equations with Neumann-type volume constraints, Du & Yang (2016) discussed AC schemes based on Fourier spectral methods for problems defined on periodic cells and Du & Yang (2017) proposed an efficient and accurate hybrid algorithm to implement the Fourier spectral methods, while Du et al. (2017b) developed discontinuous-Galerkin-based AC schemes. Moreover, AC schemes were proposed in the study by Chen et al. (2017) for nonlocal time–space models and in the study by Du et al. (2016) for the robust recovery of the nonlocal gradient of the solutions to nonlocal models. For technical reasons, much of the numerical analysis in these works, with the exception of the study by Tian & Du (2014) on conforming Galerkin finite element approximations, has focused on nonlocal operators defined for a one-dimensional spatial variable. In comparison with other AC schemes such as those based on the Galerkin finite element and spectral methods, quadrature-based collocation-type finite difference schemes not only offer optimal-order convergent approximations with simpler implementations but also preserve the discrete maximum principle and asymptotic compatibility; see details in Section 3. Furthermore, these quadrature-based schemes are closely related to other discretizations based on strong forms, such as the original mesh-free discretization developed in the study by Silling (2000). As an application of the discretization scheme, we consider in Section 4 the numerical computation of nonlocal Green’s functions Gδ = Gδ(x, y) for the nonlocal diffusion model (1.1). Nonlocal Green’s functions are important tools for nonlocal continuum models (Weckner et al., 2009; Silling et al., 2010b; Wang et al., 2016, 2017), just like their local versions, i.e., conventional Green’s functions associated with PDEs. Since analytical expressions can be expected only in simple cases, the numerical study of nonlocal Green’s functions has particular significance. We show here that while the local/conventional and nonlocal Green’s functions enjoy many similarities and are intimately connected through the local limiting process, there are also some fundamental differences among them. In particular, for a large class of nonlocal interaction kernels the associated Green’s functions Gδ = Gδ(x, y) on the continuum level for δ > 0 are measure-valued distributions, which cannot be interpreted as conventional functions defined pointwise for x and y. These differences are essential features that, if improperly handled, could present undesirable effects on their numerical approximations. This serves as another reminder of the subtlety involved in the nonlocal modeling. To better deal with these complications we discuss a splitting approach that utilizes the developed AC scheme to get effective and robust approximations of nonlocal Green’s functions. 2. Multidimensional nonlocal diffusion model For the sake of a clear illustration we let u = u(x) denote a scalar deformation field and consider (1.1) subject to a nonlocal constraint of the Dirichlet type, \begin{align}\ u=0 \;\; \textrm{on} \;\Omega_{\mathcal{I}}, \end{align} (2.1) where $$\Omega _{\mathcal{I}}$$ is a boundary layer defined as $$\Omega_{\mathcal{I}}=\{\mathbf{x}\notin \Omega|\;\textrm{dist}(\mathbf{x},\Omega)\le \delta\}.$$ For more discussions on nonlocal constraints defined on a domain with a nonzero volume we refer to Du et al. (2012, 2013). Furthermore, let us consider suitably scaled kernels so that as δ goes to zero, the local limit of the nonlocal operator $$\mathcal{L}_{\delta }$$ is exactly the Laplace operator Δ (which is denoted by $$\mathcal{L}_{0}$$). To this end, it requires that the kernel is a non-negative and nonincreasing function and has a finite second moment in z independent of the horizon parameter δ: $$\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}|\mathbf{z}|^{2}\rho_{\delta}(|\mathbf{z}| ) \, \mathrm{d}\mathbf{z}=d.$$ The limiting local model of (1.1) is then given by \begin{align} - \mathcal{L}_{0}u^{0}=f\;\; \textrm{in } \Omega\quad \mbox{ and }\quad u^{0}=0 \;\;\mbox{ on }\; \partial \Omega, \end{align} (2.2) where, without loss of generality, $$\mathcal{L}_{0}=\Delta$$ for the particular case studied here. We note that a popular choice of ρδ is a rescaled kernel given by \begin{align} \rho_{\delta}(|\mathbf{z}|)=\frac{1}{\delta^{2+d}}\rho\left(\frac{|\mathbf{z}|}{\delta}\right)\!\quad \forall\, \mathbf{z}\in \mathcal{B}_{\delta}(\boldsymbol{0}), \end{align} (2.3) where ρ = ρ(ξ) is a non-negative and nonincreasing function with compact support in [0, 1] and a normalized moment \begin{align} {\int_{0}^{1}}\rho(\xi)\xi^{1+d}\, \mathrm{d}\xi=c_{d} \end{align} (2.4) for a given constant cd, though our discussion here is not restricted to such a rescaled form. 3. Quadrature-based finite difference schemes and asymptotic compatibility AC schemes provide robust numerical approximations of nonlocal models (Tian & Du, 2013, 2014) since the convergence of such schemes is insensitive to the choices of modeling and discretization parameters, in particular with respect to either a sufficiently small horizon δ or a sufficiently refined mesh spacing (denoted by h). 3.1 The discretization scheme We now develop an AC quadrature-based finite difference discretization to the nonlocal diffusion equation (1.1) that extends a similar scheme presented first in the study by Tian & Du (2014) and also studied by Du & Tian (2014) for the one-dimensional case. For simplicity, we let $$\left \{\mathbf{x}_{\mathbf{j}}\in \Omega \cup \Omega _{\mathcal{I}}\right \}$$ be the set of nodes (grid points) of a uniform Cartesian mesh $${\mathcal{T}}_{h}$$ with mesh size h. Here, j denotes a multiindex corresponding to $$\mathbf{x}_\mathbf{j}$$ = hj. We note that the particular ordering of the nodes affects the matrix structure corresponding to the discrete system but does not affect the numerical solution. First, at any node $$\mathbf{x}_\mathbf{i}$$ ∈ Ω, the nonlocal operator can be rewritten as \begin{align} \mathcal{L}_{\delta}u(\mathbf{x}_{\mathbf{i}}) = 2\int_{\mathcal{B}_{\delta}({\boldsymbol{0}})}\dfrac{u(\mathbf{z}+\mathbf{x}_\mathbf{i})-u(\mathbf{x}_{\mathbf{i}})}{W(\mathbf{z})}W(\mathbf{z})\rho_{\delta}(|\mathbf{z} | )\, \mathrm{d}\mathbf{z}, \end{align} (3.1) where W(z) represents a weight function. Then a quadrature-based finite difference scheme of the nonlocal operator (1.2) can be given as \begin{align} \mathcal{L}_{h,\,\delta}u(\mathbf{x}_{\mathbf{i}}) = 2\int_{\mathcal{B}_{\delta}({\boldsymbol{0}})}\mathcal{I}_{h}\left(\dfrac{u(\mathbf{z}+\mathbf{x}_{\mathbf{i}})-u(\mathbf{x}_{\mathbf{i}})}{W(\mathbf{z})}\right)W(\mathbf{z})\rho_{\delta}(|\mathbf{z} | )\, \mathrm{d}\mathbf{z}, \end{align} (3.2) where $$\mathcal{I}_{h}(\cdot )$$ represents the piecewise d-multilinear interpolation operator in z associated with the mesh $${\mathcal{T}}_{h}$$, that is, $$\mathcal{I}_{h}(u)$$ is piecewise linear with respect to each component of the spatial variable. The nonlocal constraint with u = 0 is imposed at nodes in $$\Omega _{\mathcal{I}}$$. For the one-dimensional case, in the study by Tian & Du (2014) the function W(z) = |z| was used for $$z\in \mathbb{R}$$. A key ingredient in its extension to the high-dimensional case is given by \begin{align} W(\mathbf{z})=\dfrac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}, \end{align} (3.3) where the notation |⋅|1 stands for the ℓ1 norm in the d-dimensional vector space while |⋅| denotes the standard Euclidean norm. The auxiliary function W = W(z) plays an important role in ensuring asymptotic compatibility for the multidimensional case. The resulting quadrature-based finite difference scheme for solving the nonlocal diffusion problem can be more conveniently written as the following: find $$\left \{u_{h}^{\delta }(\mathbf{x}_{\mathbf{i}})\right \}$$ such that $$u_{h}^{\delta }(\mathbf{x}_{\mathbf{i}})=0$$ for $$\mathbf{x}_{\mathbf{i}}\in \Omega _{\mathcal{I}}$$ and \begin{align} -\mathcal{L}_{h,\,\delta}u_{h}^{\delta}(\mathbf{x}_{\mathbf{i}}) = \sum_{\mathbf{x}_{\mathbf{j}}\in \mathcal{B}_{\delta}(\mathbf{x }_{\mathbf{i}})}\left(u_{h}^{\delta}(\mathbf{x}_{\mathbf{j}})-u_{h}^{\delta}(\mathbf{x}_{\mathbf{i}})\right)p^{\delta}_{\mathbf{x }_{\mathbf{i}},\mathbf{x}_{\mathbf{j}}} = f(\mathbf{x}_{\mathbf{i}}) \quad \forall\; \mathbf{x}_{\mathbf{i}}\in\Omega, \end{align} (3.4) where the sum over $$\mathbf{x }_\mathbf{j}$$ is for grid points in $$\mathcal{B}_{\delta }(\mathbf{x }_{\mathbf{i}})$$ but not including $$\mathbf{x }_\mathbf{i}$$, $$p^{\delta }_{\mathbf{x}_{\mathbf{i}},\mathbf{x}_{\mathbf{j}}}={2 \beta _{\mathbf{j}\mathbf{i}}}/{W(\mathbf{x}_{\mathbf{j}}-\mathbf{x}_{\mathbf{i}})}$$ with \begin{align} \beta_{\mathbf{j}\mathbf{i}} = \int_{ \mathcal{B}_{\delta}( \boldsymbol{0} ) } \phi_{\mathbf{j}}(\mathbf{z}+\mathbf{x}_{\mathbf{i}}) W(\mathbf{z}) \rho_{\delta}(|\mathbf{z}|) \, \mathrm{d}\mathbf{z}, \end{align} (3.5) and $$\phi_\mathbf{j}$$ is the piecewise multilinear basis function satisfying $$\phi_\mathbf{j} (\mathbf{x}_\mathbf{i})$$ = 0 when $$\mathbf{i}\neq \mathbf{j}$$ and $$\phi_\mathbf{j}(\mathbf{x}_\mathbf{j})$$ = 1. Remark 3.1 The discrete scheme (3.4) is valid for general kernels ρδ = ρδ(z) with bounded second moment in z. This is another benefit of using the weight W = W(z) in (3.3). Indeed, consider a kernel ρδ(|z|) = cδ, α, d|z|−α with α ∈ [d, d + 2) for $$\mathbf{z}\in \mathcal{B}_{\delta }(\boldsymbol{0})$$; then the integral in (3.5) may be unbounded for |i −j| = 1 if we do not use such a weight, or say, if W(z) = 1 is taken instead. Remark 3.2 By the symmetry and translation invariance of the uniform mesh it can be observed that $$\beta _{\mathbf{j}\mathbf{i}}=\beta _{\delta }\left (\mathbf{x}_{\mathbf{i}}-\mathbf{x}_{\mathbf{j}}\right )$$ for some function βδ. In fact, let ϕ0 be the piecewise multilinear basis function located at the origin satisfying ϕ0$$(\mathbf{x}_\mathbf{i})$$ = 0 with $$\mathbf{x}_{\mathbf{i}}\neq{\boldsymbol{0}}$$ and ϕ0(0) = 1; then \begin{align} \beta_{\mathbf{j}\mathbf{i}} = \beta_{\delta}(\mathbf{x}_{\mathbf{i}}-\mathbf{x}_{\mathbf{j}}),\qquad \beta_{\delta}(\mathbf{x})= \int_{ \mathcal{B}_{\delta}( \boldsymbol{0} ) } \phi_{0}(\mathbf{z}+\mathbf{x}) W(\mathbf{z}) \rho_{\delta}(|\mathbf{z}|) \,\mathrm{d}\mathbf{z} . \end{align} (3.6) Moreover, it is easy to see that βδ is symmetric with respect to any axis hyperplane. Also, $$\beta_{\mathbf{j}\mathbf{i}}$$ = 0 if the support of $$\phi_\mathbf{j}$$ does not overlap with the ball Bδ$$(\mathbf{x}_\mathbf{i})$$. In comparison, we give the conventional central finite difference scheme for the local diffusion equation (2.2): \begin{align} - \mathcal{L}_{h,\,0}{u_{h}^{0}}(\mathbf{x}_{\mathbf{i}})=-\sum_{|\mathbf{x}_{\mathbf{j}} - \mathbf{x }_{\mathbf{i}}|=h}\left({u_{h}^{0}}\left(\mathbf{x}_{\mathbf{j}}\right)-{u_{h}^{0}}\left(\mathbf{x}_{\mathbf{i}}\right)\right)\frac{1}{|\mathbf{x }_{\mathbf{j}}-\mathbf{x}_{\mathbf{i}}|^{2}} = f(\mathbf{x}_{\mathbf{i}}) \quad \forall\; \mathbf{x}_{\mathbf{i}}\in\Omega, \end{align} (3.7) that is, it involves only the nearest neighbors along the axes. 3.2 Convergence analysis Due to the symmetry of the kernel we may rewrite the nonlocal operator as $$\mathcal{L}_{\delta} u(\mathbf{x})= \int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\rho_{\delta}(|\mathbf{z}|)\big(u(\mathbf{x}+\mathbf{z})+u(\mathbf{x}-\mathbf{z})-2u(\mathbf{x})\big)\, \mathrm{d}\mathbf{z}$$ and the discrete operator as $$\mathcal{L}_{h,\,\delta} u(\mathbf{x}_{\mathbf{i}})= \sum_{\mathbf{j}\neq\boldsymbol{0}}\frac{u\left(\mathbf{x}_{\mathbf{i}}+\mathbf{z}_{\mathbf{j}}\right)+u\left(\mathbf{x}_{\mathbf{i}}-\mathbf{z}_{\mathbf{j}}\right)-2u(\mathbf{x}_{\mathbf{i}})}{|\mathbf{z}_{\mathbf{j}}|^{2}}|\mathbf{z}_{\mathbf{j}}|_{1}\beta_{\delta}({\mathbf{z}_{\mathbf{j}}}) .$$ We first state a few technical lemmas concerning the nonlocal operators and the nonlocal discrete schemes that are nonlocal analogs of their well-known local versions. First, by the fact that W = W(z), the kernel ρδ, as well as the multilinear basis functions $$\phi_\mathbf{j}$$, is non-negative, we see that βδ, $$\beta_{\mathbf{j}\mathbf{i}}$$ and thus $$p^{\delta }_{\mathbf{x}_{\mathbf{j}}\mathbf{x}_{\mathbf{i}}}$$ are non-negative. This implies the M-matrix property of the coefficient matrix, just like the case corresponding to the central difference approximations of the diffusion operator. Lemma 3.3 (M-matrix). The coefficient matrix, also denoted by $$\mathcal{L}_{h,\delta }$$, corresponding to the linear system (3.4) is an M-matrix. An important stability property can be derived via the discrete maximum principle that follows from the M-matrix property of the coefficient matrix. By constructing a suitable quadratic barrier function similar to the one-dimensional case by Tao et al. (2017), as well as the local counterpart, we can get the following lemma. Lemma 3.4 (Uniform stability). For any δ > 0, $$\mathcal{L}_{h,\delta }^{-1}$$ is bounded. Furthermore, for any 0 < δ < δ0, we have \begin{align} \left\|\mathcal{L}_{h,\delta}^{-1} \right\|_{\infty}\le C(\delta_{0}) \end{align} (3.8) for some constant $$C\!\left(\delta _{0}\right )$$ independent of h and δ as h → 0. Next we consider the uniform consistency derived via the truncation error calculation. Here the uniform consistency means that the truncation error is independent of δ for small δ. We first state the following useful and vital lemma, which is a nonlocal analog of the fact that the central difference approximation to the Laplacian, as given by (3.7), is exact for quadratic polynomials. Lemma 3.5 (Quadratic exactness). For any quadratic polynomial in $$\mathbb{R}^{d}$$ given as u(x) = x ⊗ x : M, where M = (mkj) is a constant matrix, we have \begin{align} \mathcal{L}_{h,\,\delta} u(\mathbf{x}_{\mathbf{i}})=\mathcal{L}_{\delta} u(\mathbf{x}_{\mathbf{i}})=\sum_{k} m_{kk}\quad \forall\, \mathbf{i}. \end{align} (3.9) Proof. The fact that $$\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\mathbf{z}\otimes\mathbf{z}\rho_{\delta}(|\mathbf{z}| ) \,\mathrm{d}\mathbf{z}=I_{d},$$ where Id is the d × d identity matrix, yields $$\mathcal{L}_{\delta} u(\mathbf{x}_{\mathbf{i}})\equiv\sum_{k} m_{kk}.$$ Meanwhile, when u(x) = x ⊗ x : M, it is obvious that $$u\left(\mathbf{x}_{\mathbf{i}}+\mathbf{z}_{\mathbf{j}}\right)+u\left(\mathbf{x}_{\mathbf{i}}-\mathbf{z}_{\mathbf{j}}\right)-2u(\mathbf{x}_{\mathbf{i}})=\mathbf{z}_{\mathbf{j}}\otimes\mathbf{z}_{\mathbf{j}}:M.$$ Thus, $$\mathcal{L}_{h,\,\delta} u(\mathbf{x}_{\mathbf{i}}) = \sum_{\mathbf{j}\neq\boldsymbol{0}}\frac{\mathbf{z}_{\mathbf{j}}\otimes\mathbf{z}_{\mathbf{j}}:M}{\left|\mathbf{z}_{\mathbf{j}}\right|{}^{2}}\left|\mathbf{z}_{\mathbf{j}}\right|{}_{1}\beta_{\delta}\left(\mathbf{z}_{\mathbf{j}}\right)= Q:M .$$ We now need another observation on the coefficient matrix, that is, $$\beta _{\delta }\left (\mathbf{z}_{\mathbf{j}}\right )$$ is an even function; moreover, it depends only on the absolute values of the individual components of $$\mathbf{z}_\mathbf{j}$$. Thus, when doing the summation over j, it is easy to check that Qkl = 0 for $$k\neq l$$. We can further check that Qkk is independent of k, and \begin{eqnarray} \sum_{k} Q_{kk}&=&\int_{\mathcal{B}_{\delta}({\boldsymbol{0}})}\left(\sum_{\mathbf{j}\neq\boldsymbol{0}}\frac{|\mathbf{z}_{\mathbf{j}}|^{2}}{|\mathbf{z}_{\mathbf{j}}|^{2}}|\mathbf{z}_{\mathbf{j}}|_{1}\phi_{\mathbf{j}}(\mathbf{z})\right)\frac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}\rho_{\delta}(|\mathbf{z}|) \, \mathrm{d}\mathbf{z}\nonumber\\ &=&\int_{\mathcal{B}_{\delta}({\boldsymbol{0}})}|\mathbf{z}|^{2}\rho_{\delta}(|\mathbf{z}|) \,\mathrm{d}\mathbf{z}=d,\nonumber \end{eqnarray} where in the last step we have used the property that |z|1 is piecewise linear, and thus $$\sum_{\mathbf{j}\neq\boldsymbol{0}}|\mathbf{z}_{\mathbf{j}}|_{1}\phi_{\mathbf{j}}(\mathbf{z})=\sum_{\mathbf{j}}|\mathbf{z}_{\mathbf{j}}|_{1}\phi_{\mathbf{j}}(\mathbf{z})=|\mathbf{z}|_{1}.$$ Hence, we can derive that Q = Id, which gives $$\mathcal{L}_{\delta} u(\mathbf{x}_{\mathbf{i}})=\mathcal{L}_{h,\,\delta} u(\mathbf{x}_{\mathbf{i}})\equiv\sum_{k} m_{kk}.$$ This completes the proof. Before proceeding to the main theorem of this work we introduce some notation for simplicity of presentation. For any function $$u=u(\mathbf{x}): \mathbb{R}^{d}\mapsto \mathbb{R}$$ and $$\boldsymbol{\alpha }\in \mathbb{N}^{d}$$ let $$\mathbf{x}^{\boldsymbol{\alpha}}=\prod_{i=1}^{d}x_{i}^{\alpha_{i}},\quad \partial^{\boldsymbol{\alpha}}u(\mathbf{x})=\partial_{\alpha_{1}}\cdots\partial_{\alpha_{d}}u(\mathbf{x}),$$ and Dku(x) be the set containing all the kth-order derivatives in the form of ∂αu(x) with |α|1 = k, and $$\left|D^{k}u\right|{}_{\infty}=\max_{|\boldsymbol{\alpha}|_{1}=k}\sup_{\mathbf{x}}\partial^{\boldsymbol{\alpha}}u(\mathbf{x}).$$ We further use D2 to denote the Hessian and define $$H^{\mathbf{x}}(\mathbf{z})=u(\mathbf{x}+\mathbf{z})+u(\mathbf{x}-\mathbf{z})-2u(\mathbf{x}),\quad H^{\mathbf{x}}_{1}(\mathbf{z})=\mathbf{z}\otimes\mathbf{z}:D_{2}u(\mathbf{x})$$ and $$J^{\mathbf{x}}(\mathbf{z})=\frac{H^{\mathbf{x}}(\mathbf{z})}{|\mathbf{z}|^{2}}|\mathbf{z}|_{1},\quad J^{\mathbf{x}}_{1}(\mathbf{z})=\frac{H^{\mathbf{x}}_{1}(\mathbf{z})}{|\mathbf{z}|^{2}}|\mathbf{z}|_{1},$$ and let H2x(z) = Hx(z) − H1x(z) and J2x(z) = Jx(z) − J1x(z). By Taylor’s expansion and the mean value theorem we can show that \begin{align} \left\{ \begin{array}{rcl} \displaystyle H^{\mathbf{x}}_{2}(\mathbf{z})&=& \displaystyle \sum_{|\boldsymbol{\alpha}|_{1}=4}C_{\boldsymbol{\alpha}}^{1}\mathbf{z}^{\boldsymbol{\alpha}}\partial^{\boldsymbol{\alpha}}u(\mathbf{x}_{1}(\mathbf{z})),\\ \displaystyle \frac{\partial}{\partial z_{i}}H^{\mathbf{x}}_{2}(\mathbf{z})&=& \displaystyle\sum_{|\boldsymbol{\alpha}|_{1}=3}C_{\boldsymbol{\alpha}}^{2}\mathbf{z}^{\boldsymbol{\alpha}}\partial^{\boldsymbol{\alpha}+\mathbf{e}_{i}}u(\mathbf{x}_{2}(\mathbf{z})),\\ \displaystyle \frac{\partial^{2}}{\partial z_{i}\partial z_{j}}H^{\mathbf{x}}_{2}(\mathbf{z})&=& \displaystyle \sum_{|\boldsymbol{\alpha}|_{1}=2}C_{\boldsymbol{\alpha}}^{3}\mathbf{z}^{\boldsymbol{\alpha}}\partial^{\boldsymbol{\alpha}+\mathbf{e}_{i}+\mathbf{e}_{j}}u(\mathbf{x}_{3}(\mathbf{z})), \end{array} \right. \end{align} (3.10) where ei and ej are unit vectors with only one nonzero element at the ith and jth positions, respectively, and xk(z) are some points dependent on both x and z. The following is an elementary but technical result to be used later. Lemma 3.6 For J2x(z) defined above, if $$u\in C^{4}(\overline \Omega )$$, for any z and x, we have \begin{align} \big|D^{2}J^{\mathbf{x}}_{2}\big|_{\infty}\le C\big|D^{4}u\big|_{\infty}|\mathbf{z}|_{1}, \end{align} (3.11) where C is generic constant. Proof. To carry out the estimates we first note that $$\frac{\partial |\mathbf{z}|^{2}}{\partial z_{i}}=2z_{i}\quad\mbox{ and }\quad \frac{\partial |\mathbf{z}|^{2}}{\partial z_{i}\partial z_{j}}=2\delta_{ij} .$$ Now, by the definition of J2x(z), we have $$\frac{\partial^{2}}{\partial z_{i}\partial z_{j}}J^{\mathbf{x}}_{2}(\mathbf{z})=\textrm{sign}(z_{i})\frac{\partial}{\partial z_{j}}\frac{H^{\mathbf{x}}_{2}(\mathbf{z})}{|\mathbf{z}|^{2}}+\textrm{sign}(z_{j})\frac{\partial}{\partial z_{i}}\frac{H^{\mathbf{x}}_{2}(\mathbf{z})}{|\mathbf{z}|^{2}}+|\mathbf{z}|_{1}\frac{\partial^{2}}{\partial z_{i}\partial z_{j}}\frac{H^{\mathbf{x}}_{2}(\mathbf{z})}{|\mathbf{z}|^{2}},$$ where sign (⋅) is the sign function. For the first term in the above, we have from (3.10) that $$\left| \textrm{sign}(z_{i}) \frac{\partial}{\partial z_{j}}\frac{H^{\mathrm{x}}_{2}(\mathbf{z})}{|\mathbf{z}|^{2}}\right|\le C\left(\left|\frac{\mathbf{z}^{\boldsymbol{\beta}_{1}}}{|\mathbf{z}|^{2}}\right|+\left|\frac{\mathbf{z}^{\boldsymbol{\beta}_{2}}}{|\mathbf{z}|^{4}}\right|\right)|D^{4}u|_{\infty},$$ where $$\left |\boldsymbol{\beta }_{1}\right |_{1}=3$$ and $$\left |\boldsymbol{\beta }_{2}\right |_{1}=5$$. The second term satisfies a similar estimate. Moreover, by (3.10) again, we get $$\left| \frac{\partial^{2}}{\partial z_{i}\partial z_{j}}\frac{H^{\mathbf{x}}_{2}(\mathbf{z})}{{|\mathbf{z}|^{2}}}\right|\le C\left(\left|\frac{\mathbf{z}^{\boldsymbol{\alpha}_{1}}}{|\mathbf{z}|^{4}}\right|+\left|\frac{\mathbf{z}^{\boldsymbol{\alpha}_{2}}}{|\mathbf{z}|^{8}}\right|\right)|D^{4}u|_{\infty},$$ where |α1|1 = 4 and |α2|1 = 8. For any $$\boldsymbol{\alpha }\in \mathbb{N}^{d}$$ it is obvious that $$\left|\frac{\mathbf{z}^{\boldsymbol{\alpha}}}{|\mathbf{z}|^{|\boldsymbol{\alpha}|_{1}}}\right|\le 1,\quad \left|\frac{\mathbf{z}^{\boldsymbol{\alpha}}}{|\mathbf{z}|^{|\boldsymbol{\alpha}|_{1}-1}}\right|\le |\mathbf{z}|_{1}.$$ Combining the above estimates we can get the desired result (3.11). We now present the truncation error analysis for the quadrature-based difference approximation of the nonlocal model. Lemma 3.7 (Uniform nonlocal truncation error). Assume that $$u\in C^{4}(\overline{\Omega \cup \Omega _{\mathcal{I}}})$$. Then it holds that \begin{align} \max_{1\leq i\leq N_{s}} |\mathcal{L}_{h,\,\delta} u(\mathbf{x}_{\mathbf{i}})-\mathcal{L}_{\delta} u(\mathbf{x}_{\mathbf{i}})|\le C|D^{4}u|_{\infty} h^{2}\!, \end{align} (3.12) where C is a constant independent of δ and h. Proof. Without loss of generality we take xi = 0 to be the origin. Note that $$\mathcal{L}_{\delta}u(\boldsymbol{0}) =\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}J^{\boldsymbol{0}}(\mathbf{z})\frac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}\rho_{\delta}(|\mathbf{z}| ) \, \mathrm{d}\mathbf{z}$$ and $$\mathcal{L}_{h,\delta}u(\boldsymbol{0})=\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\mathcal{I}_{h} \left(J^{\boldsymbol{0}}(\mathbf{z})\right)\frac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}\rho_{\delta}(|\mathbf{z}| ) \, \mathrm{d}\mathbf{z} .$$ Simple subtraction gives us \begin{eqnarray*}\left|\mathcal{L}_{h,\delta}u(\boldsymbol{0})-\mathcal{L}_{\delta}u(\boldsymbol{0})\right|&&=\left|\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\left(\mathcal{I}_{h} \left(J^{\boldsymbol{0}}(\mathbf{z})\right)-J^{\boldsymbol{0}}(\mathbf{z})\right)\frac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}\rho_{\delta}(|\mathbf{z}| ) \, \mathrm{d}\mathbf{z}\right|,\\ && \le\left|\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\left(\mathcal{I}_{h} \left(J_{1}^{\boldsymbol{0}}(\mathbf{z})\right)-J_{1}^{\boldsymbol{0}}(\mathbf{z})\right)\frac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}\rho_{\delta}(|\mathbf{z}| ) \, \mathrm{d}\mathbf{z}\right|\\ &&\quad+\left|\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\left(\mathcal{I}_{h} \left(J_{2}^{\boldsymbol{0}}(\mathbf{z})\right)-J_{2}^{\boldsymbol{0}}(\mathbf{z})\right)\frac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}\rho_{\delta}(|\mathbf{z}| ) \,\mathrm{d}\mathbf{z}\right|\\ && :=E_{1}+E_{2}, \end{eqnarray*} where $$J_{1}^{\boldsymbol{0}}(\mathbf{z})=\frac{\mathbf{z}\otimes\mathbf{z}:D_{2}u(\boldsymbol{0})}{|\mathbf{z}|^{2}}|\mathbf{z}|_{1}\quad\mbox{ and }\quad J_{2}^{\boldsymbol{0}}(\mathbf{z})=J^{\boldsymbol{0}}(\mathbf{z})-J_{1}^{\boldsymbol{0}}(\mathbf{z}).$$ From Lemma 3.5, we know that E1 = 0. For E2, we have $$E_{2}\le Ch^{2}\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}\left| D^{2} J_{2}^{\boldsymbol{0}}\right|_{\infty}\frac{|\mathbf{z}|^{2}}{|\mathbf{z}|_{1}}\rho_{\delta}(|\mathbf{z}| ) \,\mathrm{d}\mathbf{z}\,.$$ Thus, by using Lemma 3.6 and the moment condition on the kernel that $$\int_{\mathcal{B}_{\delta}(\boldsymbol{0})}|\mathbf{z}|^{2}\rho_{\delta}(|\mathbf{z}| ) \,\mathrm{d}\mathbf{z}=d,$$ we complete the proof. Combining the stability result in Lemma 3.4 and the consistency result in Lemma 3.7, we immediately get the convergence of the quadrature-based finite difference scheme. Theorem 3.8 (Convergence to nonlocal solution). Assume that the nonlocal exact solution uδ of (1.1) is smooth enough, such as $$u^{\delta }\in C^{4}\left (\overline{\Omega \cup \Omega _{\mathcal{I}}}\right )$$, and uδ, h is the numerical solution obtained by the scheme (3.4). Then for any fixed δ there exists some constant Cδ independent of h as h → 0 such that \begin{align} \left\|{u^{\delta,h}}-I_{h}u^{\delta}\right\|_{\infty}\le Ch^{2}, \end{align} (3.13) where Ihuδ is the interpolation of the exact solution on the grid points. For any δ < δ0 the constant Cδ is uniformly bounded by some generic constant C(δ0). Next we show the asymptotic compatibility of the quadrature-based finite difference scheme to the local limiting model. We assume that the local equation (2.2) is discretized by the classical central finite difference scheme denoted by $$\mathcal{L}_{h,0}$$, whose truncation error is \begin{align} \max_{1\leq i\leq N_{s}} |\mathcal{L}_{h,\,0} u(\mathbf{x}_{\mathbf{i}})-\mathcal{L}_{0} u(\mathbf{x}_{\mathbf{i}})|=\mathcal{O}\big(h^{2}\big), \end{align} (3.14) by the standard numerical PDE analysis, if $$u\in C^{4}(\overline \Omega )$$. We now estimate the modeling error at the discrete level, that is, the error between the discrete nonlocal operator $$\mathcal{L}_{h,\delta }$$ and $$\mathcal{L}_{h,0}$$. Lemma 3.9 (Discrete model error). Assume that $$u\in C^{4}\left (\overline{\Omega \cup \Omega _{\mathcal{I}}}\right )$$. Then it holds that \begin{align} \max_{1\leq i\leq N_{s}} |\mathcal{L}_{h,\,\delta} u(\mathbf{x}_{\mathbf{i}})-\mathcal{L}_{h,0} u(\mathbf{x}_{\mathbf{i}})|=\mathcal{O}\left(\delta^{2}\right)+\mathcal{O}\left(h^{2}\right) \end{align} (3.15) when δ and h are sufficiently small. Proof. The triangle inequality gives \begin{eqnarray*} |\mathcal{L}_{h,\delta} u(\mathbf{x}_{\mathbf{i}})-\mathcal{L}_{h,0} u(\mathbf{x}_{\mathbf{i}})|&\le &|\mathcal{L}_{h,\,\delta} u(\mathbf{x}_{\mathbf{i}})-\mathcal{L}_{\delta} u(\mathbf{x}_{\mathbf{i}})|+ |\mathcal{L}_{0} u(\mathbf{x}_{\mathbf{i}})-\mathcal{L}_{\delta} u(\mathbf{x}_{\mathbf{i}})|\\ &&+ |\mathcal{L}_{0} u(\mathbf{x}_{\mathbf{i}})-\mathcal{L}_{h,\,0} u(\mathbf{x}_{\mathbf{i}})|\\ &:=& R_{1}+R_{2}+R_{3}. \end{eqnarray*} By Lemma 3.7 we have $$R_{1}= \mathcal{O}\big(h^{2}\big).$$ For $$u\in C^{4}(\overline \Omega )$$ the continuum property of nonlocal operators gives $$R_{2}=\mathcal{O}\big(\delta^{2}\big).$$ It is well known that $$R_{3}=\mathcal{O}\big(h^{2}\big).$$ The three equalities above together complete the proof. Combining the consistency in Lemma 3.9 with the stability Lemma 3.4 again, we get convergence (and asymptotic compatibility) to the local limit. Theorem 3.10 (Asymptotic compatibility). Assume that the local exact solution u0 of (2.2) is smooth enough, such as $$u^{0}\in C^{4}\left (\overline{\Omega \cup \Omega _{\mathcal{I}}}\right )$$. Then for any δ < δ0 there is some constant C independent of δ and h as h → 0 such that \begin{align} \left\|{u^{\delta,h}}-I_{h}u^{0}\right\|_{\infty}\le C\left(h^{2}+\delta^{2}\right)\!, \end{align} (3.16) where Ihu0 is the interpolation of the exact solution on the grid points. Remark 3.11 For problems on the same domain but with periodic boundary conditions (PBCs) all similar estimates also hold, which requires only that the exact solution of the local model is a $$C_{\mathrm{per}}^{4}(\overline{\Omega })$$ function. Remark 3.12 We now give a more explicit form of the algebraic system in two dimensions. We let $$u^{\delta ,h}_{j,k}$$ denote the nodal value of the nonlocal difference solution at the mesh point $$\left (x_{j}, y_{k}\right )\in \Omega =(-\pi ,\pi )^{2} \subset \mathbb{R}^{2}$$. Then let r = ⌈δ/h⌉ be the smallest integer larger than or equal to δ/h; we have $$\mathcal{L}_{h,\,\delta} u^{\delta,\,h}_{j,\,k} =\sum_{p=0}^{r}\sum_{q=0}^{r}c_{p,\,q}\left(u^{\delta,\,h}_{j+p,\,k+q}+u^{\delta,\,h}_{j-p,\,k+q}+ u^{\delta,\,h}_{j+p,\,k-q}+u^{\delta,\,h}_{j-p,\,k-q}-4u^{\delta,\,h}_{j,\,k}\right),$$ where c0, 0 = 0 and $$c_{p,\,q}= \frac{ph+qh}{p^{2}h^{2}+q^{2}h^{2}} \iint_{B_{+}(0,\, \delta)} \phi_{p,\,q}(s,t)\rho_{\delta}\left(\sqrt{s^{2}+t^{2}}\right)\frac{s^{2}+t^{2}}{s+t} \,\mathrm{d}s\, \mathrm{d}t\,,$$ with B+(0, δ) denoting the first quadrant of the disc at the origin with radius δ. One may also observe that cp, q = cq, p for any p and q. 4. Green’s functions of a nonlocal operator with integrable kernels By definition, the Green’s function of the nonlocal equation (corresponding to the particular nonlocal operator $$\mathcal{L}_{\delta }$$) solves the nonlocal equation with right-hand side given by the Dirac delta measure δy(x) = δ(x −y) of the variable x for a given y. Discussions of Green’s functions for nonlocal models such as PD have been presented in a number of earlier studies (Weckner et al., 2009; Silling et al., 2010b; Wang et al., 2016, 2017). Here we pay particular attention to the case where the nonlocal interaction kernel is given to be an integrable one. A consequence is that, unlike Green’s functions for local diffusion equations, the nonlocal analog is generically a measure function. Hence, we focus on effective constructions in this case. 4.1 Periodic case For the case of PBCs with Ω = [−π, π]d, we may impose the zero mean compatibility condition for both the right-hand side and the solution over Ω. Thus, we formally have the following modified equation of (1.1): $$-{{\mathcal{L}}_{\delta}} {G^{\delta}}_{\mathbf{y}}(\mathbf{x})={\boldsymbol{\delta}}_{\mathbf{y}}-\frac{1}{|\Omega|},$$ where |Ω| is the area of the domain. Taking the L2-inner product of the above equation with uδ(x) over x yields $$\left(-\mathcal{L}_{\delta} G^{\delta}_{\mathbf{y}},u^{\delta}\right)=u^{\delta}(\mathbf{y })-\frac{1}{|\Omega|}\int_{\Omega} u^{\delta}\, \mathrm{d}\mathbf{x} =u^{\delta}(\mathbf{y}).$$ By the nonlocal Green’s identity (Du et al., 2012, 2013), which remains valid for the periodic case, we have $$\left (-\mathcal{L}_{\delta } G^{\delta }_{\mathbf{y }},u^{\delta }\right )=\left (G^{\delta }_{\mathbf{y}},-\mathcal{L}_{\delta } u^{\delta }\right )$$. We then get $$u(\mathbf{y})=\left(f,G^{\delta}_{\mathbf{y}}\right),$$ which demonstrates how Green’s function $$G^{\delta }_{\mathbf{y}}$$ can be used to represent solutions to the nonlocal equation (1.1). Now, for the case where the kernel is integrable we denote $$C_{\delta}=\int_{\mathcal{B}_{\delta}(\bf{0})} \rho_{\delta}(|\mathbf{z}|) \, \mathrm{d}\mathbf{z}.$$ Then the nonlocal operator $$-\mathcal{L}_{\delta }$$ can be written as $$-\mathcal{L}_{\delta} u = -\rho_{\delta} \ast u + C_{\delta} u\,.$$ Let us perform the splitting $$G^{\delta }_{\mathbf{y}}=S_{1}+g_{0}$$ with S1 and g0 to be determined; we then arrive at $$-{{\mathcal{L}}_{\delta}} S_{1}(\mathbf{x})-{\rho_{\delta}} \ast g_{0} + C_{\delta} g_{0}={\boldsymbol{\delta}}_{\mathbf{y}}-\frac{1}{|\Omega|}.$$ Setting \begin{align} g_{0}=C_{\delta}^{-1}\left(\boldsymbol{\delta}_{\mathbf{y}}-\frac{1}{|\Omega|}\right) \quad \mbox{and}\quad \; g_{1}=C_{\delta}^{-1}\rho_{\delta} \ast g_{0} \end{align} (4.1) leads to \begin{align} -\mathcal{L}_{\delta} S_{1}(\mathbf{x})=\rho_{\delta} \ast g_{0} (\mathbf{x})= C_{\delta} g_{1}\, . \end{align} (4.2) Remark 4.1 Note that while g0 is a singular measure we have generically that g1 ∈ L1 for an integrable ρδ. If additional regularity of ρδ can be assumed then we also get g1 of the same regularity class as ρδ. For example, if ρδ is of the class $$C^{\infty }$$ (as a function in the whole space having compact support in the ball of radius δ) then we also have g1 and S1 in $$C^{\infty }$$. Another example, in the one-dimensional case, is that if ρδ is a bounded variation (BV) function then so are g1 and S1. In turn, if we perform a similar splitting once again as S1 = S2 + g1, we have $$-{{\mathcal{L}}_{\delta}} S_{2}(\mathbf{x})={\rho_{\delta}} \ast g_{1}.$$ Similarly, repeating this splitting technique n times we get recursively that $$g_{n}=C_{\delta}^{-1}\rho_{\delta} \ast g_{n-1},\quad -\mathcal{L}_{\delta} S_{n+1}(\mathbf{x })=\rho_{\delta} \ast g_{n},$$ and generically, the later terms obtained by the recursion are more regular than the preceding terms, similarly to the observation in Remark 4.1. In the end, we get $$G^{\delta}_{\mathbf{y}}=\sum_{n=0}^{\infty}g_{n}$$ with the following recursive relation: $$g_{0}=C_{\delta}^{-1}\left(\boldsymbol{\delta}_{\mathbf{y}}-\frac{1}{|\Omega|}\right)\!,\quad g_{n}=C_{\delta}^{-1}\rho_{\delta} \ast g_{n-1},\; n=1,\;2,\;\dots .$$ Thus, we obtain an explicit series form of $$G^{\delta }_{\mathbf{y}}(\mathbf{x})$$, which is given below. Theorem 4.2 For an integrable kernel ρδ the nonlocal Green’s function $$G^{\delta }_{\mathbf{y}}$$ for equation (1.1) with PBCs can be represented by \begin{align} G^{\delta}_{\mathbf{y}}=\sum_{n=0}^{\infty} C_{\delta}^{-(n+1)} \left(\rho_{\delta}\ast\right)^{n}\ast \left(\boldsymbol{\delta}_{\mathbf{y}}-\frac{1}{|\Omega|}\right). \end{align} (4.4) Moreover, the solution uδ of the nonlocal diffusion equation (1.1) is given by $$u^{\delta}(\mathbf{y})=\left(f,G^{\delta}_{\mathbf{y}}\right).$$ Remark 4.3 Due to the lack of elliptic smoothing of the nonlocal diffusion equation corresponding to the integrable kernels, the nonlocal Green’s function remains a composition of a singular part in the form of a Dirac delta measure and an integrable part. In this case and unlike local Green’s functions for elliptic PDEs, $$G^{\delta }_{\mathbf{y}}$$ should not be viewed as a function defined almost everywhere. Moreover, the regularity of the solution of uδ is in general the same as the right-hand side f, since the regularity $$\big (\,f,G^{\delta }_{\mathbf{y}}\big )$$ is dominated by the first term (g0, f), which is proportional to f. One could also split the solution into two parts, one given by (g0, f) while the other accounts for a more regular part. Remark 4.4 It is worth mentioning that the series expression for the Green’s function is a Neumann series expansion for the equation (I − A)G = f where $$Au = C_{\delta }^{-1}\rho _{\delta }\ast u$$ and $$f = C_{\delta }^{-1}\left (\boldsymbol{\delta }_{\mathbf{y}}-\frac{1}{|\Omega |}\right )$$. 4.2 Other cases of nonlocal constraints The approach of constructing Green’s function through a series expansion is also applicable to other nonlocal constraints. For instance, we consider the Dirichlet-type volume-constrained problems given by (1.1) and (2.1). For any y ∈ Ω let us define g0 = δ(x −y), $$\mathbf{x }\in \Omega \cup \Omega _{\mathcal{I}}$$. It is easy to see that, with gn recursively defined by $$\begin{cases} g^{n}(\mathbf{x}) =C_{\delta}^{-1}\rho_{\delta}\ast g^{n-1}(\mathbf{x}), & \mathbf{x}\in\Omega,\\ g^{n}(\mathbf{x}) =0, & \mathbf{x}\in\Omega_{\mathcal{I}},\end{cases}$$ the nonlocal Green’s function is given by $$G^{\delta}_{\mathbf{y}}(\mathbf{x})=\sum_{n=0}^{\infty} g^{n}(\mathbf{x}).$$ As in the case with PBC, we can also use the above expansion to get a series expression for a general solution, though the regularity pickup in the later terms of the series gets more delicate due to the dependence on matching the boundary data with the right-hand side and remains an interesting topic to be further studied. Remark 4.5 Other types of series expansions, in particular Fourier series expansion, of the nonlocal Green’s functions have been presented in the literature by Weckner et al. (2009), Silling et al. (2010b) and Wang et al. (2016, 2017). Such expansions are valid formally, but for integrable kernels the series do not converge pointwisely. 4.3 Discrete nonlocal Green’s function Following the discussion above on the continuum level we can present a discrete analog with the AC quadrature-based finite difference scheme. Given the uniform mesh $$\left \{\mathbf{x}_{\mathbf{j}}\right \}$$ the discrete Green’s function (matrix) $$G^{h,\delta }_{\mathbf{x}_{\mathbf{i}}}\left (\mathbf{x}_{\mathbf{j}}\right )$$ is given by $$-{{\mathcal{L}}_{h,\,\delta}} {G^{h,\,\delta}}_{{\mathbf{x}}_{\mathbf{i}}}\left(\mathbf{x }_{\mathbf{j}}\right)=\frac{1}{h^{d}} {\delta_{\mathbf{i}\mathbf{j}}}-\frac{1}{|\Omega|},$$ where δij is the Kronecker delta function, that is, δii = 1 and δij = 0 for $$\mathbf{i}\neq \mathbf{j}$$. It follows that the discrete solution to the nonlocal model with a discrete right-hand side $$\left \{f_{\mathbf{j}}\right \}$$ is given by $$u^{h,\,\delta}(\mathbf{x}_{\mathbf{i}})=\sum_{\mathbf{j}\neq\mathbf{i} } f_{\mathbf{j}} G^{h,\delta}_{\mathbf{x}_{\mathbf{i}}}\left(\mathbf{x}_{\mathbf{j}}\right)\!\quad\mbox{ for }\quad \sum_{\mathbf{j}} f_{\mathbf{j}} = 0.$$ The study of discrete nonlocal Green’s functions can be seen as a nonlocal analog of similar works on discrete Green’s functions for local PDEs (Chung & Yau, 2000). Although the nonlocal Green’s function at the continuum level is measure valued for integrable interaction kernels the discrete nonlocal Green’s function given above is well defined at all the grid points for any finite δ and h. Due to the lack of regularity, however, we do not expect strong (or pointwise) convergence of $$G^{h,\delta }_{\mathbf{y}}(\mathbf{x})$$ to $$G^{\delta }_{\mathbf{y}}(\mathbf{x})$$ as h → 0 for a given δ. Instead, we could attempt to get better approximations by the splitting technique, that is, we numerically solve for the regular part while invoking the analytical solution of singular part. Fig. 1. View largeDownload slide The function $$u^{h,\delta }_{1}$$ of the regular part of the nonlocal Green’s function. Fig. 1. View largeDownload slide The function $$u^{h,\delta }_{1}$$ of the regular part of the nonlocal Green’s function. Therefore, we obtain a hybrid approximation of $$G^{\delta }_{\mathbf{x}_{\mathbf{j}}}(\mathbf{x }_{\mathbf{i}})$$ by $$G^{h,\,\delta}_{\mathbf{x}_{\mathbf{j}}}(\mathbf{x}_{\mathbf{i}}) = C_{\delta}^{-1}\left(\boldsymbol{\delta}_{\mathbf{x}_{\mathbf{j}}}(\mathbf{x }_{\mathbf{i}})-\frac{1}{|\Omega|}\right)+u^{h,\,\delta}_{1},$$ where the approximation $$u^{h,\delta }_{1}$$ of the regular part of the nonlocal Green’s function S1 in (4.2) satisfies $$-{{\mathcal{L}}_{h,\,\delta}} {u^{h,\,\delta}}_{1} (\mathbf{x}_{\mathbf{i}})= {\rho_{\delta}} \ast g_{0}(\mathbf{x}_{\mathbf{i}}).$$ Given the higher regularity of S1 we expect to have good convergence of $$u^{h,\delta }_{1}$$ to S1 for any given δ, which is illustrated in Fig. 1 for a kernel that is a positive constant over the horizon and vanishes outside (thus piecewise continuous). Moreover, due to the AC property of the discrete scheme, we also have convergence of $$G^{h,\delta }_{\mathbf{x}_{\mathbf{j}}}(\mathbf{x}_{\mathbf{i}})$$ to the local Green’s function as both h and δ approach zero. Notice that both Theorems 3.8 and 3.10 require that the solutions to be approximated are $$C^{4}\left (\overline{\Omega \cup \Omega _{\mathcal{I}}}\right )$$. Thus, if we subtract enough terms of the Green’s function series expansion from the right-hand side, in principle we should be left with a smooth enough right-hand side that the results of both theorems are applicable. More discussions on this splitting technique are to be presented later. We thus can expect a robust approximation of the local and nonlocal Green’s functions. 4.4 More discussion on the splitting technique We next discuss a potential application of the splitting technique presented above in the numerical solution of nonlocal equations. By avoiding using derivatives, nonlocal models allow more singular solutions, which is a distinct advantage in physical modeling such as the study of material cracks. As an illustration we present a typical example showing how the splitting technique improves the performance of Fourier spectral methods for numerically solving nonlocal equations with discontinuous solutions. For simplicity, we take a one-dimensional nonlocal diffusion model (1.1) over a periodic cell Ω = [0, 2π] with a constant kernel $$\rho _{\delta }(s)=\frac{3}{\delta ^{3}}\chi _{[-\delta ,\delta ]}$$, which is of the rescaled form given in (2.3), and a discontinuous right-hand side $$f(x)=\textrm{sign}(x-\pi), \quad x\in [0,2\pi],$$ where sign(⋅) is the sign function. In this setting, the exact solution has a jump at x = π, e.g., like studied in Du & Yang (2016). One numerical solution is obtained by applying the Fourier spectral method directly. It is not surprising that the well-known Gibbs phenomena occur around x = π when standard Fourier spectral methods are directly used to approximate discontinuous solutions. On the other hand, if we use the splitting technique described above once, by first constructing the discontinuous part of the solution analytically, and then solving a new nonlocal diffusion equation with a smoother right-hand side by Fourier spectral methods, we can eliminate the Gibbs phenomena, as demonstrated in Fig. 2. Fig. 2. View largeDownload slide Numerical comparison. Left: numerical solutions; right: zoom in around $$\pi\!.$$ Fig. 2. View largeDownload slide Numerical comparison. Left: numerical solutions; right: zoom in around $$\pi\!.$$ 5. Conclusion In this work, a quadrature-based finite difference scheme for a nonlocal diffusion model in multidimensions is proposed. Our analysis is given for problems with Dirichlet volumetric constraints. The case of a Neumann-type nonlocal volumetric constraint may involve additional complications, due to the need for ghost points, which require further attention to ensure second-order accuracy. One may adapt the scheme developed here to more general models. In fact, an extension based on the idea presented in a draft version of this work has been given for nonlocal convection diffusion equations in the study by Tian et al. (2018) (where only a partial analysis of the diffusion term has been documented and one should refer to the more complete analysis presented in the current work). In deriving the error analysis our attention is focused on the case where the underlying solutions are smooth. Since nonlocal models are effective in describing physical processes with solution singularities, extending the analysis to cases with minimal regularity will be of interest. One can of course also consider extensions to systems involving vector fields such as the bond-based and state-based PD models. Our study here also illustrated an important difference between the nonlocal Green’s function and their local analog as the former may take on a possibly singular measure form. Nonlocal Green’s functions can be useful in various applications such as the mean exit time computation of jump processes and analysis of mechanical responses to point loads in nonlocal mechanics (Weckner et al., 2009; Silling et al., 2010b; Du et al., 2012; Wang et al., 2016, 2017). The splitting technique presented here offers an effective algorithm for computing singular Green’s functions. Moreover, one may further explore the applications of such techniques to the computation of singular solutions to more general nonlocal models as alluded to in Section 4.4. In comparison with the Galerkin finite element AC scheme developed in the study by Tian & Du (2014) that works for linear systems in multidimensions using unstructured meshes, we see that the quadrature-based finite difference AC schemes given here are developed only for uniform Cartesian meshes with the same mesh size in all directions. However, the schemes given here are simpler to implement as they avoid the evaluation of integrals in $$\mathbb{R}^{2d}$$. Furthermore, they lead to M-matrices and the discrete maximum principle that are in general not the case for the finite element discretization. Hence, an interesting future work is to develop hybrid schemes and particle discretization (Liu et al., 1996; Silling, 2000; Askari et al., 2008; Bessa et al., 2014) that can preserve the nice properties of the quadrature-based finite difference AC schemes while applicable to more general meshes. Funding United States National Science Foundation (DMS-1719699), Air Force Office of Scientific Research MURI center and the Army Research Office MURI Grant (W911NF-15-1-0562). References Andreu , F. , Mazón , J. M. , Rossi , J. D. & Toledo , J. ( 2010 ) Nonlocal Diffusion Problems . Mathematical Surveys and Monographs , vol. 165 . Providence, RI : AMS. Applebaum , D. ( 2004 ) Lévy Processes and Stochastic Calculus . Cambridge Studies in Advanced Mathematics , vol. 93 . Cambridge, UK : Cambridge University Press . Askari , E. , Bobaru , F. , Lehoucq , R. B. , Parks , M. L. , Silling , S. A. & Weckner , O. ( 2008 ) Peridynamics for multiscale materials modeling . J. Phys. Conf. Ser. , 125 , 12 – 78 . Google Scholar CrossRef Search ADS Bates , P. W. & Chmaj , A. ( 1999 ) An integrodifferential model for phase transitions: stationary solutions in higher space dimensions . J. Stat. Phys. , 95 , 1119 – 1139 . Google Scholar CrossRef Search ADS Bessa , M. , Foster , J. , Belytschko , T. & Liu , W. K. ( 2014 ) A meshfree unification: reproducing kernel peridynamics . Comput. Mech. , 53 , 1251 – 1264 . Google Scholar CrossRef Search ADS Bobaru , F. & Duangpanya , M. ( 2010 ) The peridynamic formulation for transient heat conduction . Int. J. Heat Mass Transf. , 53 , 4047 – 4059 . Google Scholar CrossRef Search ADS Bobaru , F. , Yang , M. , Alves , L. F. , Silling , S. A. , Askari , E. & Xu , J. ( 2009 ) Convergence, adaptive refinement, and scaling in 1D peridynamics . Int. J. Numer. Methods Eng. , 77 , 852 – 877 . Google Scholar CrossRef Search ADS Buades , A. , Coll , B. & Morel , J. M. ( 2010 ) Image denoising methods. A new nonlocal principle . SIAM Rev. , 52 , 113 – 147 . Google Scholar CrossRef Search ADS Chen , A. , Du , Q. , Li , C. & Zhou , Z. ( 2017 ) Asymptotically compatible schemes for space–time nonlocal diffusion equations . Chaos Solitons Fractals (in press) . Google Scholar CrossRef Search ADS Chen , X. & Gunzburger , M. ( 2011 ) Continuous and discontinuous finite element methods for a peridynamics model of mechanics . Comput. Methods Appl. Mech. Eng. , 200 , 1237 – 1250 . Google Scholar CrossRef Search ADS Chung , F. & Yau , S-T. ( 2000 ) Discrete Green’s functions . J. Comb. Theory A , 91 , 191 – 214 . Google Scholar CrossRef Search ADS Du , Q. , Gunzburger , M. , Lehoucq , R. B. & Zhou , K. ( 2012 ) Analysis and approximation of nonlocal diffusion problems with volume constraints . SIAM Rev. , 56 , 676 – 696 . Du , Q. , Gunzburger , M. , Lehoucq , R. B. & Zhou , K. ( 2013 ) A nonlocal vector calculus, nonlocal volume-constrained problems, and nonlocal balance laws . Math. Models Methods Appl. Sci. , 23 , 493 – 540 . Google Scholar CrossRef Search ADS Du , Q. , Ju , L. & Lu , J. ( 2017a ) A discontinuous Galerkin method for one-dimensional time-dependent nonlocal diffusion problems , Preprint . Du , Q. , Tao , Y. , Tian , X. & Yang , J. ( 2016 ) Robust a posteriori stress analysis for quadrature collocation approximations of nonlocal models via nonlocal gradients . Comput. Methods Appl. Mech. Eng. , 310 , 605 – 627 . Google Scholar CrossRef Search ADS Du , Q. & Tian , X. ( 2014 ) Asymptotically compatible schemes for peridynamics based on numerical quadratures . ASME 2014 International Mechanical Engineering Congress and Exposition (IMECE2014), vol. 1 , Paper No. IMECE2014–39620, pp . V001T01A058 . Du , Q. & Yang , J. ( 2016 ) Asymptotic compatible Fourier spectral approximations of nonlocal Allen–Cahn equations . SIAM J. Numer. Anal. , 54 , 1899 – 1919 . Google Scholar CrossRef Search ADS Du , Q. & Yang , J. ( 2017 ) Fast and accurate implementation of Fourier spectral approximations of nonlocal diffusion operators and its applications . J. Comput. Phys. , 332 , 118 – 134 . Google Scholar CrossRef Search ADS Du , Q. , Yang , J. & Zhou , Z. ( 2017b ) Analysis of a nonlocal-in-time parabolic equation . Discrete Continuous Dyn. Syst. B. , 22 , 339 – 368 . Google Scholar CrossRef Search ADS Du , Q. & Zhou , K. ( 2011 ) Mathematical analysis for the peridynamic nonlocal continuum theory . Math. Model. Numer. Anal. , 45 , 217 – 234 . Google Scholar CrossRef Search ADS Gilboa , G. & Osher , S. ( 2008 ) Nonlocal operators with applications to image processing . Multiscale Model. Simul. , 7 , 1005 – 1028 . Google Scholar CrossRef Search ADS Kilic , B. & Madenci , E. ( 2010 ) Coupling of peridynamic theory and the finite element method . J. Mech. Mater. Struct. , 5 , 707 – 733 . Google Scholar CrossRef Search ADS Liu , W. K. , Chen , Y. , Jun , S. , Chen , J. , Belytschko , T. , Pan , C. , Uras , R. & Chang , C. ( 1996 ) Overview and applications of the reproducing kernel particle methods . Arch. Comput. Methods Eng. , 3 , 3 – 80 . Google Scholar CrossRef Search ADS Lou , Y. , Zhang , X. , Osher , S. & Bertozzi , A. ( 2010 ) Image recovery via nonlocal operators . J. Sci. Comput. , 42 , 185 – 197 . Google Scholar CrossRef Search ADS Macek , R. & Silling , S. A. ( 2007 ) Peridynamics via finite element analysis . Finite Elem. Anal. Des. , 43 , 1169 – 1178 . Google Scholar CrossRef Search ADS Mengesha , T. & Du , Q. , ( 2014 ) Nonlocal constrained value problems for a linear peridynamic Navier equation . J. Elast. , 116 , 27 – 51 . Google Scholar CrossRef Search ADS Palatucci , G. , Savin , O. & Valdinoci , E. ( 2012 ) Local and global minimizers for a variational energy involving a fractional norm . Ann. Mat. Pura Appl. , 192 , 673 – 718 . Google Scholar CrossRef Search ADS Silling , S. A. ( 2000 ) Reformulation of elasticity theory for discontinuities and long-range forces . J. Mech. Phys. Solids , 48 , 175 – 209 . Google Scholar CrossRef Search ADS Silling , S. A. & Askari , E. ( 2005 ) A meshfree method based on the peridynamic model of solid mechanics . Comput. Struct. , 83 , 1526 – 1535 . Google Scholar CrossRef Search ADS Silling , S. A. & Lehoucq , R. B. ( 2010 ) Peridynamic theory of solid mechanics . Adv. Appl. Mech. , 44 , 73 – 168 . Google Scholar CrossRef Search ADS Silling , S. A. , Weckner , O. , Askari , E. & Bobaru , F. ( 2010a ) Crack nucleation in a peridynamic solid . Int. J. Fract. , 162 , 219 – 227 . Google Scholar CrossRef Search ADS Silling , S. A. , Zimmermann , M. & Abeyaratne , R. ( 2010b ) Deformation of a peridynamic bar . J. Elast. , 73 , 173 – 190 . Google Scholar CrossRef Search ADS Tadmor , E. & Tan , C. ( 2014 ) Critical thresholds in flocking hydrodynamics with non-local alignment . Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. , 372 : 20130401 . Tao , Y. , Tian , X. & Du , Q. ( 2017 ) Nonlocal diffusion and peridynamic models with Neumann type constraints and their numerical approximations . Appl. Math. Comput. , 305 , 282 – 298 . Tian , H. , Ju , L. & Du , Q. ( 2018 ) A conservative nonlocal convection-diffusion model and asymptotically compatible finite difference discretization . Comput. Methods Appl. Mech. Eng. , 320 , 46 – 67 . Google Scholar CrossRef Search ADS Tian , X. & Du , Q. ( 2013 ) Analysis and comparison of different approximations to nonlocal diffusion and linear peridynamic equations . SIAM J. Numer. Anal. , 51 , 3458 – 3482 . Google Scholar CrossRef Search ADS Tian , X. & Du . Q. ( 2014 ) Asymptotically compatible schemes and applications to robust discretization of nonlocal models . SIAM J. Numer. Anal. , 52 , 1641 – 1665 . Google Scholar CrossRef Search ADS Wang , L. , Xu , J. & Wang , J. ( 2016 ) The Green’s functions for peridynamic non-local diffusion . Proc. R. Soc. A , 472 , 20160185 . Wang , L. , Xu , J. & Wang , J. ( 2017 ) Static and dynamic Green’s functions in peridynamics . J. Elast. , 126 , 95 – 125 . Google Scholar CrossRef Search ADS Weckner , O. , Brunk , G. , Epton , M. , Silling , S. & Askari , E. ( 2009 ) Green’s functions in non-local three-dimensional linear elasticity. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. , 53 , 705 – 728 . Zhou , K. & Du , Q. ( 2010 ) Mathematical and numerical analysis of linear peridynamic models with nonlocal boundary conditions . SIAM J. Numer. Anal. , 48 , 1759 – 1780 . Google Scholar CrossRef Search ADS © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) For permissions, please e-mail: journals. permissions@oup.com

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Apr 5, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations