Log In Sign Up

Asymptotically compatible reproducing kernel collocation and meshfree integration for nonlocal diffusion

Reproducing kernel (RK) approximations are meshfree methods that construct shape functions from sets of scattered data. We present an asymptotically compatible (AC) RK collocation method for nonlocal diffusion models with Dirichlet boundary condition. The scheme is shown to be convergent to both nonlocal diffusion and its corresponding local limit as nonlocal interaction vanishes. The analysis is carried out on a special family of rectilinear Cartesian grids for linear RK method with designed kernel support. The key idea for the stability of the RK collocation scheme is to compare the collocation scheme with the standard Galerkin scheme which is stable. In addition, there is a large computational cost for assembling the stiffness matrix of the nonlocal problem because high order Gaussian quadrature is usually needed to evaluate the integral. We thus provide a remedy to the problem by introducing a quasi-discrete nonlocal diffusion operator for which no numerical quadrature is further needed after applying the RK collocation scheme. The quasi-discrete nonlocal diffusion operator combined with RK collocation is shown to be convergent to the correct local diffusion problem by taking the limits of nonlocal interaction and spatial resolution simultaneously. The theoretical results are then validated with numerical experiments. We additionally illustrate a connection between the proposed technique and an existing optimization based approach based on generalized moving least squares (GMLS).


page 1

page 2

page 3

page 4


Asymptotically compatible reproducing kernel collocation and meshfree integration for the peridynamic Navier equation

In this work, we study the reproducing kernel (RK) collocation method fo...

A Scharfetter-Gummerl stabilization scheme for HDG approximations of convection-diffusion problems

We present a Scharfetter-Gummel (SG) stabilization scheme for high-order...

Maximum Principle Preserving Finite Difference Scheme for 1-D Nonlocal-to-Local Diffusion Problems

In a recent paper (see [7]), a quasi-nonlocal coupling method was introd...

Analyses of the contour integral method for time fractional subdiffusion-normal transport equation

In this work, we theoretically and numerically discuss the time fraction...

OBMeshfree: An optimization-based meshfree solver for nonlocal diffusion and peridynamics models

We present OBMeshfree, an Optimization-Based Meshfree solver for compact...

1 Introduction

This work is motivated by the study of numerical solutions to linear nonlocal models and their local limits. Peridynamics (PD) is a nonlocal theory of continuum mechanics [32]. Unlike the classical local theory, PD models are formulated using spatial integration instead of differentiation, making them well-suited for describing discontinuities such as those present in fracture, material separation and failure. PD has been applied to hydraulic-fracture propagation problems [28], crack branching [4], damage progression in multi-layered glass [17] and others. Linear PD models also share similarities with nonlocal diffusion model [11]. Rigorous mathematical analysis and a variety of numerical methods have been developed for PD and nonlocal diffusion models [3, 7, 10, 11, 12, 13, 26, 31, 33, 36, 37]. Nonlocal models introduce a nonlocal length scale , called the horizon, which takes into account nonlocal interactions. As goes to nonlocal interactions vanish and nonlocal models recover their local equivalents, provided the limit is well-defined. It is a common practice to couple with the mesh size in engineering applications, but some standard numerical methods may converge to the wrong local limit [35]. This particular limit is of fundamental practical importance, as it leads to banded linear systems amenable to traditional preconditioning techniques.

A mathematical framework of convergence is established for PD and nonlocal diffusion models in [11, 25, 26] and asymptotically compatible (AC) discretization is introduced in [35, 36]. The AC scheme allows the numerical solution of nonlocal equations to converge to both the nonlocal solutions for a fixed and also their local limits as goes to zero, independent of the mesh size . The study of AC schemes has since then been developed for various numerical methods and model problems [5, 8, 15, 14, 20, 34, 37, 39]. Finite element methods (FEM) for nonlocal equations are studied in [35, 36] and FEM with subspaces containing piecewise linear functions are shown to be AC. However, applying FEM to nonlocal problems is computationally prohibitive because the variational formulation of nonlocal equations involves both a double integral and a geometrically costly mesh intersection calculation [7, 16]. Further, the nonlocal kernels in PD models are often singular, which adds more complexity to the computation. Finite difference methods (FDM) do not require the evaluation of a double integral but require uniform grids to obtain both AC and discrete maximum principle at the same time [14]. A meshfree discretization [33] of PD equations is widely used in engineering applications due to its simplicity. This meshfree method uses a set of particles in the domain, each with a known volume, and assumes constant fields in each nodal element. This method, however, suffers from large integration error leading to low order of convergence. Works [30, 38] have been devoted to improve the integration error and a reproducing kernel (RK) collocation approach can increase the order of convergence [29]. However, the robustness of this meshfree method needs further investigation.

The first motivation of this work is to provide a convergence analysis of RK collocation method for nonlocal diffusion models. Stability of collocation methods on integral equations is not a trivial task, due to the lack of a discrete maximum principle. A helpful view is to compare collocation schemes with Galerkin schemes [1, 2, 9], for which stability comes naturally. In this work, we use the Fourier approach [9]

and demonstrate that the Fourier symbol of the RK collocation scheme with linear interpolation order and suitable choice of

for nonlocal diffusion can be bounded below by that of the standard Galerkin scheme on Cartesian grids. Consequently, we show that the collocation scheme is stable because the standard Galerkin approximation is uniformly stable (i.e. the stability constant does not depends on ). Consistency is established using the approximation properties of the RK approximation. Therefore the proposed RK collocation method on nonlocal diffusion is AC.

Although the collocation scheme requires only a single integration to be performed for the evaluation of each element in the stiffness matrix, it is still quite expensive to do so, particularly for models with singular kernels. In practice, high-order Gauss quadrature rules are used to evaluate the integral [7, 29]. In [29], two Gauss quadrature schemes are investigated. The first one has a background mesh for integration and the second scheme places quadrature points inside the the horizon. Therefore, the second goal of this work is to develop a practical numerical method for nonlocal models. To this end, we introduce a quasi-discrete nonlocal diffusion operator which replaces the integral with a finite summation of quadrature points inside the horizon. We utilize the RK technique to calculate the quadrature weights. The quasi-discrete nonlocal diffusion operator discretized with RK collocation not only saves computational cost but also enables the use of bond breaking for fracture problems [37]. Moreover, we provide convergence analysis to the correct local limit as of the RK collocation scheme with meshfree integration. A similar technique has been proposed in [37] utilizing an optimization construction which admits interpretation as a generalized moving least squares (MLS) process. It is well known that RK and MLS shape functions are equivalent, up to a rescaling of the weighting function and for particular reproducing spaces [6]. We will show that the construction of quadrature weights using the RK technique is similarly equivalent under certain conditions to this generalized MLS approach [37], and therefore the stability proof provided here applies equally to this second class of schemes which currently lack a proof of stability. This unifies existing work in the literature using both RK [29] and MLS [37] as a framework to develop AC particle-based schemes.

This paper is organized as follows. In Section 2, we introduce the nonlocal diffusion model equations with Dirichlet boundary conditions. In Section 3, we present the RK collocation and work with the RK method with linear interpolation order with special choices of the RK support sizes. Section 4 discusses the convergence of the RK collocation method to both the nonlocal diffusion equation for fixed and the local diffusion equation as goes to zero so that the RK collocation scheme is AC. Then, a quasi-discrete nonlocal diffusion operator is developed in Section 5 and its convergence analysis is presented in Section 6. Section 7 gives numerical examples to complement our theoretical analysis. Finally, conclusions are made in Section 8.

2 Nonlocal diffusion operator and model equation

We use the following notation throughout the paper. d is a positive integer denoting spatial dimension. A generic point is expressed as . A multi-index is a collection of d nonnegative integers, and its length is . For a given , we write . Let be a bounded, open domain. The corresponding interaction domain is then defined as

and let . Following the same notations as in [11], we define the nonlocal diffusion operator , for a given , as


where is the nonlocal length and is the nonlocal diffusion kernel which is nonnegative and symmetric, i.e., . Let us consider a nonlocal diffusion problem with homogeneous Dirichlet volumetric constraint,


In this work, we study the kernels of radial type, i.e., . We denote the following scaling of the kernel


where is compactly supported in (a unit ball about ). We further assume

is a non-increasing function and has a bounded second-order moment, i.e.


The local limit of is denoted as when . We are interested in particular cases where , such that Eq. 2 goes to


We proceed to define some functional spaces. Define a space on with zero volumetric constraint on ,

The natural energy space associated with Eq. 2 is given by

The nonlocal diffusion problem as described in Eq. 2 is well-posed and uniformly stable. We give the following the result without proof and the details can be found [25].

Lemma 1

(Uniform stability) Assume for some . The bilinear form is an inner product and for any , we have

where is a constant that only depends on and .

3 RK collocation method

Since functions in the solution space are zero on the boundary layer , we can extend functions in trivially to and work on the whole space instead. We first introduce some notations. Define to be a rectilinear Cartesian grid on , namely

where and consists of discretization parameters in each dimension and denotes component-wise multiplication, i.e.,

Sometime we also write the -th component of by , which is equal to by definition. We introduce a component-wise division symbol :

Note that the grid size can be different for different . For instance, in two dimension, rectangular grids are allowed. Nonetheless, we assume the grid is quasi-uniform such that can also be written as



being a fixed vector with the maximum component to be 1. For any continuous function

, define the restriction to by


and the restriction to () as


For any sequence on , the RK interpolant operator is given as

where is the RK basis function to be introduced shortly. Denote by the trial space equipped with the RK basis on , i.e., . Let

be the interpolation projector from the space of continuous functions on to the trial space .

We proceed to recall the construction of the RK basis function. The RK approximation [24] of on is formulated as:


where is the correction function, is the nodal coefficient, and

is the kernel function defined as the tensor product of kernel functions in each dimension with support

, i.e.


where is the kernel function in the -th dimension, is the support size for and is called the window function. In this work, we use the cubic B-spline function as the window function, i.e.,


The correction function in Eq. 9 is defined as


where the vector consists of the set of monomial basis functions of order ,


is a vector containing correction function coefficients and can be obtained by satisfying the -th order polynomial reproduction condition,


Substitute Eq. 12 into Eq. 14 and obtain



where is the moment matrix and is formulated as


Each entry of the matrix is a moment given by


where , , and is the moment in the -th dimension given by


Solve the system of equations as in Eq. 15 and obtain the correction function coefficients as


Please refer to [18] for the necessary conditions of the solvability of the system Eq. 15. Finally, by substituting Eq. 19 and Eq. 12 into Eq. 9, the RK approximation of is obtained as

where is the RK basis function,

In the rest of the work, we assume the reproducing condition (14) is satisfied with , with which we call our method the linear RK approximation and the RK basis function is referred to as the linear RK basis. If the RK support size is chosen properly, it is shown that the linear RK approximation has synchronized convergence [22], a phenomenon where the convergence rates of higher-order error norms and lower-order error norms are of the same order.

Remark 1

Choose the RK support as , where is an even number, the linear RK basis function becomes a rescaling of the cubic B-spline function ([21]) and the RK approximation has synchronized convergence for and error norm ([22]). For the simplicity of presentation, we assume in the paper but the analysis also works for general even number .

Let be the parameter described in Remark 1, then the linear RK basis function can be written as ([21])


Another consequence of this choice of support size is that the one-dimensional moments up to the third order are independent of , and more precisely


for . From the one-dimensional moment, we can derive useful properties of the multi-dimensional moment which are summarized in the following lemma.

Lemma 2

The multi-dimensional moments satisfy the following properties,

  1. [label=()]

  2. and for ,

  3. for and ,


By writing out the multi-index and from Eq. 17 and Eq. 21, the desired properties follow.

Now we use the above discussed RK approximation and collocate the nonlocal diffusion equation on the grid . The RK collocation scheme is formulated as follows. Find a function such that


where is defined as

Alternatively, Eq. 22 can also be written as

where is the restriction operator given in Eq. 8.

Remark 2

With the assumption in Remark 1, the RK basis has support size of in the -th dimension. So the support of is not fully contained in but in a larger domain given as

4 Convergence analysis of the RK collocation method

In this section, we will show the convergence of the RK collocation scheme Eq. 22, which is also the method used in [29] without a convergence proof. The concern for convergence is that the numerical scheme should converge to the nonlocal problem for a fixed , and to the correct local problem as and grid size both go to zero. So the proposed RK collocation scheme is an AC scheme ([36]).

4.1 Stability of the RK collocation method

In this subsection, we provide the stability proof of our method. The key idea is to compare the RK collocation scheme with the Galerkin scheme using Fourier analysis. Similar strategies have been developed in [9].

First, define a norm in the space of sequences by


If a sequence is only defined for in a subset of , then one can always use zero extension for so that it is defined for all . Then without further explanation, is always understood as (23) with the zero extension being used. The main theorem in this subsection is now given as follows.

Theorem 1

(Stability I) For any and , we have

where is a constant that only depends on and .

The proof Theorem 1 is shown at the end of this subsection before two more lemmas are introduced. Let be the norm associated inner product, namely

For Fourier series defined on ,



The nonlocal operator defines two discrete sesquilinear forms:




The inner product in Eq. 25 is the standard inner product. Equation 25 defines a quadratic form corresponding to the Galerkin method, meanwhile, the quadratic form Eq. 26 corresponds to the collocation method. Moreover, the stiffness matrix for the collocation scheme Eq. 22 can be considered as a finite section of the infinite Toeplitz matrix induced by Eq. 26. Before applying the Fourier analysis to Eq. 25 and Eq. 26, we study the Fourier symbol of the nonlocal diffusion operator first. Take ,

is the Fourier transform of

defined by

The Fourier transform of the nonlocal diffusion operator is given as

where is the Fourier symbol of ,


More discussions on the spectral analysis of the nonlocal diffusion operator can be found in [15]. From Eq. 27, it is obvious that is real and non-negative. Now, we given a comparison of the two quadratic forms Eq. 25 and Eq. 26 using Fourier analysis.

Lemma 3

Let and be the Fourier series of the sequences respectively. Then

  1. [label=()]

  2. ,

  3. ,

  4. , for independent of and ,

and and are given by



The inverse Fourier transform of is given by

By Parseval’s identity, we arrive at

where we have used Eq. 20 and the Fourier transform of the cubic B-spline function Eq. 11 given by

Hence, we have the Fourier transform of the RK basis function as

Therefore, by Eq. 24 we obtain the Galerkin form Eq. 25 as

Next, following the same procedure, the collocation matrix may be expressed as

Therefore, the collocation form Eq. 26 is written as

This finishes the proof of 2.

With 1 and 2, it is easy to see 3 by the fact that is non-negative and .

Remark 3

defined in (29) is in fact the same as the Fourier symbol of the Galerkin scheme with basis function defined by the tensor product of linear B-splines. This shows the equivalence of some special RK collocation schemes with certain non-standard Galerkin schemes.

We need the following lemma that says the norm defined through the RK basis and discrete norm are equivalent.

Lemma 4

The following two norms are equivalent, i.e., there exist two constants independent of , such that


First, from Parsevel’s identity, we can write the norm as

Then, similar to the proof of Lemma 3 1, by replacing the nonlocal diffusion operator with the identity operator, we obtain

where is continuous and strictly positive,

Thus, is bounded above and below on . Therefore, we complete the proof.

Proof (Proof of of Theorem 1)

For all sequences , we derive via the Cauchy-Schwartz inequality and Lemma 4


Finally, for we may by definition write , and thus we have

The first line is a result of Eq. 30 and the second line is by definition of . Lemma 3 3 shows the third line and the fourth line is the stability given by Lemma 1 since .

4.2 Consistency analysis the RK collocation method

In this section, we discuss uniform consistency of the RK collocation method on the nonlocal diffusion models. The truncation error has a uniform bound independent of the nonlocal scaling parameter . Combining the stability result in Section 4.1 and the truncation error analysis to be presented shortly, we show that the RK collocation method is convergent. The asymptotic compatibility of the method is based on two important properties – the quadratic exactness of the scheme Lemma 5 and the synchronized convergence property of the linear RK approximation. Synchronized convergence means the convergence rates measured by higher-order error norms are of the same order as the ones measured by lower-order error norms. The linear RK approximation has the synchronized convergence is a long-known fact [22] and we provide the statement of it in Lemma 6 without proof.

Lemma 5

(Quadratic Exactness) For quadratic polynomials in as , , we have


Let and . From the definition of the nonlocal diffusion operator and the first-order polynomial reproduction property of the RK approximation, it is easy to see that

Then it is sufficient to show the cases when . We observe that for quadratic polynomials,


There are only two cases for . The first case is that there exist and such that . In this case