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 wellsuited for describing discontinuities such as those present in fracture, material separation and failure. PD has been applied to hydraulicfracture propagation problems [28], crack branching [4], damage progression in multilayered 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 welldefined. 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, highorder 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 quasidiscrete 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 quasidiscrete 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 particlebased 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 quasidiscrete 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 multiindex 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
(1) 
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,
(2) 
In this work, we study the kernels of radial type, i.e., . We denote the following scaling of the kernel
(3) 
where is compactly supported in (a unit ball about ). We further assume
is a nonincreasing function and has a bounded secondorder moment, i.e.
(4) 
The local limit of is denoted as when . We are interested in particular cases where , such that Eq. 2 goes to
(5) 
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 wellposed 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 componentwise multiplication, i.e.,
Sometime we also write the th component of by , which is equal to by definition. We introduce a componentwise 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 quasiuniform such that can also be written as
(6) 
with
being a fixed vector with the maximum component to be 1. For any continuous function
, define the restriction to by(7) 
and the restriction to () as
(8) 
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:
(9) 
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.(10) 
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 Bspline function as the window function, i.e.,
(11) 
The correction function in Eq. 9 is defined as
(12) 
where the vector consists of the set of monomial basis functions of order ,
(13) 
is a vector containing correction function coefficients and can be obtained by satisfying the th order polynomial reproduction condition,
(14) 
Substitute Eq. 12 into Eq. 14 and obtain
Equivalently,
(15) 
where is the moment matrix and is formulated as
(16) 
Each entry of the matrix is a moment given by
(17) 
where , , and is the moment in the th dimension given by
(18) 
Solve the system of equations as in Eq. 15 and obtain the correction function coefficients as
(19) 
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 higherorder error norms and lowerorder 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 Bspline 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])
(20) 
Another consequence of this choice of support size is that the onedimensional moments up to the third order are independent of , and more precisely
(21) 
for . From the onedimensional moment, we can derive useful properties of the multidimensional moment which are summarized in the following lemma.
Lemma 2
The multidimensional moments satisfy the following properties,

[label=()]

and for ,

for and ,
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
(22) 
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
(23) 
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 ,
(24) 
where
The nonlocal operator defines two discrete sesquilinear forms:
(25) 
and
(26) 
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 byThe Fourier transform of the nonlocal diffusion operator is given as
where is the Fourier symbol of ,
(27) 
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 nonnegative. 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

[label=()]

,

,

, for independent of and ,
and and are given by
(28) 
(29) 
Proof
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 Bspline 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
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 Bsplines. This shows the equivalence of some special RK collocation schemes with certain nonstandard 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
Proof
Proof (Proof of of Theorem 1)
For all sequences , we derive via the CauchySchwartz inequality and Lemma 4
(30)  
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 higherorder error norms are of the same order as the ones measured by lowerorder error norms. The linear RK approximation has the synchronized convergence is a longknown 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
Proof
Let and . From the definition of the nonlocal diffusion operator and the firstorder 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,
(31) 
There are only two cases for . The first case is that there exist and such that . In this case