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

Hongbin Chen¹, Changtao Sheng² and Li-Lian Wang³

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

¹College 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: hongbinchen@csuft.edu.cn (H. Chen).
²Corresponding 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: ctsheng@sufe.edu.cn (C. Sheng).
³Division 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: lilian@ntu.edu.sg (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

{comment}

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),$		(2.15)

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]).

{comment}

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}

{comment}

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].

{comment}

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)

{comment}

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)

{comment}

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)

{comment}

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)

{comment}

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. ∎

{comment}^†^†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)

{comment}

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)$ .

{comment}

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.

{comment}

•

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$ .

Refer to caption — 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]).

{comment}

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)$ .

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.

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).

{comment}

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.

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.

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$ .

{comment}

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
$-0.5$	$\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
$-0.75$	$\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
$0.5$	$\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
$1.5$	$\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.

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.

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).

{comment}

Table 3.6. The errors and time rates of semi-implicit schemes when

M=1024

$\alpha$	$N$	$\delta=0.01$		$\delta=0.1$
$\alpha$	$N$	$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$
$\alpha$	$M$	$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

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.

{comment}

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$
$\alpha$	$J$	$\mathscr{J}$	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.

{comment}

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.

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

Abstract.

Key words and phrases:

1. Introduction

2. Finite Element Method in Fourier Domain

2.1. FEM nonlocal stiffness matrix

Theorem 2.1.

Proof.

Theorem 2.2 (Case II: δ=∞\delta=\infty).

Theorem 2.3.

Proof.

2.2. FEM nonlocal stiffness matrix with δ≤h\delta\leq h

Corollary 2.1.

Proof.

Remark 2.1.

Theorem 2.4.

Proof.

Remark 2.2.

2.3. Nonlocal FEM stiffness matrix on uniform meshes

Theorem 2.5.

Proof.

Remark 2.3.

Corollary 2.2.

Proof.

Proposition 2.1.

Remark 2.4.

Remark 2.5.

3. Numerical Studies and Discussions

3.1. Conditioning of the stiffness matrix on nonuniform meshes

3.1.1. Diagonal dominance

Theorem 3.6.

Proof.

3.2. Error analysis

Lemma 3.1.

Lemma 3.2.

Lemma 3.3.

Proof.

Theorem 3.7.

Proof.

3.3. Nonlocal BVPs with various interaction kernels

3.3.1. Nonlocal BVP with λ=0\lambda=0

3.3.2. Nonlocal BVP with λ<0\lambda<0

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

3.4. Time-dependent problem

4. Concluding remarks and discussions

Appendix A Evaluation of the distance matrix

References

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

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

3.3.1. Nonlocal BVP with $\lambda=0$

3.3.2. Nonlocal BVP with $\lambda<0$

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