# 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_ADDRESS>Dipartimento di Management, Finanza e Tecnologia,
Università LUM Giuseppe Degennaro, S.S. 100 Km 18 - 70010 Casamassima (BA),
Italy<EMAIL_ADDRESS>
###### 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.

$\delta=3$.

$\delta=1$.

$\delta=0.3$.

$\delta=0.1$.

Figure 1. Comparison of the exact peridynamic eigenvalues $\lambda(\nu_{h})$ defined in (12) with respect to $\nu_{h}$ sampled at 100 equispaced points in the interval $[0,4]$ with the discrete ones computed by means of trigonometric and Chebyshev polynomials for different values of $\delta$. For the simulation we fix $N=100$. 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.

Figure 2. The first five eigenfunctions obtained by using the trigonometric
polynomials.

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

$\delta=1$.

$\delta=0.2$.

Figure 3. The comparison of the exact eigenvalues of the problem (44) in the
interval $[-1,1]$ with the discrete ones computed by means of the Chebyshev
polynomials for different choice of $N$ and $\delta$.

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

Figure 4. The comparison between the exact solution and the numeric one for
different values of $N$ at time $t=1$. On the left, the numeric solution is
computed by using the trigonometric polynomials; on the right, the numeric
solution is computed by using the Chebyshev polynomials.

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.