# On Nonlocal Computation of Eigenfrequencies of Beams Using Finite Difference and Finite Element Methods

Noël Challamel, Vincent Picandet, Issac Elishakoff, Chien Ming Wang, Bernard Collet, Thomas Michelitsch

► To cite this version:

Noël Challamel, Vincent Picandet, Issac Elishakoff, Chien Ming Wang, Bernard Collet, et al. On Nonlocal Computation of Eigenfrequencies of Beams Using Finite Difference and Finite Element Methods. International Journal of Structural Stability and Dynamics, 2015, 15 (7), 10.1142/S0219455415400088 . hal-01115607

HAL Id: hal-01115607

https://hal.science/hal-01115607

Submitted on 3 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On Nonlocal Computation of Eigenfrequencies of Beams Using Finite Difference and Finite Element Methods

Noël Challamel* and Vincent Picandet†

Université Européenne de Bretagne

University of South Brittany UBS, UBS — LIMATB

Centre de Recherche, Rue de Saint Maudé

BP92116, 56321 Lorient cedex — France

*noel.challamel@univ-ubs.fr

†vincent.picandet@univ-ubs.fr

Issac Elishakoff

Department of Ocean and Mechanical Engineering

Florida Atlantic University, Boca Raton, FL 33431-0991, USA

elishako@fau.edu

Chien Ming Wang

Engineering Science Programme and Department of Civil

and Environmental Engineering, National University of Singapore

Kent Ridge, Singapore 119260

ceewcm@nus.edu.sg

Bernard Collet‡ and Thomas Michelitsch§

Surbonne Universités, Université Pierre et Marie Curie

Paris 6, Institut Jean le Rond d'Alembert, CNRS UMR 7190

Case 162, 4 Place Jussieu, 75252 Paris Cedex 05, France

‡bernard.colletl@upmc.fr

§michel@lmm.jussieu.fr

In this paper, we show that two numerical methods, specifically the finite difference method and the finite element method applied to continuous beam dynamics problems, can be asymptotically investigated by some kind of enriched continuum approach (gradient elasticity or nonlocal elasticity). The analysis is restricted to the vibrations of elastic beams, and more specifically the computation of the natural frequencies for each numerical method. The analogy between the finite numerical approaches and the equivalent enriched continuum is demonstrated, using a continualization procedure, which converts the discrete numerical problem into a continuous one. It is shown that the finite element problem can be transformed into a system of finite difference equations. The convergence rate of the finite numerical approaches is quantified by an
equivalent Rayleigh's quotient. We also present some exact analytical solutions for a first-order finite difference method, a higher-order finite difference method or a cubic Hermitian finite element, valid for arbitrary number of nodes or segments. The comparison between the exact numerical solution and the approximated nonlocal approaches shows the efficiency of the continualization methodology. These analogies between enriched continuum and finite numerical schemes provide a new attractive framework for potential applications of enriched continua in computational mechanics.

*Keywords:* Eigenvalue problems; computational mechanics; nonlocal elasticity; gradient elasticity; finite difference methods; finite element methods; vibrations; beam mechanics; cubic Hermitian functions.

# 1. Introduction

In this paper, the vibration behavior of elastic beams will be investigated from a numerical point of view. It is known that numerical approaches such as the finite difference or the finite element methods convert a continuous vibrations problem into a discrete one. In fact, the discretization methods lead to the resolution of an algebraic problem for an initial continuous eigenvalue problem. The possibility to solve automatically the algebraic problem using a computer makes the discretization approach advantageous as compared to the initial continuous one. However, the reduction process inherent to the discretization may incur the loss of some fundamental mathematical properties of the initial continuous system. There appears to be a need to better explore the efficiency of the numerical schemes with respect to the continuous problem.

A possible way to handle these numerical discrete problems is to define a kind of an equivalent continuum that is a representative of the discrete problem. The definition of an equivalent continuum from a discrete one may be labeled as a continualization procedure. Continualization procedures are based on various approximations of the discrete operators by some continuous ones via Taylor expansion or Padé approximants.¹⁻⁷ The so-called enriched continuum equivalent to the discrete one is sometimes called a quasi-continuum.² It is generally dependent on the truncated terms in the asymptotic expansion of the difference operators. This method was pioneered by Kruskal and Zabusky¹ and initially applied to discrete wave equations, with applications to the Fermi-Pasta-Ulam model,⁸ an axial lattice with nonlinear interaction (see also the analysis of Zabusky and Kruskal within the dynamics of solitons⁹). The reader can refer to Rosenau,³ Palais¹⁰ and Maugin¹¹ for an historical perspective on the link between the Fermi-Pasta-Ulam lattice model and the continualized wave propagation equation. Kruskal and Zabusky¹ used a Taylor expansion of the second-order finite difference operator arising in the discrete lattice up to the fourth-order spatial derivative. Collins² proposed to use the inverse of the second-order finite difference operator, thereby avoiding the use of fourth-order spatial operators. Padé approximants of the finite difference operators were introduced by Rosenau³ and are shown to be efficient for capturing the wave propagation in the dynamics of axial lattice.³⁻⁷ It is worth noting that some continualization
processes have been already applied to finite difference structural problems by Cyrus and Fulton,¹²,¹³ or to finite element problems by Walz et al.¹⁴ In this paper, a similar continualization reasoning will be followed for approximating some finite numerical schemes (first-order and higher-order finite difference methods and finite element methods) in beam problems by some equivalent enriched continua. We show that the rate of convergence is strongly dependent on the order of the finite discrete scheme: higher-order finite schemes lead to higher-order enriched constitutive laws with a higher convergence rate.

The comparison of discrete methods with the equivalent continuum is not new. In fact, it was conducted by Livesley,¹⁵ Greenwood,¹⁶ Leckie and Lindberg¹⁷ for beam vibration problems using finite difference methodology for instance. Livesley¹⁵ or Leckie and Lindberg¹⁷ gave exact solutions of the vibration frequencies for the finite difference beam vibration problems. The results have been recently generalized by Zhang et al.¹⁸ for general boundary conditions. In fact, already Lagrange¹⁹ and Rayleigh²⁰ had determined the exact vibration frequencies of a string with a finite number of concentrated masses and compared the solution with the continuous system that was asymptotically obtained as a limit (see also Refs. 15 or 21 on this topic). It can be shown that the discrete string problem is equivalent to the finite difference formulation of the continuous string problem. The same analogy is also valid for elastic beams. As already discussed by Silverman,²² Hencky's bar chain²³ — which is the bending discrete system (or microstructured beam model) — is in fact equivalent to the central finite difference formulation of a continuous problem, i.e., the Euler-Bernoulli continuous beam problem. Therefore, the first-order central finite difference formulation of an Euler-Bernoulli continuous beam problem is strictly equivalent to the Hencky's microstructured chain. The performance of Finite Difference Method for solving buckling or vibration eigenvalue problems has already been evaluated in the literature (see early studies in Refs. 12, 13, 15–17, 24–31), but generally without resorting to any nonlocal mechanics perspective (except in recent papers — see Ref. 32). Furthermore, the nonlocal equivalence will be extended here for higher-order finite difference schemes. A consequence of the nonlocal equivalency principle for the modeling of discrete systems is that the finite difference system can be efficiently approached by nonlocal continuum mechanics tools. As it is known in the case of nonlocal mechanics behaviors, this result confirms the lower bound solution of such approximate Finite Difference Methods, at least for homogeneous structures (with respect to both convergence and rate of convergence arguments). We extend such a result for approximate Finite Element Methods using gradient elasticity constitutive law, which shows the upper bound solution of Finite Element results based on the work-energy formulation. Excellent monographs are available for the computation of eigenfrequencies of structural members using Finite Element Methods (see for instance Refs. 33–37). Exact solutions of the Finite Element formulation of the dynamics of Euler-Bernoulli beams have been given by Tong et al.,³⁸ Belytschko and Mindle³⁹ or Xie and Steven⁴⁰ using a cubic-based interpolation
function for the displacement. These last authors also discussed the possibility to have a lumped mass matrix or a consistent mass matrix with respect to the displacement interpolation field. Belytschko and Mindle³⁹ and Xie and Steven⁴⁰ derived the exact frequency solutions using the consistent mass matrix with cubic-based interpolation function for the displacement. Tong et al.³⁸ also obtained the asymptotic solution of the frequency parameter with respect to the size of the finite element. These numerical results are revisited in this paper using an equivalent gradient elasticity model.

## 2. Continuous Problem

The vibration problem of a continuous elastic Euler–Bernoulli beam is investigated herein. The continuous reference problem is first briefly presented. The elastic bending moment — curvature constitutive law is expressed for the Euler–Bernoulli kinematics by

$$M = EIw'', \quad (1)$$

where M is the bending moment, E the Young modulus, I the area moment of inertia, w the deflection of the beam, and the prime denotes differentiation with respect to x. The equation of motion including the inertia forces is given by

$$M'' = -\mu\ddot{w}, \quad (2)$$

where the superdot denotes the differentiation with respect to time, and μ is the mass per unit length of the beam material. By combining Eqs. (1) and (2), one obtains the fourth-order partial differential equation of motion:

$$EIw'' + \mu\ddot{w} = 0, \quad (3)$$

By assuming harmonic motion with ω as the angular frequency of vibration, the free vibration problem of the Euler–Bernoulli beam is governed by the linear differential equation:

$$EIw'' + \mu\omega^2 w = 0, \quad (4)$$

which can be equivalently handled in a weak format using the principle of virtual work, i.e.,

$$\int_{0}^{L} \{EIw''\delta w'' - \mu\omega^2 w\delta w\} dx = 0. \quad (5)$$

The Rayleigh’s quotient $R[w]$ can be introduced for the computation of $\omega^2$:

$$R[w] = \omega^2 = \frac{\int_{0}^{L} EIw''^2 dx}{\int_{0}^{L} \mu w^2 dx}. \quad (6)$$
For simply supported boundary conditions, the exact natural kth vibration mode is the trigonometric solution given by:

$$w(x) = W \sin\left(\frac{k\pi x}{L}\right). \qquad (7)$$

The substitution of Eq. (7) in the Rayleigh's quotient Eq. (6) leads to the well-known solution of the continuous Euler-Bernoulli beam problem:

$$\omega_{k,\infty}^2 = \frac{EI}{\mu} \left(\frac{k\pi}{L}\right)^4, \qquad (8)$$

where the subscript k refers to the kth mode considered, and the infinite character refers to the continuous reference problem which possesses an infinite number of degrees-of-freedom.

### 3. First-Order Finite Difference Method

The first-order finite difference formulation of this vibration problem will now be presented, and exactly solved for the simply supported boundary condition. The constitutive law Eq. (1) is now written in this finite difference formulation as:

$$M_i = EI \frac{w_{i-1} - 2w_i + w_{i+1}}{a^2}. \qquad (9)$$

a is the uniform grid spacing also defined by $a = L/n$ where n is the number of grid spacing (and $n + 1$ is the number of uniformly spaced grid points). The equilibrium Eq. (2) is expressed with the second-order finite difference equation:

$$\frac{M_{i-1} - 2M_i + M_{i+1}}{a^2} = -\mu\ddot{w}_i. \qquad (10)$$

By using again both the equilibrium equation and the constitutive law expressed in the finite difference formulation, one obtains the discretized equation of motion

$$EI \frac{w_{i+2} - 4w_{i+1} + 6w_i - 4w_{i-1} + w_{i-2}}{a^4} + \mu\ddot{w}_i = 0. \qquad (11)$$

By assuming harmonic motion, the free vibration equation of motion is given by

$$EI \frac{w_{i+2} - 4w_{i+1} + 6w_i - 4w_{i-1} + w_{i-2}}{a^4} - \mu\omega^2 w_i = 0. \qquad (12)$$

An exact solution to this problem can now be investigated. The same methodology can be followed for the vibration equation, as detailed for instance in Ref. 31, from the fourth-order linear finite difference Eq. (12) restricted to the vibration terms:

$$w_{i+2} - 4w_{i+1} + 6w_i - 4w_{i-1} + w_{i-2} - \frac{\Omega^2}{n^4}w_i = 0 \quad \text{with } \Omega^2 = \omega^2 \frac{\mu L^4}{EI}. \qquad (13)$$
The characteristic equation is obtained by replacing $w_i = A\lambda^i$ in Eq. (13) which leads to:

$$ \left(\frac{1}{\lambda} + \lambda\right)^2 - 4\left(\frac{1}{\lambda} + \lambda\right) + 4 - \frac{\Omega^2}{n^4} = 0, \quad (14) $$

where A is a constant.

Equation (14) is symmetrical with respect to interchanging $\lambda$ and $1/\lambda$, and admits four solutions written as (see also Ref. 31):

$$ \begin{gathered} \lambda_{1,2} = \cos \phi \pm j \sin \phi \quad \text{and} \quad \lambda_{3,4} = 2 - \cos \phi \pm \sqrt{(2 - \cos \phi)^2 - 1} \\ \text{with} \quad \phi = \arccos \left(1 - \frac{\Omega}{2n^2}\right) \quad \text{and} \quad j^2 = -1. \end{gathered} \quad (15) $$

For the simply supported discrete beam system, the natural vibration modes are obtained from the trigonometric shape function $w_i = W \sin(\phi i)$, thus leading to the natural vibration frequency $\phi n = k\pi$ as obtained in Refs. 17 or 31:

$$ \Omega_{k,n} = 4n^2 \sin^2 \left( \frac{k\pi}{2n} \right) = (k\pi)^2 \left( 1 - \frac{(k\pi)^2}{12n^2} \right) + o \left( \frac{1}{n^4} \right). \quad (16) $$

We then obtain for the square of the natural frequencies:

$$ \frac{\omega_{\text{discrete}}^2}{\omega_E^2} = 1 - \frac{(k\pi)^2}{6n^2} + o\left(\frac{1}{n^4}\right), \quad (17) $$

where $\omega_E^2 = \frac{EI}{\mu} \left(\frac{k\pi}{L}\right)^4$ is the natural frequency parameter of the continuous beam.

The discrete equations are extended to an equivalent continuum via a con- tinualization method. The following relation between the discrete and the equivalent continuous system $w_i = w(x = ia)$ holds for a sufficiently smooth deflection function as:

$$ w(x + a) = \sum_{k=0}^{\infty} \frac{a^k \partial_x^k}{k!} w(x) = e^{a \partial_x} w(x) \quad \text{with} \quad \partial_x = \frac{\partial}{\partial x}. \quad (18) $$

The pseudo-differential operators can be introduced as:

$$ w_{i-1} + w_{i+1} - 2w_i = [e^{a\partial_x} + e^{-a\partial_x} - 2]w(x) = 4\sinh^2\left(\frac{a}{2}\partial_x\right)w(x). \quad (19) $$

The pseudo-differential operator can be efficiently approximated by the Padé's approximant (see for instance Refs. 3-7):

$$ \frac{4}{a^2} \sinh^2\left(\frac{a}{2} \partial_x\right) = \frac{\partial_x^2}{1 - l_c^2 \partial_x^2} + \dots \quad \text{with} \quad l_c^2 = \frac{a^2}{12}. \quad (20) $$

By using such a pseudo-differential operator, the constitutive law Eq. (9) using Eqs. (19) and (20) may be continually approximated as:

$$ M - l_c^2 M'' = EIw'' \quad \text{with} \quad l_c^2 = \frac{a^2}{12}. \quad (21) $$
One recognizes an Eringen's type differential equation<sup>41</sup> applied at the beam scale, thereby leading to a nonlocal bending moment — curvature constitutive law.

By using the same methodology, the equilibrium equations [Eq. (10)] may be also continualized as:

$$M'' = -\mu\ddot{w} + \mu l_c^2 \dot{w}'' \quad \text{with} \quad l_c^2 = \frac{a^2}{12}. \qquad (22)$$

In view of Eqs. (21) and (22) and neglecting the terms in $l_c^4$, one obtains the modified nonlocal bending wave equation as:

$$\mu(1 - 2l_c^2 \partial_x^2) \partial_t^2 w + EI \partial_x^4 w = 0 \quad \text{with} \quad l_c^2 = \frac{a^2}{12}. \qquad (23)$$

Equation (23) can be also directly obtained from the expansion of the discrete operator in Eq. (12) as:

$$\frac{w_{i+2} - 4w_{i+1} + 6w_i - 4w_{i-1} + w_{i-2}}{a^4} = \frac{[e^{2a\partial_x} - 4e^{a\partial_x} + 6 - 4e^{-a\partial_x} + e^{-2a\partial_x}]}{a^4} w(x) \\ = \frac{16}{a^4} \sinh^4\left(\frac{a}{2}\partial_x\right) w(x) = \partial_x^4\left(1 + \frac{a^2}{6}\partial_x^2\right) w(x) + \dots = \frac{\partial_x^4}{1 - \frac{a^2}{6}\partial_x^2} w(x) + \dots. \qquad (24)$$

The continualized problem governed by Eq. (23) is also equivalent to the following nonlocal model, where the nonlocal operator is applied only to the constitutive law and the equilibrium equations remain local:

$$M - 2l_c^2 M'' = EIw'' \quad \text{and} \quad M'' = -\mu\ddot{w} \quad \text{with} \quad l_c^2 = \frac{a^2}{12}. \qquad (25)$$

Within this point of view, it is worth mentioning that a factor 2 affects the length scale calibration in the nonlocal law.

The Rayleigh's quotient for the computation of $\omega^2$ of the nonlocal problem can be presented as:

$$R[w] = \frac{\int_{0}^{L} EIw''^2 dx}{\int_{0}^{L} \mu w^2 + 2\mu l_c^2 w'^2 dx} \le \frac{\int_{0}^{L} EIw''^2 dx}{\int_{0}^{L} \mu w^2 dx}. \qquad (26)$$

For the simply supported beam problem, the substitution of Eq. (7) in the Rayleigh's quotient Eq. (26) leads to the nonlocal solution of the continuous Euler-Bernoulli beam problem:

$$\omega_{k,n}^2 = \frac{\omega_{k,\infty}^2}{1 + \left(\frac{knl_c}{L}\right)^2} \quad \text{with} \quad \omega_{k,\infty}^2 = \frac{EI}{\mu} \left(\frac{k\pi}{L}\right)^4 \quad \text{and} \quad l_c^2 = \frac{a^2}{6}. \qquad (27)$$

Equation (27) is consistent with the asymptotic expansion of the exact solution given by Eq. (17), and it also coincides with the asymptotic expansion given by Leckie and Lindberg.<sup>17</sup> We note that there is a factor 2 between the equivalent length scale for the buckling problem and the one of the dynamics problem<sup>32</sup>:

$$l_{c,dynamics}^2 = \frac{a^2}{6} = 2l_{c,statics}^2. \qquad (28)$$
## 4. Higher-Order Finite Difference Method

The same continuous problem will now be handled using a higher-order finite difference scheme. An improved finite difference analysis may be based on the introduction of the second-order central difference for the expressions of the first and the second derivatives of the displacement (see Refs. 16, 28–30). The higher-order finite difference formulation of the constitutive equation [Eq. (1)] now reads as

$$M_i = EI \frac{-w_{i-2} + 16w_{i-1} - 30w_i + 16w_{i+1} - w_{i+2}}{12a^2} \quad (29)$$

The higher-order finite difference equilibrium equation [Eq. (2)] is presented as:

$$\frac{-M_{i-2} + 16M_{i-1} - 30M_i + 16M_{i+1} - M_{i+2}}{12a^2} = -\mu\dot{w}_i. \quad (30)$$

The consideration of harmonic motion and in view of Eqs. (29) and (30) lead to the eight-order finite difference equation:

$$EI \frac{w_{i-4} - 32w_{i-3} + 316w_{i-2} - 992w_{i-1} + 1414w_i - 992w_{i+1} + 316w_{i+2} - 32w_{i+3} + w_{i+4}}{144a^4} - \mu\omega^2 w_i = 0. \quad (31)$$

It is also possible to get an exact analytical solution to this numerical scheme, by solving the linear finite difference equation. The characteristic equation is obtained by replacing $w_i = A\lambda^i$ in Eq. (31) which leads, with $\Omega^2 = \omega^2 \frac{\mu L_i^4}{EI}$ to:

$$\begin{aligned} & (\lambda^{-4} - 32\lambda^{-3} + 316\lambda^{-2} - 992\lambda^{-1} + 1414 - 992\lambda + 316\lambda^2 - 32\lambda^3 + \lambda^4) - 144 \frac{\Omega^2}{n^4} \\ & = 0, \end{aligned} \quad (32)$$

which can also be presented as:

$$\left[ \left( \frac{1}{\lambda} + \lambda \right) - 2 \right]^2 \left[ \left( \frac{1}{\lambda} + \lambda \right) - 14 \right]^2 - 144 \frac{\Omega^2}{n^4} = 0. \quad (33)$$

It is possible to factorize in the following way:

$$\left\{ \left[ \left( \frac{1}{\lambda} + \lambda \right) - 2 \right] \left[ \left( \frac{1}{\lambda} + \lambda \right) - 14 \right] - 12 \frac{\Omega}{n^2} \right\} \times \left\{ \left[ \left( \frac{1}{\lambda} + \lambda \right) - 2 \right] \left[ \left( \frac{1}{\lambda} + \lambda \right) - 14 \right] + 12 \frac{\Omega}{n^2} \right\} = 0. \quad (34)$$

This equation admits eight solutions for $\lambda$. By solving the first term, one gets

$$\left(\frac{1}{\lambda} + \lambda\right)^2 - 16\left(\frac{1}{\lambda} + \lambda\right) + 28 - 12\frac{\Omega}{n^2} = 0$$

or equivalently

$$\lambda^2 - 8 \pm 2\sqrt{3\left(3 + \frac{\Omega}{n^2}\right)} \lambda + 1 = 0. \quad (35)$$
The first group of solutions are obtained as:

$$ 
\begin{aligned}
\lambda_{1,2} &= 4 + \sqrt{3\left(3 + \frac{\Omega}{n^2}\right)} \pm \sqrt{15 + 8\sqrt{3\left(3 + \frac{\Omega}{n^2}\right)} + 3\left(3 + \frac{\Omega}{n^2}\right)}, \\
\lambda_{3,4} &= 4 - \sqrt{3\left(3 + \frac{\Omega}{n^2}\right)} \pm \sqrt{15 - 8\sqrt{3\left(3 + \frac{\Omega}{n^2}\right)} + 3\left(3 + \frac{\Omega}{n^2}\right)}.
\end{aligned} 
\quad (36) $$

The second group of solutions comes from the second-order polynomial equation:

$$ \left(\frac{1}{\lambda} + \lambda\right)^2 - 16\left(\frac{1}{\lambda} + \lambda\right) + 28 + 12\frac{\Omega}{n^2} = 0 $$

or equivalently

$$ \lambda^2 - 8 \pm 2\sqrt{3\left(3 - \frac{\Omega}{n^2}\right)} \lambda + 1 = 0, \quad (37) $$

which also possesses four solutions:

$$ \lambda_{5,6} = 4 + \sqrt{3\left(3 - \frac{\Omega}{n^2}\right)} \pm \sqrt{15 + 8\sqrt{3\left(3 - \frac{\Omega}{n^2}\right)} + 3\left(3 - \frac{\Omega}{n^2}\right)}, \quad (38) $$

$$ \lambda_{7,8} = 4 - \sqrt{3\left(3 - \frac{\Omega}{n^2}\right)} \pm \sqrt{15 - 8\sqrt{3\left(3 - \frac{\Omega}{n^2}\right)} + 3\left(3 - \frac{\Omega}{n^2}\right)}. $$

It is also possible to rewrite the set of solutions $\lambda_{3,4}$ using the trigonometric functions:

$$ \lambda_{3,4} = \cos \phi \pm j \sin \phi \quad \text{with} \quad \phi = \arccos \left[ 4 - \sqrt{3\left(3 + \frac{\Omega}{n^2}\right)} \right] \quad \text{and} \quad j^2 = -1. \quad (39) $$

For the simply supported discrete system, the natural vibration modes are obtained from the trigonometric shape function $w_i = W \sin(\phi i)$, thus leading to the natural vibration frequency $\theta n = k\pi$, which can also be expressed by:

$$ \frac{\Omega}{n^2} = \frac{7 - 8 \cos \phi + \cos^2 \phi}{3} \quad \text{with} \quad \phi = \arccos \left[ 4 - \sqrt{3\left(3 + \frac{\Omega}{n^2}\right)} \right], \quad (40) $$

which has hitherto not been reported to the best of the authors' knowledge. It is possible to show from Eq. (40) by using an asymptotic expansion that:

$$ \frac{\Omega_{k,n}^2}{\Omega_{k,\infty}^2} = \frac{\omega_{k,n}^2}{\omega_{k,\infty}^2} = 1 - \frac{(k\pi)^4}{45n^4} + o\left(\frac{1}{n^6}\right). \quad (41) $$
The higher-order pseudo-differential operator in Eq. (29) or in Eq. (30) can also be
expanded as:

$$
\begin{align}
& \frac{-w_{i-2} + 16w_{i-1} - 30w_i + 16w_{i+1} - w_{i+2}}{12a^2} \nonumber \\
&= \frac{[-e^{-2a\partial_x} + 16e^{-a\partial_x} - 30 + 16e^{a\partial_x} - e^{2a\partial_x}]}{12a^2} w(x) \nonumber \\
&= \partial_x^2 \left(1 - \frac{a^4 \partial_x^4}{90}\right) w(x) + \dots = \frac{\partial_x^2}{1 + \frac{a^4 \partial_x^4}{90}} w(x) + \dots \tag{42}
\end{align}
$$

The constitutive law Eq. (29) can then be continualized such as:

$$
M + l_c^4 M^{(4)} = EIw'' \quad \text{with} \quad l_c^4 = \frac{a^4}{90}. \qquad (43)
$$

The equilibrium equations [Eq. (30)] can be continualized as well:

$$
M'' = -\mu\ddot{w} - \mu l_c^4 \dot{w}^{(4)} \quad \text{with} \quad l_c^4 = \frac{a^4}{90}. \qquad (44)
$$

By combining Eq. (43) with Eq. (44), one obtains the nonlocal bending wave equation as

$$
\mu(1 + 2l_c^4 \partial_x^4) \partial_t^2 w + EI \partial_x^4 w = 0 \quad \text{with} \quad l_c^4 = \frac{a^4}{90}. \qquad (45)
$$

Equation (45) could have been obtained directly from the asymptotic expansion:

$$
\begin{align}
& \frac{w_{i-4} - 32w_{i-3} + 316w_{i-2} - 992w_{i-1} + 1414w_i - 992w_{i+1} + 316w_{i+2} - 32w_{i+3} + w_{i+4}}{144a^4} \nonumber \\
&= \partial_x^4 \left(1 - \frac{a^4}{45}\partial_x^4\right) w(x) + \dots = \frac{\partial_x^4}{1 + \frac{a^4}{45}\partial_x^4} w(x) + \dots \tag{46}
\end{align}
$$

Equation (46) is also equivalent to the nonlocal model:

$$
M + 2l_c^4 M^{(4)} = EIw'' \quad \text{and} \quad M'' = -\mu\ddot{w} \quad \text{with} \quad l_c^4 = \frac{a^4}{90}. \qquad (47)
$$

As observed for the first-order finite difference method, a factor 2 also affects the
length scale calibration in the nonlocal law for the higher-order finite difference
method. The Rayleigh's quotient for the computation of $\omega^2$ of the new nonlocal
problem can be presented as:

$$
R = \frac{\int_{0}^{L} EIw''^2 dx}{\int_{0}^{L} \mu w^2 + 2\mu l_c^4 w''^2 dx} \le \frac{\int_{0}^{L} EIw''^2 dx}{\int_{0}^{L} \mu w^2 dx}. \quad (48)
$$

For the simply supported beam problem, the substitution of Eq. (7) in the Rayleigh’s quotient Eq. (48) leads to the nonlocal solution of the continuous Euler–Bernoulli beam problem:

$$
\omega_{k,n}^2 = \frac{\omega_{k,\infty}^2}{1 + \left(\frac{k\pi l_c}{L}\right)^4} \quad \text{with} \quad \omega_{k,\infty}^2 = \frac{EI}{\mu} \left(\frac{k\pi}{L}\right)^4 \quad \text{and} \quad l_c^4 = \frac{a^4}{45}. \qquad (49)
$$
We also note that there is a factor 2 between the equivalent length scale for the (static) buckling problem and the one of the dynamics problem:

$$l_{c,dynamics}^4 = \frac{a^4}{45} = 2l_{c,statics}^4. \qquad (50)$$

## 5. Finite Element Method

The Finite Element Method is now applied to the beam vibrations problem, using Hermitian cubic-based interpolation functions. The Hermitian cubic functions can be used for the interpolation function of the displacement field:

$$w = w_{i-1}(1 - 3\xi^2 + 2\xi^3) + w_i\xi^2(3 - 2\xi) + \theta_{i-1}a\xi(1 - \xi)^2 - \theta_i a\xi^2(1 - \xi), \quad (51)$$

where $\xi = x/a$. By substituting this Hermitian cubic function into the Rayleigh's quotient, one obtains

$$R[w, \theta] = \frac{\sum_{i=1}^{n} \int_{0}^{1} \frac{EI}{a^4} \left( \frac{d^2 w}{d\xi^2} \right)^2 d\xi}{\sum_{i=1}^{n} \int_{0}^{1} \mu w^2 d\xi}, \qquad (52)$$

which is now calculated for the cubic-based Hermitian interpolation function:

$$R[w, \theta] = \frac{\sum_{i=1}^{n} \frac{EI}{3a^4} \left[ (6w_i - 6w_{i-1})^2 - (6w_i - 6w_{i-1})(6a\theta_{i-1} + 6a\theta_i) + 12a^2(\theta_{i-1}^2 + \theta_{i-1}\theta_i + \theta_i^2) \right]}{\sum_{i=1}^{n} \frac{\mu}{35} \left[ (13w_{i-1}^2 + 9w_{i-1}w_i + 13w_i^2) + a^2 \left( \frac{\theta_{i-1}^2}{3} - \frac{\theta_{i-1}\theta_i}{2} + \frac{\theta_i^2}{3} \right) + \frac{11}{3}(w_{i-1}a\theta_{i-1} - w_i a\theta_i) + \frac{13}{6}(w_i a\theta_{i-1} - w_{i-1} a\theta_i) \right]} \quad (53)$$

By taking the stationarity conditions of the Rayleigh's quotient $R[w, \theta] = 0$ for the two-variable field $(w_i, \theta_i)$, we obtain the coupled system of finite difference equations:

$$\left\{
\begin{aligned}
& 4(w_{i+1} - 2w_i + w_{i-1}) - 2(a\theta_{i+1} - a\theta_{i-1}) \\
& \quad + \frac{\mu\omega^2 a^4}{210EI} \left[ 9w_{i+1} + 52w_i + 9w_{i-1} - \frac{13}{6}(a\theta_{i+1} - a\theta_{i-1}) \right] = 0, \\[2ex]
& 6(w_{i+1} - w_{i-1}) - 2(a\theta_{i+1} + 4a\theta_i + a\theta_{i-1}) \\
& \quad + \frac{\mu\omega^2 a^4}{70EI} \left[ \frac{13}{6}(w_{i+1} - w_{i-1}) + \frac{-3a\theta_{i+1} + 8a\theta_i - 3a\theta_{i-1}}{6} \right] = 0,
\end{aligned}
\right.
\qquad (54)$$

which would have been equivalently obtained by using the weak formulation Eq. (5) of the problem with the Hermitian cubic functions. It is possible to define the
following finite difference operators

$$
\begin{align}
\delta_0 &= \frac{1}{6} [e^{-a\partial_x} + 4 + e^{a\partial_x}]; \quad \delta_1 = \frac{1}{2a} [-e^{-a\partial_x} + e^{a\partial_x}] \quad \text{and} \tag{55} \\
\delta_2 &= \frac{1}{a^2} [e^{-a\partial_x} - 2 + e^{a\partial_x}]. \nonumber
\end{align}
$$

The finite difference system associated with the Finite Element Method can then be
presented using the finite difference operators:

$$
\left\{
\begin{aligned}
& \left( 4\delta_2 + \frac{3\mu\omega^2 a^4}{70EI} \delta_2 + \frac{\mu\omega^2 a^2}{3EI} \right) w - \left( 4 + \frac{13\mu\omega^2 a^4}{630EI} \right) \delta_1 \theta = 0, \\
& \left( 12 + \frac{13\mu\omega^2 a^4}{210EI} \right) \delta_1 w + \left[ -12\delta_0 - \frac{3\mu\omega^2 a^4}{70EI} \delta_0 + \frac{\mu\omega^2 a^4}{21EI} \right] \theta = 0.
\end{aligned}
\right.
\tag{56}
$$

The finite difference equation is then obtained in a single format as:

$$
\left[
\begin{aligned}
& \left( 4\delta_2 + \frac{3\mu\omega^2 a^4}{70EI} \delta_2 + \frac{\mu\omega^2 a^2}{3EI} \right) \left( -12\delta_0 - \frac{3\mu\omega^2 a^4}{70EI} \delta_0 + \frac{\mu\omega^2 a^4}{21EI} \right) \\
& \qquad + \left( 4 + \frac{13\mu\omega^2 a^4}{630EI} \right) \left( 12 + \frac{13\mu\omega^2 a^4}{210EI} \right) \delta_1^2 \\
& \quad - w = 0,
\end{aligned}
\right] \tag{57}
$$

which can also be written as

$$
\begin{equation}
\begin{aligned}
& \frac{\mu^2 \omega^4 a^6}{63(EI)^2} w_i + \frac{\mu \omega^2 a^4}{21EI} \left(4 + \frac{3\mu \omega^2 a^4}{70EI}\right) \frac{w_{i+1} - 2w_i + w_{i-1}}{a^2} \\
& \quad + \frac{\mu \omega^2 a^2}{3EI} \left(-12 - \frac{3\mu \omega^2}{70EI} a^4\right) \frac{w_{i+1} + 4w_i + w_{i-1}}{6} \\
& \quad + \left(4 + \frac{3\mu \omega^2 a^4}{70EI}\right) \left(-12 - \frac{3\mu \omega^2}{70EI} a^4\right) \frac{w_{i-2} + 2w_{i-1} - 6w_i + 2w_{i+1} + w_{i+2}}{6a^2} \\
& \quad + 3 \left(4 + \frac{13\mu \omega^2 a^4}{630EI}\right)^2 \frac{w_{i-2} - 2w_i + w_{i+2}}{4a^2} = 0.
\end{aligned}
\tag{58}
\end{equation}
$$

The characteristic equation is obtained by replacing $w_i = A\lambda^i$ in Eq. (58) which
leads, with $\Omega^2 = \omega^2 \frac{\mu L^4}{EI}$ to:

$$
\left(48 + \frac{4\Omega^2}{35n^4} + \frac{\Omega^4}{6300n^8}\right) \left(\lambda + \frac{1}{\lambda}\right)^2 + \left(-192 - \frac{296\Omega^2}{35n^4} - \frac{2\Omega^4}{175n^8}\right) \left(\lambda + \frac{1}{\lambda}\right) \\
+ 192 - \frac{1104\Omega^2}{35n^4} + \frac{13\Omega^4}{315n^8} = 0, \tag{59}
$$

that admits the following four solutions

$$
\lambda_{1,2} = \cos \phi \pm j \sin \phi \quad \text{and} \quad \lambda_{3,4} = \frac{192 + \frac{296\Omega^2}{35n^4} + \frac{2\Omega^4}{175n^8}}{2\left(48 + \frac{4\Omega^2}{35n^4} + \frac{\Omega^4}{6300n^8}\right)} \\
- \cos \phi \pm \sqrt{\left[ \frac{192 + \frac{296\Omega^2}{35n^4} + \frac{2\Omega^4}{175n^8}}{2\left(48 + \frac{4\Omega^2}{35n^4} + \frac{\Omega^4}{6300n^8}\right)} - \cos \phi \right]^2 - 1} \quad \text{with}
$$
$$ \phi = \arccos \left[ \frac{192 + \frac{296\Omega^2}{35n^4} + \frac{2\Omega^4}{175n^8} - \sqrt{\left(192 + \frac{296\Omega^2}{35n^4} + \frac{2\Omega^4}{175n^8}\right)^2 - 4\left(48 + \frac{4\Omega^2}{35n^4} + \frac{\Omega^4}{6300n^8}\right)\left(192 - \frac{1104\Omega^2}{35n^4} + \frac{13\Omega^4}{315n^8}\right)}}{4\left(48 + \frac{4\Omega^2}{35n^4} + \frac{\Omega^4}{6300n^8}\right)} \right]. \quad (60) $$

A second-order polynomial equation in $\Omega^2$ is obtained from Eq. (60):

$$ \frac{\Omega^4}{n^8} \left( \frac{4\cos^2\phi}{1575} - \frac{16\cos\phi}{175} + \frac{52}{315} \right) + \frac{\Omega^2}{n^4} \left( \frac{64\cos^2\phi}{35} - \frac{2368\cos\phi}{35} - \frac{4416}{35} \right) + 768\cos^2\phi - 1536\cos\phi + 768 = 0. \quad (61) $$

This equation has been first obtained by Belytschko and Mindle³⁹ and then used again by Xie and Steven.⁴⁰ For the simply supported discrete beam system, the natural vibration modes are obtained from the trigonometric shape function $w_i = W \sin(\phi i)$, thus leading to the natural vibration frequency $\phi n = k\pi$, which is also expressed by:

$$ \frac{\Omega^2}{n^4} = \frac{\frac{-64\cos^2\phi + 2368\cos\phi + 4416}{35}}{2\left(\frac{4\cos^2\phi}{1575} - \frac{16\cos\phi}{175} + \frac{52}{315}\right)} \frac{\sqrt{\left(\frac{-64\cos^2\phi + 2368\cos\phi + 4416}{35}\right)^2 - 4\left(\frac{4\cos^2\phi}{1575} - \frac{16\cos\phi}{175} + \frac{52}{315}\right)\left(768\cos^2\phi - 1536\cos\phi + 768\right)}}{2\left(\frac{4\cos^2\phi}{1575} - \frac{16\cos\phi}{175} + \frac{52}{315}\right)} \quad (62) $$

An asymptotic expansion shows that:

$$ \Omega_{k,n} = (k\pi)^2 \left(1 + \frac{(k\pi)^4}{1440n^4}\right) + o\left(\frac{1}{n^6}\right). \quad (63) $$

We then obtain for the square of the natural frequencies:

$$ \frac{\omega_{\text{discrete}}^2}{\omega_E^2} = 1 + \frac{(k\pi)^4}{720n^4} + o\left(\frac{1}{n^6}\right), \quad (64) $$

where $\omega_E^2 = \frac{EI}{\mu} (\frac{k\pi}{L})^4$ is the natural frequency parameter of the continuous beam. This result Eq. (64) with the factor 1/720 has been already obtained by Tong et al.³⁸

Now, by using a continualization procedure, an asymptotic expansion of each difference operator in Eq. (57) gives:

$$ \left\{ \left[ \left( 4 + \frac{(a\partial_x)^2}{3} + \frac{1}{90}(a\partial_x)^4 + \frac{1}{5040}(a\partial_x)^6 \right) \partial_x^2 + \frac{3\mu\omega^2 a^4}{70EI} \partial_x^2 + \frac{3\mu\omega^2 a^6}{840EI} \partial_x^4 + \frac{\mu\omega^2 a^2}{3EI} \right] \times \left[ -12 - 2a^2 \partial_x^2 - \frac{1}{6}(a\partial_x)^4 - \frac{1}{180}(a\partial_x)^6 - \frac{3\mu\omega^2 a^4}{70EI} - \frac{\mu\omega^2 a^6}{140EI} \partial_x^2 + \frac{\mu\omega^2 a^4}{21EI} \right] + \left( 48 + \frac{52\mu\omega^2 a^4}{105EI} \right) \left( 1 + \frac{(a\partial_x)^2}{3} + \frac{2}{45}(a\partial_x)^4 + \frac{1}{315}(a\partial_x)^6 \right) \partial_x^2 + o(a^8) \right\} w = 0, \quad (65) $$
which can be efficiently approximated by the following sixth-order differential equation, when collecting the terms up to the fourth-order in $a^4$:

$$ \left\{ \partial_x^4 - \frac{\mu\omega^2}{EI} + a^2 \left( \frac{\partial_x^6}{6} - \frac{\mu\omega^2\partial_x^2}{6EI} \right) + a^4 \frac{\partial_x^8}{80} - \frac{29\mu\omega^2\partial_x^4}{2520EI} + \frac{1}{2520} \left[ \frac{\mu\omega^2}{EI} \right]^2 \right\} w = 0. \quad (66) $$

The differential equation [Eq. (66)] can also be factorized as

$$ \begin{aligned} & \partial_x^4 - \frac{\mu\omega^2}{EI} + a^2 \left( \frac{\partial_x^6}{6} - \frac{\mu\omega^2\partial_x^2}{6EI} \right) + a^4 \frac{\partial_x^8}{80} - \frac{29\mu\omega^2\partial_x^4}{2520EI} + \frac{1}{2520} \left[ \frac{\mu\omega^2}{EI} \right]^2 \\ & = \left[ 1 + \frac{(a\partial_x)^2}{6} + \frac{(a\partial_x)^4}{90} - \frac{\mu\omega^2 a^4}{2520EI} \right] \left[ \partial_x^4 - \frac{\mu\omega^2}{EI} + \frac{a^4}{720} \partial_x^8 \right] + o(a^6), \end{aligned} \quad (67) $$

which means that the cubic-based finite element model can be equivalently reduced to the eight-order differential equation:

$$ EI \frac{a^4}{720} w^{(8)} + EIw^{(4)} - \mu\omega^2 w = 0. \quad (68) $$

Walz et al.¹⁴ also obtained an eight-order differential equation for the continualized bending problem which was investigated by the Hermitian-based Finite Element model, with the correct coefficient 1/720 but with a different sign. For the finite element model considered herein, the associated Rayleigh's quotient can then be expressed by:

$$ R = \frac{\int_{0}^{L} EI(w''^2 + l_c^4 w^{(4)^2}) dx}{\int_{0}^{L} \mu w^2 dx} \geq \frac{\int_{0}^{L} EIw''^2 dx}{\int_{0}^{L} \mu w^2 dx} \quad \text{with } l_c^4 = \frac{a^4}{720} \quad (69) $$

leading to the gradient elasticity solution Eq. (68), associated with the gradient elasticity constitutive law:

$$ M = EI(w'' + l_c^4 w^{(6)}) \quad \text{with} \quad l_c^4 = \frac{a^4}{720}. \quad (70) $$

Considering again a simply supported beam, and introducing a sinusoidal shape function $w(x) = W \sin(\pi x/L)$ as a test function into the Rayleigh's quotient leads to the natural frequencies from this approximated gradient elasticity solution. Equation (69) shows that the Finite Element model gives an upper bound of the “local” problem asymptotically found for $n$ tending towards infinite.

$$ \frac{\omega_{k,n}^2}{\omega_{k,\infty}^2} = 1 + \frac{(k\pi)^4}{720n^4}. \quad (71) $$

Figures 1 and 2 show the comparison between the exact numerical solution (first-order finite difference method, higher-order finite difference method, and Hermitian cubic-based finite element method). The closeness of the approximated nonlocal
<figure><img src="image_2.png" /><figcaption>Fig. 1. Comparison between finite difference methods, finite element methods and continualization procedure for fundamental frequency with respect to the size of discretized problem.</figcaption></figure>

approaches with respect to the exact discretized problem shows the efficiency of the continualization methodology. The upper bound status of the finite element approach and the lower bound status of the finite difference approach are confirmed in Fig. 1. Moreover, as shown by Fig. 2, the finite element solution based on the cubic

<figure><img src="image_3.png" /><figcaption>Fig. 2. Error function with respect to size of discretized problem for each numerical method; Comparison between exact solution and approximated nonlocal or gradient elasticity solution.</figcaption></figure>
Hermitian interpolation function appears to be the most efficient numerical approach for this problem in comparison to both the first-order and the higher-order finite difference approaches. Of course, the efficiency of the finite element method strongly depends on the adopted displacement interpolation field.

## 6. Concluding Remarks

The vibration behavior of elastic beams is studied using standard numerical methods such as finite difference or finite element methods. Some exact solutions of the finite-dimensional discretized problems are presented for archetypal boundary conditions. The finite numerical approaches are then approximated by some enriched nonlocal or gradient elastic problems. The analogy between the finite numerical approaches and the equivalent enriched continuum is demonstrated, using a continualization procedure, which converts the discrete numerical problem into a continuous nonlocal or gradient one. For different orders of finite difference (FD) schemes, different nonlocal length scales are obtained, with respect to the step size $a$ of the numerical scheme, i.e. $l_c = a/\sqrt{6} \approx 0.408a$ for first-order FD scheme as shown in Eq. (27), and $l_c = a/\sqrt{45} \approx a/\sqrt{6.708} \approx 0.386a$ for higher-order FD scheme as shown in Eq. (49). In both cases, the nonlocal length scale of the dynamics problem is larger than the one of the statics problem. The finite difference schemes leads to a lower bound of the “local” continuous problem.

For the cubic Hermitian finite element, the discretized system is efficiently approximated by a gradient elasticity law. In this case, the gradient elasticity length scale is the same than the one of the statics problem. The finite element formulation leads to an upper bound of the “local” continuous problem. The comparison between the exact numerical solution and the approximated nonlocal approaches shows the efficiency of the continualization methodology. These analogies between enriched continuum and finite numerical schemes provide a new attractive framework for potential applications of enriched continua in computational mechanics.

## References

1. M. D. Kruskal and N. J. Zabusky, Stroboscopic perturbation for treating a class of nonlinear wave equations, *J. Math. Phys.* **5** (1964) 231–244.

2. M. A. Collins, A quasicontinuum approximation for solitons in an atomic chain, *Chem. Phys. Lett.* **77**(2) (1981) 342–347.

3. P. Rosenau, Dynamics of nonlinear mass-spring chains near the continuum limit, *Phys. Lett. A* **118**(5) (1986) 222–227.

4. P. Rosenau, Dynamics of dense lattices, *Phys. Rev. B* **36**(11) (1987) 5868–5876.

5. J. A. D. Wattis, Quasi-continuum approximations to lattice equations arising from the discrete non-linear telegraph equation, *J. Phys. A* **33** (2000) 5925–5944.

6. P. G. Kevrekidis, I. G. Kevrekidis, A. R. Bishop and E. S. Titi, Continuum approach to discreteness, *Phys. Rev. E* **65**(046613) (2002) 1–13.

7. I. V. Andrianov, J. Awrejcewicz and D. Weichert, Improved continuous models for discrete media, *Math. Prob. Eng.* **986242** (2010) 1–35.
8. E. Fermi, J. Pasta and S. Ulam, Studies of nonlinear problems. Los Alamos Report LA-1940 (1955); Published later in *Collected Papers of Enrico Fermi*, E. Segré (eds.) (University of Chicago Press, 1965).

9. N. J. Zabusky and M. D. Kruskal, Interaction of solitons in a collisionless plasma and the recurrence of initial states, *Phys. Rev. Lett.* **15** (1965) 240–243.

10. R. S. Palais, The symmetries of solitons, *Bull. Am. Math. Soc.* **34**(4) (1997) 339–403.

11. G. A. Maugin, Solitons in elastic solids (1938)–(2010), *Mech. Res. Comm.* **38** (2011) 341–349.

12. N. J. Cyrus and R. E. Fulton, Finite difference accuracy in structural analysis, *Proc. ASCE J. Struct. Div.* **92**(576) (1966) 459–471.

13. N. J. Cyrus and R. E. Fulton, Accuracy study of finite difference methods, NASA Technical Note, NASA-TN D-4372, NASA Langley Research Center (1968), pp. 1–29.

14. J. E. Walz, R. E. Fulton and N. J. Cyrus, Accuracy and convergence of finite element approximations, Report AFFDL-TR-68-150, Air Force Flight Dynamics Laboratory (1968) 996–1027.

15. R. K. Livesley, The equivalence of continuous and discrete mass distributions in certain vibration problems, *Q. J. Mech. Appl. Math.* **8**(3) (1955) 353–360.

16. D. T. Greenwood, The use of higher-order finite differences in beam vibrations, NASA TND-964 (1961).

17. F. A. Leckie and G. M. Lindberg, The effect of lumped parameters on beam frequencies, *Aeronaut. Quart.* **14** (1963) 224–240.

18. Z. Zhang, C. M. Wang, N. Challamel and I. Elishakoff, Obtaining Eringen's length scale coefficient for vibrating nonlocal beams via continualization method, *J. Sound Vibr.* **333** (2014) 4977–4990.

19. J. L. Lagrange, *Mécanique Analytique* (Paris, 1788).

20. L. Rayleigh, *The Theory of Sound*, 2nd edn. (London, 1894).

21. F. Gantmacher, *Lectures in Analytical Mechanics* (Mir Publishers, Moscow, 1970).

22. I. K. Silverman, Discussion on the paper of "Salvadori M. G., Numerical computation of buckling loads by finite differences," *ASCE 116* (1951) 625–626.

23. H. Hencky, Über die angenäherte Lösung von stabilitätsproblemen im raum mittels der elastischen gelenkkette, *Der Eisenbau* **11** (1920) 437–452.

24. M. G. Salvadori, Numerical computation of buckling loads by finite differences, *Trans. ASCE 116* (1951) 590–636.

25. C. T. Wang, Discussion on the paper of Salvadori M. G., Numerical computation of buckling loads by finite differences, *ASCE 116* (1951) 629–631.

26. C. T. Wang, *Applied Elasticity* (McGraw-Hill, New York, 1953).

27. P. Seide, Accuracy of some numerical methods for column buckling, *J. Eng. Mech.*, *ASCE 101*(5) (1975) 549–560.

28. K. B. Subrahmanyam and K. R. V. Kaza, An improved finite-difference analysis of uncoupled vibrations of tapered cantilever beams, NASA Technical Memorandum 83495, September (1983), pp. 1–46.

29. K. B. Subrahmanyam and A. W. Leissa, An improved finite difference analysis of uncoupled vibrations of cantilevered beams, *J. Sound Vibr.* **98**(1) (1985) 1–11.

30. I. Elishakoff and R. Santoro, Error in the finite difference based probabilistic dynamic analysis: Analytical evaluation, *J. Sound Vibr.* **281**(3–5) (2005) 1195–1206.

31. R. Santoro and I. Elishakoff, Accuracy of the finite difference method in stochastic setting, *J. Sound Vibr.* **291**(1–2) (2006) 275–284.

32. N. Challamel, C. M. Wang and I. Elishakoff, Discrete systems behave as nonlocal structural elements: Bending, buckling and vibration analysis, *Eur. J. Mech. A-Solid.* **44** (2014) 125–135.
33. G. Strang and G. J. Fix, *An Analysis of the Finite Element Method* (Prentice-Hall, Englewood Cliffs, NJ, 1973).

34. T. J. R. Hughes, *The Finite Element Method: Linear Static and Dynamic Finite Element analysis* (Prentice Hall, Englewood Cliffs, 1987).

35. J. N. Reddy, *An Introduction to the Finite Element Method*, 2nd edn. (McGrawHill, New-York, 1992).

36. K. D. Hjelmstad, *Fundamental of Structural Mechanics*, 2nd Edn. (Springer, New York, 2005).

37. P. Wriggers, *Nonlinear Finite Element Methods* (Springer, New York, 2008).

38. P. Tong, T. H. H. Pian and L. L. Bucciarelli, Mode shapes and frequencies by finite element method using consistent and lumped masses, *Comput. Struct.* **1** (1971) 623–638.

39. T. Belytschko and W. L. Mindle, Flexural wave propagation behaviour of lumped mass approximation, *Comput. Struct.* **12** (1980) 805–812.

40. Y. M. Xie and G. P. Steven, Explicit formulas for correcting finite-element predictions of natural frequencies, *Comm. Num. Meth. Eng.* **9** (1993) 671–680.

41. A. C. Eringen, On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves, *J. Appl. Phys.* **54** (1983) 4703–4710.