# FEM on nonuniform meshes for nonlocal Laplacian: Semi-analytic
Implementation in One Dimension

Hongbin Chen1, Changtao Sheng2 and Li-Lian Wang3

###### Abstract.

In this paper, we compute stiffness matrix of the nonlocal Laplacian
discretized by the piecewise linear finite element on nonuniform meshes, and
implement the FEM in the Fourier transformed domain. We derive useful integral
expressions of the entries that allow us to explicitly or semi-analytically
evaluate the entries for various interaction kernels. Moreover, the limiting
cases of the nonlocal stiffness matrix when the interactional radius
$\delta\rightarrow 0$ or $\delta\rightarrow\infty$ automatically lead to
integer and fractional FEM stiffness matrices, respectively, and the FEM
discretisation is intrinsically compatible. We conduct ample numerical
experiments to study and predict some of its properties and test on different
types of nonlocal problems. To the best of our knowledge, such a semi-analytic
approach has not been explored in literature even in the one-dimensional case.

###### Key words and phrases:

Finite element method, nonlocal operator, nonuniform meshes, nonlocal
stiffness matrix, nonlocal Allen-Cahn equation

1College of Computer Science and Mathematics, Central South University of
Forestry and Technology, Changsha, Hunan, 410004, China. The research of the
author is partially supported by the Natural Science Foundation of Hunan
Province (No: 2022JJ30996). Email<EMAIL_ADDRESS>(H. Chen).
2Corresponding author. School of Mathematics, Shanghai University of Finance
and Economics, Shanghai, 200433, China. The research of the author is
partially supported by the National Natural Science Foundation of China (Nos.
12201385 and 12271365) and the Fundamental Research Funds for the Central
Universities (No: 2021110474). Email<EMAIL_ADDRESS>(C. Sheng).
3Division of Mathematical Sciences, School of Physical and Mathematical
Sciences, Nanyang Technological University (NTU), 637371, Singapore. The
research of the author is partially supported by Singapore MOE AcRF Tier 2
Grants: MOE2018-T2-1-059 and MOE2017-T2-2-144. Email<EMAIL_ADDRESS>(L.-L.
Wang).

##  1\. Introduction

We consider the following nonlocal constrained value problem on a nonzero
measure volume:

$\begin{cases}\mathcal{L}_{\delta}u(x)=f(x),\quad&\text{on}\
\Omega=(a,b),\\\\[4.0pt] u(x)=0,\quad&\text{on}\
\Omega_{\mathcal{\delta}}=[a-\delta,a]\cup[b,b+\delta],\end{cases}$ (1.1)

where the nonlocal operator is defined as

$\mathcal{L}_{\delta}u(x)=\frac{1}{2}\int_{-\delta}^{\delta}(u(x+s)-u(x))\rho_{\delta}(|s|)\rm{d}s,$
(1.2)

where $\delta$ is called a horizon parameter and the kernel
$\rho_{\delta}=\rho_{\delta}(s)$ is a nonnegative, radial function (i.e.,
$\rho_{\delta}(s)=\rho_{\delta}(|s|)$ for any $s$), compactly supported in
$|s|\leq\delta$, and has a bounded second moment. It is known that nonlocal
models (1.1) provides an alternative way to modeling many phenomena and
complex systems in diverse fields, for instance, the materials science with
crack nucleation and growth, fracture, and failure of composites [5], the
peridynamics model for mechanics [9, 28, 36], and nonlocal heat conduction
[12]. We also refer to the up-to-date review paper [12, 15, 18] and the
monograph [14] for more detail.

There has been much interest in developing numerical algorithms for nonlocal
models, particularly for PDEs involving the integral fractional Laplacian,
which can be viewed as a limiting case of the nonlocal operator [17, 38],
i.e., $\delta\rightarrow\infty.$ The recent works on PDEs involving integral
fractional Laplacian operator particularly include finite difference methods
[16, 22, 23], finite element methods [2, 3], and spectral methods [31, 33,
35]. The interested readers may also refer to [7, 18, 27] for three up-to-date
review papers. It should be noted that finite element method in frequency
domain space for integral fractional Laplacian is provided in [10, 29, 34],
where the entries of the stiffness matrix can provide exact form in one
dimension, and can be simply expressed as one-dimensional integral on a finite
interval in two dimension. This fact encourages us to try to extend this
method for more general problems involving the nonlocal operator (1.2).

Besides, extensive numerical results have been conducted on various kinds
nonlocal models. Chen and Gunzburger [9] provided a systematic study of
algorithmic and implementation issues using continuous and discontinuous
Galerkin finite element methods for peridynamics model. Tian and Du [37]
introduce a finite element and quadrature-based finite difference scheme,
compare the similarities and differences of the nonlocal stiffness matrices
for the two methods, and study fundamental numerical analysis issues. Wang and
Tian [32] develop a fast and faithful collocation method for a two-dimensional
nonlocal diffusion model. Liu, Cheng and Wang [28] introduced and analyzed a
fast Galerkin and adaptive Galerkin finite element methods to solve a steady-
state peridynamic model. Additionally, other numerical methods for nonlocal
models have been proposed, including the finite difference method [40], the
finite element methods [24, 6], the collocation methods [25, 21, 30, 8], and
so on. The development of the aforementioned methods relies on constructing
specific quadrature formulas to discretize the underlying nonlocal operators,
which unavoidably introduces truncation errors. However, these truncation
errors often hinder our ability to analyze the attractive structure of the
nonlocal stiffness matrices and affect the numerical stability of our solver
on some extreme meshes, e.g., graded and geometric meshes. Recently, Li, Liu,
and Wang [26] developed an efficient Hermite spectral-Galerkin method for
nonlocal diffusion equations in unbounded domains, where the derivation
therein was built upon the convolution property of the Hermite functions.

In this paper, we provide a semi-analytic means for computing the piecewise
linear FEM stiffness matrix for nonlocal problem (1.1). Then, we can obtain an
explicit form of the one-dimensional integral for the entries of the stiffness
matrix without generating any truncation errors. This one-dimensional integral
can even be evaluated exactly for most of the kernel functions, such as the
fractional kernel $\rho_{\delta}(|s|)=C_{\delta,\alpha}/s^{1+\alpha}$, with
the normalized constant $C_{\delta,\alpha}$ and the index $\alpha\in[-1,2)$.
Note that such an integral representation is derived from (i) implementation
of FEM in the Fourier transformed domain, and (ii) use of some analytic
formulas in the handbook [1, 11] to evaluation of the inner integral. Benefit
from this explicit expression, we show the two limiting cases of the nonlocal
stiffness matrix: (a) when $\delta\rightarrow 0$, it reduces to the stiffness
matrix associated with the classical Laplacian operator $-\Delta$, see (2.33),
and (b) when $\delta\rightarrow\infty$, it becomes the stiffness matrix
corresponding to the integral fractional Laplacian operator
$(-\Delta)^{\frac{\alpha}{2}}$. As a byproduct, we numerically predicted the
condition number of the nonlocal stiffness matrix $\bm{{S_{\delta}}}$ (even
for some extreme mesh) behaves like
$\text{Cond}(\bm{S_{\delta}})\leq\min\\{\left(h_{\max}/h_{\min}\right)^{\alpha-1}N^{\alpha}\delta^{\alpha-2},\,\,\left(h_{\max}/h_{\min}\right)^{2}N^{2}\\}$,
where $h_{\max}$ and $h_{\min}$ denote the maximum and minimum of the mesh
size, and $N$ denotes the number of meshes. Our approach is suitable for wider
range of nonlocal interaction kernel on the uniform and nonuniform meshes.
Numerical experiments show that our scheme is asymptotically compatible scheme
[37]. Note that we only focus on one dimensional problems in this work to
better illustrate the idea. The multidimensional generalizations are possible
and will be carried out in separate works.

The rest of this paper is organized as follows. In section 2, we present a
piecewise linear finite element method on nonuniform and uniform meshes for
model problem (1.1), and derive the explicit integral form for the entries of
the nonlocal stiffness matrix. In section 3, ample numerical results for
nonlocal problems with various both regular and singular kernels, including
accuracy testing, limiting behaviors, eigenvalue problem and nonlocal Allen-
Cahn equation, are presented to show the efficiency and accuracy. Finally, we
give some concluding remarks and discussions in the last section.

##  2\. Finite Element Method in Fourier Domain

Following the same notation as in [36], we introduce the energy space

$\mathcal{V}(\Omega)=\Big{\\{}u\in
L^{2}(\Omega\cup\Omega_{\mathcal{\delta}}):\int_{\Omega}\int_{-\delta}^{\delta}(u(x+s)-u(x))^{2}\rho_{\delta}(|s|){\rm
d}s{\rm d}x<\infty\Big{\\}},$ (2.1)

equipped with the semi-norm

$|u|_{\mathcal{V}}=\Big{\\{}\int_{\Omega}\int_{-\delta}^{\delta}(u(x+s)-u(x))^{2}\rho_{\delta}(|s|){\rm
d}s{\rm d}x\Big{\\}}^{1/2},$ (2.2)

and the norm

$\|u\|_{\mathcal{V}}=(|u|_{\mathcal{V}}^{2}+\|u\|^{2})^{1/2},\quad\|u\|^{2}=\int_{\Omega}|u(x)|^{2}{\rm
d}x.$ (2.3)

We further define the constrained energy space as

$\mathcal{V}_{0}(\Omega)=\big{\\{}u\in\mathcal{V}(\Omega)\ \big{|}\,u=0\ {\rm
in}\,\,\Omega_{\delta}\big{\\}}.$ (2.4)

The bilinear form on $L^{2}(\Omega)\times L^{2}(\Omega)$ and the inner product
on $L^{2}(\Omega)$ is, respectively, denoted by

$\begin{split}&\mathcal{A}_{\delta}(u,v)=\frac{1}{2}\int_{\Omega}\int_{-\delta}^{\delta}(u(x+s)-u(x))(v(x+s)-v(x))\rho_{\delta}(|s|){\rm
d}s{\rm d}x,\\\ &(u,v)_{\Omega}=\int_{\Omega}u(x)v(x)\,{\rm d}x.\end{split}$
(2.5)

Then, the weak formulation of (1.1) is to find $u\in\mathcal{V}_{0}(\Omega)$
such that

$\mathcal{A}_{\delta}(u,v)=(f,v)_{\Omega},\quad\forall
v\in\mathcal{V}_{0}(\Omega).$ (2.6)

In this section, we compute the nonlocal FEM stiffness matrix. To fix the
idea, we focus on the piecewise linear FEM with the partition:

$\Omega_{h}:=\\{x_{j}:a=x_{0}<x_{1}<\cdots<x_{N+1}=b\\},$ (2.7)

and the piecewise linear finite element basis:

$\phi_{j}(x)=\begin{dcases}\frac{x-x_{j-1}}{h_{j}},\quad&x\in[x_{j-1},x_{j}],\\\
\frac{x_{j+1}-x}{h_{j+1}},\quad&x\in(x_{j},x_{j+1}],\\\
0,\quad&\text{elsewhere
on}\;\;\Omega\cup\overline{\Omega}_{\delta},\end{dcases}\;h_{j}=x_{j}-x_{j-1},$
(2.8)

for $1\leq j\leq N.$ Correspondingly, we have the FE approximation space:

$\mathcal{V}_{h}^{0}:={\rm span}\big{\\{}\phi_{j}\,:\,1\leq j\leq N\big{\\}}.$

Then the FE approximation to (1.1) is to find $u_{h}\in\mathcal{V}_{h}^{0}$
such that

$\begin{split}\mathcal{A}_{\delta}(u_{h},v_{h})&=\frac{1}{2}\int_{\Omega}\int_{-\delta}^{\delta}(u_{h}(x+s)-u_{h}(x))(v_{h}(x+s)-v_{h}(x))\rho_{\delta}(|s|)\,\rm{d}s\rm{d}x\\\
&=(I_{h}f,v_{h})_{\Omega},\quad\forall\,\upsilon_{h}\in\mathcal{V}_{h}^{0},\end{split}$
(2.9)

where $I_{h}f$ is the standard piecewise linear interpolation of $f$ on
$\Omega_{h}$. We write

$u_{h}(x)=\sum\limits_{j=1}^{N}u_{j}\phi_{j}(x)\in\mathcal{V}_{h}^{0},$

and obtain the following linear system of (2.9):

$\bm{S}_{\delta}\bm{u}=\bm{f},$

where $\bm{S}_{\delta}$ is the nonlocal FEM stiffness matrix with entries
given by

$(\bm{S}_{\delta})_{jk}=(\bm{S}_{\delta})_{kj}=\frac{1}{2}\int_{\Omega}\int_{-\delta}^{\delta}(\phi_{j}(x+s)-\phi_{j}(x))(\phi_{k}(x+s)-\phi_{k}(x))\rho_{\delta}(|s|)\,{\rm
d}s{\rm d}x,$ (2.10)

and

$\bm{u}=(u_{1},u_{2},\cdots,u_{N})^{\top},\;\;\;\bm{f}=(f_{1},f_{2},\cdots,f_{N})^{\top},\;\;\;f_{k}=(I_{h}f,\phi_{k})_{\Omega}.$

###  2.1. FEM nonlocal stiffness matrix

Through the implementation in the Fourier transformed domain, we can derive
the following exact integral formulas for computing the stiffness matrix.

######  Theorem 2.1.

Let

$d_{j}^{k}=|x_{j}-x_{k}|,\quad\bm{c}_{j}:=\Big{(}\frac{1}{h_{j}},-\frac{1}{h_{j}}-\frac{1}{h_{j+1}},\frac{1}{h_{j+1}}\Big{)},$
(2.11)

and define

$g(z):=g(z;s):=-\frac{1}{12}(|z+s|^{3}-2|z|^{3}+|z-s|^{3}).$ (2.12)

Then the stiffness matrix can be evaluated by the component-wise integral

$\begin{split}\bm{S}_{\delta}=&\int_{0}^{\delta}\bm{I}(s)\,\rho_{\delta}(s)\,{\rm
d}s,\end{split}$ (2.13)

where the matrix-valued function $\bm{I}(s)\in{\mathbb{R}}^{N\times N}$ with
the entries given by

$\begin{split}\bm{I}_{jk}(s)=\bm{I}_{kj}(s)=\bm{c}_{j}g(\bm{D}_{j}^{k})\bm{c}_{k}^{\top},\quad\bm{D}_{j}^{k}:=\begin{pmatrix}d_{j-1}^{k-1}&d_{j-1}^{k}&d_{j-1}^{k+1}\\\\[2.0pt]
d_{j}^{k-1}&d_{j}^{k}&d_{j}^{k+1}\\\\[2.0pt]
d_{j+1}^{k-1}&d_{j+1}^{k}&d_{j+1}^{k+1}\end{pmatrix}.\end{split}$ (2.14)

Here, $g(\bm{D}_{j}^{k})$ means the function $g$ performs component-wisely
upon the matrix $\bm{D}_{j}^{k}$.

###### Proof.

Let $\tilde{\phi}_{j}(x)$ be the zero extension of $\phi_{j}(x)$ on
$\Omega\cup\overline{\Omega}_{\delta}$ to $\mathbb{R}$. Then by [10, Lemma
2.1], its Fourier transform is

$\displaystyle\mathscr{F}[\tilde{\phi}_{j}(x)](\xi)$
$\displaystyle=-\frac{1}{\sqrt{2\pi}}\frac{\beta_{j}e^{-{\rm
i}x_{j-1}\xi}-(\beta_{j}+\beta_{j+1})e^{-{\rm i}x_{j}\xi}+\beta_{j+1}e^{-{\rm
i}x_{j+1}\xi}}{\xi^{2}}$ (2.15)
$\displaystyle=-\frac{1}{\sqrt{2\pi}\xi^{2}}\bm{c}_{j}\cdot\bm{e}_{j}(\xi),$

where

$\beta_{j}=\frac{1}{h_{j}},\quad\bm{e}_{j}(\xi):=\big{(}e^{-{\rm
i}x_{j-1}\xi},e^{-{\rm i}x_{j}\xi},e^{-{\rm i}x_{j+1}\xi}\big{)}.$

Using Parseval identity and the translation property:
${\mathscr{F}}[u(x+s)](\xi)=e^{{\rm i}\xi s}{\mathscr{F}}[u(x)](\xi)$, we
derive from (2.8), (2.10) and (2.15) that

$\begin{split}\mathcal{A}_{\delta}(\phi_{j},\phi_{k})&=\frac{1}{2}\int_{\Omega}\int_{-\delta}^{\delta}(\phi_{j}(x+s)-\phi_{j}(x))(\phi_{k}(x+s)-\phi_{k}(x))\rho_{\delta}(|s|){\rm
d}s{\rm d}x\\\
&=\frac{1}{2}\int_{{\mathbb{R}}}\int_{-\delta}^{\delta}(\tilde{\phi}_{j}(x+s)-\tilde{\phi}_{j}(x))(\tilde{\phi}_{k}(x+s)-\tilde{\phi}_{k}(x))\rho_{\delta}(|s|){\rm
d}s{\rm d}x\\\
&=\frac{1}{2}\int_{-\delta}^{\delta}\rho_{\delta}(|s|)\int_{{\mathbb{R}}}{\mathscr{F}}[\tilde{\phi}_{j}(x+s)-\tilde{\phi}_{j}(x)](\xi)\overline{{\mathscr{F}}[\tilde{\phi}_{k}(x+s)-\tilde{\phi}_{k}(x)](\xi)}{\rm
d}\xi{\rm d}s\\\
&=\frac{1}{2}\int_{-\delta}^{\delta}\rho_{\delta}(|s|)\int_{{\mathbb{R}}}|e^{{\rm
i}\xi
s}-1|^{2}{\mathscr{F}}[\tilde{\phi}_{j}(x)](\xi)\overline{{\mathscr{F}}[\tilde{\phi}_{k}(x)](\xi)}{\rm
d}\xi{\rm d}s\\\
&=\int_{-\delta}^{\delta}\rho_{\delta}(|s|)\int_{{\mathbb{R}}}(1-\cos(\xi
s)){\mathscr{F}}[\tilde{\phi}_{j}(x)](\xi)\overline{{\mathscr{F}}[\tilde{\phi}_{k}(x)](\xi)}{\rm
d}\xi{\rm d}s\\\
&=2\int_{0}^{\delta}\rho_{\delta}(s)\int_{{\mathbb{R}}}(1-\cos(\xi
s)){\mathscr{F}}[\tilde{\phi}_{j}(x)](\xi)\overline{{\mathscr{F}}[\tilde{\phi}_{k}(x)](\xi)}{\rm
d}\xi{\rm d}s.\end{split}$ (2.16)

Using (2.15), we obtain from direct calculation that

$\begin{split}&\int_{{\mathbb{R}}}(1-\cos(\xi
s)){\mathscr{F}}[\tilde{\phi}_{j}(x)](\xi)\overline{{\mathscr{F}}[\tilde{\phi}_{k}(x)](\xi)}{\rm
d}\xi\\\
\quad&=\frac{1}{2\pi}\bm{c}_{j}\Big{(}\int_{\mathbb{R}}{\xi}^{-4}(1-\cos(\xi
s))\bm{e}_{j}(\xi)\cdot\bm{e}_{k}(-\xi){\rm
d}\xi\Big{)}\bm{c}_{k}^{\top}=\frac{1}{\pi}\int_{0}^{\infty}{\xi}^{-4}f_{jk}(\xi)\,{\rm
d}\xi,\end{split}$ (2.17)

where

$f_{jk}(\xi)=\bm{c}_{j}\bm{F}_{jk}(\xi)\bm{c}_{k}^{\top},\quad\bm{F}_{jk}(\xi)=(1-\cos(\xi
s))\cos(\bm{D}_{j}^{k}\xi).$ (2.18)

Using the trigonometric identities

$\begin{split}(1-\cos(\xi
s))\cos(d_{j}^{k}\xi)&=\cos(d_{j}^{k}\xi)-\frac{1}{2}\big{(}\cos((d_{j}^{k}+s)\xi)+\cos((d_{j}^{k}-s)\xi)\big{)}\\\
&=\frac{1}{2}\big{(}2\cos(d_{j}^{k}\xi)-\cos(|d_{j}^{k}+s|\xi)-\cos(|d_{j}^{k}-s|\xi)\big{)},\end{split}$
(2.19)

we obtain from (2.14) and (2.18) readily that

$\begin{split}f^{\prime\prime\prime}_{jk}(\xi)=\bm{c}_{j}\bm{F}^{\prime\prime\prime}_{jk}(\xi)\bm{c}_{k}^{\top}\quad{\rm
and}\quad\bm{F}^{\prime\prime\prime}_{jk}(\xi)=\frac{1}{2}g_{1}(\bm{D}_{j}^{k};\xi),\end{split}$
(2.20)

where we denote

$g_{1}(z;\xi):=2|z|^{3}\sin(|z|\xi)-|z+s|^{3}\sin(|z+s|\xi)-|z-s|^{3}\sin(|z-s|\xi).$

One verifies readily that $\partial_{\xi}^{k}g_{1}(z;0)=f_{jk}^{(k)}(0)=0$ for
$k=0,1,2,3.$ Thus, we derive from integration by parts that

$\begin{split}\int_{0}^{\infty}{\xi}^{-4}f_{jk}(\xi)\,{\rm
d}\xi&=\frac{1}{3}\int_{0}^{\infty}{\xi}^{-3}f_{jk}^{\prime}(\xi)\,{\rm
d}\xi=\frac{1}{6}\int_{0}^{\infty}{\xi}^{-2}f_{jk}^{\prime\prime}(\xi)\,{\rm
d}\xi\\\
&=\frac{1}{6}\int_{0}^{\infty}{\xi}^{-1}f_{jk}^{\prime\prime\prime}(\xi)\,{\rm
d}\xi.\end{split}$ (2.21)

Recall the sine integral formula (see [11, P. 427]):

$\int_{0}^{\infty}\frac{\sin(ax)}{x}\,{\rm
d}x=\frac{\pi}{2}\operatorname{sign}(a)=\frac{\pi}{2}\begin{cases}-1,\;\;&{\rm
if}\;\;a<0,\\\ 0,\;\;&{\rm if}\;\;a=0,\\\ 1,\;\;&{\rm if}\;\;a>0.\end{cases}$
(2.22)

Inserting (2.20) into (2.21), we find from (2.22) immediately that

$\int_{0}^{\infty}{\xi}^{-4}f_{jk}(\xi)\,{\rm
d}\xi=\frac{1}{6}\int_{0}^{\infty}{\xi}^{-1}f_{jk}^{\prime\prime\prime}(\xi)\,{\rm
d}\xi=\frac{\pi}{2}\bm{c}_{j}g(\bm{D}_{j}^{k})\bm{c}_{k}^{\top}.$ (2.23)

Then we obtain from (2.16)-(2.17) and (2.23) that

$(\bm{S}_{\delta})_{jk}=\frac{2}{\pi}\int_{0}^{\delta}\rho_{\delta}(s)\Big{(}\int_{0}^{\infty}{\xi}^{-4}f_{jk}(\xi)\,{\rm
d}\xi\Big{)}{\rm
d}s=\bm{c}_{j}\Big{(}\int_{0}^{\delta}g(\bm{D}_{j}^{k})\,\rho_{\delta}(s)\,{\rm
d}s\Big{)}\bm{c}_{k}^{\top}.$ (2.24)

This completes the proof. ∎

Theorem 2.1 is valid for the interaction radius $\delta<\infty.$ The bandwidth
of the stiffness matrix increases as $\delta/h$ increases, which becomes full
as $\delta\geq b-a.$ It is known that the nonlocal operator (1.2) tend to
integral fractional Laplacian as $\delta\to\infty$ when the corresponding
kernel function is chosen to satisfy

$\rho_{\delta}(s)=\rho_{\infty}(s)\mathbf{\chi}_{(0,\delta)},$ (2.25)

where $\mathbf{\chi}_{(0,\delta)}$ denotes the characteristic (indicator)
function of $(0,\delta)$ and

$\rho_{\infty}(s)=\frac{C_{\alpha}}{s^{1+\alpha}}\,\;\text{with}\,\;C_{\alpha}=\frac{2^{\alpha-1}\alpha\Gamma(\frac{1+\alpha}{2})}{\sqrt{\pi}\Gamma(1-\frac{\alpha}{2})},$
(2.26)

is the fractional kernel (see [15, 17, 38]).

######  Theorem 2.2 (Case II: $\delta=\infty$).

Let $\delta=\infty$, $\alpha\in(0,2)$, and

$\rho_{\infty}(s)=\frac{C_{\alpha}}{s^{1+\alpha}},\,\;\text{with}\,\;C_{\alpha}=\frac{2^{\alpha-1}\alpha\Gamma(\frac{1+\alpha}{2})}{\sqrt{\pi}\Gamma(1-\frac{\alpha}{2})},$
(2.27)

then the entries of the stiffness matrix $\bm{S}_{\infty}$ in (2.13) reduces
to

$(\bm{S}_{\infty})_{jk}=\frac{\bm{c}_{j}\,(\bm{D}^{k}_{j})^{3-\alpha}\,\bm{c}_{k}^{\top}}{2\Gamma(4-\alpha)\cos(\alpha\pi/2)},$
(2.28)

where the notation $(\bm{D}_{j}^{k})^{3-\alpha}$ performs component-wisely
upon the matrix $\bm{D}_{j}^{k}$, and $\bm{c}_{j}$ and $\bm{D}_{j}^{k}$ are
defined in (2.12) and (2.14), respectively.

With the interchange in the order of integration in the above proof, we can
derive the nonlocal stiffness matrix for $\delta\to\infty.$

######  Theorem 2.3.

Let $\alpha\in(0,2)$, and choose the kernel function (2.25), then the entries
of the stiffness matrix $\bm{S}_{\infty}$ in (2.13) reduces to

$(\bm{S}_{\infty})_{jk}=\begin{dcases}\widehat{C}_{\alpha}\bm{c}_{j}\,(\bm{D}^{k}_{j})^{3-\alpha}\,\bm{c}_{k}^{\top},\quad\alpha\not=1,\\\
\widehat{C}_{\alpha}\bm{c}_{j}\,(\bm{D}^{k}_{j})^{2}\ln{\bm{D}^{k}_{j}}\,\bm{c}_{k}^{\top},\quad\alpha=1,\\\
\end{dcases}\,\;\text{with}\,\;\widehat{C}_{\alpha}=\frac{1}{2\Gamma(4-\alpha)\cos(\alpha\pi/2)},$
(2.29)

where the notation $(\bm{D}_{j}^{k})^{3-\alpha}$ and
$(\bm{D}^{k}_{j})^{2}\ln{\bm{D}^{k}_{j}}$ performs component-wisely upon the
matrix $\bm{D}_{j}^{k}$, and $\bm{c}_{j}$ and $\bm{D}_{j}^{k}$ are defined in
(2.12) and (2.14), respectively.

###### Proof.

Note that when $\delta\to\infty$, we have $\rho_{\delta}(s)=\rho_{\infty}(s)$
is independent of $\delta$. Recall the integral formula (cf. [19, (3.12)]):
for any $\xi\in\mathbb{R}$, we have

$\int_{\mathbb{R}}\frac{1-\cos(\xi s)}{|s|^{1+\alpha}}\,{\rm
d}s=C_{\alpha}^{-1}|\xi|^{\alpha}.$ (2.30)

We change the order of integration in the last equality of (2.16), and derive
from (2.30) that

$\begin{split}(\bm{S}_{\infty})_{jk}&=2C_{\alpha}\int_{0}^{\infty}s^{-1-\alpha}\int_{{\mathbb{R}}}(1-\cos(\xi
s)){\mathscr{F}}[\tilde{\phi}_{j}(x)](\xi)\overline{{\mathscr{F}}[\tilde{\phi}_{k}(x)](\xi)}\,{\rm
d}\xi{\rm d}s\\\
&=C_{\alpha}\int_{{\mathbb{R}}}\Big{(}\int_{{\mathbb{R}}}\frac{1-\cos(\xi
s)}{|s|^{1+\alpha}}\,{\rm
d}s\Big{)}{\mathscr{F}}[\tilde{\phi}_{j}(x)](\xi)\overline{{\mathscr{F}}[\tilde{\phi}_{k}(x)](\xi)}\,{\rm
d}\xi\\\
&=\int_{{\mathbb{R}}}|\xi|^{\alpha}{\mathscr{F}}[\tilde{\phi}_{j}(x)](\xi)\overline{{\mathscr{F}}[\tilde{\phi}_{k}(x)](\xi)}\,{\rm
d}\xi,\end{split}$ (2.31)

which is identical to the FEM stiffness matrix derived in [10, Theorem 1] with
the formula given in (2.30). ∎

The above explicit formula allows for the exact evaluation or accurate
computation of the stiffness matrix for general interaction kernel function
$\rho_{\delta}(s),$ and study the intrinsic properties of this matrix.

###  2.2. FEM nonlocal stiffness matrix with $\delta\leq h$

To study the compatibility to the usual FEM stiffness matrix, when $\delta\to
0$, we consider the special case with

$\delta\leq h=\min\limits_{1\leq j\leq
N}\\{h_{j}\\};\quad\rho_{\delta}(s)=\frac{C_{\delta,\alpha}}{s^{1+\alpha}},\quad
C_{\delta,\alpha}:=\frac{2-\alpha}{\delta^{2-\alpha}},\;\;\;\alpha\in[-1,2),$
(2.32)

where $C_{\delta,\alpha}$ is chosen so that it satisfies the second moment
condition: $\int_{0}^{\delta}s^{2}\rho_{\delta}(s){{\rm d}}s=1.$ We derive
from the formulas in Theorem 2.1 straightforwardly the following interesting
identity, which leads to the usual FEM stiffness matrix when $\delta\to 0.$

######  Corollary 2.1.

Under (2.32), we have the identity

$\bm{S}_{\delta}=\bm{S}_{0}-c_{\alpha}\delta\bm{S}_{0}^{2},\quad
c_{\alpha}:=\frac{2-\alpha}{6(3-\alpha)},$ (2.33)

where $\bm{S}_{0}$ is the usual FEM stiffness matrix for
$-u^{\prime\prime}(x)=f(x)$ on $(a,b)$ with $u(a)=u(b)=0,$ that is,

$\bm{S}_{0}={\rm
diag}\Big{(}\\!\\!-\frac{1}{h_{j}},\frac{1}{h_{j}}+\frac{1}{h_{j+1}},-\frac{1}{h_{j+1}}\Big{)}.$

Consequently, we have $\lim\limits_{\delta\to
0^{+}}\bm{S}_{\delta}=\bm{S}_{0}.$

###### Proof.

As $s\leq\delta\leq h,$ we can simplify $g(z;s)$ in (2.12) and evaluate (2.14)
directly to derive

$\begin{split}\bm{I}_{jk}(s)&=\bm{c}_{j}g(\bm{D}_{j}^{k})\bm{c}_{k}^{\top}\\\
&=\begin{dcases}-2\Big{(}\frac{1}{h_{j}^{2}}+\Big{(}\frac{1}{h_{j}}+\frac{1}{h_{j+1}}\Big{)}^{2}+\frac{1}{h_{j+1}^{2}}\Big{)}s^{3}+12\Big{(}\frac{1}{h_{j}}+\frac{1}{h_{j+1}}\Big{)}s^{2},\;\;&j=k,\\\\[1.0pt]
\frac{2}{h_{j+1}}\Big{(}\frac{1}{h_{j}}+\frac{2}{h_{j+1}}+\frac{1}{h_{j+2}}\Big{)}s^{3}-\frac{12}{h_{j+1}}s^{2},\;\;&j=k-1,\\\\[1.0pt]
\frac{2}{h_{j}}\Big{(}\frac{1}{h_{j-1}}+\frac{2}{h_{j}}+\frac{1}{h_{j+1}}\Big{)}s^{3}-\frac{12}{h_{j}}s^{2},\;\;&j=k+1,\\\\[1.0pt]
-\frac{2}{h_{j+1}h_{j+2}}s^{3},\;\;&j=k-2,\\\\[1.0pt]
-\frac{2}{h_{j-1}h_{j}}s^{3},\;\;&j=k+2,\\\\[1.0pt] 0,\;\;&{\rm
otherwise}.\end{dcases}\end{split}$

Direct calculation using (2.13) shows that $\bm{I}(s)$ is a penta-diagonal
symmetric matrix with the nonzero entries on the main diagonal and two upper
off-diagonals given by

Together with (2.13), we can obtain the explicit form of stiffness matrix for
given nonlocal interaction kernels, e.g., the following nonlocal interaction
kernels

$\rho_{\delta}(s)=\frac{C_{\delta,\alpha}}{s^{1+\alpha}}\,\,\text{with}\,\,C_{\delta,\alpha}=\frac{2-\alpha}{\delta^{2-\alpha}},\quad
s>0,\quad\alpha\in[-1,2),$ (2.34)

where $C_{\delta,\alpha}$ is chosen to satisfy the normalization

$\int_{0}^{\delta}s^{2}\rho_{\delta}(s){\rm d}s=1.$

Then we have Thus for given $\rho_{\delta}(s)$ in (2.32), we find

$\begin{split}(\bm{S}_{\delta})_{jk}&=\int_{0}^{\delta}\bm{I}_{jk}(s)\,\rho_{\delta}(s)\,{\rm
d}s\\\
&=\begin{dcases}-c_{\alpha}\Big{(}\frac{1}{h_{j}^{2}}+\Big{(}\frac{1}{h_{j}}+\frac{1}{h_{j+1}}\Big{)}^{2}+\frac{1}{h_{j+1}^{2}}\Big{)}\delta+\frac{1}{h_{j}}+\frac{1}{h_{j+1}},\;\;&j=k,\\\\[1.0pt]
c_{\alpha}\Big{(}\frac{1}{h_{j}h_{j+1}}+\frac{2}{h^{2}_{j+1}}+\frac{1}{h_{j+1}h_{j+2}}\Big{)}\delta-\frac{1}{h_{j+1}},\;\;&j=k-1,\\\\[1.0pt]
c_{\alpha}\Big{(}\frac{1}{h_{j-1}h_{j}}+\frac{2}{h^{2}_{j}}+\frac{1}{h_{j}h_{j+1}}\Big{)}\delta-\frac{1}{h_{j}},\;\;&j=k+1,\\\\[1.0pt]
-c_{\alpha}\frac{1}{h_{j+1}h_{j+2}}\delta,\;\;&j=k-2,\\\\[1.0pt]
-c_{\alpha}\frac{1}{h_{j-1}h_{j}}\delta,\;\;&j=k+2,\\\\[1.0pt] 0,\;\;&{\rm
otherwise},\end{dcases}\\\
&=(\bm{S}_{0})_{jk}-c_{\alpha}\delta(\bm{S}_{0}^{2})_{jk},\end{split}$ (2.35)

by direct calculaton, we can obtain the identity (2.33). ∎

######  Remark 2.1.

The compatibility of the discretisation is essential for the numerical
solution of nonlocal models. With the evaluation of the stiffness matrix by
the analytic formulas, the FE scheme is automatically compatible [37, 38].

######  Theorem 2.4.

(Limiting case $\delta\rightarrow\infty$) Let

$\rho_{\delta}(s)=\frac{C_{\alpha}}{s^{1+\alpha}},\,\;\text{with}\,\;C_{\alpha}=\frac{2^{\alpha-1}\alpha\Gamma(\frac{1+\alpha}{2})}{\sqrt{\pi}\Gamma(1-\frac{\alpha}{2})},\;\;\;\alpha\in(0,2),$
(2.36)

then the stiffness matrix $\bm{S}_{\delta}$ in (2.13) reduces to

$\begin{split}\bm{S}_{\delta}=\end{split}$ (2.37)

###### Proof.

In this case, $\delta\gg diam(\Omega)$, then we have

$\begin{split}(\bm{S}_{\delta})_{jk}&=(\bm{S}_{\delta})_{kj}=\int_{0}^{\delta}\bm{c}_{j}g(\bm{D}_{j}^{k})\bm{c}_{k}^{\top}\,\rho_{\delta}(s)\,{\rm
d}s=\widehat{C}_{\alpha}\bm{c}_{j}(\bm{D}_{j}^{k})^{3-2\alpha}\bm{c}_{k}^{\top},\end{split}$
(2.38)

where

$\widehat{C}_{\alpha}=C_{\alpha}\int_{0}^{\delta}g(ones(N))\,s^{-(1+\alpha)}\,{\rm
d}s=\frac{1}{2\Gamma(4-2\alpha)\cos(\alpha\pi)}.$ (2.39)

This completes the proof. ∎

$\begin{cases}\mathcal{L}_{\infty}u(x)=f(x),\quad&\text{on}\
\Omega=(a,b),\\\\[4.0pt] u(x)=0,\quad&\text{on}\ \Omega^{c},\end{cases}$
(2.40)

######  Remark 2.2.

When $\delta\to\infty$, the stiffness matrix $\bm{S}_{\delta}$ in (2.13)
recovers to the FEM stiffness matrix for the integral fractional Laplacian in
[10, Theorem 1].

###  2.3. Nonlocal FEM stiffness matrix on uniform meshes

As a special case of Theorem 2.1, the nonlocal stiffness matrix of FEM on
uniform meshes reduces to a symmetric Toeplitz matrix.

######  Theorem 2.5.

Given a uniform partition $\\{x_{j}=a+jh\\}_{j=0}^{N+1}$ with $h=(b-a)/(N+1)$,
the nonlocal stiffness matrix $\bm{S}_{\delta}$ is a symmetric Toeplitz matrix
of the form

$\bm{S}_{\delta}=\begin{bmatrix}t_{0}&t_{1}&t_{2}&\cdots&t_{N-3}&t_{N-2}&t_{N-1}\\\\[1.0pt]
t_{1}&t_{0}&t_{1}&\ddots&\cdots&t_{N-3}&t_{N-2}\\\\[-1.0pt]
t_{2}&t_{1}&t_{0}&\ddots&\ddots&\vdots&t_{N-3}\\\\[0.0pt]
\vdots&\ddots&\ddots&\ddots&\ddots&\ddots&\vdots\\\\[2.0pt]
t_{N-3}&\vdots&\ddots&\ddots&\hskip 8.0ptt_{0}&t_{1}&t_{2}\\\\[-2.0pt]
t_{N-2}&t_{N-3}&\cdots&\ddots&\hskip 8.0ptt_{1}&t_{0}&t_{1}\\\\[3.0pt]
t_{N-1}&t_{N-2}&t_{N-3}&\cdots&\hskip 8.0ptt_{2}&t_{1}&t_{0}\end{bmatrix},$
(2.41)

i.e., $\\{(\bm{S}_{\delta})_{jk}=t_{|j-k|}\\}_{j,k=1}^{N}.$ Here, the
generating vector $\bm{t}=(t_{0},t_{1},\cdots,t_{N-1})$ can be evaluated by

$t_{p}=\begin{dcases}h^{2}\int_{0}^{\frac{\delta}{h}}{\mathbb{I}}_{p}(\tau)\,\rho_{\delta}(\tau
h)\,{\rm d}\tau,\quad&0\leq p<\delta/h+2,\\\
0,\quad&p\geq\delta/h+2,\end{dcases}$ (2.42)

where $p=|j-k|$ and

$\begin{split}\mathbb{I}_{p}(\tau):=-\frac{1}{12}\sum\limits_{i=-2}^{2}\eta_{i}\big{\\{}|p+i+\tau|^{3}-2|p+i|^{3}+|p+i-\tau|^{3}\big{\\}},\end{split}$
(2.43)

with the constants $\eta_{0}=6$, $\eta_{\pm 1}=-4$ and $\eta_{\pm 2}=1.$

###### Proof.

In this case, the entries in (2.14) can be simplified into

$\begin{split}{\bm{I}}_{jk}(s)&={\bm{I}}_{kj}(s)=\bm{c}_{j}g(\bm{D}_{j}^{k})\bm{c}_{k}^{\top}\\\
&=\frac{1}{h^{2}}\big{(}1,-2,1\big{)}\begin{pmatrix}g(|j-k|h)&g(|j-k-1|h)&g(|j-k-2|h)\\\\[2.0pt]
g(|j-k+1|h)&g(|j-k|h)&g(|j-k-1|h)\\\\[2.0pt]
g(|j-k+2|h)&g(|j-k+1|h)&g(|j-k|h)\end{pmatrix}\big{(}1,-2,1\big{)}^{\top}\\\
&=\frac{1}{h^{2}}\sum\limits_{i=-2}^{2}\eta_{i}g(|j-k+i|h)=\frac{1}{h^{2}}\sum\limits_{i=-2}^{2}\eta_{i}g((j-k+i)h)\\\
&=\frac{1}{h^{2}}\sum\limits_{i=-2}^{2}\eta_{i}g((|j-k|+i)h),\end{split}$

where we observed from the definition (2.12) that $g$ is an even function in
$z$ for fixed $s\geq 0,$ and also noted that the symmetric properties of
coefficients $\\{\eta_{i}\\}$ and summation in $i.$ Thus we can denote

$\begin{split}\bm{I}_{jk}(s)&=\frac{1}{h^{2}}\sum\limits_{i=-2}^{2}\eta_{i}g((p+i)h):={\mathbb{I}}_{p}(s)={\mathbb{I}}_{|j-k|}(s),\end{split}$

which implies the entries on each diagonal are the same, i.e., $\bm{I}(s)$ is
a symmetric Toeplitz matrix, so is $\bm{S}_{\delta}.$ The expression
(2.51)-(2.52) is a direct consequence of (2.12)-(2.13).

If $p\geq\delta/h+2,$ we find from (2.52) that for $\tau\in[0,\delta/h],$

$\begin{split}{\mathbb{I}}_{p}(\tau)&=-\frac{1}{12}\sum\limits_{i=-2}^{2}\eta_{i}\big{\\{}|p+i+\tau|^{3}-2|p+i|^{3}+|p+i-\tau|^{3}\big{\\}}\\\
&=-\frac{1}{12}\sum\limits_{i=-2}^{2}\eta_{i}\big{\\{}(p+i+\tau)^{3}-2(p+i)^{3}+(p+i-\tau)^{3}\big{\\}}=0.\end{split}$

This completes the proof. ∎

######  Remark 2.3.

_For any $\tau\geq 0$ and fixed integer $p\geq 0,$ ${\mathbb{I}}_{p}(\tau)$ in
(2.52) has the following piecewise representation:_

  * (i)

_If_ $p=0,1,$ then we have

${\mathbb{I}}_{0}(\tau)=\begin{cases}12(2-\tau)\tau^{2},\quad&\tau\in[0,1),\\\\[1.0pt]
16+4(\tau-2)^{3},\;&\tau\in[1,2),\\\\[1.0pt] 16,\quad&\tau\geq 2,\end{cases}$
(2.44)

_and_

${\mathbb{I}}_{1}(\tau)=\begin{cases}-4(3-2\tau)\tau^{2},\quad&\tau\in[0,1),\\\\[1.0pt]
14-42\tau+30\tau^{2}-6\tau^{3},\;\;&\tau\in[1,2),\\\\[1.0pt]
4+2(\tau-3)^{3},\quad&\tau\in[2,3),\\\\[1.0pt] 4,\quad&\tau\geq 3.\end{cases}$
(2.45)

  * (ii)

_If_ $p\geq 2,$ then we have

${\mathbb{I}}_{p}(\tau)=\begin{cases}2(p-2-\tau)^{3},\quad&\tau\in[p-2,\,p-1),\\\\[1.0pt]
2(p-2-\tau)^{3}-8(p-1-\tau)^{3},\;\;&\tau\in[p-1,\,p),\\\\[1.0pt]
8(p+1-\tau)^{3}-2(p+2-\tau)^{3},\;&\tau\in[p,\,p+1),\\\\[1.0pt]
-2(p+2-\tau)^{3},\quad&\tau\in[p+1,\,p+2),\\\\[1.0pt] 0,\quad&{\rm
otherwise}.\end{cases}$ (2.46)

Case II ($\delta\to\infty$): Consider the following fractional Laplacian
kernel

$\rho_{\infty}(s)=C_{\alpha}s^{-(1+2\alpha)},\quad s>0,\quad\alpha\in(0,2),$

with
$C_{\alpha}=\frac{4^{\alpha}\alpha\Gamma(\frac{1}{2}+\alpha)}{\sqrt{\pi}\Gamma(1-\alpha)}$
and $\Gamma(\cdot)$ in the Gamma function. Then we have

$\begin{split}(\bm{S}_{\infty})_{jk}&=(\bm{S}_{\infty})_{kj}=\bm{c}_{j}\Big{(}\int_{0}^{\infty}g(\bm{D}_{j}^{k})\,\rho_{\infty}(s)\,{\rm
d}s\Big{)}\bm{c}_{k}^{T}\\\
&=\widehat{C}_{\alpha}\bm{c}_{j}(\bm{D}_{j}^{k})^{3-2\alpha}\bm{c}_{k}^{T},\quad
1\leq j,\,k\leq N.\end{split}$ (2.47)

where

$\widehat{C}_{\alpha}=C_{\alpha}\int_{0}^{\infty}g(ones(N))\,s^{-(1+2\alpha)}\,{\rm
d}s=\frac{1}{2\Gamma(4-2\alpha)\cos(\alpha\pi)}.$ (2.48)

Correspondingly, we define piecewise linear FEM basis associated with uniform
partition

$\psi_{j}(x)=\begin{cases}\frac{x-x_{j-1}}{h},\quad&{\rm
if}\;\;x\in(x_{j-1},x_{j}],\\\\[2.0pt] \frac{x_{j+1}-x}{h},\quad&{\rm
if}\;\;x\in(x_{j},x_{j+1}],\\\\[2.0pt] 0,\quad&{\rm
elsewhere\,\,on\,\,\mathbb{R}}.\end{cases}$ (2.49)

According to (2.13), we can easily obtain the entries of the nonlocal
stiffness matrix on uniform mesh, which enjoys the Toeplitz structure.

######  Corollary 2.2.

For uniform partition $\\{x_{j}=a+jh\\}$, the nonlocal stiffness matrix
$\bm{S}_{\delta}$ is a symmetric Toeplitz matrix given by

$\bm{S}_{\delta}=h^{2}\begin{bmatrix}t_{0}&t_{1}&t_{2}&\cdots&t_{N-3}&t_{N-2}&t_{N-1}\\\\[1.0pt]
t_{1}&t_{0}&t_{1}&\ddots&\cdots&t_{N-3}&t_{N-2}\\\\[-1.0pt]
t_{2}&t_{1}&t_{0}&\ddots&\ddots&\vdots&t_{N-3}\\\\[0.0pt]
\vdots&\ddots&\ddots&\ddots&\ddots&\ddots&\vdots\\\\[2.0pt]
t_{N-3}&\vdots&\ddots&\ddots&\hskip 8.0ptt_{0}&t_{1}&t_{2}\\\\[-2.0pt]
t_{N-2}&t_{N-3}&\cdots&\ddots&\hskip 8.0ptt_{1}&t_{0}&t_{1}\\\\[3.0pt]
t_{N-1}&t_{N-2}&t_{N-3}&\cdots&\hskip 8.0ptt_{2}&t_{1}&t_{0}\end{bmatrix},$
(2.50)

is generated by the vector $(t_{0},t_{1},\cdots,t_{N-1})$ in the first row or
column of $\bm{S}_{\delta}$ with

$\begin{split}t_{p}=\int_{0}^{\frac{\delta}{h}}{\mathbb{I}}_{p}(\tau)\,\rho_{\delta}(\tau
h)\,{\rm
d}\tau=\frac{1}{h}\int_{0}^{\delta}{\mathbb{I}}_{p}\big{(}\frac{s}{h}\big{)}\,\rho_{\delta}(s)\,{\rm
d}s,\end{split}$ (2.51)

where

$\begin{split}{\mathbb{I}}_{p}(\tau):=-\frac{1}{12}\sum\limits_{i=-2}^{2}c_{i}\big{\\{}|p+i+\tau|^{3}-2|p+i|^{3}+|p+i-\tau|^{3}\big{\\}},\end{split}$
(2.52)

and the constants $c_{\pm 2}=1,c_{\pm 1}=-4,c_{0}=6$.

###### Proof.

One can follow the same procedure as in the proof of Theorem 2.1 to obtain the
desired result. ∎

††margin:  Can simplify uniform grid case! Whether can be directly derived
from Theorem 2.1?

######  Proposition 2.1.

For any $\tau\geq 0$ and fixed integer $p\geq 0,$ ${\mathbb{I}}_{p}(\tau)$ in
(2.52) has the following piecewise representation:

  * (i)

If $p=0,1,$ then we have

${\mathbb{I}}_{0}(\tau)=\begin{cases}12(2-\tau)\tau^{2},\quad&\tau\in[0,1),\\\\[1.0pt]
16+4(\tau-2)^{3},\;&\tau\in[1,2),\\\\[1.0pt] 16,\quad&\tau\geq 2,\end{cases}$
(2.53)

and

${\mathbb{I}}_{1}(\tau)=\begin{cases}-4(3-2\tau)\tau^{2},\quad&\tau\in[0,1),\\\\[1.0pt]
14-42\tau+30\tau^{2}-6\tau^{3},\;\;&\tau\in[1,2),\\\\[1.0pt]
4+2(\tau-3)^{3},\quad&\tau\in[2,3),\\\\[1.0pt] 4,\quad&\tau\geq 3.\end{cases}$
(2.54)

  * (ii)

If $p\geq 2,$ then we have

${\mathbb{I}}_{p}(\tau)=\begin{cases}2(p-2-\tau)^{3},\quad&\tau\in[p-2,\,p-1),\\\\[1.0pt]
2(p-2-\tau)^{3}-8(p-1-\tau)^{3},\;\;&\tau\in[p-1,\,p),\\\\[1.0pt]
8(p+1-\tau)^{3}-2(p+2-\tau)^{3},\;&\tau\in[p,\,p+1),\\\\[1.0pt]
-2(p+2-\tau)^{3},\quad&\tau\in[p+1,\,p+2),\\\\[1.0pt] 0,\quad&{\rm
otherwise}.\end{cases}$ (2.55)

######  Remark 2.4.

We observe from (2.51) that the bandwidth of the above matrix (2.50) is
directly determined by the ratio of $\delta/h$, which plays an important role
in studying the sparsity properties of the stiffness matrix. A particular case
is when $\delta\leq h$ and the fractional interaction kernels (2.36), we can
derive the explicit form of the entries

$\begin{split}(\bm{S}_{\delta})_{kj}&=(\bm{S}_{\delta})_{jk}=\frac{h^{2}}{12}\int_{0}^{\delta/h}{\mathbb{I}}_{p}(\tau)\,\rho_{\delta}(\tau
h)\,{\rm
d}\tau=\frac{1}{h}\begin{cases}2-\frac{2-\alpha}{3-\alpha}\frac{\delta}{h},\;\;&{\rm
if}\;\;p=0,\\\\[1.0pt]
-1+\frac{2(2-\alpha)}{3(3-\alpha)}\frac{\delta}{h},\;\;&{\rm
if}\;\;p=1,\\\\[1.0pt]
-\frac{2-\alpha}{6(3-\alpha))}\frac{\delta}{h},\;\;&{\rm
if}\;\;p=2,\\\\[1.0pt] 0,\;\;&{\rm otherwise},\end{cases}\end{split}$ (2.56)

which implies that as $\delta\rightarrow 0$ the matrix $\bm{S}_{\delta}$ in
(2.56) reduces to the usual (tri-diagonal) FEM stiffness matrix
$\bm{S}_{0}=\frac{1}{h}{\rm diag}(-1,2,-1)$.

######  Remark 2.5.

Denote $r=\frac{\delta}{h}$, then we can obtain the explicit form of stiffness
matrix for given fractional interaction kernels (2.36)

$\begin{split}(\bm{S}_{\delta})_{kj}&=(\bm{S}_{\delta})_{jk}=\frac{h^{2}}{12}\int_{0}^{r}{\mathbb{I}}_{|k-j|}(\tau)\,\rho_{\delta}(\tau
h)\,{\rm
d}\tau={\color[rgb]{0,0,1}\frac{1}{h}\begin{cases}2-\frac{2-\alpha}{3-\alpha}\frac{\delta}{h},\;\;&{\rm
if}\;\;j=k,\\\\[1.0pt] 0,\;\;&{\rm otherwise}.\end{cases}}\end{split}$ (2.57)

##  3\. Numerical Studies and Discussions

In this section, we first conduct a numerical study of the conditioning of the
stiffness matrix on three sets of typical non-uniform grids used in practice
to resolve singular solutions. We then numerically solve several types of
nonlocal problems by using the FEM nonlocal matrix to show the efficiency and
accuracy of the proposed scheme.

  * •

Two sets of typical non-uniform grids: graded & Geometric (can show the
advantage of exact computation)!

  * •

Numerical results include (i) accuracy test; (ii) the use of non-uniform grids

###  3.1. Conditioning of the stiffness matrix on nonuniform meshes

An issue of interest is the condition number estimation for the FEM nonlocal
stiffness matrix. It is a well-studied topic in the uniform meshes case, but
much less known in this nonuniform meshes setting. In general, the mesh
geometry affects not only the approximation error of the finite element
solution but also the spectral properties of the corresponding stiffness
matrix.

Extensive research has been conducted on the condition number of the nonlocal
stiffness matrix on uniform grids (cf. [4, 42]), which the entries of the
nonlocal stiffness matrix can be computed exactly based on the definition of
hypersingular integrals form therein. With the aid of the analytical
expressions, the following condition number bounds for uniform mesh were
reported in [42, p.1772]:

$\text{Cond}(\bm{S_{\delta}})\leq\min\\{N^{\alpha}\delta^{\alpha-2},\,\,N^{2}\\}.$
(3.1)

However, the study on the condition number of the nonlocal stiffness matrix on
the non-uniform mesh is still open. Inspired by the result of the condition
number for FEM stiffness on nonuniform mesh in [20, (26)], we predict the
condition number of FEM stiffness matrix (2.13) on nonuniform meshes
associated with fractional kernel (2.32) is

$\text{\rm
Cond}(\bm{S_{\delta}})\leq\min\\{\left(h_{\max}/h_{\min}\right)^{\alpha-1}N^{\alpha}\delta^{\alpha-2},\,\,\left(h_{\max}/h_{\min}\right)^{2}N^{2}\\},$
(3.2)

where $h_{\max}=\max_{1\leq j\leq N}h_{j}$ and $h_{\min}=\min_{1\leq j\leq
N}h_{j}$. It is clear that the above prediction can reduce to the result with
uniform mesh (3.1) when $h_{\max}=h_{\min}$.

We will explore this numerically on graded, geometric and Shishkin meshes in
order to demonstrate the prediction or conjecture. In the following numerical
experiment, we will choose the mesh division as follows:

  * (i)

Graded meshes: A graded meshes on the interval $(a,b)$ is defined as follows

$x_{j}=\begin{cases}a+\dfrac{b-a}{2}\Big{(}\dfrac{2j}{N}\Big{)}^{\gamma},&j=0,1,\cdots,N/2-1,\\\\[6.0pt]
b-\dfrac{b-a}{2}\Big{(}2-\dfrac{2j}{N}\Big{)}^{\gamma},\;\;&j=N/2,N/2+1,\cdots,N.\end{cases}$
(3.3)

From (3.2), the condition number of FEM stiffness matrix on the above graded
meshes is

$\text{Cond}(\bm{S_{\delta}})\sim N^{\gamma(\alpha-1)+1}.$ (3.4)

  * (ii)

Geometric meshes: A geometric meshes on the interval $(a,b)$ is defined
according to

$x_{j}=\begin{cases}a+q^{N-j}\dfrac{b-a}{2},&\quad
j=1,2,\cdots,N-1,\\\\[4.0pt] b-q^{j-N}\dfrac{b-a}{2},&\quad
j=N,N+1,\cdots,2N-1,\end{cases}$ (3.5)

where $q$ ($0<q<1$, $q$ is independent of $j$) is subdivision ratio.

Using (3.2), the condition number of FEM stiffness matrix on the above
geometric meshes is

$\text{Cond}(\bm{S_{\delta}})\sim q^{1-N}N^{\alpha}.$ (3.6)

  * (iii)

Shishkin meshes: A Shishkin meshes is piecewise uniform mesh over $[a,b]$ with
$2M+N$ elements. For $\eta\in(0,\frac{1}{2})$, we refer to the mesh size over
$[0,\eta]$ and $[1-\eta,\,1]$ as $h=\frac{\eta}{M}(b-a)$ and to the mesh size
over $[\eta,\,1-\eta]$ as $H=\frac{1-2\eta}{N}(b-a)$.

According to (3.2), the condition number of FEM stiffness matrix on the above
Shishkin meshes is

$\text{Cond}(\bm{S_{\delta}})\sim\Big{(}\frac{M(1-2\eta)}{N\eta}\Big{)}^{\alpha-1}N^{\alpha}.$
(3.7)

In Figure 3.1 (a)-(c), we plot the condition number of the stiffness matrix
against various $N$, fixed horizon $\delta=0.3$ and different $\alpha$. We
observe a good agreement between the numerical results and (3.2), which
indicates a clear dependence of the condition number on three parameters: the
mesh ratio $h_{\max}/h_{\min}$, the size of the nonlocality horizon parameter
$\delta$, and index of the fractional kernel function $\alpha$.

(a) $\text{Cond}(\bm{S_{\delta}})\sim N^{\gamma(\alpha-1)+1}$.

(b) $\text{Cond}(\bm{S_{\delta}})\sim q^{1-N}N^{\alpha}$.

(c)
$\text{Cond}(\bm{S_{\delta}})\sim(\frac{2(1-2\eta)}{\eta})^{\alpha-1}N^{\alpha}$.

(d) $\lambda_{\min}\sim N^{-1}$.

Figure 3.1. Conditioning and the smallest eigenvalue of the stiffness matrix
$\bm{S}_{\delta}$ for
$\rho_{\delta}(s)=\frac{2-\alpha}{\delta^{2-\alpha}}s^{-(1+\alpha)}$ on
different meshes with $\delta=0.3$. Condition number: (a) Graded meshes with
$\gamma=2$, (b) Geometric meshes with $q=0.999$, (c) Shishkin meshes with
$M=2N$ and $\eta=0.2$, (d) The smallest eigenvalue on graded meshes with
$\gamma=2$.

Meanwhile, we plot the first smallest eigenvalue of the stiffness matrix
against various $N$ and different $\alpha$ in Figure 3.1 (d), where we
observed that the smallest eigenvalue of $\bm{S}_{\delta}$ satisfy

$\lambda_{\min}(\bm{S}_{\delta})=ch_{\max}=cN^{-1},\quad\alpha\in[1,2].$

This result is consistent with the classical usual Laplacian (i.e.,
$\delta=0$) and the integral fractional Laplacian (i.e., $\delta\to\infty$)
(cf. [10]).

####  3.1.1. Diagonal dominance

The study of diagonal dominance of the stiffness matrix has been a research
subject of a longstanding interest in numerical linear algebra and numerical
analysis. The diagonal dominance property of the stiffness matrix associated
with fractional Laplacian operators has been proved in [29], which can be
viewed as a limiting case of this paper (cf. [15, 17, 38]). In what follows,
we restrict our attention to the uniform meshes case and discuss the intrinsic
property of the stiffness matrix as $h\to 0$, that is, we extended the
following diagonal dominance property of the stiffness matrix when $\delta\leq
h$.

######  Theorem 3.6.

For $\delta\leq h$, the FEM stiffness matrix $\bm{S}_{\delta}$ associated with
fractional kernel function (2.36) on uniform meshes is diagonally dominant and
off-diagonal entries are negative.

###### Proof.

For $\delta\leq h$, let $r=\frac{\delta}{h}$, then $0<r<1$. From (2.56), we
have

$\displaystyle(\bm{S}_{\delta})_{kj}$
$\displaystyle=(\bm{S}_{\delta})_{jk}=\frac{1}{h}\begin{cases}2-\dfrac{2-\alpha}{3-\alpha}r,\qquad\qquad\text{if}\
j=k,\\\ -1+\dfrac{2(2-\alpha)}{3(3-\alpha)}r,\qquad\text{if}\ |j-k|=1,\\\
-\dfrac{2-\alpha}{6(3-\alpha)}r,\qquad\qquad\text{if}\ |j-k|=2,\\\
0,\qquad\qquad\qquad\qquad\ \text{otherwise.}\end{cases}$ (3.8)

Define

$d_{k}:=|(\bm{S}_{\delta})_{kk}|-\sum_{k\neq
j,\,j=1}^{N-1}|(\bm{S}_{\delta})_{kj}|.$

For $\alpha\in[-1,2)$, we can easily obtain $0<\frac{2-\alpha}{3-\alpha}<1$
and

$(\bm{S}_{\delta})_{kk}>0,\quad(\bm{S}_{\delta})_{kj}<0\quad\text{for}\quad|j-k|=1,2.$

It knows that $\bm{S}_{\delta}$ is a symmetric Toeplitz matrix from Corollary
2.1 and $\bm{S}_{\delta}$ is a five-diagonal matrix, then we have

$d_{1}=d_{N-1},\quad d_{2}=d_{N-2},\quad d_{3}=\cdots=d_{N-3}.$

Next, we will consider $d_{1},d_{2},d_{3}$, respectively.

$\displaystyle d_{1}$
$\displaystyle=|(\bm{S}_{\delta})_{11}|-|(\bm{S}_{\delta})_{12}|-|(\bm{S}_{\delta})_{13}|=(\bm{S}_{\delta})_{11}+(\bm{S}_{\delta})_{12}+(\bm{S}_{\delta})_{13}=\frac{1}{h}\bigg{[}1-\frac{2-\alpha}{2(3-\alpha)}r\bigg{]}>0,$
$\displaystyle d_{2}$
$\displaystyle=|(\bm{S}_{\delta})_{22}|-2|(\bm{S}_{\delta})_{21}|-|(\bm{S}_{\delta})_{24}|=(\bm{S}_{\delta})_{22}+2(\bm{S}_{\delta})_{21}+(\bm{S}_{\delta})_{24}=\frac{1}{h}\frac{2-\alpha}{6(3-\alpha)}r>0,$
$\displaystyle d_{3}$
$\displaystyle=|(\bm{S}_{\delta})_{33}|-2|(\bm{S}_{\delta})_{31}|-2|(\bm{S}_{\delta})_{32}|=(\bm{S}_{\delta})_{33}+2(\bm{S}_{\delta})_{31}+2(\bm{S}_{\delta})_{32}=0.$

That is $d_{1}=d_{N-1}>0,\ d_{2}=d_{N-2}>0,\ d_{3}=\cdots=d_{N-3}=0$, so the
stiffness matrix $\bm{S}_{\delta}$ is diagonally dominant and off-diagonal
entries are negative. ∎

###  3.2. Error analysis

In this part, we will give the error analysis for the piecewise linear FEM in
analogous to results widely known for local problems.

Following the same notation as in [36], we introduce the energy space

$\mathcal{V}(\Omega)=\Big{\\{}u\in
L^{2}(\Omega\cup\Omega_{\mathcal{\delta}}):\int_{\Omega}\int_{-\delta}^{\delta}(u(x+s)-u(x))^{2}\rho_{\delta}(|s|){\rm
d}s{\rm d}x<\infty\Big{\\}},$ (3.9)

equipped with the semi-norm

$|u|_{\mathcal{V}}=\Big{\\{}\int_{\Omega}\int_{-\delta}^{\delta}(u(x+s)-u(x))^{2}\rho_{\delta}(|s|){\rm
d}s{\rm d}x\Big{\\}}^{1/2},$ (3.10)

and the norm

$\|u\|_{\mathcal{V}}=(|u|_{\mathcal{V}}^{2}+\|u\|^{2})^{1/2},\quad\|u\|^{2}=\int_{\Omega}|u(x)|^{2}{\rm
d}x.$ (3.11)

The following bound plays an important part in the error analysis.

######  Lemma 3.1.

(see [37], Lemma 3.1) For $\alpha\in(0,2]$ and a kernel $\rho_{\delta}$
satisfying $|s|^{\alpha}\rho_{\delta}(|s|)\in L^{1}(\mathbb{R})$, we have a
constant $C$ depending only on $\Omega$ such that

$|u|^{2}_{\mathcal{V}}\leq
C\Big{(}\int_{\Omega}|s|^{\alpha}\rho_{\delta}(|s|){\rm
d}s\Big{)}\|u\|^{2}_{H^{\alpha/2}(\Omega)}.$ (3.12)

The following regularity result is already derived in [12].

######  Lemma 3.2.

(see [12], Theorem 6.2) Let $m$ be a non-negative integer and $0<r<1$. Suppose
that $u\in\mathcal{V}(\Omega)\bigcap H^{m+t-r}(\Omega)$, where $r\leq t\leq
1$. Then, there exists a constant $C>0$ with value independent of $h$ such
that, for sufficiently small $h$,

$\|u-u_{h}\|_{H^{r}(\Omega)}\leq Ch^{m+t-r}\|u\|_{H^{m+t}(\Omega)}.$ (3.13)

In particular, if $m=1$, then second-order convergence with respect to the
$L^{2}(\Omega)$ norm is obtain for piecewise linear elements if $u$ is
sufficiently smooth, i.e.,

$\|u-u_{h}\|_{L^{2}(\Omega)}\leq Ch^{2}\|u\|_{H^{2}(\Omega)}.$ (3.14)

For $u\in\mathcal{V}(\Omega)$, let us denote
$P_{h}u\in\mathcal{V}^{h}(\Omega)$ be the Ritz-projection of $u$ given by

$\mathcal{A}_{\delta}(P_{h}u,\upsilon_{h})=\mathcal{A}_{\delta}(u,\upsilon_{h}),\quad\forall\upsilon_{h}\in\mathcal{V}^{h}(\Omega).$

With the aid of Lemma 3.1, we can derive the error estimate of the Ritz-
projection $P_{h}u$ as follows.

######  Lemma 3.3.

(see [12], Theorem 6.2) Let $m$ be a non-negative integer and $0<\mu<1$.
Suppose that $u\in\mathcal{V}(\Omega)\bigcap H^{m+t-\mu}(\Omega)$, where
$\mu\leq t\leq 1$. Then, there exists a constant $C>0$ with value independent
of $h$ such that, for sufficiently small $h$,

$|u-P_{h}u|_{\mathcal{V}}\leq\min\limits_{\upsilon_{h}\in\mathcal{V}_{h}}|u-\upsilon_{h}|_{\mathcal{V}}\leq
Ch^{m+t-\mu}\|u\|_{H^{m+t}(\Omega)}.$ (3.15)

In particular, if $m=1$, then second-order convergence with respect to the
$L^{2}(\Omega)$ norm is obtain for piecewise linear elements if $u$ is
sufficiently smooth, i.e.,

$|u-P_{h}u|_{\mathcal{V}}\leq\min\limits_{\upsilon_{h}\in\mathcal{V}_{h}}|u-\upsilon_{h}|_{\mathcal{V}}\leq
Ch^{2}\|u\|_{H^{2}(\Omega)}.$ (3.16)

###### Proof.

Using the definition of the Ritz-projection $P_{h}u$, we have

$\mathcal{A}_{\delta}(u-P_{h}u,\upsilon_{h})=0,\quad\forall\upsilon_{h}\in\mathcal{V}^{h}(\Omega),$

Using the orthogonality and Cauchy-Schwartz inequality, we have

$\begin{split}\mathcal{A}_{\delta}(u-P_{h}u,u-P_{h}u)&=\mathcal{A}_{\delta}(u-P_{h}u,u-\upsilon_{h}+\upsilon_{h}-P_{h}u)=\mathcal{A}_{\delta}(u-P_{h}u,u-\upsilon_{h})\\\
&\leq\mathcal{A}_{\delta}(u-P_{h}u,u-P_{h}u)^{1/2}\mathcal{A}_{\delta}(u-\upsilon_{h},u-\upsilon_{h})^{1/2},\quad\forall\upsilon_{h}\in\mathcal{V}^{h}(\Omega),\end{split}$

which leads to

$\mathcal{A}_{\delta}(u-P_{h}u,u-P_{h}u)\leq\min\limits_{\upsilon_{h}\in\mathcal{S}_{h}}\mathcal{A}_{\delta}(u-\upsilon_{h},u-\upsilon_{h}).$

Using Lemma 3.1, we have

$|u-P_{h}u|_{\mathcal{V}}\leq\min\limits_{\upsilon_{h}\in\mathcal{V}_{h}}|u-\upsilon_{h}|_{\mathcal{V}}.$

This completes the proof. ∎

######  Theorem 3.7.

Let the solution $u$ of the nonlocal problem (1.1) be sufficiently smooth,
then for the piecewise linear elements solutions $u_{h}$, we have

$|u-u_{h}|_{\mathcal{V}}\leq{\color[rgb]{0,0,1}Ch^{2}\|u\|_{H^{2}(\Omega)}.}$
(3.17)

###### Proof.

Decompose the error $e_{h}=u-u_{h}$ as $e_{h}=\eta_{h}-\xi_{h}$, where

$\eta_{h}=u-P_{h}u\quad\text{and}\quad\xi_{h}=u_{h}-P_{h}u.$

Subtracting (2.6) to (2.9), we have

$\mathcal{A}_{\delta}(u-u_{h},\upsilon_{h})=\mathcal{A}_{\delta}(\eta_{h}-\xi_{h},\upsilon_{h})=0,\quad\forall\upsilon_{h}\in\mathcal{S}^{h}(\Omega).$
(3.18)

Choosing $\upsilon_{h}=\xi_{h}$ and applying the Cauchy-Schwartz inequality,
we obtain

$\mathcal{A}_{\delta}(\xi_{h},\xi_{h})\leq\mathcal{A}_{\delta}(\eta_{h},\eta_{h}).$

Thus,

$|\xi_{h}|_{\mathcal{V}}\leq|\eta_{h}|_{\mathcal{V}}.$

Using the triangle inequality and (3.14), we finish the proof. ∎

###  3.3. Nonlocal BVPs with various interaction kernels

In this subsection, we focus on numerical approximation the following nonlocal
problems

$\begin{cases}{\mathcal{L}}_{\delta}u(x)=\lambda
u(x)+f(x),\;\;\;&\text{on}\;\Omega,\\\
u(x)=0,&\text{on}\;\Omega_{\delta},\end{cases}$ (3.19)

which can be considered as nonlocal BVP if $\lambda=0$; or the nonlocal
Helmholtz problem if $\lambda<0$ is given; or eigenvalue problem if $\lambda$
is unknown and $f(x)=0$.

####  3.3.1. Nonlocal BVP with $\lambda=0$

In this subsection, we will consider the accuracy test and limiting behaviours
for nonlocal BVPs (3.19) with $\lambda=0$.

Example 1. (Accuracy test of smooth solutions). We first test the accuracy of
the proposed method for the nonlocal BVP (3.19) with $\lambda=0$ and involving
a fractional kernel function (2.32) with different values of the parameters
$\alpha$.

In the simulation, we use the FEM on uniform meshes, and take the exact
solution $u(x)=x^{2}(1-x)^{2}$ and $\Omega=(0,1)$ (cf. [9, 36]). In Figure
3.2, we plot the numerical errors against various $h$ with different
parameters $\alpha=0.1,0.3,0.5,0.7,0.9$. We observe from Figure 3.2 that the
convergence order of the piecewise linear FEM on uniform meshes is $O(h^{2})$
for fixed horizon $h\leq\delta\leq\frac{b-a}{2}$, while for fixed horizon
$0\leq\delta\leq h$, the convergence rate is $O(h)$.

(a) $h\leq\delta\leq\frac{b-a}{2}$.

(b) $0<\delta<h$.

Figure 3.2. The numerical errors for
$\rho_{\delta}(s)=\frac{2-\alpha}{\delta^{2-\alpha}}s^{-(1+\alpha)}$ with
different $\alpha\in(-1,0)$.

Example 2. (Accuracy test of discontinuous solution). Next, we consider the
FEM on uniform and graded meshes for the nonlocal BVP (3.19) with $\lambda=0$
and the discontinuous exact solution on $\Omega=(0,1)$ given by (cf. [9, 36])

$u(x)=\begin{cases}2x^{2},\quad&x<0.5,\\\\[1.0pt]
(1-x)^{2},\quad&\text{otherwise}.\end{cases}$ (3.20)

In this computation, we take the constant box potential
$\rho_{\delta}(s)=\frac{3}{\delta^{3}}$, and choose the following graded
meshes

$x_{j}=\begin{cases}\dfrac{b+a}{2}-\dfrac{b-a}{2}\Big{(}1-\dfrac{2j}{N}\Big{)}^{\gamma},&j=0,1,\cdots,\dfrac{N}{2}-1,\\\\[6.0pt]
\dfrac{b+a}{2}+\dfrac{b-a}{2}\Big{(}\dfrac{2j}{N}-1\Big{)}^{\gamma},\;\;&j=\dfrac{N}{2},\dfrac{N}{2}+1,\cdots,N,\end{cases}$
(3.21)

where $\gamma>1$ is a suitable grading exponent. Different from (3.3), the
above graded meshes clustering near the discontinuous points $x=0.5$, which is
beneficial for capturing the discontinuous nature of the solution.

In Figure 3.3(a)-(b), we plot the numerical errors of the FEM on uniform
meshes with various $\delta$, which shows that the convergence order is
$O(h^{0.5})$ in $L^{\infty}$-norm and $O(h)$ in $L^{2}$-norm. As a comparison,
we plot numerical errors of the FEM on graded meshes with various $\delta$ in
Figure 3.3(c)-(d). We observe that, for fixed grading exponent
$1\leq\gamma\leq 2$, the $L^{2}$-norm errors are roughly of
$O(h^{\frac{\gamma}{2}})$. For fixed grading exponent $\gamma\geq 2$, i.e.,
the $L^{2}$-norm errors are roughly of $O(h^{\min(2,\gamma-1)})$. The
convergence order of FEM on graded meshes is significantly higher than that of
FEM on uniform meshes when there are discontinuities in the solution.

(a) $L^{\infty}$-norm errors.

(b) $L^{2}$-norm errors.

(c) $1\leq\gamma<2$.

(d) $\gamma\geq 2$.

Figure 3.3. The numerical errors of discontinuous solution. (a)-(b): FEM on
uniform meshes, (c)-(d): FEM on graded meshes in $L^{2}$-norm.

(a) uniform meshes.

(b) graded meshes with $\gamma=2$.

Figure 3.4.  Analytical and numerical solutions for
$\rho_{\delta}(s)=\frac{3}{\delta^{3}}$.

It should be pointed out that no instability has been observed in any of our
computations. The piecewise linear FEM we use is not only convergent to the
nonlocal constrained value problem with fixed $\delta$ but also convergent to
the local differential equation with a fixed ratio between $\delta$ and $h$.
Thus, our scheme is an asymptotically compatible scheme [37]. In addition, we
find from Figure 3.4 that the numerical solutions of FEM on graded meshes have
a better approximation effect and stability than the FEM on uniform meshes at
the discontinuity point.

Example 3. (Limiting behaviours). It is interesting to numerically study the
local limit (i.e., $\delta\to 0$) and the nonlocal interactions approach
infinite (i.e., $\delta\to\infty$) of the nonlocal BVP (3.19) with
$\lambda=0$. More specifically, we compare solutions of the nonlocal BVP
(3.19) with $\lambda=0$ to the ODE (resp. fractional Laplacian equation)
equation with special interest in observing behaviors in the local limit
$\delta\to 0$ (resp. fractional limit $\delta\to\infty$).

Case 1 (local limit): To study the local limit (i.e., $\delta\to 0$), we adopt
the piecewise linear FEM for solving the nonlocal BVP (3.19) with $\lambda=0$,
source function $f(x)=2$ and kernel function (2.32). The exact solution is
unknown, but the following ODE

$\begin{cases}-u^{\prime\prime}(x)=2,\;&\text{for}\ x\in(-1,1),\\\\[4.0pt]
u(-1)=u(1)=0,\end{cases}$ (3.22)

has the exact solution $u(x)=1-x^{2}$.

In our computation, we use the graded meshes (3.3). We compare solutions of
the nonlocal model to the ODE (3.22) with special interest in observing
behaviors in the local limit $\delta\to 0$. Figure 3.10 shows the comparison
of numerical solutions for the nonlocal model and the solution of the
governing equation given by (3.22) for $\gamma=2$ and $\gamma=3$,
respectively. It can be clearly observed that as $\delta\to 0$, the solution
of the nonlocal diffusion model converges to the solution of the local one
(3.22).

(a) $\gamma=2$.

(b) $\gamma=3$.

(c) $\gamma=2$.

(d) $\gamma=3$.

Figure 3.5. The evolution process of the numerical solutions on graded meshes
as $\delta\rightarrow 0$. Top: $\rho_{\delta}(s)=\frac{3}{\delta^{3}}$,
Bottom: $\rho_{\delta}(s)=\frac{2-\alpha}{\delta^{2-\alpha}}s^{-(1+\alpha)}$
and $\alpha=-0.25$.

Figure 3.6. The evolution process of the numerical solutions as
$\delta\rightarrow 0$ with
$\rho_{\delta}(s)=\frac{2-\alpha}{\delta^{2-\alpha}}s^{-(1+\alpha)}$ and
$\alpha=-0.75$ on graded meshes. Left: $\gamma=2$, Right: $\gamma=3$.

Case 2 (global limit): Now, we turn to the global limit (i.e.,
$\delta\rightarrow\infty$), that is, consider the nonlocal models with the
following fractional Laplacian kernel (2.32). It is known that the solution of
the nonlocal problem with a fractional Laplacian kernel exhibits singularities
near the boundary of $\Omega$. We take $\Omega=(-1,1)$ and consider a
piecewise linear FEM on graded meshes (3.3) for solving the nonlocal BVP
(3.19) with $\lambda=0$ and $f(x)=1$. The exact solution is unknown, but the
following fractional Poisson equation

$\begin{cases}(-\Delta)^{\alpha}u(x)=1,\quad&x\in\Omega=(-1,1),\\\\[4.0pt]
u(x)=0,\quad&x\in\Omega^{c}=\mathbb{R}\backslash\bar{\Omega},\end{cases}$
(3.23)

where

$(-\Delta)^{\alpha}u(x)=C_{\alpha}\,{\rm
p.v.}\\!\int_{\mathbb{R}}\frac{u(x)-u(y)}{|x-y|^{1+2\alpha}}{\rm d}y,\;\;\;\;$
(3.24)

with “p.v.” stands for the principle value. The equation (3.23) admits a exact
solution $u(x)=(1-x^{2})^{\alpha}_{+}/{\Gamma(2\alpha+1)}$.

As we can see from Figure 3.7, the solution of the nonlocal BVP converges to
the solution of the fractional Laplacian equation (3.23) as the nonlocal
interactions become infinite, i.e., $\delta\to\infty$. Thus the conclusions in
[17] are verified numerically.

(a) $\gamma=2$.

(b) $\gamma=3$.

Figure 3.7. The evolution process of the numerical solutions on graded meshes
with $\gamma=1.1$ as $\delta\rightarrow\infty$.

####  3.3.2. Nonlocal BVP with $\lambda<0$

In this subsection, we will consider the dynamics of nonlocal BVPs (3.19) with
$\lambda<0$.

Example 4. (The dynamics of nonlocal Helmholtz problem). In this example, we
will consider the FEM on uniform meshes for the following nonlocal Helmholtz
problem

$\begin{cases}\mathcal{L}_{\delta}u(x)-k^{2}n(x)u(x)=f(x),\;&\text{on}\;\Omega=(-L,L),\\\
u(x)=0,&\text{on}\;\Omega_{\delta}=[-L-\delta,L+\delta].\end{cases}$ (3.25)

where $n(x)$ is a given function and $k^{2}$ is a given constant. Here, we
focus on the dynamics of the solution with the following two cases:

Case 1: $n(x)$ is continuous function

$n(x)=\sin{x},\quad x\in(-L,L).$ (3.26)

Case 2: $n(x)$ is piecewise constant

$n(x)=\begin{cases}0.5,\;&x\in(-L,0),\\\\[1.0pt] 1,\;&x\in(0,L).\end{cases}$
(3.27)

For the convenience of description, we take $f(x)=k^{2}$ and the kernel
function (2.32), i.e.,
$\rho_{\delta}(s)=\frac{2-\alpha}{\delta^{2-\alpha}}s^{-(1+\alpha)}$ with
fixed $\alpha=0.5$.

In Figure 3.8, we depict the profile of the numerical solution for Case 1 on
the computation domain $\Omega=[-4\pi,4\pi]$ with different $k^{2}$ and
$\delta$. We observe from the top of Figure 3.8 that oscillation in regions
with $\sin(x)<0$ exponential growth and decay in regions with $\sin(x)>0$.
This observation is consistent with the one described in [39, p.81] for the
local model. Moreover, we find from the bottom of Figure 3.8 that as $k^{2}$
decreases, the oscillation frequency within the oscillatory region amplifies.

(a) $\delta=0.1,\,k^{2}=1000/3$ (b) $\delta=0.001,\,k^{2}=1000/3$ (c)
$\delta=0.0001,\,k^{2}=1000/3$ (d) $k^{2}=10/3,\,\delta=0.001$

(e) $k^{2}=1000/3,\,\delta=0.001$ (f) $k^{2}=10000/3,\,\delta=0.001$

Figure 3.8. The profiles of numerical solutions for Case 1. Top: versus
different $\delta$ with fixed $k^{2}=1000/3$; Bottom: versus different $k^{2}$
with fixed $\delta=0.001$.

In Figure 3.9, we plot the profile of the numerical solution for Case 2 on the
computation domain $\Omega=[-100,100]$ with $k^{2}=2$. It can be seen that the
amplitude of the solution is affected by magnitude of the piecewise constant
function $n(x)$, that is, the amplitude is larger within the interval
$(-100,0)$, while it is relatively smaller within the interval $(0,100)$. We
also observe that as $\delta$ decreases, the solution converges to the
solution of the local Helmholtz problem, i.e., $\delta=0$.

(a) $\delta=0.1,\,k^{2}=2$ (b) $\delta=0.01,\,k^{2}=2$ (c)
$\delta=0,\,k^{2}=2$

Figure 3.9. The profiles of numerical solution for Case 2 with $\alpha=0.5$
and $k^{2}=2$.

####  3.3.3. Nonlocal BVP with $f(x)=0$ and $\lambda$ is unknown

One fundamental issue in the study of the nonlocal operator is its eigenvalue
problem, i.e., (3.19) with $f(x)=0$ and $\lambda$ is unknown. Since the
spectrum of nonlocal eigenvalue problems depends on the horizon size $\delta$
and the parameters $\alpha$, few analytical results of the eigenvalues or
eigenfunctions can be found in the literature, which prompts us to numerical
study the nonlocal eigenvalue problem.

Example 5. (Eigenvalue problem) Consider the eigenvalue problem of nonlocal
problem with a fractional kernel function (2.32) with different values of the
parameters $\alpha$.

$\begin{cases}\mathcal{L}_{\delta}u(x)=\lambda u(x),\quad&\text{on}\
\Omega=(-1,1),\\\\[4.0pt] u(x)=0,\quad&\text{on}\
\Omega_{\mathcal{\delta}}=[-1-\delta,-1]\cup[1,1+\delta].\end{cases}$ (3.28)

Using FEM discretization, we can obtain the following standard matrix
eigenvalue problem

$\bm{S}_{\delta}\bm{u}=\lambda\bm{M}\bm{u},$ (3.29)

where $\bm{M}$ is the usual FEM mass matrix.

When $\delta\to 0$, $\mathcal{L}_{\delta}u(x)$ is approximating Laplacian
operator $u_{xx}(x)$. As is well known, the Laplacian operator eigenvalue
problem with Dirichlet boundary conditions is

$\begin{cases}-u_{xx}(x)=\lambda
u(x),\;&\text{on}\,\,\Omega=(-1,1),\\\\[4.0pt] u(-1)=u(1)=0.\end{cases}$
(3.30)

Its eigenvalues and eigenfunction expressions are

$\lambda_{k}=\frac{k^{2}\pi^{2}}{4},\quad u_{k}(x)=\sin(\frac{k\pi}{2}(1+x)).$

Table 3.1. The eigenvalues of FEM on graded meshes when $N=1024$ . |  | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\lambda_{4}$ | $\lambda_{5}$
---|---|---|---|---|---|---
| $\delta=0$ | 2.46740110 | 9.86960440 | 22.20660990 | 39.47841760 | 61.68502751
$\alpha$ | $\delta=0.1$ | $\lambda^{1}_{\delta}$ | $\lambda^{2}_{\delta}$ | $\lambda^{3}_{\delta}$ | $\lambda^{4}_{\delta}$ | $\lambda^{5}_{\delta}$
| $\delta/2^{1}$ | 2.41759862 | 9.66390067 | 21.71945040 | 38.55295196 | 60.11616374
| $\delta/2^{2}$ | 2.44244850 | 9.76813715 | 21.97209805 | 39.04605489 | 60.97843581
0.4 | $\delta/2^{3}$ | 2.45491393 | 9.81923809 | 22.09172241 | 39.27027858 | 61.35199331
| $\delta/2^{5}$ | 2.46115693 | 9.84452507 | 22.14980400 | 39.37647857 | 61.52385393
| $\delta/2^{5}$ | 2.46428142 | 9.85710596 | 22.17843509 | 39.42816470 | 61.60622225
| $\delta/2^{6}$ | 2.46584534 | 9.86339180 | 22.19270072 | 39.45381799 | 61.64691606
| $\delta/2^{1}$ | 2.42751075 | 9.70453885 | 21.81459327 | 38.73025554 | 60.41328634
| $\delta/2^{2}$ | 2.44745354 | 9.78841357 | 22.01868176 | 39.13126116 | 61.11637305
0.8 | $\delta/2^{3}$ | 2.45742959 | 9.82936677 | 22.11476158 | 39.31185528 | 61.41820329
| $\delta/2^{4}$ | 2.46241914 | 9.84959215 | 22.16127618 | 39.39704646 | 61.55634551
| $\delta/2^{5}$ | 2.46491475 | 9.85964614 | 22.18417827 | 39.43844094 | 61.62241698
| $\delta/2^{6}$ | 2.46616355 | 9.86466801 | 22.19558573 | 39.45897877 | 61.65504623
| $\delta/2^{1}$ | 2.43870151 | 9.75058814 | 21.92302394 | 38.93499308 | 60.75718811
| $\delta/2^{2}$ | 2.45296532 | 9.81079299 | 22.07028410 | 39.22610053 | 61.27079425
1.2 | $\delta/2^{3}$ | 2.46007728 | 9.84004411 | 22.13911423 | 39.35596052 | 61.48875749
| $\delta/2^{4}$ | 2.46362860 | 9.85445490 | 22.17231554 | 39.41690918 | 61.58786858
| $\delta/2^{5}$ | 2.46540354 | 9.86161116 | 22.18863941 | 39.44646660 | 61.63515483
| $\delta/2^{6}$ | 2.46629141 | 9.86518432 | 22.19676697 | 39.46112478 | 61.65849527
| $\delta/2^{1}$ | 2.41843091 | 9.67125122 | 21.75106140 | 38.64553780 | 60.33751981
| $\delta/2^{2}$ | 2.42618623 | 9.70412710 | 21.83198239 | 38.80666199 | 60.62389058
1.6 | $\delta/2^{3}$ | 2.43003864 | 9.72000928 | 21.86949303 | 38.87775885 | 60.74384622
| $\delta/2^{4}$ | 2.43195879 | 9.72781211 | 21.88751329 | 38.91094086 | 60.79800753
| $\delta/2^{5}$ | 2.43291768 | 9.73168259 | 21.89635993 | 38.92700088 | 60.82378576
| $\delta/2^{6}$ | 2.43339712 | 9.73361430 | 21.90076261 | 38.93496166 | 60.83650391

We will test the accuracy of the FEM discretization for eigenvalue of the
nonlocal operator (1.2) with a fractional kernel function of the form

$\rho_{\delta}(s)=\frac{2-\alpha}{\delta^{2-\alpha}}s^{-(1+\alpha)},\quad
s>0,\quad\alpha\in[-1,2),$

with different values of the parameters $\alpha$, which lead to different
types of singularities.

Table 3.2. The eigenvalues on uniform meshes when $N=1024$. $\alpha$ | $\delta$ | 0 | $10^{-4}$ | $10^{-3}$ | $10^{-2}$ | $10^{-1}$ | 1
---|---|---|---|---|---|---|---
$-0.5$ | $\lambda_{1}$ | 2.4674 | 2.4673 | 2.4668 | 2.4564 | 2.3384 | 1.4074
$\lambda_{2}$ | 9.8606 | 9.8694 | 9.8671 | 9.8251 | 9.3226 | 4.4369
$-0.75$ | $\lambda_{1}$ | 2.4674 | 2.4673 | 2.4668 | 2.4558 | 2.3320 | 1.3706
$\lambda_{2}$ | 9.8696 | 9.8694 | 9.8670 | 9.8227 | 9.2957 | 4.2549
$0.5$ | $\lambda_{1}$ | 2.4674 | 2.4673 | 2.4668 | 2.4564 | 2.3384 | 1.4074
$\lambda_{2}$ | 9.8696 | 9.8694 | 9.8675 | 9.8274 | 9.4741 | 5.4788
$1.5$ | $\lambda_{1}$ | 2.4674 | 2.4673 | 2.4668 | 2.4564 | 2.3384 | 1.4074
$\lambda_{2}$ | 9.8696 | 9.8694 | 9.8671 | 9.8251 | 9.3226 | 4.4369
Table 3.3. The errors of eigenvalues on graded meshes when $N=1024$ . $\alpha$ | $\delta=0.1$ | $|\lambda_{1}-\lambda^{1}_{\delta}|$ | $|\lambda_{2}-\lambda^{2}_{\delta}|$ | $|\lambda_{3}-\lambda^{3}_{\delta}|$ | $|\lambda_{4}-\lambda^{4}_{\delta}|$ | $|\lambda_{5}-\lambda^{5}_{\delta}|$
---|---|---|---|---|---|---
| $\delta/2^{1}$ | 4.98024805e-02 | 2.05703733e-01 | 4.87159503e-01 | 9.26512343e-01 | 1.56886377e+00
| $\delta/2^{2}$ | 2.49526043e-02 | 1.01467251e-01 | 2.34511852e-01 | 4.32362714e-01 | 7.06591693e-01
0.4 | $\delta/2^{3}$ | 1.24871709e-02 | 5.03663062e-02 | 1.14887493e-01 | 2.08139023e-01 | 3.33034195e-01
| $\delta/2^{4}$ | 6.24417136e-03 | 2.50793326e-02 | 5.68059049e-02 | 1.01939032e-01 | 1.61173579e-01
| $\delta/2^{5}$ | 3.11968427e-03 | 1.24984377e-02 | 2.81748121e-02 | 5.02529059e-02 | 7.88052573e-02
| $\delta/2^{6}$ | 1.55576053e-03 | 6.21259956e-03 | 1.39091786e-02 | 2.45996104e-02 | 3.81114443e-02
| $\delta/2^{1}$ | 3.98903502e-02 | 1.65065555e-01 | 3.92016628e-01 | 7.48162060e-01 | 1.27174116e+00
| $\delta/2^{2}$ | 1.99475595e-02 | 8.11908268e-02 | 1.87928141e-01 | 3.47156447e-01 | 5.68654457e-01
0.8 | $\delta/2^{3}$ | 9.97151249e-03 | 4.02376321e-02 | 9.18483274e-02 | 1.66562326e-01 | 2.66824215e-01
| $\delta/2^{4}$ | 4.98195687e-03 | 2.00122486e-02 | 4.53337268e-02 | 8.13711451e-02 | 1.28681992e-01
| $\delta/2^{5}$ | 2.48634990e-03 | 9.95825886e-03 | 2.24316343e-02 | 3.99766639e-02 | 6.26105263e-02
| $\delta/2^{6}$ | 1.23755302e-03 | 4.93639234e-03 | 1.10241763e-02 | 1.94388314e-02 | 2.99812714e-02
| $\delta/2^{1}$ | 2.86995852e-02 | 1.19016259e-01 | 2.83585958e-01 | 5.43424529e-01 | 9.27839399e-01
| $\delta/2^{2}$ | 1.44357763e-02 | 5.88114115e-02 | 1.36325797e-01 | 2.52317070e-01 | 4.14233255e-01
1.2 | $\delta/2^{3}$ | 7.32381792e-03 | 2.95602930e-02 | 6.74956675e-02 | 1.22457082e-01 | 1.96270021e-01
| $\delta/2^{4}$ | 3.77250484e-03 | 1.51494988e-02 | 3.42943660e-02 | 6.15084242e-02 | 9.71589232e-02
| $\delta/2^{5}$ | 1.99755832e-03 | 7.99324535e-03 | 1.79704905e-02 | 3.19510088e-02 | 4.98726797e-02
| $\delta/2^{6}$ | 1.10969274e-03 | 4.42007585e-03 | 9.84292928e-03 | 1.72928194e-02 | 2.65322374e-02
| $\delta/2^{1}$ | 4.89701911e-02 | 1.98353178e-01 | 4.55548502e-01 | 8.32879809e-01 | 1.34750769e+00
| $\delta/2^{2}$ | 4.12148682e-02 | 1.65477295e-01 | 3.74627512e-01 | 6.71755614e-01 | 1.06113693e+00
1.6 | $\delta/2^{3}$ | 3.73624593e-02 | 1.49595119e-01 | 3.37116875e-01 | 6.00658756e-01 | 9.41181290e-01
| $\delta/2^{4}$ | 3.54423138e-02 | 1.41792288e-01 | 3.19096612e-01 | 5.67476741e-01 | 8.87019975e-01
| $\delta/2^{5}$ | 3.44834248e-02 | 1.37921811e-01 | 3.10249976e-01 | 5.51416724e-01 | 8.61241744e-01
| $\delta/2^{6}$ | 3.40039801e-02 | 1.35990099e-01 | 3.05847296e-01 | 5.43455943e-01 | 8.48523597e-01
Table 3.4. The convergence rates of FEM on graded meshes when $N=1024$. $\alpha$ | $\delta=0.1$ | $Rate_{1}$ | $Rate_{2}$ | $Rate_{3}$ | $Rate_{4}$ | $Rate_{5}$
---|---|---|---|---|---|---
| $\delta/2^{2}$ | 0.9970 | 1.0196 | 1.0547 | 1.0996 | 1.1508
| $\delta/2^{3}$ | 0.9987 | 1.0105 | 1.0294 | 1.0547 | 1.0852
0.4 | $\delta/2^{4}$ | 0.9999 | 1.0060 | 1.0161 | 1.0298 | 1.0471
| $\delta/2^{5}$ | 1.0011 | 1.0048 | 1.0116 | 1.0204 | 1.0323
| $\delta/2^{6}$ | 1.0038 | 1.0085 | 1.0184 | 1.0306 | 1.0481
| $\delta/2^{2}$ | 0.9998 | 1.0237 | 1.0607 | 1.1078 | 1.1612
| $\delta/2^{3}$ | 1.0003 | 1.0128 | 1.0329 | 1.0595 | 1.0917
0.8 | $\delta/2^{4}$ | 1.0011 | 1.0077 | 1.0187 | 1.0335 | 1.0521
| $\delta/2^{5}$ | 1.0027 | 1.0069 | 1.0151 | 1.0254 | 1.0393
| $\delta/2^{6}$ | 1.0065 | 1.0124 | 1.0249 | 1.0402 | 1.0623
| $\delta/2^{2}$ | 0.9914 | 1.0170 | 1.0567 | 1.1068 | 1.1634
| $\delta/2^{3}$ | 0.9790 | 0.9924 | 1.0142 | 1.0430 | 1.0776
1.2 | $\delta/2^{4}$ | 0.9571 | 0.9644 | 0.9768 | 0.9934 | 1.0144
| $\delta/2^{5}$ | 0.9173 | 0.9224 | 0.9323 | 0.9449 | 0.9621
| $\delta/2^{6}$ | 0.8481 | 0.8547 | 0.8685 | 0.8857 | 0.9105
| $\delta/2^{2}$ | 0.2487 | 0.2614 | 0.2821 | 0.3102 | 0.3447
| $\delta/2^{3}$ | 0.1416 | 0.1456 | 0.1522 | 0.1614 | 0.1731
1.6 | $\delta/2^{4}$ | 0.0761 | 0.0773 | 0.0793 | 0.0820 | 0.0855
| $\delta/2^{5}$ | 0.0396 | 0.0399 | 0.0406 | 0.0414 | 0.0425
| $\delta/2^{6}$ | 0.0202 | 0.0203 | 0.0206 | 0.0210 | 0.0215

Tables 3.1 provide numerical eigenvalues of FEM on graded meshes. A large
amount of experimental results can observe three phenomena: 1) When $\delta\to
0$, the eigenvalues of nonlocal operator tend to the corresponding eigenvalues
of local Laplacian operator (3.30); 2) The eigenvalues of local Laplacian
operator (3.30) are the upper bound of the corresponding eigenvalues of
nonlocal operator; 3) The smaller the horizon parameter $\delta$, the larger
the eigenvalues of nonlocal operator.

(1)$\delta=0.1$.

(2)$\delta=0.3$.

(3)$\delta=0.5$.

(4)$\delta=0.7$.

Figure 3.10. The different states of on uniform meshes for different
$\delta$.(The first row is the ground states, and the second row is the first
excited state)

(1)$\delta=0.1$.

(2)$\delta=0.3$.

(3)$\delta=0.5$.

(4)$\delta=0.7$.

Figure 3.11. The different states of on non-uniform meshes for different
$\delta$.(The first row is the ground states, and the second row is the first
excited state)

As shown in the Figure 3.10-3.11, whether using a uniform grid or a non-
uniform grid, it becomes evident that as $\delta\to 0$, the eigenvalues and
eigenfunctions of nonlocal operator (1.2) all convergence to those of the
classical Laplace operator.

(1)$\delta=0.01$.

(2)$\delta=0.1$.

(3)$\delta=0.5$.

(1)$\delta=0.01$.

(2)$\delta=0.1$.

(3)$\delta=0.5$.

Figure 3.12. The first eigenvalue for different $\delta$ on two
meshes.($\lambda_{\delta}^{1}$ represents the first numerical eigenvalue,
$\lambda_{1}$ represents the first eigenvalue of the exact solution)

According to the Figure (3.12), we observe that as $\delta$ increases, the
error between the numerical solution of the eigenvalue and the exact solution
becomes progressively larger. However, the calculated numerical solutions are
consistently smaller than the exact solutions, allowing us to establish a
lower bound for the eigenvalue obtained through this method. We can get the
following conclusions ($\delta\geq 0.001$???):

$\lambda_{\delta}^{1}\leq\lambda_{1}.$

Table 3.5. The errors and rates when $\alpha=0.5$ and $N=128$. $\alpha$ | $\delta=0.1$ | $|\lambda_{1}-\lambda^{1}_{\delta}|$ | $Rate$ | $|\lambda_{2}-\lambda^{2}_{\delta}|$ | $Rate$
---|---|---|---|---|---
| $\delta/2^{1}$ | 4.03e-02 |  | 1.68e-01 |
| $\delta/2^{2}$ | 1.65e-02 | 1.288 | 6.76e-02 | 1.313
0.25 | $\delta/2^{3}$ | 7.13e-03 | 1.210 | 2.85e-02 | 1.244
| $\delta/2^{4}$ | 3.29e-03 | 1.118 | 1.24e-02 | 1.201
| $\delta/2^{5}$ | 1.51e-03 | 1.125 | 4.92e-03 | 1.336
| $\delta/2^{1}$ | 3.17e-02 |  | 1.32e-01 |
| $\delta/2^{2}$ | 1.31e-02 | 1.280 | 5.34e-02 | 1.310
0.75 | $\delta/2^{3}$ | 5.36e-03 | 1.285 | 2.13e-02 | 1.330
| $\delta/2^{4}$ | 2.07e-03 | 1.370 | 7.46e-03 | 1.510
| $\delta/2^{5}$ | 5.33e-04 | 1.960 | 9.74e-04 | 2.937
| $\delta/2^{1}$ | 3.65e-02 |  | 1.50e-01 |
| $\delta/2^{2}$ | 2.40e-02 | 0.606 | 9.66e-02 | 0.634
1.25 | $\delta/2^{3}$ | 1.85e-02 | 0.375 | 7.36e-02 | 0.392
| $\delta/2^{4}$ | 1.61e-02 | 0.206 | 6.33e-02 | 0.219
| $\delta/2^{5}$ | 1.49e-02 | 0.109 | 5.83e-02 | 0.117
| $\delta/2^{1}$ | 5.51e-01 |  | 2.21e+00 |
| $\delta/2^{2}$ | 5.46e-01 | - | 2.19e+00 | -
1.75 | $\delta/2^{3}$ | 5.43e-01 | - | 2.17e+00 | -
| $\delta/2^{4}$ | 5.43e-01 | - | 2.17e+00 | -
| $\delta/2^{5}$ | 5.42e-01 | - | 2.17e+00 | -

###  3.4. Time-dependent problem

Consider the following nonlocal Allen-Cahn equation:

$\begin{cases}u_{t}(x,t)+\epsilon^{2}\mathcal{L}_{\delta}u(x,t)+f(u)=0,\quad&\text{on}\
\Omega=(a,b),\quad 0<t\leq T,\\\\[4.0pt] u(x,t)=0,\quad&\text{on}\
\Omega_{\mathcal{\delta}}=[a-\delta,a]\cup[b,b+\delta],\quad 0<t\leq
T,\\\\[4.0pt] u(x,0)=u_{0}(x),\quad&\text{on}\ \Omega=(a,b).\end{cases}$
(3.31)

where $u_{0}(x)$ is the initial data, $f(u)$ is given function on
${\mathbb{R}}$ satisfies $f(u)=F^{\prime}(u)$ with
$F(u)=\frac{1}{4}(u^{2}-1)^{2}$. In the above, $\mathcal{L}_{\delta}$ is the
nonlocal operator (1.2) and $\epsilon$ is a constant.

Let $N$ be a positive integer and $\tau=T/N$ be a time step size. We can get
the time grid $t_{n}=n\tau,0\leq n\leq N$ and use
$u_{h}^{n}\in\mathcal{V}_{h}^{0}$ to approximate $u_{h}(x,t_{n})$,
$n=0,1,\cdots,N$. Using backward Euler method to approximate $u_{t}(x,t_{n})$,
we can get

$u_{t}(x,t_{n})\approx\delta_{t}u_{h}^{n}=\frac{1}{\tau}(u^{n}_{h}-u^{n-1}_{h}),\quad
x\in\Omega,\quad 1\leq n\leq N.$

Here, we employ the fully discrete scheme with the semi-implicit scheme in
time and FEM discretization in space for solving the problem (3.31): find
$u_{h}^{n}\in H^{1}(0,T;\mathcal{V}_{h}^{0}(\Omega))$, such that

$\begin{cases}(\delta_{t}u_{h}^{n},v_{h})+\epsilon^{2}\mathcal{A}_{\delta}\left(u_{h}^{n},v_{h}\right)+(f(u_{h}^{n-1}),v_{h})=0,~{}\forall\,v_{h}\in\mathcal{V}_{h}^{0}(\Omega),\quad
1\leq n\leq N,\\\
(u_{h}^{0},v_{h})=(u_{0}(x),v_{h}),\quad\forall\,v_{h}\in\mathcal{V}_{h}^{0}(\Omega).\end{cases}$
(3.32)

Then, it is easy to obtain the matrix form of (3.32) as follows:

$(\bm{M}+\tau\epsilon^{2}\bm{S_{\delta}})\bm{U^{n}}=\bm{M}\bm{U^{n-1}}-\tau\bm{F^{n-1}},$
(3.33)

where $\bm{M}$ is the usual (tridiagonal) FEM mass matrix, the stiffness
matrix $\bm{S_{\delta}}$ is defined in (2.13), and $\bm{F^{n-1}}$ is defined
as

$\bm{F}^{n-1}=\int_{\Omega}f(u_{h}^{n-1})\bm{\Phi}\,{\rm
d}x,\;\;\;\bm{\Phi}=\big{(}\varphi_{1}(x),\varphi_{2}(x),\cdots,\varphi_{N-1}(x)\big{)}^{\top}.$

In this example, we adopt the proposed fully discrete scheme (3.32) to
simulate the phase evolution behavior and energy dissipation. We take
$\Omega=(-1,1)$, $\epsilon=0.01$, the initial date $u_{0}(x)=e^{-100x^{2}}$,
and take the fractional kernel function of the form (2.32).

Table 3.6. The errors and time rates of semi-implicit schemes when $M=1024$. $\alpha$ | $N$ | $\delta=0.01$ | $\delta=0.1$
---|---|---|---
$Error$ | $Rate_{t}$ | $Error$ | $Rate_{t}$
| $2^{3}$ | 2.45e-03 |  | 2.55e-03 |
| $2^{4}$ | 1.23e-03 | 0.995 | 1.28e-03 | 0.993
0.25 | $2^{5}$ | 6.16e-04 | 0.998 | 6.43e-04 | 0.997
| $2^{6}$ | 3.08e-04 | 0.999 | 3.22e-04 | 0.999
| $2^{7}$ | 1.54e-04 | 0.999 | 1.61e-04 | 0.999
| $2^{3}$ | 2.45e-03 |  | 2.54e-03 |
| $2^{4}$ | 1.23e-03 | 0.995 | 1.27e-03 | 0.994
0.75 | $2^{5}$ | 6.15e-04 | 0.998 | 6.38e-04 | 0.997
| $2^{6}$ | 3.08e-04 | 0.998 | 3.20e-04 | 0.999
| $2^{7}$ | 1.54e-04 | 0.999 | 1.60e-04 | 0.999
| $2^{3}$ | 2.45e-03 |  | 2.52e-03 |
| $2^{4}$ | 1.23e-03 | 0.995 | 1.26e-03 | 0.994
1.25 | $2^{5}$ | 6.16e-04 | 0.998 | 6.33e-04 | 0.997
| $2^{6}$ | 3.08e-04 | 0.999 | 3.17e-04 | 0.999
| $2^{7}$ | 1.54e-04 | 0.999 | 1.59e-04 | 0.999
| $2^{3}$ | 2.57e-03 |  | 2.60e-03 |
| $2^{4}$ | 1.29e-03 | 0.995 | 1.30e-03 | 0.994
1.75 | $2^{5}$ | 6.45e-04 | 0.997 | 6.53e-04 | 0.997
| $2^{6}$ | 3.23e-04 | 0.999 | 3.27e-04 | 0.999
| $2^{7}$ | 1.61e-04 | 0.999 | 1.63e-04 | 0.999
Table 3.7. The errors and space rates of semi-implicit schemes when $N=1024$. $\alpha$ | $M$ | $\delta=0.01$ | $\delta=0.1$
---|---|---|---
$Error$ | $Rate_{s}$ | $CPU(s)$ | $Error$ | $Rate_{s}$ | $CPU(s)$
| $2^{6}$ | 9.57e-03 |  | 0.145 | 9.90e-03 |  | 0.154
| $2^{7}$ | 2.38e-03 | 2.006 | 0.409 | 2.42e-03 | 2.034 | 0.359
0.25 | $2^{8}$ | 6.02e-04 | 1.985 | 1.869 | 6.05e-04 | 1.999 | 1.925
| $2^{9}$ | 1.51e-04 | 1.998 | 8.405 | 1.51e-04 | 2.000 | 8.535
| $2^{10}$ | 3.77e-05 | 1.999 | 47.98 | 3.78e-05 | 2.000 | 49.09
| $2^{6}$ | 9.56e-03 |  | 0.142 | 9.87e-03 |  | 0.142
| $2^{7}$ | 2.38e-03 | 2.006 | 0.360 | 2.41e-03 | 2.032 | 0.380
0.75 | $2^{8}$ | 6.01e-04 | 1.986 | 1.974 | 6.04e-04 | 1.997 | 1.884
| $2^{9}$ | 1.50e-04 | 1.997 | 8.415 | 1.51e-04 | 2.000 | 8.461
| $2^{10}$ | 3.77e-05 | 1.998 | 47.17 | 3.78e-05 | 1.999 | 47.11
| $2^{6}$ | 9.55e-03 |  | 0.142 | 9.81e-03 |  | 0.140
| $2^{7}$ | 2.38e-03 | 2.006 | 0.380 | 2.40e-03 | 2.029 | 0.350
1.25 | $2^{8}$ | 5.99e-04 | 1.988 | 1.884 | 6.04e-04 | 1.993 | 1.952
| $2^{9}$ | 1.50e-04 | 1.996 | 8.461 | 1.51e-04 | 2.000 | 8.400
| $2^{10}$ | 3.76e-05 | 1.997 | 47.11 | 3.78e-05 | 1.999 | 47.10
| $2^{6}$ | 9.64e-03 |  | 0.200 | 9.78e-03 |  | 0.211
| $2^{7}$ | 2.38e-03 | 2.017 | 0.438 | 2.40e-03 | 2.027 | 0.544
1.75 | $2^{8}$ | 6.00e-04 | 1.989 | 2.565 | 6.04e-04 | 1.991 | 2.262
| $2^{9}$ | 1.51e-04 | 1.997 | 12.73 | 1.51e-04 | 1.999 | 13.02
| $2^{10}$ | 3.76e-05 | 1.997 | 63.23 | 3.79e-05 | 1.997 | 56.96

(a) $\delta=0.01$.

(b) $\delta=0.1$.

(c) $\delta=0.01$.

(d) $\delta=0.1$.

Figure 3.13. Convergence order of FEM with different $\delta$. (a)-(b) are the
time convergence order, (c)-(d) are the spatial convergence order.

In Figure 3.13, we plot the numerical errors and the corresponding convergence
rates of the proposed method max errors for different $\alpha$ and $\delta$.
As we can see the temporal convergence rate is $O(\tau)$, and the space
convergence rate is $O(h^{2})$. Moreover, we display the waveform of the
numerical solution at the different $\delta$ and time, see Figures 3.14. The
dates indicate that the $\delta$ affects the propagation velocity of the
solitary wave. For a smaller $\delta$, the propagation of the soliton becomes
slower, thus indicating the presence of quantum sub-diffusion. Finally, we
plot the evolution of maximum value at various times with different $\alpha$
and $\delta$ in Figure 3.15, which shows that the maximum principle is
preserved numerically. We observe from Figure 3.15 (d) that the maximum value
of the steady state is increased as $\alpha$ increases.

(a) $\delta=0.02$ (b) $\delta=0.2$ (c) $\delta=2$ (d) $t=10$

Figure 3.14.  The waveforms of the numerical solution with the $\alpha=1.25$
on nonuniform meshes. (a)-(c) are the snapshots of $u(x)$ are taken at
$t=1,5,10$. (d) is the snapshots of $u(x)$ are taken at $t=10$ for different
$\delta$.

(a) $\alpha=-0.25$.

(b) $\alpha=0.25$.

(c) $\alpha=0.75$.

(d) $\delta=0.1$.

Figure 3.15. Evolution of maximum value as $\delta$ increases on nonuniform
meshes. (a)-(c) are the maximum value at $\alpha=-0.25,0.25,0.75$. (d) is the
maximum value at $\delta=0.1$ for different $\alpha$.

we consider the following two-grid system based on finite element method,
which includes the space coarse mesh $H$ in $S_{H}$ and the space fine mesh
$h$ in $S_{h}$. Now, we use the two-grid method and give its fully discrete
scheme.
Step 1: We choose the coarse meshes which the step size is
$H_{j}(j=1,2,\cdots,J)$ in $S_{H}$, then we have

$a=x^{H}_{0}<x^{H}_{1}<\cdots<x^{H}_{J-1}<x^{H}_{J}=b,\quad
H_{j}=x_{j}^{H}-x_{j-1}^{H}.$

If $u_{H}^{n}\in H^{1}(0,T;\mathcal{S}_{H}^{0}(\Omega))$, then we can building
the fully implicit finite element scheme of the problem (3.31) such follows:

$\begin{cases}(\delta_{t}u_{H}^{n},v_{H})+\epsilon^{2}\mathcal{A}_{\delta}\left(u_{H}^{n},v_{H}\right)+(f(u_{H}^{n}),v_{H})=0,\,\,\forall
v_{H}\in\mathcal{S}_{H}^{0}(\Omega),\,\,0<t\leq T,\\\
(u_{H}^{0},v_{H})=(\psi(x),v_{H}),~{}\forall
v_{H}\in\mathcal{S}_{H}^{0}(\Omega).\end{cases}$ (3.34)

Step 2: We choose the fine meshes which the step size is
$h_{j}(j=1,2,\cdots,\mathscr{J})$. Then we have

$a=x^{h}_{0}<x^{h}_{1}<\cdots<x^{h}_{J-1}<x^{h}_{J}=b,\quad
h_{j}=x_{j}^{h}-x_{j-1}^{h}.$

In the fine meshes, the nonlinear term is calculated by $u_{H}^{n}$, then we
can building the semi-implicit finite element scheme as follows:

$\begin{cases}(\delta_{t}u_{h}^{n},v_{h})+\epsilon^{2}\mathcal{A}_{\delta}\left(u_{h}^{n},v_{h}\right)+(f(u_{H}^{n}),v_{h})\\\
\qquad\qquad+(f^{\prime}(u_{H}^{n})(u_{h}^{n}-u_{H}^{n}),v_{h})=0,~{}\forall
v_{h}\in\mathcal{S}_{h}^{0}(\Omega),\quad 0<t\leq T,\\\
(u_{h}^{0},v_{h})=(\psi(x),v_{h}),\quad\forall
v_{h}\in\mathcal{S}_{h}^{0}(\Omega).\end{cases}$ (3.35)

Finally, we can give the next matrix equations for completing (3.34)-(3.35)
which are:
Step 1: The matrix equation of (3.34) is

$(\bm{A}_{H}\bm{+}\tau\epsilon^{2}\bm{S}_{H})\bm{U}_{H}^{n}=\bm{A}_{H}\bm{U}_{H}^{n-1}-\tau\bm{F}^{n}_{H},$
(3.36)

where $\bm{A}_{H}$ is the usual (tridiagonal) FEM stiffness matrix which step
size is $H_{j}$ on uniform meshes, the stiffness matrix $\bm{S}_{H}$ is
defined in the above. $\bm{F}^{n}_{H}$ are defined as follows:

$\bm{F}_{H}^{n}=\int_{\Omega}f(u_{H}^{n})\bm{\Phi}dx.$

In this step, we use (3.37) to complete (3.36). It is

$(\bm{A}_{H}\bm{+}\tau\epsilon^{2}\bm{S}_{H})\bm{U}_{H}^{n,(k+1)}=\bm{A}_{H}\bm{U}_{H}^{n-1}-\tau\bm{F}_{H}^{n,(k)}.$
(3.37)

where $\bm{A}_{H}$ is the usual (tri-diagonal) FEM stiffness matrix $\bm{A}$
in coarse meshes, the stiffness matrix $\bm{S}$ is defined in the above. And
$\bm{F}_{H}^{n,(k)}$ are defined as follows:

$\bm{F}_{H}^{n,(k)}=\int_{I}f(u_{H}^{n,(k)})\bm{\Phi}dx,$

where $k$ is the iteration steps ($k=0,1,2,\cdots$).
Step 2: Based on the solutions $u_{H}^{n}$ in Step 1 and the basis functions
$\\{\varphi_{j}(x)\\}_{j=1}^{J-1}$ of the finite element method in the coarse
meshes, we have

$\displaystyle
u^{n}_{H}(x)=\sum_{j=1}^{J-1}u^{n}_{H,j}\varphi_{j}(x)=\frac{x^{H}_{j+1}-x}{x^{H}_{j+1}-x^{H}_{j}}u^{n}_{H,j}+\frac{x-x^{H}_{j}}{x^{H}_{j+1}-x^{H}_{j}}u^{n}_{H,j+1},$
$\displaystyle x\in(x^{H}_{j},x^{H}_{j+1}),~{}~{}j=0,1,2,\cdots,J-1.$ (3.38)

Step 3: Utilizing the solutions obtained from Step 2, the matrix equation of
(3.35) is

$(\bm{A}_{h}\bm{+}\tau\epsilon^{2}\bm{S}_{h}+\tau\bm{G}^{n})\bm{U}^{n}_{h}=(\bm{A}_{h}+\tau\bm{G}^{n})\bm{U}^{n-1}_{h}-\tau\bm{F}^{n}_{H},$
(3.39)

where $\bm{A}_{h}$ is the usual (tri-diagonal) FEM stiffness matrix $\bm{A}$
in the fine meshes, the stiffness matrix $\bm{S}_{h}$ is defined in the above.
$\bm{F}^{n}_{H}$ is defined as follows:

$\bm{F}_{H}^{n}=\int_{\Omega}f(u_{H}^{n})\bm{\Phi}\,{\rm d}x.$

Table 3.8. The errors and space rates of semi-implicit schemes when $\delta=0.01$. $\alpha$ | $J$ | $\mathscr{J}$ | $N=128$ | $N=1024$
---|---|---|---|---
Max error | Rate | CPU(s) | Max error | Rate | CPU(s)
| $2^{3}$ | $2^{6}$ | 6.17e-01 |  | 0.269 | 6.25e-01 |  | 2.083
| $2^{4}$ | $2^{8}$ | 3.04e-02 | 2.171 | 1.312 | 3.08e-02 | 2.171 | 10.42
0.25 | $2^{5}$ | $2^{10}$ | 1.22e-03 | 2.320 | 7.847 | 1.22e-03 | 2.328 | 69.42
| $2^{6}$ | $2^{12}$ | 5.96e-05 | 2.176 | 73.88 | 5.96e-05 | 2.179 | 664.88
| $2^{7}$ | $2^{14}$ | 3.73e-06 | 2.000 | 1363.3 | 3.73e-06 | 2.000 | 12720.3

##  4\. Concluding remarks and discussions

Different from the implementation of FEM in the physical space, we evaluated
the nonlocal stiffness matrix of piecewise linear FEM for the nonlocal
operator in the Fourier-transformed domain, and derived the explicit integral
form of the entries on non-uniform meshes. This allowed us to adopt the
different meshes for different situations, and enabled us to resort to the
Gauss quadrature for one dimension integral accurately. Plenty of numerical
examples for solving the nonlocal models are presented to demonstrate the
accuracy and efficiency of the proposed method.

The current study is largely based on a simple one-dimensional linear model
for the sake of offering insight without being impeded by tedious
calculations. More careful studies on the stability and convergence analysis
of piecewise linear FEM on non-uniform meshes, and one may extend the study to
two-dimensional rectangular elements, but it is more involved, which we shall
report in a future work.

Declarations

Conflict of interest The authors declare that they have no conflict of
interest.

##  Appendix A Evaluation of the distance matrix

For $j,k=1,2,\cdots,N-1$, tedious but straightforward calculations show that
$\bm{D}_{j}^{k}(s)$ in (2.14) are

  * (i)

For $j=k,$

$\begin{split}\bm{D}_{j}^{j}(s)=\begin{pmatrix}g(0)&g(h_{j})&g(h_{j}+h_{j+1})\\\\[2.0pt]
g(h_{j})&g(0)&g(h_{j})\\\\[2.0pt]
g(h_{j}+h_{j+1})&g(h_{j})&g(0)\end{pmatrix},\end{split}$ (A.1)

  * (ii)

For $j=k+1,$

$\begin{split}\bm{D}_{j}^{j-1}(s)=\begin{pmatrix}g(h_{j-1})&g(0)&g(h_{j})\\\\[2.0pt]
g(h_{j-1}+h_{j})&g(h_{j})&g(0)\\\\[2.0pt]
g(h_{j-1}+h_{j}+h_{j+1})&g(h_{j}+h_{j+1})&g(h_{j+1})\end{pmatrix},\end{split}$
(A.2)

  * (iii)

For $j=k-1,$

$\begin{split}\bm{D}_{j}^{j+1}(s)=\begin{pmatrix}g(h_{j})&g(h_{j}+h_{j+1})&g(h_{j}+h_{j+1}+h_{j+2})\\\\[2.0pt]
g(0)&g(h_{j+1})&g(h_{j+1}+h_{j+2})\\\\[2.0pt]
g(h_{j+1})&g(0)&g(h_{j+2})\end{pmatrix},\end{split}$ (A.3)

  * (iv)

For $j=k+2,$

$\begin{split}\bm{D}_{j}^{j-2}(s)=\begin{pmatrix}g(h_{j-2}+h_{j-1})&g(h_{j-1})&g(0)\\\\[2.0pt]
g(h_{j-2}+h_{j-1}+h_{j})&g(h_{j-1}+h_{j})&g(h_{j})\\\\[2.0pt]
g(h_{j-2}+h_{j-1}+h_{j}+h_{j+1})&g(h_{j-1}+h_{j}+h_{j+1})&g(h_{j}+h_{j+1})\end{pmatrix},\end{split}$
(A.4)

  * (v)

For $j=k-2,$

$\begin{split}\bm{D}_{j}^{j+2}(s)=\begin{pmatrix}g(h_{j}+h_{j+1})&g(h_{j}+h_{j+1}+h_{j+2})&g(h_{j}+h_{j+1}+h_{j+2}+h_{j+3})\\\\[2.0pt]
g(h_{j+1})&g(h_{j+1}+h_{j+2})&g(h_{j+1}+h_{j+2}+h_{j+3})\\\\[2.0pt]
g(0)&g(h_{j+2})&g(h_{j+2}+h_{j+3})\end{pmatrix},\end{split}$ (A.5)

  * (vi)

For $j=k+m,(m\geq 3)$

$\begin{split}\bm{D}_{j}^{j-m}(s)=\begin{pmatrix}g(\sum\limits_{i=1}^{m}h_{j-i})&g(\sum\limits_{i=1}^{m-1}h_{j-i})&g(\sum\limits_{i=1}^{m-2}h_{j-i})\\\\[2.0pt]
g(\sum\limits_{i=0}^{m}h_{j-i})&g(\sum\limits_{i=0}^{m-1}h_{j-i})&g(\sum\limits_{i=0}^{m-2}h_{j-i})\\\\[2.0pt]
g(\sum\limits_{i=-1}^{m}h_{j-i})&g(\sum\limits_{i=-1}^{m-1}h_{j-i})&g(\sum\limits_{i=-1}^{m-2}h_{j-i})\end{pmatrix},\end{split}$
(A.6)

  * (vii)

For $j=k-m,(m\geq 3)$

$\begin{split}\bm{D}_{j}^{j+m}(s)=\begin{pmatrix}g(\sum\limits_{i=0}^{m-1}h_{j+i})&g(\sum\limits_{i=0}^{m}h_{j+i})&g(\sum\limits_{i=0}^{m+1}h_{j+i})\\\\[2.0pt]
g(\sum\limits_{i=1}^{m-1}h_{j+i})&g(\sum\limits_{i=1}^{m}h_{j+i})&g(\sum\limits_{i=1}^{m+1}h_{j+i})\\\\[2.0pt]
g(\sum\limits_{i=2}^{m-1}h_{j+i})&g(\sum\limits_{i=2}^{m}h_{j+i})&g(\sum\limits_{i=2}^{m+1}h_{j+i})\end{pmatrix}.\end{split}$
(A.7)

## References

  * [1] M. Abramovitz and I. A. Stegun,  Handbook of Mathematical Functions,  Dover, New York, 1972.
  * [2] G. Acosta and J. P. Borthagaray,  A fractional Laplace equation: regularity of solutions and finite element approximations,  SIAM J. Numer. Anal., 55(2017)472–495.
  * [3] M. Ainsworth and C. Glusa,  Aspects of an adaptive finite element method for the fractional Laplacian: a priori and a posteriori error estimates, efficient implementation and multigrid solver,  Comput. Methods Appl. Mech. Engrg., 327(2017)4-35.
  * [4] B. Aksoylu and Z. Unlu,  Conditioning analysis of nonlocal integral operator in fractional Sobolev spaces,  SIAM J. Numer. Anal., 52(2014)653-677.
  * [5] E. Askari, F. Bobaru, R. B. Lehoucq, M. L. Parks, S. A. Silling and O. Weckner,  Peridynamics for multiscale materials modeling,  Journal of Physics: Conference Series, 2008.
  * [6] E. Aulisa, G. Capodaglio, A. Chierici, M. D’Elia, Marta,  Efficient quadrature rules for finite element discretizations of nonlocal equations,  Numer. Methods Partial Differential Equations, 38(2022)1767–1793.
  * [7] A. Bonito, J. P. Borthagaray, R. H. Nochetto, E. Otárola and A. J. Salgado,  Numerical methods for fractional diffusion,  Comput. Vis. Sci., 19 (2018)19-46.
  * [8] R. Cao, M. Chen, Y. Qi, J. Shi, X. Yin,  Analysis of (shifted) piecewise quadratic polynomial collocation for nonlocal diffusion model,  Appl. Numer. Math. 185(2023)120–140.
  * [9] X. Chen, M. Gunzburger,  Continuous and discontinuous finite element methods for a peridynamics model of mechanics,  Comput. Methods Appl. Mech. Engrg., 200(2011)1237-1250.
  * [10] H. Chen, C. Sheng, and L.-L Wang,  On explicit form of the FEM stiffness matrix for the integral fractional Laplacian on non-uniform meshes,  Appl. Math. Lett., 113(2021)106864.
  * [11] I. S. Gradshteyn and I. M. Ryzhik,  Table of Integrals, Series and Products,  Seven Edition, Elsevier, 2007.
  * [12] Q. Du, M. Gunzburger, R. Lehoucq and K. Zhou,  Analysis and approximation of nonlocal diffusion problems with volume constraints,  SIAM Review, 54(2012)667-696.
  * [13] Q. Du, M. Gunzburger, R. Lehoucq and K. Zhou,  A posteriori error analysis of finite element method for linear nonlocal diffusion and peridynamic models,  Math. Comp., 82(2013)1889-1922.
  * [14] Q. Du,  Nonlocal Modeling, Analysis, and Computation,  Society for Industrial and Applied Mathematics, Philadelphia, 2019.
  * [15] Q. Du, X. Tian and Z. Zhou  Nonlocal diffusion models with consistent local and fractional limits,  A³N²M: Approximation, Applications, and Analysis of Nonlocal, Nonlinear Models, The IMA Volumes in Mathematics and its Applications, Springer, Cham. 165(2023)175–213.
  * [16] S. Duo, H. van Wyk and Y. Zhang,  A novel and accurate finite difference method for the fractional Laplacian and the fractional poisson problem,  J. Comput. Phys., 355(2018)233-252.
  * [17] M. D’Elia, M. Gunzburger,  The fractional Laplacian operator on bounded domains as a special case of the nonlocal diffusion operator,  Comput. Math. Appl., 66(2013)1245– 1260.
  * [18] M. D’Elia, Q. Du, C. Glusa, M. Gunzburger, X. Tian, and Z. Zhou,  Numerical methods for nonlocal and fractional models,  Acta Numer., 29(2020)1-124.
  * [19] E. Di Nezza, G. Palatucci, and E. Valdinoci,  Hitchhiker’s guide to the fractional Sobolev spaces,  Bull. Sci. Math. 136(2012)521–573.
  * [20] I. Fried,  Condition of finite element matrices generated from nonuniform meshes,  AIAA J. 10(1972)219-221.
  * [21] Q. Guan, M. Gunzburger, X. Zhang,  Collocation method for one dimensional nonlocal diffusion equations,  Numer. Methods Partial Differential Equations, 38(2022)1618–1635.
  * [22] Z. Hao, Z. Zhang and R. Du,  Fractional centered difference scheme for high-dimensional integral fractional Laplacian,  J. Comput. Phys., 424(2021)109851.
  * [23] Y. Huang, A. Oberman,  Numerical methods for the fractional Laplacian: a finite difference-quadrature approach,  SIAM J. Numer. Anal., 52(2014)3056-3084.
  * [24] M. Klar, G. Capodaglio, M. D’Elia, C. Glusa, M. Gunzburger, Max, and C. Vollmann,  A scalable domain decomposition method for FEM discretizations of nonlocal equations of integrable and fractional type,  Comput. Math. Appl. 151(2023)434-448.
  * [25] Y. Leng, X. Tian, N. Trask, J. T. Foster,  Asymptotically compatible reproducing kernel collocation and meshfree integration for nonlocal diffusion,  SIAM J. Numer. Anal. 59(2021)88–118.
  * [26] H. Li, R. Liu and L.-L. Wang,  Efficient Hermite spectral-Galerkin methods for nonlocal diffusion equations in unbounded domains,  Numer. Math. Theo. Meth. Appl., 15 (2022)1009-1040.
  * [27] A. Lischke, G. Pang, M. Gulian, and et al.,  What is the fractional Laplacian? A comparative review with new results,  J. Comput. Phys., 404(2020)109009.
  * [28] Z. Liu, A. Cheng and H. Wang,  An hp-Galerkin method with fast solution for linear peridynamic models in one dimension,  Comp. Math. Appl., 73(2017)1546-1565.
  * [29] H. Liu, C. Sheng, L.-L. Wang, and H. Yuan,  On diagonal dominance of FEM stiffness matrix of fractional Laplacian and maximum principle preserving schemes for fractional Allen-Cahn equation,  J. Sci. Comput., 86(2021)19.
  * [30] J. Lu, Y. Nie,  A reduced-order fast reproducing kernel collocation method for nonlocal models with inhomogeneous volume constraints,  Comput. Math. Appl. 121(2022)52–61.
  * [31] Z. Mao, J. Shen,  Hermite spectral methods for fractional PDEs in unbounded domains,  SIAM J. Sci. Comput., 39(2017)A1928-A1950.
  * [32] H. Wang, H. Tian,  A fast and faithful collocation method with efficient matrix assembly for a two-dimensional nonlocal diffusion model,  Comput. Methods Appl. Mech. Engrg., 273(2014)19-36.
  * [33] C. Sheng, J. Shen, T. Tang, L. L. Wang and H. Yuan,  Fast Fourier-like mapped Chebyshev spectral-Galerkin methods for PDEs with integral fractional Laplacian in unbounded domains,  SIAM J. Numer. Anal., 58(2020)2435-2464.
  * [34] C. Sheng, L. L. Wang, H. Chen and H. Li,  Fast implementation of FEM for integral fractional Laplacian on rectangular meshes,  Commun. Comput. Phys., 35(2024).
  * [35] T. Tang, L.L. Wang, H. Yuan and T. Zhou,  Rational spectral methods for PDEs involving fractional Laplacian in unbounded domains,  SIAM J. Sci. Comput., 42 (2020)A585-A611.
  * [36] X. Tian and Q. Du,  Analysis and comparison of different approximations to nonlocal diffusion and linear peridynamic equations,  SIAM J. Numer. Anal., 51(2013)3458-3482.
  * [37] X. Tian and Q. Du,  Asymptotically compatible schemes for robust discretization of nonlocal models and their local limits,  SIAM J. Numer. Anal., 52(2014)1641-1665.
  * [38] X. Tian, Q. Du and M. Gunzburger,  Asymptotically compatible schemes for the approximation of fractional Laplacian and related nonlocal diffusion problems on bounded domains,  Adv. Comput. Math., 42(2016)1363–1380.
  * [39] L. Trefethen, Á. Birkisson, and T. A. Driscoll,  Exploring ODEs.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2018.
  * [40] Q. Ye and X. Tian,  Monotone meshfree methods for linear elliptic equations in non-divergence form via nonlocal relaxation,  J. Sci. Comput. 96 (2023)33.
  * [41] C. Zheng, J. Hu, Q. Du, and J. Zhang,  Numerical solution of the nonlocal diffusion equation on the real line,  SIAM J. Sci. Comput., 39(2017)A1951-A1968.
  * [42] K. Zhou and Q. Du,  Mathematical and numerical analysis of linear peridynamic models with nonlocal boundary conditions,  SIAM J. Numer. Anal., 48(2010)1759-1780.