A conforming DG method for linear nonlocal models with integrable kernels

TL;DR

This paper introduces a new discontinuous Galerkin method for solving linear nonlocal models with integrable kernels, achieving optimal convergence in 1D and near-optimal in 2D, ensuring efficiency and accuracy.

Contribution

It proposes a novel conforming DG method tailored for nonlocal problems with integrable kernels, demonstrating asymptotic compatibility and optimal convergence rates.

Findings

Method is asymptotically compatible.

Achieves optimal convergence in 1D.

Attains nearly optimal convergence in 2D.

Abstract

Numerical solution of nonlocal constrained value problems with integrable kernels are considered. These nonlocal problems arise in nonlocal mechanics and nonlocal diffusion. The structure of the true solution to the problem is analyzed first. The analysis leads naturally to a new kind of discontinuous Galerkin method that efficiently solve the problem numerically. This method is shown to be asymptotically compatible. Moreover, it has optimal convergence rate for one dimensional case under very weak assumptions, and almost optimal convergence rate for two dimensional case under mild assumptions.

Tables

Table 1. Table 1: Example 1: Errors in L 2 superscript 𝐿 2 L^{2} and H 1 superscript 𝐻 1 H^{1} norms, corresponding convergence rates and spectral condition numbers on uniform meshes, Cond is abbreviation for the spectral condition number

$\delta /h$	4	8	16	32	64	128
$\Vert u-u_h\Vert$	7.45e-4	1.86e-4	4.66e-5	1.16e-5	2.91e-6	7.28e-7
Rate	–	2.0000	2.0000	2.0000	2.0000	2.0000
$\Vert u^{\prime }-u_h^{\prime }\Vert$	5.77e-2	2.89e-2	1.44e-2	7.22e-3	3.61e-3	1.80e-3
Rate	–	1.0000	1.0000	1.0000	1.0000	1.0000
Cond	5.8875	6.7743	6.9926	7.0294	7.0282	7.0225

Table 2. Table 2: Example 1: Errors in L 2 superscript 𝐿 2 L^{2} and H 1 superscript 𝐻 1 H^{1} norms, corresponding convergence rates and spectral condition numbers on non-uniform meshes ( 5.3 )

$\delta /h$	4	8	16	32	64
$\Vert u-u_h\Vert$	8.04e-4	2.50e-4	6.00e-5	1.24e-5	3.14e-6
Rate	–	1.6827	2.0629	2.2710	1.9846
$\Vert u^{\prime }-u_h^{\prime }\Vert$	5.84e-2	2.95e-2	1.48e-2	7.30e-3	3.65e-3
Rate	–	0.9848	0.9965	1.0179	0.9983
Cond	5.9830	6.7507	7.0829	7.0786	7.0646

Table 3. Table 3: Example 2: Errors in L 2 superscript 𝐿 2 L^{2} and H 1 superscript 𝐻 1 H^{1} norms, corresponding convergence rates and spectral condition numbers on uniform meshes

$\delta /h$	4	8	16	32	64
$\Vert u-u_h\Vert$	7.54e-3	1.79e-3	4.38e-4	1.08e-4	2.69e-5
Rate	–	2.0717	2.0349	2.0170	2.0084
$\Vert u^{\prime }-u_h^{\prime }\Vert$	4.90e-1	2.41e-1	1.19e-1	5.94e-2	2.97e-2
Rate	–	1.0259	1.0114	1.0053	1.0026
Cond	5.8875	6.7743	6.9926	7.0294	7.0282

Table 4. Table 4: Example 2: Errors in L 2 superscript 𝐿 2 L^{2} and H 1 superscript 𝐻 1 H^{1} norms, and convergence rates on non-uniform meshes ( 5.3 )

$\delta /h$	4	8	16	32	64
$\Vert u-u_h\Vert$	9.34e-3	3.41e-3	1.24e-3	4.26e-4	1.59e-4
Rate	–	1.4519	1.4611	1.5411	1.4222
$\Vert u^{\prime }-u_h^{\prime }\Vert$	3.95e-1	2.66e-1	2.13e-1	1.42e-1	1.21e-1
Rate	–	0.5712	0.3204	0.5815	0.3432

Table 5. Table 5: Example 2: Errors in L 2 superscript 𝐿 2 L^{2} and H 1 superscript 𝐻 1 H^{1} norms, and convergence rates on non-uniform meshes ( 5.3 ) with δ 𝛿 \delta and 1 − δ 1 𝛿 1-\delta as grids

$\delta /h$	4	8	16	32	64
$\Vert u-u_h\Vert$	6.98e-3	1.72e-3	4.10e-4	1.01e-4	2.49e-5
Rate	–	2.0229	2.0661	2.0243	2.0180
$\Vert u^{\prime }-u_h^{\prime }\Vert$	4.17e-1	2.14e-1	1.01e-1	4.74e-2	2.42e-2
Rate	–	0.9646	1.0781	1.0941	0.9671

Table 6. Table 6: Example 3: Errors in L 2 superscript 𝐿 2 L^{2} and H 1 superscript 𝐻 1 H^{1} norms, and convergence rates on uniform meshes

$\delta /h$	4	8	16	32	64
$\Vert u-u_h\Vert$	1.18e-2	2.74e-3	6.61e-4	1.62e-4	4.03e-5
Rate	–	2.1019	2.0525	2.0246	2.0116
$\Vert u^{\prime }-u_h^{\prime }\Vert$	7.20e-1	3.60e-1	1.79e-1	8.90e-2	4.44e-2
Rate	–	1.0018	1.0089	1.0056	1.0030

Table 7. Table 7: Example 3: Errors in L 2 superscript 𝐿 2 L^{2} and H 1 superscript 𝐻 1 H^{1} norms, and convergence rates on nonuniform meshes ( 5.3 )

$\delta /h$	4	8	16	32	64
$\Vert u-u_h\Vert$	1.11e-2	2.61e-3	6.08e-4	1.51e-4	3.76e-5
Rate	–	2.0891	2.0993	2.0143	1.9985
$\Vert u^{\prime }-u_h^{\prime }\Vert$	6.70e-1	3.32e-1	1.51e-1	7.47e-2	3.67e-2
Rate	–	1.0128	1.1331	1.0189	1.0261

Equations

$$\mathcal{L}u({\bf x})=\int_{\mathbb{R}^{d}}{\big{(}}u({\bf y})-u({\bf x}){\big{)}}\gamma({\bf x},{\bf y})d{\bf y}\quad\forall{\bf x}\in\Omega,$$

$$\left\{\begin{array}[]{rl}-\mathcal{L}u({\bf x})&=b({\bf x})\>\mbox{on}\>\Omega,\\ u({\bf x})&=g({\bf x})\>\mbox{on}\>\Omega_{I},\end{array}\right.$$

$$\gamma (\mathbf{x},\mathbf{y})\ge {\gamma }_0\forall \mathbf{y}\in B_{\delta /2}(\mathbf{x}),\gamma (\mathbf{x},\mathbf{y})=0\forall \mathbf{y}\in (\Omega \cup {\Omega }_I)\setminus B_{\delta }(\mathbf{x}),$$

$$\gamma (\mathbf{x},\mathbf{y})=\tilde{\gamma }(|\mathbf{y}-\mathbf{x}|)\ge 0,{\int }_{{\mathbb{R}}^d}\tilde{\gamma }(|\mathbf{z}|)d\mathbf{z}=1.$$

$$\overline{b}({\bf x})=\left\{\begin{array}[]{ll}b({\bf x}),&{\bf x}\in\Omega,\\ g({\bf x}),&{\bf x}\in\Omega_{I},\end{array}\right.$$

$$\overline{u}(\mathbf{x})=u(\mathbf{x})-\overline{b}(\mathbf{x}).$$

$$-\mathcal{L}\overline{u}(\mathbf{x})$$

$$={\int }_{B_{\delta }(\mathbf{x})\cap \Omega }b(\mathbf{y})\gamma (\mathbf{y}-\mathbf{x})d\mathbf{y}+{\int }_{B_{\delta }(\mathbf{x})\cap {\Omega }_I}g(\mathbf{y})\gamma (\mathbf{y}-\mathbf{x})d\mathbf{y},\forall \mathbf{x}\in \Omega .$$

$$\left\{\begin{array}[]{rl}-\mathcal{L}\overline{u}({\bf x})&=f({\bf x}),\>\>\mbox{on}\>\Omega,\\ \overline{u}({\bf x})&=0,\qquad\mbox{on}\>\Omega_{I}.\end{array}\right.$$

$$\left\{\begin{array}[]{rl}-\mathcal{L}u({\bf x})&=b({\bf x}),\>\>\mbox{on}\>\Omega,\\ u({\bf x})&=0,\quad\>\>\>\mbox{on}\>\Omega_{I},\end{array}\right.$$

$$u(\mathbf{x})\in L^2(\Omega \cup {\Omega }_I).$$

$$u(\mathbf{x})=b(\mathbf{x})+{\int }_{B_{\delta }(\mathbf{x})}u(\mathbf{y})\gamma (\mathbf{y}-\mathbf{x})d\mathbf{y},\forall \mathbf{x}\in \Omega .$$

$${\int }_{B_{\delta }(\mathbf{x})}u(\mathbf{y})\gamma (\mathbf{y}-\mathbf{x})d\mathbf{y}$$

$$u(\mathbf{x})=b(\mathbf{x})+{\int }_{B_{\delta }(0)}u(\mathbf{x}+\mathbf{s})\gamma (\mathbf{s})d\mathbf{s}.$$

$$\frac{\partial u(\mathbf{x})}{\partial \mathbf{t}}=\frac{\partial b(\mathbf{x})}{\partial \mathbf{t}}+{\int }_{B_{\delta }(0)}\frac{\partial u(\mathbf{x}+\mathbf{s})}{\partial \mathbf{t}}\gamma (\mathbf{s})d\mathbf{s}.$$

$${\int }_{B_{\delta }(0)}\frac{\partial u(\mathbf{x}+\mathbf{s})}{\partial \mathbf{t}}\gamma (\mathbf{s})d\mathbf{s}$$

$$={\int }_{\partial B_{\delta }(0)}u(\mathbf{x}+\mathbf{s})\gamma (\mathbf{s})\mathbf{t}\cdot {\mathbf{n}}_{\mathbf{s}}d{\mathbf{S}}_{\mathbf{s}}-{\int }_{B_{\delta }(0)}u(\mathbf{x}+\mathbf{s})\mathbf{t}\cdot \nabla \gamma (\mathbf{s})d\mathbf{s}$$

$$=\tilde{\gamma }(\delta ){\int }_{\partial B_{\delta }(0)}u(\mathbf{x}+\mathbf{s})\mathbf{t}\cdot {\mathbf{n}}_{\mathbf{s}}d{\mathbf{S}}_{\mathbf{s}}-{\int }_{B_{\delta }(0)}u(\mathbf{x}+\mathbf{s}){\tilde{\gamma }}^{\prime }(|\mathbf{s}|)\mathbf{t}\cdot {\mathbf{n}}_{\mathbf{s}}d\mathbf{s}$$

$$=-{\int }_{B_{\delta }(0)}u(\mathbf{x}+\mathbf{s}){\tilde{\gamma }}^{\prime }(|\mathbf{s}|)\mathbf{t}\cdot {\mathbf{n}}_{\mathbf{s}}d\mathbf{s}.$$

$$\frac{\partial u(\mathbf{x})}{\partial \mathbf{t}}=\frac{\partial b(\mathbf{x})}{\partial \mathbf{t}}-{\int }_{B_{\delta }(0)}u(\mathbf{x}+\mathbf{s}){\tilde{\gamma }}^{\prime }(|\mathbf{s}|)\mathbf{t}\cdot {\mathbf{n}}_{\mathbf{s}}d\mathbf{s}.$$

$$\frac{\partial u({\mathbf{x}}_0)}{\partial \mathbf{t}}$$

$$-{\int }_{B_{\delta }(0)}u({\mathbf{x}}_0+\mathbf{s}){\tilde{\gamma }}^{\prime }(|\mathbf{s}|)\mathbf{t}\cdot {\mathbf{n}}_{\mathbf{s}}d\mathbf{s},$$

$${\int }_{\partial B_{\delta }(0)}u(\mathbf{x}+\mathbf{s})\mathbf{t}\cdot {\mathbf{n}}_{\mathbf{s}}d{\mathbf{S}}_{\mathbf{s}}$$

$$\frac{{\partial }^{2}u(\mathbf{x})}{\partial {\mathbf{t}}_1\partial {\mathbf{t}}_2}$$

$$\underset{h\to 0+}{\lim }u(\mathbf{x}+h\mathbf{t})=u(\mathbf{x}-),\forall \mathbf{t}\cdot {\mathbf{n}}_{\mathbf{x}}<0.$$

$$(u,v)_{\||}:={\Big{(}}\int_{\Omega\cup\Omega_{I}}\int_{\Omega\cup\Omega_{I}}{\big{(}}u({\bf y})-u({\bf x}){\big{)}}{\big{(}}v({\bf y})-v({\bf x}){\big{)}}\gamma({\bf x},{\bf y})d{\bf y}d{\bf x}{\Big{)}},$$

$$\Vert |u\Vert |:=(u,u{)}_{\Vert |}^{1/2},$$

$$V(\Omega \cup {\Omega }_I):=\{u\in L^2(\Omega \cup {\Omega }_I):\Vert |u\Vert |<\infty \},$$

$$V_{c,0}(\Omega \cup {\Omega }_I):=\{u\in V(\Omega \cup {\Omega }_I):u(\mathbf{x})=0\text{mbox}on{\Omega }_I\},$$

$$L_{c,0}^2(\Omega \cup {\Omega }_I):=\{u\in L^2(\Omega \cup {\Omega }_I):u(\mathbf{x})=0\text{mbox}on{\Omega }_I\}.$$

$$V_{c,0}^{pc}(\Omega \cup {\Omega }_I)=\{u\in V_{c,0}(\Omega \cup {\Omega }_I):u{|}_{\Omega }\in C(\Omega )\},$$

$$V_{c,0}^{pc,h}\subset V_{c,0}^{pc}(\Omega \cup {\Omega }_I).$$

$$\widehat{\phi}_{j}({\bf x})=\left\{\begin{array}[]{ll}\phi_{j}({\bf x})|_{\Omega_{h}},&{\bf x}\in\Omega_{h},\\ 0,&{\bf x}\in(\Omega\cup\Omega_{I})\setminus\Omega_{h}.\end{array}\right.$$

Full text

A conforming DG method for linear nonlocal models with integrable kernels

Qiang Du111Department of Applied Physics and Applied Mathematics, Columbia University, NY 10027, USA (email: qd2125@columbia.edu)  and Xiaobo Yin222Corresponding author, School of Mathematics and Statistics & Hubei Key Laboratory of Mathematical Sciences, Central China Normal University, Wuhan 430079, China (email: yinxb@mail.ccnu.edu.cn)

Abstract.   Numerical solution of nonlocal constrained value problems with integrable kernels are considered. These nonlocal problems arise in nonlocal mechanics and nonlocal diffusion. The structure of the true solution to the problem is analyzed first. The analysis leads naturally to a new kind of discontinuous Galerkin method that efficiently solve the problem numerically. This method is shown to be asymptotically compatible. Moreover, it has optimal convergence rate for one dimensional case under very weak assumptions, and almost optimal convergence rate for two dimensional case under mild assumptions.

Keywords. nonlocal diffusion; peridynamic model; nonlocal model; integrable kernel; discontinuous Galerkin; finite element; convergence analysis; condition number

**2010 Mathematics Subject Classification 82C21, 65R20, 74S05, 46N20, 45A05 **

1 Introduction

Nonlocal models have generated much interests in recent years [14]. For example, the peridynamic (PD) model proposed by Silling [26] is an integral-type nonlocal theory of continuum mechanics which provides an alternative setup to classical continuum mechanics based on PDEs. Since PD avoids the explicit use of spatial derivatives, it is especially effective for modeling problems involving discontinuities or other singularities in the deformation [6, 19, 22, 24, 29, 28]. The nonlocal diffusion (ND) model, described in [16] is another example of integro-differential equations. More recently, mathematical analysis of nonlocal models is also paid more attention, which could be found in [2, 3, 4, 9, 15, 17, 20, 21]. Meanwhile, to simulate nonlocal models, various numerical approximations have been proposed and studied, including finite difference, finite element, meshfree method, quadrature and particle-based methods [7, 10, 18, 22, 23, 25, 27, 31, 32, 33, 37].

Let $\Omega\subset\mathbb{R}^{d}$ denote a bounded, open and convex domain with Lipschitz continuous boundary. For $u({\bf x}):\Omega\rightarrow\mathbb{R}$, the nonlocal operator $\mathcal{L}$ on $u({\bf x})$ is defined as

[TABLE]

where the nonnegative symmetric mapping $\gamma({\bf x},{\bf y}):\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ is called a kernel. The operator $\mathcal{L}$ is regarded nonlocal since the value of $\mathcal{L}u$ at a point $\bf x$ involves information about $u$ at points ${\bf y}\neq{\bf x}$. In this article, we consider the following nonlocal Dirichlet volume-constrained diffusion problem

[TABLE]

where $\Omega_{I}=\{{\bf y}\in\mathbb{R}^{d}\setminus\Omega:\,\mbox{dist}({\bf y},\partial\Omega)<\delta\}$ with $\delta$ a constant called horizon parameter, $b({\bf x})\in L^{2}(\Omega)$ and $g({\bf x})\in L^{2}(\Omega_{I})$ are given functions. Volume constraints are natural extensions, to the nonlocal case, of boundary conditions for differential equations. Nonlocal versions of Neumann and Robin boundary conditions are also naturally defined [16].

We assume that the nonnegative symmetric kernel $\gamma({\bf x},{\bf y})$ satisfies, there exists a positive constant $\gamma_{0}$, for all ${\bf x}\in\Omega\cup\Omega_{I}$,

[TABLE]

where $B_{\delta}({\bf x}):=\{{\bf y}\in\Omega\cup\Omega_{I}:|{\bf y}-{\bf x}|\leq\delta\}$. Obviously, (1.3) implies that although interactions are nonlocal, they are limited to a ball of radius $\delta$. A few important classes of kernels are considered in [16]. Of particular interests here is a special choice of $\gamma$ being a radial function of ${\bf x}-{\bf y}$ (which also makes the kernel translation invariant). Such a case has also been studied earlier in [4] where

[TABLE]

This condition on $\gamma$ implies that $\mathcal{L}$ is a bounded mapping from $L^{2}(\mathbb{R}^{d})$ to $L^{2}(\mathbb{R}^{d})$. As we will discuss here, even though for smooth enough $b({\bf x})$ in $\Omega$, unlike the classical local PDE boundary value problems, the solution $u({\bf x})$ is possibly discontinuous across $\partial\Omega$ which makes the numerical solution to the corresponding nonlocal problem challenging.

An intuitive idea to overcome the possible loss of continuity is to use discontinuous Galerkin (DG) methods. The latter are in fact conforming, which is in stark contrast to DG methods for the discretization of second order elliptic partial differential equations for which they are nonconforming [5]. While nonconforming DG has also been studied recently for nonlocal models [12], if the structure of the solution could be studied carefully, a well designed conforming DG method could be a more competitive option as it could lessen the cost of computation. In this work, we propose a new kind of conforming DG method to approximate the nonlocal problem (1.2) with kernels satisfying (1.3) and (1.4) where the key is to adopt a hybrid version of continuous elements with DG at the boundary.

The paper is organized as follows. In Section 2, the structure of the solution for the given problem is analyzed, which is a generalization of the results in [30]. We also convert the original inhomogeneous problem (1.2) into the homogeneous problem (2.3) whose solution is smoother, so we just need to discuss a smoother homogeneous problem (2.4). Based on the structure of the solution, in Section 3 we propose a new DG method which solves the problem (2.4) efficiently. Convergence analysis and condition number estimates along with asymptotic compatibility for the method are given in Section 4. Results of numerical experiments are reported in Section 5.

2 The structure of the solution

To design more efficient numerical discretization, we first present some regularity studies on the nonlocal constrained value problem. We recall first some one dimensional results presented in [30]: when using peridynamic theory to model the elasticity on $\mathbb{R}=(-\infty,\infty)$, the displacement field $u$ has the same smoothness as the body force field $b$. In addition, any discontinuity in the kernel $\gamma$ (or in one of its derivatives) has some additional effect on the smoothness of $u$. For a peridynamic bar, suppose that $b$ has a discontinuity in its $N$th derivative at some $x=x_{b}$, and $\gamma$ has a discontinuity in its $L$th derivative at some $x=x_{c}$, then $u$ has a discontinuity in its $(N+nL+n)$th derivative at $x=x_{b}+nx_{c},n=1,2,\cdots$. This propagation of discontinuities is illustrated schematically in [30, Figure 3]. Roughly speaking, the smoothness of $u$ increases as one moves away from the location where the solution is discontinuous due to the discontinuity of the body force field $b$. These types of step-wise improved regularity associated with a finite horizon parameter have also been observed for nonlocal initial value problems of nonlocal-in-time dynamic systems in [13].

2.1 The structure of the solution for general dimensional cases on bounded domains

Recall for the 1D case in an unbounded domain, the regularity of the solution for a nonlocal problem is affected by both the right hand side function and the kernel function, assuming good behavior of the solution at infinities. In this subsection we consider the effect of these two sources on the regularity of the solution for general multidimensional constrained value problem on a bounded domain.

First, let us present a result to reduce the problem (1.2) which we are concerned with to be a problem with a homogeneous nonlocal constraint. Denoted by

[TABLE]

and

[TABLE]

Then the nonlocal operator $-\mathcal{L}$ on $\overline{u}({\bf x})$ is

[TABLE]

Thus, $\overline{u}({\bf x})$ is the solution of the following homogeneous nonlocal problem

[TABLE]

Due to $\overline{b}({\bf y})\in L^{2}(\Omega\cup\Omega_{I})$, $\gamma({\bf s})\in L^{2}(\mathbb{R}^{d})$, and the fact that convolution of functions in dual $L^{p}(\mathbb{R}^{d})$-spaces is continuous, we know that $f({\bf x})\in C(\Omega)$.

The problem (1.2) is equivalent to the problem (2.3) which has a homogeneous nonlocal constraint. That is, we just need to study the following homogeneous nonlocal problem:

[TABLE]

with $b({\bf x})\in C(\Omega)$, $\gamma({\bf x},{\bf y})$ satisfying (1.3) and (1.4). We will show that if $\gamma({\bf x},{\bf y})$ satisfies some mild assumptions, the results recalled earlier for the one dimensional case can be generalized to multidimensional case on a bounded domain.

2.2 Discontinuities due to the right hand side function

Theorem 2.1**.**

If $\gamma({\bf x},{\bf y})$ satisfies (1.3) and (1.4), and $b({\bf x})\in C(\Omega)$, the solution of (2.4) is continuous in $\Omega$, i.e., $u({\bf x})\in C(\Omega)$.

Proof.

Since $\gamma({\bf x},{\bf y})$ satisfies (1.3) and (1.4), we easily see that

[TABLE]

From (2.4), we have

[TABLE]

Since $u({\bf y})\in L^{2}(\Omega\cup\Omega_{I})$ and $\gamma({\bf s})\in L^{2}(\mathbb{R}^{d})$, we have that

[TABLE]

is continuous for ${\bf x}$ in $\Omega$. Together with the condition $b({\bf x})\in C(\Omega)$, we complete the proof. ∎

It is obvious from (2.5) that the discontinuity at a point $\bf x$ of the right hand side $b$ will lead to the discontinuity at the same point of the solution $u$. With an additional assumption on the kernel, we can bootstrap a higher order regularity result as follows.

Theorem 2.2**.**

Suppose that $\gamma({\bf x},{\bf y})$ satisfies (1.3) and (1.4). If $\tilde{\gamma}(r)\in C^{1}(0,\delta)$, $\tilde{\gamma}(\delta)=0$ and $b({\bf x})\in C^{1}(\Omega)$, then $u({\bf x})\in C^{1}(\Omega)$.

Proof.

From (2.5), we know

[TABLE]

For any unit vector ${\bf t}$, we take the corresponding directional derivative for (2.6), so

[TABLE]

Since $\tilde{\gamma}(r)\in C^{1}(0,\delta)$ and $\tilde{\gamma}(\delta)=0$, we have

[TABLE]

Thus,

[TABLE]

Since the convolution of functions in dual $L^{p}(\mathbb{R}^{d})$-spaces is continuous, the second term in the right hand side of (2.7) is continuous with respect to ${\bf x}$. Together with the condition $b({\bf x})\in C^{1}(\Omega)$, we complete the proof. ∎

From (2.7) we get that under the assumptions of Theorem 2.2, if the first derivative of the right hand side $b$ is discontinuous at some point $\bf x$, the first derivative of the solution $u$ will be discontinuous at the same point. However, if we may have both the condition $\tilde{\gamma}(\delta)>0$ and the condition $b({\bf x})\notin C^{1}(\Omega)$, the conclusion $u({\bf x})\in C^{1}(\Omega)$ may still hold (see Example 1, where in fact $u({\bf x})\in C^{\infty}(\Omega)$). This is not contradicting as the sum of two discontinuous function could be continuous, and even infinitely differentiable.

2.3 Propagation of discontinuities due to the kernel

For the convenience of discussion, let us denote by $\Omega_{1}=\{{\bf x}\in\Omega:\mbox{dist}(x,\partial\Omega)>\delta\}$ and $\Omega_{2}=\Omega\setminus\overline{\Omega_{1}}$. In fact, the significance of the result Theorem 2.1 is twofold. First, it indicates the smoothness of $u({\bf x})$ is the same as $b({\bf x})$ for general multidimensional case on a bounded domain. Second, it reveals the possible propagation of discontinuities due to the kernel. Although $b({\bf x})\in C(\Omega)$, it might be discontinuous across $\partial\Omega$. So does $u({\bf x})$, and this discontinuity will be propagated to those points on $\partial\Omega_{1}$, which are $\delta$ distance from $\partial\Omega$ onto one order higher derivatives by Theorem 2.1, which is consistent with the conclusion for 1D case on $\mathbb{R}$ with $N=0$ and $L=0$. Similar argument can be given for Theorem 2.2 which is consistent with the conclusion for 1D case on $\mathbb{R}$ with $N=0$ and $L=1$. This bootstrap procedure could be repeated again and again, and the corresponding results for general $N$ and $L$ then follow.

Let us now emphasize on the importance and necessity for the smoothness of the kernel function. For instance, in Theorem 2.2 $\tilde{\gamma}(\delta)=0$ is required such that $\tilde{\gamma}$ has a discontinuity in its first (but not zeroth) derivative at $x=\delta$. If $\tilde{\gamma}(\delta)>0$, then from the proof of the Theorem 2.2, we see that $u({\bf x})\in C^{1}(\Omega)$ may not hold. In fact, for all ${\bf x}_{0}\in\partial\Omega_{1}$ and any unit vector ${\bf t}$, since

[TABLE]

if $\tilde{\gamma}(\delta)>0$ (that is $\tilde{\gamma}$ has a discontinuity in its zeroth derivative) and $b$ (thus $u$) is discontinuous across $\partial\Omega$, then the term

[TABLE]

is likely to be discontinuous across ${\bf x}_{0}$. If so, $\frac{\partial u}{\partial{\bf t}}$ would be discontinuous at ${\bf x}_{0}$ ($u$ has a discontinuity in its first derivative). This situation might happen, as illustrated in Example 2. Using a similar argument we could prove the following theorem.

Theorem 2.3**.**

Suppose that $\gamma({\bf x},{\bf y})$ satisfies (1.3) and (1.4). If the following two conditions hold:

(i) $b({\bf x})\in C^{1}(\Omega)$, $b\in C^{2}(\Omega_{1})$, $b\in C^{2}(\Omega_{2})$,

*(ii) $\tilde{\gamma}(r)\in C^{1}(0,\delta)$ and $\tilde{\gamma}(\delta)=0$.

Then $u\in C^{1}(\Omega)$, $u\in C^{2}(\Omega_{1})$ and $u\in C^{2}(\Omega_{2})$.*

Proof.

The conditions of Theorem 2.2 hold due to the conditions (i) and (ii), so $u\in C^{1}(\Omega)$. Furthermore, for any two unit vectors ${\bf t}_{1}$ and ${\bf t}_{2}$, take a directional derivative for (2.7),

[TABLE]

Applying this equality in $\Omega_{1}$ and $\Omega_{2}$ will lead to the conclusion $u\in C^{2}(\Omega_{1})$ and $u\in C^{2}(\Omega_{2})$, respectively. ∎

To get the optimal convergence order, we always need the regularity $u\in H^{2}(\Omega)$, the following corollary give a sufficient condition to guarantee this property.

Corollary 2.1**.**

Suppose that $\gamma({\bf x},{\bf y})$ satisfies (1.3) and (1.4). If the following two conditions hold:

(i) $b({\bf x})\in C^{1}(\Omega)$, $b\in H^{2}(\Omega_{1})$, $b\in H^{2}(\Omega_{2})$,

*(ii) $\tilde{\gamma}(r)\in C^{1}(0,\delta)$ and $\tilde{\gamma}(\delta)=0$.

Then $u\in H^{2}(\Omega)$.*

Proof.

Using the density of $C^{2}(\Omega_{i})$ in $H^{2}(\Omega_{i})$ ($i=1,2$) and Theorem 2.3 we could get the result $u\in H^{2}(\Omega_{i})$. Since $u\in C^{1}(\Omega)$ is proven, the result $u\in H^{2}(\Omega)$ holds. ∎

3 A new DG method for nonlocal models with integrable kernels

Here and after, for a function $\phi({\bf x})$ we denote $\lim\limits_{h\rightarrow 0-}\phi({\bf x}+h{\bf n_{x}})$ by $\phi({\bf x}-)$ in the case of no ambiguity. Under the condition of Theorem 2.1, we know that for given ${\bf x}\in\partial\Omega$,

[TABLE]

However $u({\bf x}-)$ does not need to be zero, that is $u({\bf x})$ is possibly discontinuous across $\partial\Omega$. Thus, it might be inefficient to use continuous FEM on the whole domain $\Omega\cup\Omega_{I}$. Moreover, since we do not specify the value of the right hand side function on $\Omega_{I}$ a priori, we have no control on the amount of the jump across $\partial\Omega$ where the solution is likely to be discontinuous. Thus, we propose a suitable DG method by adopting a hybrid version of DG (across the domain boundary) and continuous elements (in the interior domain).

As in [16] the nonlocal energy inner product, the nonlocal energy norm, nonlocal energy space, and nonlocal volume constrained energy space are defined by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

respectively. Similar to the definition $V_{c,0}(\Omega\cup\Omega_{I})$, the subspace of $L^{2}(\Omega\cup\Omega_{I})$ is defined as follows:

[TABLE]

The authors in [16] prove that if the kernel function $\gamma({\bf x},{\bf y})$ satisfies (1.3) and (1.4), then $V_{c,0}(\Omega\cup\Omega_{I})$ is equivalent to $L^{2}_{c,0}(\Omega\cup\Omega_{I})$. Denote by

[TABLE]

where the superscripts $pc$ means partly continuous (continuous in $\Omega$).

For a given triangulation $\mathcal{T}_{h}$ of $\Omega\cup\Omega_{I}$ that simultaneously triangulates $\Omega$, let $\Omega_{h}=\mathcal{T}_{h}\cap\overline{\Omega}$. Next, let $V^{pc,h}_{c,0}$ consist of those functions in $V_{c,0}^{pc}(\Omega\cup\Omega_{I})$ that are piecewise linear. Since $\Omega$ is convex, this conforming property is satisfied, that is,

[TABLE]

We assume that $\mathcal{T}_{h}$ is shape-regular and quasi-uniform [8] as $h\rightarrow 0$, where $h$ denotes the diameter of the largest element in $\mathcal{T}_{h}$. For a fixed $\mathcal{T}_{h}$, the set of inner nodes of $\Omega_{h}$, i.e., all nodes in $\Omega_{h}\setminus\partial\Omega$, is denoted by $N\!I=\{{\bf x}_{j}:j=1,2,\cdots,m\}$, with piecewise linear basis functions defined on $\mathcal{T}_{h}$ being $\phi_{j}({\bf x}),\>j=1,2,\cdots,m$. The set of all nodes in $\Omega_{h}\cap\partial\Omega$ is denoted by $N\!B=\{{\bf x}_{m+j}:j=1,2,\cdots,n\}$ with piecewise linear basis functions defined on $\mathcal{T}_{h}$ being $\phi_{m+j}({\bf x}),\>j=1,2,\cdots,n$. The basis functions for the space $V_{c,0}^{pc,h}$ are as follows: for $j=1,2,\cdots,m+n$,

[TABLE]

Throughout the paper, the generic constant $C$ is always independent of the finite element mesh parameter $h$.

Since we know the structure of the true solution and the space it belongs to, we could design a DG method for its approximation. First, variational form in $V_{c,0}^{pc}(\Omega\cup\Omega_{I})$ finds $u({\bf x})\in V_{c,0}^{pc}(\Omega\cup\Omega_{I})$, such that

[TABLE]

The finite dimensional approximation for (3.2) finds $u_{h}({\bf x})\in V^{pc,h}_{c,0}$, such that

[TABLE]

Set $u_{h}({\bf x})=\sum\limits_{j=1}^{m+n}u_{j}\widehat{\phi}_{j}({\bf x})$, ${\bf u}=(u_{1},u_{2},\cdots,u_{m+n})^{T}$. Denote by

[TABLE]

and

[TABLE]

with

[TABLE]

Set $w_{h}=\widehat{\phi}_{i}$, $i=1,2,\cdots,m+n$, the algebraic system of (3.3) is

[TABLE]

with

[TABLE]

In the process to solve for $u_{h}({\bf x})$ we use the finite element space $V^{pc,h}_{c,0}$ which is continuous in $\Omega_{h}$, however, discontinuous across $\partial\Omega_{h}$, thus we regard it a conforming but hybrid version of DG and continuous FEM. This method possesses some advantages as follows:

(i) The method leads to a linear algebraic system with the coefficient matrix $A$ in (3.4) that is symmetric and positive definite, just as in the case using either the conforming DG or continuous FEM, thus many efficient solvers suitable to such linear systems could still be used.

(ii) The method is asymptotically compatible: as shown in [32], as long as the finite element space contains continuous piecewise linear functions (which is the case for our hybrid algorithm), the Galerkin finite element approximation is always asymptotically compatible, and thus offers robust numerical discretizations to problems involving nonlocal interactions.

(iii) The method has optimal convergence rate provided that the solution is smooth on $\Omega$, that is $O(h^{2})$ ($O(h)$) for error in $L^{2}$ ($H^{1}$) norm provided that the true solution $u\in H^{2}(\Omega)$. This result is in sharp contrast to the assumption given that in [16] where to insure the optimal convergence rate the true solution is required to be in $H^{2}(\Omega\cup\Omega_{I})$ which generally not holds for the problem (2.4). Furthermore, it has optimal convergence rate for 1D case under very weak assumptions, and nearly optimal convergence rate for two dimensional (2D) case under some mild assumptions, as shown in the next section.

(iv) The method, in comparison with the direct use of DG in all discrete elements, uses a smaller degree of freedoms. For example, the degree of freedoms is $n+1$ versus $2n$ for a mesh with $n+1$ nodes in 1D case, and $(n+1)^{2}$ versus $6n^{2}$ for a uniform triangulation with $n^{2}$ nodes in 2D case.

4 Theoretical Analysis

We now provide further theoretical analysis on the new DG approximations. Given what has already been discussed in (ii) of the above section, the asymptotic compatibility is assured and we thus focus on the case where the problems remain strictly nonlocal.

4.1 Convergence

The following convergence result describes the best approximation property of the finite-dimensional Ritz-Galerkin solution.

Theorem 4.1**.**

If $\gamma({\bf x},{\bf y})$ satisfies (1.3) and (1.4), $b({\bf x})\in C(\Omega)$, $u({\bf x})$ is the solution of (2.4), $u_{h}({\bf x})$ is the solution of (3.3). We define

[TABLE]

Then we have

[TABLE]

Consequently,

[TABLE]

Proof.

Since $V^{pc,h}_{c,0}\subset V^{pc}_{c,0}(\Omega\cup\Omega_{I})$ as in (3.1), then for all $w_{h}\in V_{c,0}^{pc,h}$,

[TABLE]

together with (3.3), we have

[TABLE]

Using the nonlocal Green’s first identity [17], we have

[TABLE]

Then we get the following estimate

[TABLE]

and then

[TABLE]

By the equivalence between $V_{c,0}(\Omega\cup\Omega_{I})$ and $L^{2}_{c,0}(\Omega\cup\Omega_{I})$, we complete the proof. ∎

Let us note that due to the use of norm equivalence in the above proof, generally speaking, the constant $C$ in the lemma could depend on the nonlocal space and thus the nonlocal kernel. One may not infer that this constant remains uniformly bounded when we consider the local limit of the nonlocal problem. Fortunately, as alluded to earlier, with the asymptotic compatibility already established in [32], we hereby only focus on the strict nonlocal case.

We now combine the theory of the interpolation error estimate and (4.2) to give the convergence rate estimate.

Theorem 4.2**.**

If $\gamma({\bf x},{\bf y})$ satisfies (1.3) and (1.4), $b({\bf x})\in C(\Omega)$, $u({\bf x})$ is the solution of (2.4), $u_{h}({\bf x})$ is the solution of (3.3). Suppose that $u\in H^{t}(\Omega)$ holds, there exists a constant $C$ such that, for sufficiently small $h$,

[TABLE]

with $s=\min(t,2)>d/2$. If $s>1$

[TABLE]

Moreover, if the following two conditions hold:

(i) $b({\bf x})\in C^{1}(\Omega)$, $b\in H^{2}(\Omega_{1})$, $b\in H^{2}(\Omega_{2})$;

*(ii) $\gamma({\bf x},{\bf y})$ is a radial function such that $\tilde{\gamma}(r)\in C^{1}(0,\delta)$ and $\tilde{\gamma}(\delta)=0$.

Then $u\in H^{2}(\Omega)$, thus $s=2$.*

Proof.

Denote by $\mathcal{I}_{h}u$ the Lagrange interpolant from $C(\Omega_{h})$ to $V_{c,0}^{pc,h}|_{\Omega_{h}}$

[TABLE]

then $w_{h}\in V^{pc,h}_{c,0}$, and

[TABLE]

with $s=\min(t,2)$. Combination of (4.2) and (4.5) leads to (4.3).

Using the inverse estimate for finite element space, we have

[TABLE]

This is the desired result (4.4).

The conditions (i) and (ii) lead to $u\in H^{2}(\Omega)$ due to Corollary 2.1. ∎

We recall by Theorem 6.2 in [16] that, when continuous FEM is used to approximate the nonlocal problem (2.4), the approximation $u_{n}$ has an error estimate of the form $\|u-u_{n}\|_{\Omega\cup\Omega_{I}}\leq Ch^{s}\|u\|_{s,\Omega\cup\Omega_{I}}$. Since the solution of nonlocal problem (2.4) could be discontinuous across $\partial\Omega$, we see that $u\in H^{s}(\Omega\cup\Omega_{I})$ does not hold for $s\geq 1/2$, let alone for $s=2$. For 1D case, the best to expect is $s=1/2-\epsilon$ for arbitrary small positive $\epsilon$ if continuous FEM is used. Theorem 4.2 improves the convergence rate from $1/2-\epsilon$ to $3/2-\epsilon$ for this case since we have the regularity of $H^{3/2-\epsilon}(\Omega)$. This convergence rate is still not optimal. If the points on $\partial\Omega_{1}$ are selected as the mesh grids, the optimal convergence rate could be obtained.

Theorem 4.3**.**

For the 1D case, assume $u\in C(0,1)$, $u\in H^{2}(0,\delta)$, $u\in H^{2}(\delta,1-\delta)$, $u\in H^{2}(1-\delta,1)$. If $\delta$ and $1-\delta$ are all selected as the mesh grids, then there exists a constant $C$ such that, for sufficiently small $h$,

[TABLE]

Proof.

Using the interpolation error estimate in three intervals $(0,\delta)$, $(\delta,1-\delta)$ and $(1-\delta,1)$, respectively and add them together, we get

[TABLE]

Together with (4.2) and the inverse estimate we get (4.6). ∎

For 2D case, under mild assumptions on the smoothness of the boundary, the regularity of the solution, and the conformity between the mesh and the boundary, we have the almost optimal convergence rate, that is optimal up to a factor $|\log h|^{1/2}$.

Theorem 4.4**.**

For 2D case, assume $u\in C(\Omega)$, $u\in H^{2}(\Omega_{1})$, $u\in H^{2}(\Omega_{2})$. If $\partial\Omega$ is of class $C^{2}$, and $N\!B\subset\partial\Omega_{1}$, then there exists a constant $C$ such that, for sufficiently small $h$,

[TABLE]

Proof.

We cite from [11] the results for the linear interpolation error estimate for the interface problem, that is

[TABLE]

Together with (4.2) and the inverse estimate we complete the proof. ∎

Remark 4.1**.**

Due to the structure of the solution for the problem (2.4) we have discussed, the solution $u$ is often discontinuous across $\partial\Omega$. It may cause the discontinuity of the first derivative across $\partial\Omega_{1}$ if $\tilde{\gamma}(r)$ has a discontinuity at $r=\delta$ (that is $\tilde{\gamma}(\delta)>0$), or the discontinuity of the second derivative across $\partial\Omega_{1}$ if $\tilde{\gamma}(r)$ is continuous at $r=\delta$ (that is $\tilde{\gamma}(\delta)=0$) but $\tilde{\gamma}^{\prime}(r)$ has a discontinuity at $r=\delta$ (that is $\tilde{\gamma}^{\prime}(\delta-)\neq 0$). In the next section, we will discuss two kinds of kernels representing the above two cases, respectively.

4.2 Condition number estimate

The condition number of the stiffness matrix is an indicator of the sensitivity of the discrete solution with respect to the data and the performance of iterative solvers such as the conjugate-gradient method. For the DG method we propose in Section 3, consider the $(m+n)\times(m+n)$ stiffness matrix $A$ defined in (3.4). We have the following condition number estimate whose proof is standard and is given for completeness. Similar discussions can be found in earlier studies[1, 36].

Theorem 4.5**.**

For the stiffness matrix $A$ defined in (3.4) associated with the kernel $\gamma({\bf x},{\bf y})$ that satisfies (1.3) and (1.4), there exists a constant $C$ such that

[TABLE]

Proof.

For the given finite element nodal basis, there exist two generic constants $c_{2}\geq c_{1}>0$ such that

[TABLE]

where $\{w_{j}\}$, $j=1,2,\cdots,m+n$, are components of the vector ${\bf w}$. Since the space $V_{c,0}(\Omega\cup\Omega_{I})$ is equivalent to the space $L^{2}_{c,0}(\Omega\cup\Omega_{I})$, we get the theorem immediately. ∎

Remark 4.2**.**

We note again that the constant $C$ may depend on the kernel, as demonstrated in [36], hence the result is only meaningful for nonlocal problems with a fixed kernel that satisfies the assumptions (1.3) and (1.4).

5 Numerical results

We now report results of numerical experiments which substantiate the theoretical analysis in Section 4. For 1D case, the problem (2.4) becomes the following form

[TABLE]

We solve the nonlocal problem first on uniform meshes and take $\delta$ to be a constant multiple of $h$ and reduce $h$ to check convergence and condition number properties of the proposed DG method, and then solve the problem on non-uniform meshes which are obtained by random disturbance to uniform meshes. Here we choose two popular examples of kernel functions representing two cases as discussed in Remark 4.1.

5.1 Constant kernel function

We first consider the following kernel function

[TABLE]

Obviously $\gamma$ defined as (5.2) is discontinuous at points $\pm\delta$, if $b$ is discontinuous at $x=0$ or $x=1$, the solution of (5.1) will probably be (however, not necessarily) discontinuous in its first derivative at $x=\delta$ or $x=1-\delta$. In fact, in Example 1 the smoothness pick-up is beyond first order, that is, although $b$ is discontinuous at $x=1$, $u$ is infinitely continuously differentiable at $x=1-\delta$. While in Example 2 the smoothness pick-up is only first order and could not be improved, that is, since $b$ is discontinuous at $x=0$ and $x=1$, the first derivative of $u$ is is discontinuous at $x=\delta$ and $x=1-\delta$.

Example 1. In order to get simpler benchmark solutions, we calculate the right-hand side of (5.1) based on an exact solution $u(x)=x^{2},\>x\in\Omega=(0,1)$ and $u(x)=0,\>x\in\Omega_{I}=(-\delta,0)\cup(1,1+\delta)$, with kernel function (5.2). This naturally leads to a $\delta$-dependent right-hand side $b(x)=b_{\delta}(x)$, see Figure 1 for the plots of $u(x)$ and $b(x)$. The DG method we proposed in Section 3 is used to discretize it with $\delta=0.4$.

We first use the proposed DG method on uniform meshes and conclude from Table 1 that convergence rates for errors in $L^{2}$ and $H^{1}$ norms are all optimal. The spectral condition number of the stiffness matrix is almost constant when the mesh size decreases, indicating the insensitivity of the discrete solution regardless how small $h$ is.

We then use a kind of non-uniform meshes obtained by random disturbance to uniform meshes. To be specific, for a fixed $m$, let $h=\delta/m$, the non-uniform mesh is obtained by adding a random vector ${\bf\varepsilon}\in\mathbb{R}^{m-1}$ (which obeys the uniform distribution on $[-0.1h,0.1h]$) to $x_{i}$ to reach $x_{i}+\varepsilon_{i},\>i=1,2,\cdots,m-1$. Together with $x_{0}$ and $x_{m}$ we get the new mesh grids

[TABLE]

We have done over twenty tests, and the convergence rates and the spectral condition numbers are all similar. Thus, instead of listing all of them, we just select one test to verify our theoretical analysis. Similar actions and presentations are made also in later examples. It is seen from Table 2 that the errors in $L^{2}$ and $H^{1}$ norms and convergence rates are comparable with that in uniform meshes case. This is consistent with the theoretical result in Theorem 4.2 since $u\in C^{\infty}(\Omega)$, thus $s=2$. The spectral condition numbers of the stiffness matrices behave similar as in the uniform meshes case, too.

Example 2. We consider (5.1) with kernel function (5.2) and $b(x)=e^{x}$. The DG method we proposed in Section 3 is used to discretize it with $\delta=0.4$.

Since the exact solution $u(x)$ is not known for this problem we compute errors using the solution on finer meshes as approximation of the true solution. We first use the proposed DG method on uniform meshes. The right hand side function $b(x)$ and the approximation $u_{h}$ with $h=0.00625$ are plotted in Figure 2. Although $b(x)$ is in $C^{\infty}(0,1)$, it has two discontinuous points $x=0$ and $x=1$, which causes the discontinuity in the first derivative of $u$ at $x=\delta=0.4$ and $x=1-\delta=0.6$. From Table 3 it is seen that the method has optimal convergence rates for errors in $L^{2}$ and $H^{1}$ norms. The spectral condition numbers of the stiffness matrices for the method behave similarly as Example 1. Since the results for the spectral condition numbers of the stiffness matrices are all similar in the rest of the numerical tests, we no longer list them to avoid repetitions.

Next, we use the non-uniform meshes (5.3) to solve the problem. In this example the true solution $u(x)\in H^{1.5-\epsilon}(\Omega)$ for arbitrary small positive $\epsilon$. So for general non-uniform meshes, we expect the convergence rates for errors in $L^{2}$ ($H^{1}$) norm to be at most $1.5$ ($0.5$). This fact is verified in Table 4. However, since we know the two discontinuous points of $u^{\prime}(x)$, if they are selected as mesh grids, optimal convergence rates could be recovered. To be specific, random disturbances are added to the original mesh grids except $\delta$ and $1-\delta$. The corresponding results are shown in Table 5. That is the convergence rates for errors in $L^{2}$ and $H^{1}$ norms are $2$ and $1$, respectively, which is consistent with the theoretical result (4.6). We then re-examine the computation on the meshes used in Table 4. The convergence rates, however, are not necessarily similar, which is different with Example 1. To be specific, if the perturbation for the points $\delta$ and $1-\delta$ is large in absolute value, such as $0.1h$, the convergence rates remain similar to Table 4. On the other hand, if the perturbation for the two points is small or close to zero, the convergence rates of the test are similar to Table 5.

5.2 Non-constant kernel function

We then consider the kernel function

[TABLE]

Obviously the first derivative of $\gamma$ is discontinuous at points $\pm\delta$, if $b$ is discontinuous at $x=0$ or $x=1$, thus $u$ will likely be discontinuous in its second derivative at $x=\delta$ or $x=1-\delta$. This is the case in Example 3 where the regularity pick-up is only second order and could not be further improved. That is, since $b$ is discontinuous at $x=0$ and $x=1$, the second derivative of $u$ is is discontinuous at $x=\delta$ and $x=1-\delta$.

Example 3. We consider (5.1) with kernel function (5.4) and $b(x)=0.01e^{6x}$. The proposed DG method in Section 3 is used to discretize it with $\delta=0.4$.

As in Example 2, since we do not know the exact solution $u(x)$, errors are computed using the solution on finer meshes as approximation of the true solution. We first implement the proposed DG method on uniform meshes. The approximation $u_{h}(x)$ with $h=0.00625$ are plotted in Figure 3. To see the discontinuity of the second derivative, the first derivative of $u_{h}(x)$ is also plotted in Figure 3. Although $b(x)\in C^{\infty}(0,1)$, however, $b(x)$ has two discontinuous points $x=0$ and $x=1$, which causes the discontinuity for the second derivative of $u$ at $x=\delta=0.4$ and $x=1-\delta=0.6$. From Table 6 it is seen that optimal convergence rates for errors in $L^{2}$ and $H^{1}$ norms are achieved.

Next, we use non-uniform meshes (5.3) to solve the problem. In this example the true solution satisfies $u(x)\in H^{2.5-\epsilon}(\Omega)$ for arbitrary small positive $\epsilon$. So for general non-uniform meshes the theoretical convergence rates for errors in $L^{2}$ and $H^{1}$ norms are all optimal. This is indeed verified in Table 7.

Remark 5.1**.**

We have considered the homogeneous problem (2.4) with a right hand side function $b({\bf x})\in C(\Omega)$. The inhomogeneous problem (1.2) could be converted to a homogeneous one like (2.4) via the transforms (2.1) and (2.2).

6 Conclusion

We propose a new kind of DG method in this paper to numerically solve the nonlocal models with integrable kernels. The existed references tell us that if the right hand side function and the volume constraint function are in $L^{2}(\Omega)$ and $L^{2}(\Omega_{I})$, respectively, then the true solution of that nonlocal model also belongs to $L^{2}(\Omega\cup\Omega_{I})$. Such a general result makes the numerical approximation difficult to operate, or easy to operate but not so efficiently. To make the approximation easier and more efficient simultaneously, we first convert the original nonhomogeneous problem with right hand side function in $L^{2}(\Omega)$ to a homogeneous problem with right hand side function continuous in $\Omega$. Then we analyze the structure of true solution of the homogeneous problem, especially for higher dimensional cases. The main result is, this kind of problem often encounters the discontinuity across the boundary $\partial\Omega$, thus possibly causes the discontinuity of first or second derivative (perhaps higher order derivatives, depending on the smoothness of the kernels) across $\partial\Omega_{1}$. Based on this observation, an appropriate DG method is proposed which has some good properties, such as, the matrix of the algebraic system is symmetrical positive definite and has almost constant spectral condition number independent of the mesh size, the method is asymptotically compatible and uses less degrees of freedom compared with direct use of DG method. Moreover, it has optimal convergence rate for 1D case under very weak assumption, and the almost optimal convergence rate for 2D case under mild assumption. This is the essential improvement over the existed theory for standard approximation like continuous FEM.

The error and condition number estimates for the method are proven in Section 4 for any dimensional case. However, the numerical experiments are implemented just for 1D case in Section 5. This is because the implementation of a FEM for a nonlocal problem involves calculation of double integral which is rather complicated for higher dimensional cases. To be more specific, since the kernels we discussed are supported on a ball of radius $\delta$, when the inner integral of the double integral is written by a sum of some components (also integrals) over the elements which overlap with the support of the kernel, the integral region for some components is strictly contained within the support of the kernel. For such special integrals, some specifically designed quadrature rule should be used to obtain good accuracy. Authors in [34, 35] discuss this issue in details and use the quadrature rule for multiscale implementation for PD models. Interested readers are referred to them and references cited therein.