Computation of Eigenvalues for Nonlocal Models by Spectral Methods

Luciano Lopez and Sabrina Francesca Pellegrino Dipartimento di Matematica, Università degli Studi di Bari Aldo Moro, via E. Orabona 4, 70125 Bari, Italy [email protected] Dipartimento di Management, Finanza e Tecnologia, Università LUM Giuseppe Degennaro, S.S. 100 Km 18 - 70010 Casamassima (BA), Italy [email protected]

Abstract.

The purpose of this work is to study spectral methods to approximate the eigenvalues of nonlocal integral operators. Indeed, even if the spatial domain is an interval, it is very challenging to obtain closed analytical expressions for the eigenpairs of peridynamic operators. Our approach is based on the weak formulation of eigenvalue problem and in order to compute the eigenvalues we consider an orthogonal basis consisting of a set of Fourier trigonometric or Chebyshev polynomials. We show the order of convergence for eigenvalues and eigenfunctions in $L^{2}$ -norm, and finally, we perform some numerical simulations to compare the two proposed methods.

Keywords. Eigenvalues computation, Spectral methods, Fourier trigonometric polynomials, Chebyshev polynomials, Nonlocal peridynamics.

1. Introduction

Eigenvalues and eigenfunctions are used widely in science and engineering. They have many applications, particularly in physics, where they are usually related to vibrations. Indeed, eigenfunctions can provide the direction of spread of data, while eigenvalues represent the intensity of spread in the direction of the corresponding eigenfunction.

Eigenvalue problems occur naturally in the vibration analysis of mechanical structures with many degrees of freedom and in the context of non linear and nonlocal wave equations. The eigenvalues are used to determine the natural frequencies of vibration, and the relative eigenfunctions determine the shape of these vibrational modes. The orthogonality properties of the eigenfunctions allow a differential equation to be decoupled into a certain number of independent directions whose linear combination help to better understand the system under consideration.

In particular, eigenvalue problems can be achieved when someone looks for the solutions of non linear problems by using the method of separation of variables.

The eigenvalue problem of complex structures is often solved using finite element analysis or spectral techniques.

Moreover, they have applications in other areas of applied mathematics such as elastic theory, deformation, finance, quantum mechanics, image recognition and dimensional reduction.

In this work, we restrict our attention to the applications of the eigenvalue problems to nonlocal peridynamic framework.

Peridynamics is a nonlocal theory for elasticity and dynamic fracture analysis consisting in a second order in time partial integro-differential equation introduced by Silling in [30]. It is suitable to solve problems with spontaneous formation of discontinuities or singularities as it avoids the use of spatial partial derivatives. It is a nonlocal theory as the interactions between material points have a long-range: this means that a material particle interacts with other units in its neighborhood directly across finite distance called horizon (see [31, 27, 10, 5, 12, 26, 22, 18]).

The horizon represents the locality of the model, and one can prove that when the horizon goes to zero, the peridynamic equation approaches the wave equation (see [23, 3, 4, 8]). Additionally, the peridynamic operator has been used to approximate the fractional Laplacian operator (see [11]).

In peridynamic setting, the eigenvalue analysis predicts the theoretical buckling strength of a linear elastic structure and the eigenvalue problem is correlated to the static structural analysis of the system under consideration (see [16, 18, 17]). It can be useful to explore the deformation mechanism and the engineering properties of an elastic material. Moreover, it contributes to the determination of the preload in the nonlinear analysis can be utilized to investigate the effect of crack length, crack orientation and plate thickness. Indeed, the eigenfunctions represent the buckling mode that the material under consideration performs.

Several numerical methods are proposed to study the Laplacian eigenvalue problem. In particular the weak Galerkin finite element method and its accelerated version such as the two-grid and two-space methods are highly flexible and efficient methods used for the computation of the Laplacian eigenvalues, (see [35, 34]).

Our work is inspired by some studies of weak eigenvalues formulation for some local and nonlocal peridynamic operators (see [2, 1, 13]).

Applying the nonlinear mechanical properties in practice is a difficult issue, and, as a result, finite element method is frequently adopted for the spatial discretization. Additionally, for the simple geometric nonlinearities, analytical solution can be obtained.

For radially symmetric kernels, the evaluation of eigenvalues of linear nonlocal peridynamic operators under periodic boundary conditions can be used to obtain a spectral approximation of such operators. Following this idea, an hybrid spectral discretization obtained by combing truncated series expansion and high order Runge-Kutta ODEs solvers is proposed in [13].

In [1], the authors define local periodic and antiperiodic operators to enforce boundary conditions to the linear one dimensional peridynamic model, derive an explicit expression of their eigenvalues in terms of the horizon and study the conditioning and error analysis of such operators. For the numerical tests they use the Nyström method with the trapezoidal rule for the discretization.

An important issue in peridynamics is finding the solution of time dependent initial and boundary value problems. In [29, 33], the authors proposed coupled meshless finite point and finite element methods to discretize the peridynamic equation. Recently, spectral methods based on the discrete Fourier transform and convolution properties have been proposed to efficiently approximate the solution, (see [9, 25, 19, 20, 24, 21]). A different approach can be developed by using eigenvalues, (see [16]). Indeed, the computation of the eigenvalues of nonlocal operators can be used to solve periodic nonlocal time dependent problems. In particular, in [2], the authors proposed efficient and accurate spectral solvers based on the hypergeometric representation of the Fourier multipliers to evaluate the eigenvalues of the linear peridynamic integral operator.

In this paper, we focus on the weak formulation of a class of nonlocal peridynamic models and we study the spectrum properties of the involved operator. We show that the solution of the problem can be obtained as linear combination of a certain number of vibrational modes. We consider two different orthogonal basis to compute the eigenvalues and eigenfunctions of the peridynamic operator. The first one consists of Fourier trigonometric polynomials and is suitable for problems with periodic boundary conditions, while the second one does not require the periodicity condition and is given by the Chebyshev polynomials.

Let $\Omega$ be a one-dimensional bounded domain, then, the general form of the nonlocal peridynamic operator is the following

\mathcal{L}w(x)=\int_{\Omega\cap B_{\delta}(x)}f(x^{\prime}-x,w(x^{\prime})-w(x))dx^{\prime},\qquad x\in\Omega,

(1)

where $B_{\delta}(x)=\{x^{\prime}\in\Omega\,:\,|x-x^{\prime}|\leq\delta\}$ , $\delta>0$ is the horizon, $f$ is the pairwise force function, which describes the interaction between two material points, and the unknown $w$ represents the displacement field.

In what follows we restrict our attention to the case in which the evolution of the material body is given by a linear pairwise force function of convolution type in separable form

f(x^{\prime}-x,w(x^{\prime})-w(x))=C(x^{\prime}-x)\left(w(x^{\prime})-w(x)\right),

(2)

where the function $C\in L^{\infty}(\Omega)\cap C^{1}(\overline{\Omega})$ is a bounded positive even function, i.e. $C(x^{\prime}-x)=C(x-x^{\prime})$ , called micromodulus function, such that $C(x^{\prime}-x)\equiv 0$ , for all $x$ , $x^{\prime}\in\Omega$ such that $|x-x^{\prime}|>\delta$ . Examples of micromodulus functions for linear microelastic materials can be found in [32].

Thus the peridynamic operator becomes

\mathcal{L}w(x)=\int_{\Omega\cap B_{\delta}(x)}C(x^{\prime}-x)\left(w(x^{\prime})-w(x)\right)dx^{\prime}.

(3)

In the present paper we aim at studying the following eigenvalue problem

-\mathcal{L}w(x)=\lambda w(x),

(4)

with periodic boundary conditions on a one-dimensional bounded domain $\Omega$ , by means of spectral Fourier and Chebyshev approximation, where $\mathcal{L}$ is defined in (3).

The paper is organized as follows. Section 2 recalls the main theoretical results about the existence of eigenvalues for the variational formulation of the problem (4). Section 3 is devoted to the description of a numerical method for solving the eigenvalue problem by using the Fourier trigonometric and the Chebyshev polynomials. Moreover, this section contains some estimates about the distance between the continuous eigenpairs and their discrete approximation. Section 4 collects the numerical simulations and finally, Section 5 concludes the paper.

2. Main theoretical results

We recall the main theoretical results about the functional space in which we will work and the characterization of eigenvalues in terms of the min-max principle. To do so, we start by introducing the nonlocal integration by part rule.

Lemma 1.

The following equality holds

\int_{\Omega}\mathcal{L}w(x)v(x)dx=-\frac{1}{2}\int_{\Omega}\int_{\Omega}C(x^{\prime}-x)\left(w(x^{\prime})-w(x)\right)\left(v(x^{\prime})-v(x)\right)dx^{\prime}dx.

(5)

Proof.

Using the symmetry of $C$ and by relabeling the variable, we find

\int_{\Omega}\mathcal{L}w(x)v(x)dx=\int_{\Omega}\int_{\Omega}C(x^{\prime}-x)\left(w(x)-w(x^{\prime})\right)v(x^{\prime})dxdx^{\prime}

Moreover, by applying Fubini-Tonelli’s theorem there holds

\int_{\Omega}\mathcal{L}w(x)v(x)dx=-\int_{\Omega}\int_{\Omega}C(x^{\prime}-x)\left(w(x^{\prime})-w(x)\right)v(x^{\prime})dx^{\prime}dx.

Therefore,

	$\displaystyle\int_{\Omega}\mathcal{L}w(x)v(x)dx$	$\displaystyle=\frac{1}{2}\int_{\Omega}\mathcal{L}w(x)v(x)dx+\frac{1}{2}\int_{\Omega}\mathcal{L}w(x)v(x)dx$
		$\displaystyle=-\frac{1}{2}\int_{\Omega}\int_{\Omega}C(x^{\prime}-x)\left(w(x^{\prime})-w(x)\right)\left(v(x^{\prime})-v(x)\right)dx^{\prime}dx.$

∎

In what follows, we use $(\cdot,\cdot)$ and $\left\|\cdot\right\|$ to denote the inner product and the norm of the Hilbert space $L^{2}(\Omega)$ , respectively, namely

\left(w,v\right)=\int_{\Omega}w(x)v(x)dx,\quad\left\|w\right\|^{2}=\left(w,w\right),\qquad w,\,v\in L^{2}(\Omega).

Let $V\subset L^{2}(\Omega)$ be a suitable Banach space, then, using (5), we define the energetic extension of the peridynamic operator $\mathcal{L}:V\to V^{\prime}$ as follows

	$\displaystyle\left(\mathcal{L}w,v\right)$	$\displaystyle=\int_{\Omega}\mathcal{L}w(x)v(x)dx$		(6)
		$\displaystyle=-\frac{1}{2}\int_{\Omega}\int_{\Omega}C(x^{\prime}-x)\left(w(x^{\prime})-w(x)\right)\left(v(x^{\prime})-v(x)\right)dx^{\prime}dx,$

where $V^{\prime}$ is the dual space of $V$ .

The choice of the function space $V$ depends on the degree of singularity of the micromodulus function $C$ at $x^{\prime}-x=0$ . In case of strong singularity, one need to work with fractional Sobolev spaces, otherwise, we can work in the framework of $L^{p}$ spaces, for $p\geq 2$ (see [14]).

In particular, in what follows we consider the functional space

V=\{v\in L^{2}(\Omega):|v|_{V}^{2}<\infty\},

(7)

where $|\cdot|_{V}$ is the seminorm

|v|_{V}^{2}=\int_{\Omega}\int_{\Omega}C(x^{\prime}-x)|v(x^{\prime})-v(x)|^{2}dx^{\prime}dx,\qquad v\in L^{2}(\Omega).

(8)

We have that $V$ is an Hilbert space endowed with the norm $\left\|\cdot\right\|_{V}=\left\|\cdot\right\|_{L^{2}(\Omega)}+|\cdot|_{V}$ . Additionally, due to the assumptions on the micromodulus function, the $V$ -norm is equivalent to the $L^{2}(\Omega)$ -norm, in fact, for every $v\in V$ it holds

\left\|v\right\|_{L^{2}(\Omega)}\leq\left\|v\right\|_{V}\leq(1+M)\left\|v\right\|_{L^{2}(\Omega)},

for some positive constant $M$ .

As a consequence, the following theorem holds.

Theorem 1.

The embedding $V\hookrightarrow L^{2}(\Omega)$ is continuous.

In this context, the eigenvalue problem (4) has the following variational formulation: find $\lambda\in\mathbb{R}$ and $w\in V$ , $w\neq 0$ such that

\left(\mathcal{L}w,v\right)=\lambda(w,v),\qquad v\in V.

(9)

It is obtained by multiplying by a test function the equation (4) and by using the nonlocal integration by part rule (5).

Existence and uniqueness of weak solution $w\in V$ of the problem (4) is a consequence of the Lax-Milgram theorem (see [15]).

Remark 1.

In this paper we focus on the weak formulation (9) of (4), as it fits very well with the spectral methods we proposed. As a consequence, in what follows, unless different specified, when we talk about exact eigenvalues, we always refer to the weak ones. A classical eigenvalue is also a weak eigenvalue, but the reverse is not true in general. However, in our setting, using classical results of functional analysis for elliptic operators, one can prove the convergence of the weak problem (9) to the classical one (4).

A useful characterization of eigenvalues in terms of the Rayleigh quotient can be derived by using the so-called min-max principle for elliptic operators, which can be applied to self-adjoint operators in a Hilbert space. (see [15]).

Let $\lambda_{1}\leq\lambda_{2}\leq\cdots\leq\lambda_{k}\leq\cdots$ be the eigenvalues of $-\mathcal{L}$ defined in (3). The corresponding eigenfunctions $\{w_{k}\}_{k}$ form a complete orthogonal basis in $L^{2}(\Omega)$ .

Definition 1.

Let $w\in L^{2}(\Omega)$ , with $w\neq 0$ . Then, the Rayleigh quotient associated to the peridynamic operator $\mathcal{L}$ is defined as follows

R(w)=\frac{\left(\mathcal{L}w,w\right)}{\left\|w\right\|^{2}}.

(10)

We report here a well-known result on the Rayleigh quotient.

Theorem 2.

The $k$ -th eigenvalue $\lambda_{k}$ of $\mathcal{L}$ satisfies the following condition

\lambda_{k}=\min_{W_{k}}\max_{\stackrel{{\scriptstyle w\in W_{k}}}{{w\neq 0}}}R(w),

(11)

where $W_{k}$ is a $k$ -dimensional subspace.

3. Spectral approximation

In this section, we will consider two methods to approximate the eigenvalues of the nonlocal peridynamic operator. They deal with the weak formulation of the corresponding eigenvalue problem and use, as basis for the approximation, the Fourier trigonometric and the Chebyshev polynomials.

For each method, we will discuss its approximation properties and establish some convergence results for the nonlocal peridynamic problem.

As pointed out in [2], in the case of periodic domains, the classical eigenvalues of the peridynamic operator $\mathcal{L}$ are given by its Fourier multipliers defined through the Fourier transform

\lambda(\nu)=\int_{B_{\delta}(0)}C(y)\left(1-cos(\nu y)\right)dy.

(12)

In order to approximate the eigenvalues in (12) one can use a quadrature formula to compute numerically the integral, and the order of accuracy of the method provides the discretization error in the space variable. Instead, in this work, we compute an approximation of the eigenvalues of the peridynamic operator $\mathcal{L}$ by means of the weak formulation (9) of the problem (4). As we will show in the next section, our approach can benefit from the properties of the fast Fourier Transform (FFT), so that the dominant cost of the computation is those of a spectral solver.

More precisely, to obtain a numerical method for the problem (9), we replace the infinite dimensional space $V$ by a finite dimensional subspace $V_{N}$ for some integer $N>0$ . So, the problem becomes to find $\lambda_{N}\in\mathbb{R}$ and $w_{N}\in V_{N}$ , $w_{N}\neq 0$ such that

\left(\mathcal{L}w_{N},v_{N}\right)=\lambda_{N}(w_{N},v_{N}),\quad\text{for all $v_{N}\in V_{N}$}.

(13)

As we will see later, since the problem (13) is finite dimensional, it can be rewritten as a matrix eigenvalue problem. In particular, different approximations correspond to different choices of the subspace $V_{N}$ .

Let $\Omega=[a,b]$ be the computational domain. We fix $N>0$ an even integer number and let $\{\psi_{k}(x)\}_{|k|\leq N/2}$ be an orthogonal basis in the space $V\subset L^{2}([a,b])$ , so that $V_{N}=Span\{\psi_{-N/2}(x),\dots,\psi_{N/2}(x)\}$ .

Then, we recall that the variational formulation of the eigenvalue problem (4) is equivalent to finding a function $w$ , such that

\int_{a}^{b}\left(-\int_{\Omega\cap B_{\delta}(x)}C(x-x^{\prime})\left(w(x^{\prime})-w(x)\right)dx^{\prime}\right)\psi_{k}(x)\,dx=\int_{a}^{b}\lambda w(x)\psi_{k}(x)\,dx.

(14)

We seek $w_{N}\in V_{N}$ an approximation of $w(x)$ in form of linear combination of the basis $\{\psi_{k}(x)\}_{|k|\leq N/2}$ , that is

w_{N}(x)=\sum_{|h|\leq N/2}\tilde{w}_{h}\psi_{h}(x),

(15)

where $\tilde{w}_{h}$ are the coefficients of the linear combination.

Substituting (15) into (14), we get

		$\displaystyle\int_{a}^{b}\left(-\int_{\Omega\cap B_{\delta}(x)}C(x-x^{\prime})\left(\sum_{\|h\|\leq N/2}\tilde{w}_{h}(\psi_{h}(x^{\prime})-\psi_{h}(x))\right)dx^{\prime}\right)\psi_{k}(x)\,dx$		(16)
		$\displaystyle\quad=\int_{a}^{b}\lambda_{N}\sum_{\|h\|\leq N/2}\tilde{w}_{h}\psi_{h}(x)\psi_{k}(x)\,dx$

for all $k=-N/2,\cdots,N/2$ .

Notice that we can rewrite the left hand side of (16) as follows

	$\displaystyle\int_{a}^{b}\left(-\int_{\Omega\cap B_{\delta}(x)}C(x-x^{\prime})\left(\sum_{\|h\|\leq N/2}\tilde{w}_{h}(\psi_{h}(x^{\prime})-\psi_{h}(x))\right)dx^{\prime}\right)\psi_{k}(x)\,dx$
	$\displaystyle\quad=\sum_{\|h\|\leq N/2}\tilde{w}_{h}\left(-\int_{a}^{b}\psi_{k}(x)\int_{\Omega\cap B_{\delta}(x)}C(x-x^{\prime})(\psi_{h}(x^{\prime})-\psi_{h}(x))dx^{\prime}dx\right),$

while the right hand side of (16) becomes

\int_{a}^{b}\lambda_{N}\sum_{|h|\leq N/2}\tilde{w}_{h}\psi_{h}(x)\psi_{k}(x)\,dx=\lambda_{N}\sum_{|h|\leq N/2}\tilde{w}_{h}\left(\psi_{h}(x),\psi_{k}(x)\right).

Hence, in matrix notation, equation (16) gives the following eigenvalue problem

A\tilde{w}=\lambda_{N}M\tilde{w},

(17)

where $A=(a_{hk})_{h,k}$ and $M=(m_{hk})_{h,k}$ , with

	$\displaystyle a_{hk}$	$\displaystyle=-\int_{a}^{b}\psi_{k}(x)\int_{\Omega\cap B_{\delta}(x)}C(x-x^{\prime})(\psi_{h}(x^{\prime})-\psi_{h}(x))\,dx^{\prime}dx$
		$\displaystyle=-\int_{a}^{b}\psi_{k}(x)\int_{\Omega\cap B_{\delta}(x)}C(x-x^{\prime})\psi_{h}(x^{\prime})dx^{\prime}dx+\beta\int_{a}^{b}\psi_{k}(x)\psi_{h}(x)dx,$

m_{hk}=\int_{a}^{b}\psi_{h}(x)\psi_{k}(x)\,dx,

and $\beta=\int_{B_{\delta}(0)}C(x)dx$ .

Therefore, the approximate problem (13) is equivalent to the matrix eigenvalue problem (17).

We can observe that, due to the orthogonality of the basis $\{\psi_{k}(x)\}_{|k|\leq N/2}$ , the matrix $M$ is diagonal with non-zero diagonal entries.

Additionally, we can order the approximate eigenvalues of (13) as follows

\lambda_{N}^{1}\leq\lambda_{N}^{2}\leq\cdots\leq\lambda_{N}^{N}.

They fulfill an analogue characterization of the continuous eigenvalues:

\lambda_{N}^{k}=\min_{W_{k,N}\subset V_{N}}\quad\max_{\stackrel{{\scriptstyle w_{N}\in W_{k,N}}}{{w_{N}\neq 0}}}\frac{\left\|w_{N}\right\|_{W_{k,N}}^{2}}{\left\|w_{N}\right\|_{L^{2}(\Omega)}^{2}},

where $W_{k,N}$ is a $k$ -dimensional subspace of $V_{N}$ .

Since $V_{N}\subset V$ , as a consequence of Theorem 2, the following lemma holds.

Lemma 2.

Let $V_{N}\subset V$ . If $\lambda_{k}$ and $\lambda_{N}^{k}$ are the $k$ -th smallest eigenvalues of (9) and (13), respectively, then $\lambda_{k}\leq\lambda_{N}^{k}$ .

This lemma will be useful to prove the error estimate of the distance between the continuous and the discrete eigenvalues.

3.1. Convergence of the spectral method

This section is devoted to the proof of the following convergence result.

Theorem 3.

If $(\lambda_{N}^{k},w_{N}^{k})$ is the $k$ -th smallest eigenpair corresponding to the discrete problem (13). Then, there exists a positive constant $M=M(\Omega)$ independent of $N$ such that

\lambda_{k}\leq\lambda_{N}^{k}\leq\lambda_{k}\left(1+\frac{M}{N^{2}}\right),

(18)

where $\lambda_{k}$ is the $k$ -th smallest eigenvalue of the continuous problem (9).

In order to prove the theorem, we need to introduce some definitions and to recall some regularity results.

Let start by introducing the projection operator with respect to the $V$ -norm and the interpolant operator.

Definition 2.

The projection operator is the function $P_{N}:V\to V_{N}$ which associates to every $v\in V$ the only function $P_{N}v\in V_{N}$ satisfying the following orthogonality condition

\left(\mathcal{L}(v-P_{N}v),v_{N}\right)=0,\quad\text{for all $v_{N}\in V_{N}$},

(19)

or equivalently

\left\|v-P_{N}v\right\|_{V}=\inf_{v_{N}\in V_{N}}\left\|v-v_{N}\right\|_{V}

(20)

Definition 3.

Let $\{x_{j}\}_{j=1}^{N}$ be a partition of $\Omega$ . The interpolant operator is the function $I_{N}:V\to V_{N}$ which associates to every $v\in V$ the interpolant of degree $N$ defined as follows

I_{N}v=\sum_{j=1}^{N}v(x_{j})\psi_{j},

(21)

where $\{\psi_{j}\}_{j=1}^{N}$ is the orthogonal basis in $V$ we fixed for the approximation.

Remark 2.

Condition (20) ensures that

\left\|v-P_{N}v\right\|_{V}\leq\left\|v-I_{N}v\right\|_{V},\quad v\in V.

(22)

We now prove an interpolation estimate.

Lemma 3.

Let $I_{N}$ be the interpolant operator defined in (21). Then, there exists a positive constant $M=M(\Omega)$ independent of $N$ such that

\left\|v-I_{N}v\right\|_{V}\leq\frac{M}{N^{2}}\left\|v\right\|_{L^{2}(\Omega)}.

(23)

Proof.

We have

	$\displaystyle\left\\|v-I_{N}v\right\\|_{V}^{2}$	$\displaystyle\leq\int_{\Omega}\int_{\Omega}C(x^{\prime}-x)\left((v-I_{N}v)(x^{\prime})-(v-I_{N}v)(x)\right)^{2}dx^{\prime}dx$
		$\displaystyle\leq M\left\\|v-I_{N}v\right\\|_{L^{2}(\Omega)}^{2}\leq\frac{M}{N^{4}}\left\\|v\right\\|_{L^{2}(\Omega)}^{2},$

where the last inequality is a consequence of the spectral properties of the interpolant operator (see[7]). ∎

We recall the regularity result, (see [6]): given a function $f\in L^{2}(\Omega)$ , let $w\in V$ be the solution of the boundary-valued problem

\begin{cases}-\mathcal{L}w=f&\text{in $\Omega$}\\ w=0&\text{on $\partial\Omega$}.\end{cases}

Then, the following regularity estimate holds

\left\|w\right\|_{L^{2}(\Omega)}\leq M\left\|f\right\|_{L^{2}(\Omega)},

(24)

for some positive constant $M=M(\Omega)$ .

The next lemma studies the $L^{2}$ -convergence of the projection operator.

Lemma 4.

Let $v\in V$ be an eigenfunction of the problem (9). Then, there is a positive constant $M$ independent of $N$ such that

\left\|v-P_{N}v\right\|_{L^{2}(\Omega)}\leq\frac{M}{N^{2}}.

(25)

Proof.

Let $w\in V$ be the weak solution of the problem

-\mathcal{L}w=v-P_{N}v.

Thus, for every $\tilde{v}\in V$ , there holds

\left(\mathcal{L}w,\tilde{v}\right)=(v-P_{N}v,\tilde{v}).

(26)

If we choose $\tilde{v}=v-P_{N}v$ , since $I_{N}w$ and $v-P_{N}v$ are orthogonal, thanks to (26), we obtain

	$\displaystyle\left\\|v-P_{N}v\right\\|_{L^{2}(\Omega)}^{2}$	$\displaystyle=\left(\mathcal{L}w,v-P_{N}v\right)=\left(\mathcal{L}(w-I_{N}w),v-P_{N}v\right)$		(27)
		$\displaystyle\leq\left\\|w-I_{N}w\right\\|_{V}\left\\|v-P_{N}v\right\\|_{V}.$

The interpolant estimate (23) and the regularity estimate (24) give us

\left\|w-I_{N}w\right\|_{V}\leq\frac{M}{N^{2}}\left\|w\right\|_{L^{2}(\Omega)}\leq\frac{M}{N^{2}}\left\|v-P_{N}v\right\|_{L^{2}(\Omega)}.

(28)

Finally, using (25), we get the claim. ∎

We are now ready to prove Theorem 3.

Theorem 3.

Let $w_{k}$ be the eigenfunction of the continuous problem (9) corresponding to the $k$ -th smallest eigenvalue $\lambda_{k}$ . We consider the $k$ -dimensional subspaces

E_{k}=Span\{w_{1},\cdots,w_{k}\},\quad\tilde{E}_{k}=P_{N}E_{k}.

If $v\in E_{k}$ , with $v\neq 0$ , then

\left\|P_{N}v\right\|\geq\left\|v\right\|-\left\|v-P_{N}v\right\|\geq\left(1-\frac{M}{N^{2}}\right)\left\|v\right\|>0.

(29)

Therefore, $P_{N}$ is an isomorphism, and as a consequence the spaces $E_{k}$ and $\tilde{E}_{k}$ have the same dimension.

Thanks to Theorem 2, we find

	$\displaystyle\lambda_{k}\leq\lambda_{N}^{k}$	$\displaystyle\leq\max_{\stackrel{{\scriptstyle v_{N}\in\tilde{E}_{k}}}{{v_{N}\neq 0}}}\frac{\left(\mathcal{L}v_{N},v_{N}\right)}{\left\\|v_{N}\right\\|^{2}}=\max_{\stackrel{{\scriptstyle v\in E_{k}}}{{v\neq 0}}}\frac{\left(\mathcal{L}(P_{N}v),P_{N}v\right)}{\left\\|P_{N}v\right\\|^{2}}$
		$\displaystyle\leq\max_{\stackrel{{\scriptstyle v\in E_{k}}}{{v\neq 0}}}\frac{\left(\mathcal{L}v,v\right)}{\left\\|P_{N}v\right\\|^{2}}\leq\underbrace{\max_{\stackrel{{\scriptstyle v\in E_{k}}}{{v\neq 0}}}\frac{\left(\mathcal{L}v,v\right)}{\left\\|v\right\\|^{2}}}_{{}=\lambda_{k}}\max_{\stackrel{{\scriptstyle v\in E_{k}}}{{v\neq 0}}}\frac{\left\\|v\right\\|^{2}}{\left\\|P_{N}v\right\\|^{2}}.$

Finally, using (29), we get

\lambda_{k}\leq\lambda_{N}^{k}\leq\lambda_{k}\left(1+\frac{M}{N^{2}}\right),

and this concludes the proof. ∎

We can also derive a convergence result for the eigenfunctions.

Theorem 4.

Let $(\lambda_{k},w_{k})$ and $(\lambda_{N}^{k},w_{N}^{k})$ be the $k$ -th smallest eigenpair corresponding to the continuous problem (9) and to the discrete problem (13), respectively. Then, there exists a positive constant $M=M(\Omega)$ independent of $N$ such that

\left\|w_{k}-w_{N}^{k}\right\|_{L^{2}(\Omega)}\leq\frac{M}{N^{2}}.

(30)

Proof.

For simplicity we assume that $\lambda_{k}$ and $\lambda_{N}^{k}$ have algebraic multiplicity equals to one. The proof of the general case is similar.

We define

v_{N}^{k}=\left(P_{N}w_{k},w_{N}^{k}\right)w_{N}^{k},

(31)

and

\theta_{N}^{k}=\max_{i\neq k}\frac{\lambda_{k}}{|\lambda_{k}-\lambda_{N}^{i}|}.

(32)

We have

	$\displaystyle\left\\|w_{k}-w_{N}^{k}\right\\|_{L^{2}(\Omega)}$	$\displaystyle\leq\left\\|w_{k}-P_{N}w_{k}\right\\|_{L^{2}(\Omega)}+\left\\|P_{N}w_{k}-v_{N}^{k}\right\\|_{L^{2}(\Omega)}$		(33)
		$\displaystyle+\left\\|v_{N}^{k}-w_{N}^{k}\right\\|_{L^{2}(\Omega)}.$

Let us estimate the term $\left\|v_{N}^{k}-w_{N}^{k}\right\|_{L^{2}(\Omega)}$ . As long as we normalize the discrete and the continuous eigenfunctions, by (31) we have

\left\|v_{N}^{k}\right\|_{L^{2}(\Omega)}=|\left(P_{N}w_{k},w_{N}^{k}\right)|,

and

v_{N}^{k}-w_{N}^{k}=\left(1-\left(P_{N}w_{k},w_{N}^{k}\right)\right)w_{N}^{k},

so that

\left\|v_{N}^{k}-w_{N}^{k}\right\|_{L^{2}(\Omega)}=|1-\left(P_{N}w_{k},w_{N}^{k}\right)|.

If we choose the sign of $w_{N}^{k}$ such that the inner product $\left(P_{N}w_{k},w_{N}^{k}\right)\geq 0$ , then we obtain

	$\displaystyle\left\\|v_{N}^{k}-w_{N}^{k}\right\\|_{L^{2}(\Omega)}$	$\displaystyle=\|1-\left(P_{N}w_{k},w_{N}^{k}\right)\|$
		$\displaystyle=\bigl{\|}1-\|\left(P_{N}w_{k},w_{N}^{k}\right)\|\bigr{\|}$
		$\displaystyle=\bigl{\|}\left\\|w_{k}\right\\|_{L^{2}(\Omega)}-\left\\|v_{N}^{k}\right\\|_{L^{2}(\Omega)}\bigr{\|}$
		$\displaystyle\leq\left\\|w_{k}-v_{N}^{k}\right\\|_{L^{2}(\Omega)}.$

Therefore,

	$\displaystyle\left\\|v_{N}^{k}-w_{N}^{k}\right\\|_{L^{2}(\Omega)}$	$\displaystyle\leq\left\\|w_{k}-v_{N}^{k}\right\\|_{L^{2}(\Omega)}$		(34)
		$\displaystyle\leq\left\\|w_{k}-P_{N}w_{k}\right\\|_{L^{2}(\Omega)}+\left\\|P_{N}w_{k}-v_{N}^{k}\right\\|_{L^{2}(\Omega)}.$

Thanks to (34), equation (33) becomes

\left\|w_{k}-w_{N}^{k}\right\|_{L^{2}(\Omega)}\leq 2\left\|w_{k}-P_{N}w_{k}\right\|_{L^{2}(\Omega)}+2\left\|P_{N}w_{k}-v_{N}^{k}\right\|_{L^{2}(\Omega)}.

(35)

Lemma 4 implies that there exists a positive constant $M=M(\Omega)$ independent of $N$ such that

\left\|w_{k}-P_{N}w_{k}\right\|_{L^{2}(\Omega)}\leq\frac{M}{N^{2}}.

(36)

Finally, let $i\neq k$ . Since

	$\displaystyle\left(P_{N}w_{k},w_{N}^{i}\right)$	$\displaystyle=\frac{1}{\lambda_{N}^{i}}\left(\mathcal{L}(P_{N}w_{k}),w_{N}^{i}\right)=\frac{1}{\lambda_{N}^{i}}\left(\mathcal{L}w_{k},w_{N}^{i}\right)$
		$\displaystyle=\frac{\lambda_{k}}{\lambda_{N}^{i}}\left(w_{k},w_{N}^{i}\right),$

we find

	$\displaystyle\left(w_{k}-P_{N}v_{k},v_{N}^{i}\right)$	$\displaystyle=\left(w_{k},w_{N}^{i}\right)-\left(P_{N}w_{k},w_{N}^{i}\right)$
		$\displaystyle=\left(\frac{\lambda_{N}^{i}-\lambda_{k}}{\lambda_{k}}\right)\left(P_{N}w_{k},w_{N}^{i}\right).$

So, using (32)

\bigl{|}\left(P_{N}w_{k},w_{N}^{i}\right)\bigr{|}\leq\theta_{N}^{k}\bigl{|}\left(w_{k}-P_{N}v_{k},v_{N}^{i}\right)\bigr{|}.

As a consequence,

$\displaystyle\left\\|P_{n}w_{k}-v_{N}^{k}\right\\|_{L^{2}(\Omega)}^{2}$	$\displaystyle=\sum_{i\neq k}\left(P_{N}w_{k},w_{N}^{i}\right)^{2}$	(37)
	$\displaystyle\leq(\theta_{N}^{k})^{2}\sum_{i\neq k}\left(w_{k}-P_{N}w_{k},w_{N}^{i}\right)^{2}$
	$\displaystyle\leq(\theta_{N}^{k})^{2}\left\\|w_{k}-P_{N}w_{k}\right\\|_{L^{2}(\Omega)}^{2}\leq\frac{M}{N^{4}}.$

Collecting (36) and (37) into (35), we get the claim. ∎

In the following sections, we study the problem by considering the method of Fourier trigonometric and the Chebyshev polynomials as basis. The first choice is suitable in presence of periodic boundary condition, while the second one can be applied to non-periodic problem.

Remark 3.

This work focuses on the one-dimensional peridynamic model, however, the variational approach can be easily extended to higher dimensional problems, by using other functional spaces and test functions. For instance, for the bi-dimensional case, we can implement the 2-D fast Fourier Transform as done in [24].

3.2. Spectral approximation by Fourier trigonometric polynomials

Let $\Omega=[0,2\pi]$ be the computational domain. Assume that the problem is subject to periodic boundary conditions $w(0)=w(2\pi)$ .

Let $N>0$ be an even number and $\{\psi_{k}(x)\}_{|k|\leq N/2}$ be an orthogonal basis of Fourier trigonometric polynomials in $L^{2}([0,2\pi])$ such that

\psi_{k}(x)=e^{ikx},\quad x\in\Omega,

(38)

with $\psi_{k}(0)=\psi_{k}(2\pi)$ , for all $k=-N/2,\cdots,N/2$ .

We seek an approximation of $w(x)$ in form of real-valued Fourier trigonometric polynomials

w_{N}(x)=\sum_{|h|\leq N/2}\tilde{w}_{h}\psi_{h}(x).

(39)

Then, using the same computation as before, namely, substituting (39) and (38) into (14), we obtain the eigenvalue problem defined in (17), where matrix values of $A$ and $M$ have the following expression

	$\displaystyle a_{hk}$	$\displaystyle=-\int_{0}^{2\pi}\psi_{k}(x)\int_{\Omega\cap B_{\delta}(x)}C(x-x^{\prime})(\psi_{h}(x^{\prime})-\psi_{h}(x))\,dx^{\prime}dx$
		$\displaystyle=-\int_{0}^{2\pi}\psi_{k}(x)\psi_{h}(x)dx\,\int_{\Omega\cap B_{\delta}(0)}C(x^{\prime})\cos(hx^{\prime})dx^{\prime}+2\pi\beta\delta_{h,-k}$
		$\displaystyle=-2\pi\delta_{h,-k}\int_{\Omega\cap B_{\delta}(0)}C(x^{\prime})\cos(hx^{\prime})dx^{\prime}+2\pi\beta\delta_{h,-k},$

m_{hk}=\int_{0}^{2\pi}\psi_{h}(x)\psi_{k}(x)\,dx=2\pi\delta_{h,-k},

and $\delta_{h,-k}$ is the Kronecker delta.

More in details, the matrices $A$ and $M$ have the following representation

A=\begin{bmatrix}0&&*\\ &\iddots&\\ *&&0\end{bmatrix},\qquad M=\begin{bmatrix}0&&*\\ &\iddots&\\ *&&0\end{bmatrix}.

We can observe that, in this case the matrices $A$ and $M$ are anti-diagonal and the product $M^{-1}A$ is a diagonal definite-positive matrix. Moreover, the integral $\gamma_{h}=\int_{\Omega\cap B_{\delta}(0)}C(x^{\prime})\cos(hx^{\prime})dx^{\prime}$ which appear in the definition of $a_{hk}$ can be accurately computed by means of the Fast Fourier transform (FFT) as $C$ is an even real-valued function.

Additionally, when $h=0$ , such integral is equal to $\beta$ , therefore the matrix $A$ has a zero value on its diagonal, so $\lambda=0$ is an eigenvalue for the peridynamic operator.

Since the matrices $A$ and $M$ are anti-diagonal in this setting, and since the product of two anti-diagonal matrices gives a diagonal matrix, a complete derivation of the spectrum of the peridynamic operator is possible in this case. Indeed we have

\lambda_{h}=\beta-\gamma_{h}^{N},\quad h=-N/2,\cdots,N/2,

where $\gamma_{h}^{N}$ is the discrete Fourier transform of $\gamma_{h}$ and is less than $\beta$ . They corresponds to the diagonal values of the matrix $M^{-1}A$ . Moreover, the corresponding eigenfunctions are given by

w_{N}^{h}(x)=\psi_{h}(x)=e^{ihx},\qquad x\in\Omega,\quad h=-N/2,\dots,N/2.

Because the discrete sine and cosine functions of different frequencies are mutually orthogonal, we find $N$ linearly independent eigenfunctions. Note that, except for $\lambda_{0}$ , the eigenvalues have multiplicity higher than one because of the symmetry between $h$ and $-h$ .

We now compare the discrete eigenvalues with those of the continuous problem, found in (12). The following theorem establishes the estimate for the distance between the continuous eigenvalue and its Fourier trigonometric approximation. This is a trivial consequence of the estimates of the distance between a function and its Fourier approximation.

Theorem 5 (see [7]).

Let $\lambda_{k}$ and $\lambda_{N}^{k}$ be the $k$ -th smallest eigenvalues of the continuous and discrete eigenvalue problem respectively. If the micromodulus function $C\in L^{2}([0,2\pi])\cap\mathcal{C}^{1}([0,2\pi])$ , then there exists a positive constant $M$ independent of $N$ such that

|\lambda_{k}-\lambda_{N}^{k}|\leq\frac{M}{N^{2}}\left\|C^{\prime}\right\|_{L^{2}([0,2\pi])}.

(40)

3.3. Spectral approximation by Chebyshev polynomials

We consider the shifted Chebyshev polynomials $T_{k}(x)$ in the computational domain $\Omega=[0,2\pi]$ , for $k=0,\cdots,N$ , defined as

T_{k}(x)=\cos\left(k\arccos\left(\frac{x-\pi}{\pi}\right)\right).

They form an orthogonal basis in the weighted space $L^{2}_{\omega}([0,2\pi])$ , where the weighted function is $\omega(x)=\left(\sqrt{1-\left(\frac{x-\pi}{\pi}\right)^{2}}\right)^{-1}$ .

Indeed

\int_{0}^{2\pi}T_{h}(x)T_{k}(x)\omega(x)dx=\frac{\pi^{2}}{2}c_{h}\delta_{h,k},

where

c_{h}=\begin{cases}2,&h=0,\\ 1,&\text{otherwise}.\end{cases}

In this case, the equation (14) becomes

	$\displaystyle\int_{0}^{2\pi}$	$\displaystyle\left(-\int_{B_{\delta}(x)}C(x-x^{\prime})\left(w(x^{\prime})-w(x)\right)dx^{\prime}\right)T_{k}(x)\omega(x)dx$		(41)
		$\displaystyle\qquad=\int_{0}^{2\pi}\lambda w(x)T_{k}(x)\omega(x)dx,$

and we approximate $w(x)$ as linear combination of the Chebyshev polynomials:

w_{N}(x)=\sum_{h=0}^{N}\hat{w}_{h}T_{h}(x).

(42)

If we substitute equation (42) into (41), we find again the eigenvalue problem

A\hat{w}=\lambda M\hat{w},

with matrices $A=(a_{hk})_{h,k}$ and $M=(m_{hk})_{h,k}$ defined as follows

a_{hk}=-\int_{0}^{2\pi}\omega(x)T_{k}(x)\int_{\Omega\cap B_{\delta}(x)}C(x-x^{\prime})T_{h}(x^{\prime})dx^{\prime}dx+\frac{\pi^{2}}{2}\beta c_{h}\delta_{h,k},

m_{hk}=\int_{0}^{2\pi}T_{h}(x)T_{k}(x)\omega(x)\,dx=\frac{\pi^{2}}{2}c_{h}\delta_{h,k}.

If we denote by

\gamma_{h}(x)=\int_{\Omega\cap B_{\delta}(x)}C(x-x^{\prime})T_{h}(x^{\prime})dx^{\prime},

(43)

we can observe that the integral $\int_{0}^{2\pi}\omega(x)\gamma_{h}(x)T_{k}(x)dx$ can be efficiently computed by using the Chebyshev trasform.

Additionally, using the Chebyshev polynomials for the approximation of the peridynamic operator, we have that the matrix $A$ is nondiagonal in general.

Remark 4.

The advantage of using the Chebyshev polynomials relies in the fact that they can be implemented to solve problem without imposing periodic boundary conditions. However, in order to make a comparison between this method and the one based on the Fourier trigonometric polynomials, we can apply the Chebyshev polynomials to periodic problems. To do this, we need to require that $T_{k}(0)=T_{k}(2\pi)$ , which implies that $k$ must be even.

Thus, we have to restrict our attention to the following orthogonal sub-basis $\{T_{k}(x)\}_{k}$ , for $k=0,\dots,N$ , $k$ even.

4. Numerical simulations

In this section we validate and compare the methods proposed in Section 3 and collect some simulations in order to investigate the properties of the solutions of the eigenvalue problem (4).

To validate the results, we consider the following one-dimensional benchmark problem. We consider a bar in the spatial domain $[0,2\pi]$ and we take the micromodulus function $C(x)=4e^{-x^{2}}/\sqrt{\pi}$ , so in this case $\beta\approx 4$ .

When we deal with the Fourier trigonometric polynomials, we discretize the computational domain by using the uniform mesh $x_{j}=2\pi j/N$ for $j=0,\cdots,N$ ; while, when we consider the Chebyshev polynomials, we take the non-uniform mesh given by the Guass-Chebyshev-Lobatto nodes $x_{j}=\pi+\pi\left(-\cos(\pi j/N)\right)$ , for $j=0,\cdots,N$ . This choice of nodes is suitable as it avoids the Gibb’s phenomenon at the boundaries.

Figure 1 shows a comparison between the exact eigenvalues of the peridynamic operator with the discrete ones computed by means of the Fourier trigonometric and Chebyshev polynomials for various choice of the horizon $\delta$ . The convergence of the approximated eigenvalues is evaluated by computing the relative error in the discrete $L^{2}(\Omega)$ norm:

E_{L^{2}}=\frac{\sum_{k}\bigl{|}\lambda_{N}^{k}-\lambda_{k}\bigr{|}^{2}}{\sum_{k}\bigl{|}\lambda_{N}^{k}\bigr{|}^{2}},

where $\lambda_{k}$ represents the $k$ -th smallest eigenvalue of the continuous problem (9), and $\lambda_{N}^{k}$ denotes $k$ -th smallest eigenvalue of the discrete problem (13) computed using both the Fourier trigonometric and the Chebyshev polynomials.

Table 1 shows the relative error and the convergence rate for the two methods for different values of $N$ and for $\delta=3$ . We can observe that the order of convergence is in accordance with the theoretical results.

Trigonometric approximation			Chebyshev approximation
$N$	$E_{L^{2}}$	convergence rate	$N$	$E_{L^{2}}$	convergence rate
$20$	$4.8016\times 10^{-3}$	$-$	$20$	$8.8896\times 10^{-3}$	$-$
$40$	$1.1021\times 10^{-3}$	$2.0465$	$40$	$8.6192\times 10^{-4}$	$2.0377$
$80$	$2.5438\times 10^{-4}$	$2.0616$	$80$	$2.9249\times 10^{-4}$	$2.0485$
$160$	$5.5865\times 10^{-5}$	$2.0943$	$160$	$5.1822\times 10^{-5}$	$2.0867$
$320$	$1.1000\times 10^{-5}$	$2.1481$	$320$	$1.1802\times 10^{-5}$	$2.1267$

Table 1. The relative

L^{2}

-error as function of

N

corresponding to the trigonometric and the Chebyshev approximation.

Table 2 depicts the eigenvalues computed by the two methods for $N=20$ and $\delta=3$ . We can notice that they are all positive and bounded by the value of $\beta$ .

Moreover, we find that using the Fourier trigonometric polynomials as basis for the approximation, except for $\lambda_{0}\approx 0$ , the algebraic multiplicity of the other eigenvalues is equal to two. This is correlated with the fact that trigonometric polynomials are characterized by two different frequencies.

$\lambda_{n}^{Trig}$ , $n=0,\cdots,20$	$\lambda_{n}^{Cheby}$ , $n=0,\cdots,20$
0.0000	0.1905
0.8846	0.7084
0.8846	1.4173
2.5285	2.1567
2.5285	2.8007
3.5782	3.2864
3.5787	3.6105
3.9267	3.8042
3.9267	3.9090
3.9921	3.9607
3.9921	3.9841
3.9995	3.9940
3.9995	3.9978
3.9999	3.9993
3.9999	3.9997
3.9999	3.9999
3.9999	3.9999
3.9999	3.9999
3.9999	3.9999
3.9999	3.9999
3.9999	3.9999

Table 2. The peridynamic eigenvalues computed by using trigonometric and Chebyshev polynomials

Figure 2 shows the first five eigenfunctions obtained by using the trigonometric polynomials for the approximation of the peridynamic operator.

In addition, we make a simulation considering the peridynamic eigenvalue problem subject to antiperiodic boundary condition on the spatial domain $\Omega=[-1,1]$ :

\begin{cases}-\mathcal{L}w(x)=\lambda w(x),\qquad x\in[-1,1],\\ w(-1)=-1,\,w(1)=1.\end{cases}

(44)

The well-posedness of (44) is achieved in [36, 28].For the simulation we use the same micromodulus function as before and we compute the eigenvalues by using the Chebyshev polynomials. In Figure 3, we display a comparison between the exact eigenvalues and the discrete ones for different values of $N$ and $\delta$ . Here, the exact eigenvalues are computed by using the proposed numerical method by using $N=200$ .

4.1. The solution of the peridynamic problem

Let $[0,T]$ , for some $T>0$ , be the time domain under investigation. Consider a microelastic bar with constant mass density $\rho>0$ occupying a region $\Omega\subset{\mathbb{R}}$ . Then, the peridynamic model describes the dynamics of the body and its equation is given by

\rho\partial_{tt}^{2}{u}(x,t)=\mathcal{L}u(x,t)=\int_{\Omega\cap B_{\delta}(x)}C(x^{\prime}-x)\left(u(x^{\prime},t)-u(x,t)\right)dx^{\prime},

(45)

for $x\in\Omega$ , and $t\in[0,T]$ , with initial conditions

u(x,0)=u_{0}(x),\quad\partial_{t}u(x,0)=v_{0}(x),\qquad x\in\Omega,

(46)

where $u$ is the displacement field.

We assume that the bar is subjected to the uniform initial displacement $u_{0}(x)=e^{-x^{2}}$ , $v_{0}(x)=0$ and we fix $\delta=3$ as the size of the horizon.

Moreover, we assume that the constant density of the body is $\rho(x)=1$ and the interaction between two material particle is given by the micromodulus function $C$ defined as in the previous section.

If the spectrum of the nonlocal operator $-\mathcal{L}$ is known, one can formally express the solution of (45)-(46) as

u(x,t)=\sum_{k=0}^{\infty}w_{k}(x)\left(a_{k}\cos(\sqrt{\lambda_{k}}t)+b_{k}\sin(\sqrt{\lambda_{k}}t)\right),\qquad x\in\Omega,\quad t\in[0,T],

(47)

where $\lambda_{k}$ and $w_{k}(x)$ represent respectively the $k$ -th eigenvalue and eigenfunction of the peridynamic operator $-\mathcal{L}$ on the domain $\Omega$ , and the coefficients $a_{k}$ and $b_{k}$ are computed from the initial conditions (46) by

	$\displaystyle a_{k}$	$\displaystyle=\frac{\int_{\Omega}u_{0}(x)w_{k}(x)dx}{\int_{\Omega}w_{k}^{2}(x)dx},$		(48)
	$\displaystyle b_{k}$	$\displaystyle=\frac{1}{\sqrt{\lambda_{k}}}\frac{\int_{\Omega}v_{0}(x)w_{k}(x)dx}{\int_{\Omega}w_{k}^{2}(x)dx}.$

Thus, the problem of finding the solution of (45)-(46) reduces to study an eigenvalue problem.

In Figure 4 we make a comparison between the exact solution of the peridynamic problem (45)-(46) with the numerical solution obtained by using both the Fourier trigonometric polynomials and the Chebyshev polynomials as basis to solve the eigenvalue problem (17). We can observe in both cases a good agreement as the total number of mesh-points $N$ increases.

5. Conclusions

In this paper we propose two spectral methods to numerically compute the eigenvalues of the one-dimensional nonlocal linear peridynamic operator. They deal with the weak formulation of peridynamic model and make use of both the Fourier trigonometric and the Chebyshev polynomials. The first ones are suitable in presence of periodic boundary conditions, while the others can be applied even to non periodic problems. We prove the convergence of the discrete eigenpairs to the continuous ones in the $L^{2}$ -norm and our simulations confirm the theoretical results.

Acknowledgements

This paper has been partially supported by GNCS of Istituto Nazionale di Alta Matematica, by PRIN 2017 “Discontinuous dynamical systems: theory, numerics and applications” coordinated by Nicola Guglielmi and by Regione Puglia, “Programma POR Puglia 2014/2020-Asse X-Azione 10.4 Research for Innovation-REFIN - (D1AB726C)”. The authors thank the anonymous referees for critical remarks, suggestions and for the careful reading of the manuscript.

References

[1] B. Aksoylu and A. Kaya. Conditioning and error analysis of nonlocal operators with local boundary conditions. Journal of Computational and Applied Mathematics, 335:1–19, 2018.
[2] B. Alali and N. Albin. Fourier Spectral Methods for Nonlocal Models. Journal of Peridynamics and Nonlocal Modeling, 2:317–335, 2020.
[3] R. Alebrahim. Peridynamic modeling of Lamb wave propagation in bimaterial plates. Composite Structures, 214:12 – 22, 2019.
[4] R. Alebrahim, P. Packo, M. Zaccariotto, and U. Galvanetto. Wave propagation improvement in two-dimensional bond-based peridynamics model. In Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 2021.
[5] F. Bobaru, M. Yang, S. Alves, F.and Silling, E. Askari, and J. Xu. Convergence, adaptive refinement, and slaning in 1D peridynamics. Int. J. Numer. Mech. Eng., 77:852 – 877, 2009.
[6] L. Caffarelli and L. Silvestre. The Evans-Krylov theorem for nonlocal fully nonlinear equations. Annals of Mathematics, 174:1163 – 1187, 2011.
[7] C. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang. Spectral Methods: Fundamentals in Single Domains. Springer-Verlag Berlin Heidelberg, 2006.
[8] J. Chen, Y. Tian, and X. Cui. Free and forced vibration analysis of peridynamic finite bar. International Journal of Applied Mechanics, 10(1):1850003, 2018.
[9] G. M. Coclite, A. Fanizzi, L. Lopez, F. Maddalena, and S. F. Pellegrino. Numerical methods for the nonlocal wave equation of the peridynamics. Applied Numerical Mathematics, 155:119 – 139, 2020.
[10] M. D’Elia, Q. Du, and M. Gunzburger. Recent Progress in Mathematical and Computational Aspects of Peridynamics. Springer International Publishing, 2017.
[11] M. D’Elia and M. Gunzburger. The fractional Laplacian operator on bounded domains as a special case of the nonlocal diffusion operator. Computers & Mathematics with Applications, 66(7):1245 – 1260, 2013.
[12] C. Diyaroglu, E. Oterkus, S. Oterkus, and E. Madenci. Peridynamics for bending of beams and plates with transverse shear deformation. International Journal of Solids and Structures, 69-70:152 – 168, 2015.
[13] Q. Du and J. Yang. Fast and accurate implementation of Fourier spectral approximations of nonlocal diffusion operators and its applications. Journal of Computational Physics, 332:118–134, 2017.
[14] E. Emmrich and D. Puhst. Survey of existence results in nonlinear peridynamics in comparison with local elastodynamics. Comput. Methods Appl. Math., 15(4):483–496, 2015.
[15] L. C. Evans. Partial Differential Equations. American Mathematical Society, 1998.
[16] J. Heo, Z. Yang, W. Xia, S. Oterkus, and E. Oterkus. Buckling analysis of cracked plates using peridynamics. Ocean Engineering, page 107817, 2020.
[17] J. Heo, Z. Yang, W. Xia, S. Oterkus, and E. Oterkus. Free vibration analysis of cracked plates using peridynamics. Ships and Offshore Structures, 15(sup1):S220–S229, 2020.
[18] A. Jafari, M. Ezzati, and A.A. Atai. Static and free vibration analysis of Timoshenko beam based on combined peridynamic-classical theory besides FEM formulation. Computers & Structures, 213:72–81, 2019.
[19] S. Jafarzadeh, A. Larios, and F. Bobaru. Efficient solutions for nonlocal diffusion problems via boundary-adapted spectral methods. Journal of Peridynamics and Nonlocal Modeling, 2:85 – 110, 2020.
[20] S. Jafarzadeh, F. Mousavi, A. Larios, and F. Bobaru. A general and fast convolution-based method for peridynamics: applications to elasticity and brittle fracture. arXiv preprint arXiv:2105.06055v1, 2021.
[21] S. Jafarzadeh, L. Wang, A. Larios, and F. Bobaru. A fast convolution-based method for peridynamic transient diffusion in arbitrary domains. Computer Methods in Applied Mechanics and Engineering, 375:113633, 2021.
[22] B. Kilic and E. Madenci. An adaptive dynamic relaxation method for quasi-static simulations using the peridynamic theory. Theoretical and Applied Fracture Mechanics, 53(3):194 – 204, 2010.
[23] R. Lipton. Dynamic brittle fracture as a small horizon limit of peridynamics. J. Elasticity, 117:21–50, 2014.
[24] L. Lopez and S. F. Pellegrino. A space-time discretization of a nonlinear peridynamic model on a 2D lamina. Computers and Mathematics with Applications. In Press, 2021.
[25] L. Lopez and S. F. Pellegrino. A spectral method with volume penalization for a nonlinear peridynamic model. International Journal for Numerical Methods in Engineering, 122(3):707–725, 2021.
[26] E. Madenci and E. Oterkus. Peridynamic theory and its applications. Springer New York, 2013.
[27] E. Madenci and S. Oterkus. Ordinary state-based peridynamics for thermoviscoelastic deformation. Engineering Fracture Mechanics, 175:31 – 45, 2017.
[28] R. Servadei and E. Valdinoci. Variational methods for non-local operators of elliptic type. Discrete and Continuous Dynamical Systems, 33(5):2105–2137, 2013.
[29] A. Shojaei, T. Mudric, M. Zaccariotto, and U. Galvanetto. A coupled meshless finite point/Peridynamic method for 2D dynamic fracture analysis. International Journal of Mechanical Sciences, 119:419 – 431, 2016.
[30] S.A. Silling. Reformulation of elasticity theory for discontinuities and long-range forces. J. Mech. Phys. Solids, 48(17–18):175–209, 2000.
[31] S.A. Silling and F. Bobaru. Peridynamic modeling of membranes and fibers. International Journal of Non-Linear Mechanics, 40(2):395 – 409, 2005. Special Issue in Honour of C.O. Horgan.
[32] O. Weckner and R. Abeyaratne. The effect of long-range forces on the dynamics of a bar. Journal of the Mechanics and Physics of Solids, 53(3):705 – 728, 2005.
[33] M. Zaccariotto, T. Mudric, D. Tomasi, A. Shojaei, and U. Galvanetto. Coupling of FEM meshes with Peridynamic grids. Computer Methods in Applied Mechanics and Engineering, 330:471 – 497, 2018.
[34] Q. Zhai, H. Xie, R. Zhang, and Z. Zhang. Acceleration of Weak Galerkin Methods for the Laplacian Eigenvalue Problem. Journal of Scientific Computing, 79:914–934, 2019.
[35] Q. Zhai, H. Xie, R. Zhang, and Z. Zhang. The weak Galerkin method for elliptic eigenvalue problems. Communications in Computational Physics, 26(1):160–191, 2019.
[36] K. Zhou and Q. Du. Mathematical and numerical analysis of linear peridynamic models with nonlocal boundary conditions. SIAM Journal on Numerical Analysis, 48(5):1759–1780, 2010.

	$\displaystyle\left\\|w_{k}-w_{N}^{k}\right\\|_{L^{2}(\Omega)}$	$\displaystyle\leq\left\\|w_{k}-P_{N}w_{k}\right\\|_{L^{2}(\Omega)}+\left\\|P_{N}w_{k}-v_{N}^{k}\right\\|_{L^{2}(\Omega)}$		(33)
		$\displaystyle+\left\\|v_{N}^{k}-w_{N}^{k}\right\\|_{L^{2}(\Omega)}.$

	$\displaystyle\left\\|v_{N}^{k}-w_{N}^{k}\right\\|_{L^{2}(\Omega)}$	$\displaystyle=\|1-\left(P_{N}w_{k},w_{N}^{k}\right)\|$
		$\displaystyle=\bigl{\|}1-\|\left(P_{N}w_{k},w_{N}^{k}\right)\|\bigr{\|}$
		$\displaystyle=\bigl{\|}\left\\|w_{k}\right\\|_{L^{2}(\Omega)}-\left\\|v_{N}^{k}\right\\|_{L^{2}(\Omega)}\bigr{\|}$
		$\displaystyle\leq\left\\|w_{k}-v_{N}^{k}\right\\|_{L^{2}(\Omega)}.$

	$\displaystyle\left\\|v_{N}^{k}-w_{N}^{k}\right\\|_{L^{2}(\Omega)}$	$\displaystyle\leq\left\\|w_{k}-v_{N}^{k}\right\\|_{L^{2}(\Omega)}$		(34)
		$\displaystyle\leq\left\\|w_{k}-P_{N}w_{k}\right\\|_{L^{2}(\Omega)}+\left\\|P_{N}w_{k}-v_{N}^{k}\right\\|_{L^{2}(\Omega)}.$