# Stability of the Rao–Nakra sandwich beam with a dissipation of fractional
derivative type: theoretical and numerical study

K. Ammari  LR Analysis and Control of PDEs, LR 22ES03, Department of
Mathematics, Faculty of Sciences of Monastir, University of Monastir, 5019
Monastir, Tunisia<EMAIL_ADDRESS>,  V. Komornik  Département de
mathématique, Université de Strasbourg, 7 rue René Descartes, 67084 Strasbourg
Cedex, France<EMAIL_ADDRESS>,  M. Sepúlveda  DIM and CI2MA,
Universidad del Concepciòn, Concepciòn, Chile<EMAIL_ADDRESS>and  O. Vera
Departamento de Matemática, Universidad de Tarapacá, Av. 18 de Septiembre
2222, Arica, Chile<EMAIL_ADDRESS>
###### Abstract.

This paper is devoted to the solution and stability of a one-dimensional model
depicting Rao–Nakra sandwich beams, incorporating damping terms characterized
by fractional derivative types within the domain, specifically a generalized
Caputo derivative with exponential weight. To address existence, uniqueness,
stability, and numerical results, fractional derivatives are substituted by
diffusion equations relative to a new independent variable, $\xi$, resulting
in an augmented model with a dissipative semigroup operator. Polynomial decay
of energy is achieved, with a decay rate depending on the fractional
derivative parameters. Both the polynomial decay and its dependency on the
parameters of the generalized Caputo derivative are numerically validated. To
this end, an energy-conserving finite difference numerical scheme is employed.

###### Key words and phrases:

Rao–Nakra sandwich beam, fractional derivative, numerical polynomial stability

###### 2010 Mathematics Subject Classification:

35A01, 35A02, 35M33, 93D20, 34B05, 34D05, 34H05, 35B40, 74K20, 74F05

The author was supported by the following grants: NSFC No. 11871348, CAPES:
No. 88881.520205/2020-01, MATH AMSUD: 21-MATH-03.

The author was supported by the following grants: Fondecyt-ANID project
1220869, and ANID-Chile through Centro de Modelamiento Matemático (FB210005)

The author is partially financed by project Fondecyt 1191137 and UTA MAYOR
2022-2023, 4764-22.

###### Contents

  1. 1 Introduction
  2. 2 Augmented Model
  3. 3 Setting of the Semigroup
  4. 4 Non-uniform stabilization
    1. 4.1 Polynomial stability for $\eta=0.$
    2. 4.2 Polynomial stability for $\eta>0$
  5. 5 Numerical study
    1. 5.1 Linear equations of Motion
    2. 5.2 Time discretization
    3. 5.3 Approximation of fractional derivatives
      1. 5.3.1 Approximation using a truncation of the Grünwald-–Letnikov derivative.
      2. 5.3.2 Approximation using the Mbodje augmented model
    4. 5.4 Decay of the discrete energy using the Mbodje augmented model
    5. 5.5 Numerical examples

##  1\. Introduction

In 1999, Liu, Trogdon and Yong [28] developed the following general three-
layer laminated beam–plate model:

(1.6) $\displaystyle\left\\{\begin{array}[]{l}\varrho_{1}h_{1}u_{\bf
tt}-E_{1}h_{1}u_{xx}-\tau=0,\\\ \varrho_{3}h_{3}v_{\bf
tt}-E_{3}h_{3}u_{xx}+\tau=0,\\\ \varrho hw_{\bf
tt}+EIw_{xxxx}-G_{1}h_{1}(w_{x}+\phi_{1})_{x}-G_{3}h_{3}(w_{x}+\phi_{3})_{x}-h_{2}\tau_{x}=0,\\\
\varrho_{1}I_{1}\phi_{1{\bf
tt}}-E_{1}I_{1}\phi_{1xx}-\frac{h_{1}}{2}\tau+G_{1}h_{1}(w_{x}+\phi_{1})=0,\\\
\varrho_{3}I_{3}\phi_{3{\bf
tt}}-E_{3}I_{3}\phi_{3xx}-\frac{h_{3}}{2}\tau+G_{3}h_{3}(w_{x}+\phi_{3})=0.\end{array}\right.$

The physical parameters $h_{i},$ $\varrho_{i},$ $E_{i},$ $G_{i},$ $I_{i}$ $>0$
are the thickness, density, Young’s modulus, shear modulus and moments of
inertia of the $i$-th layer for $i=1,\,2,\,3,$ from bottom to top,
respectively. In addition, $\varrho h=\varrho_{1}h_{1}+\varrho_{3}h_{3}$ and
$EI=E_{1}I_{1}+E_{3}I_{3}.$

The Rao–Nakra sandwich beam model [45] is the following:

(1.13) $\displaystyle\left\\{\begin{array}[]{l}\rho_{1}u_{\bf tt}(x,\,{\bf
t})-\vartheta u_{xx}(x,\,{\bf t})-k(-u(x,\,{\bf t})+v(x,\,{\bf t})+\gamma
w_{x}(x,\,{\bf t}))=0,\\\ \rho_{2}v_{\bf tt}(x,\,{\bf t})-\chi v_{xx}(x,\,{\bf
t})+k(-u(x,\,{\bf t})+v(x,\,{\bf t})+\gamma w_{x}(x,\,{\bf t}))=0,\\\
\rho_{3}w_{\bf tt}(x,\,{\bf t})+\zeta w_{xxxx}(x,\,{\bf t})-k(-u(x,\,{\bf
t})+v(x,\,{\bf t})+\gamma w_{x}(x,\,{\bf t}))_{x}=0,\\\ u(x,\,0)=u_{0}(x),\
u_{\bf t}(x,\,0)=u_{1}(x),\\\ v(x,\,0)=v_{0}(x),\ v_{\bf
t}(x,\,0)=v_{1}(x),\\\ w(x,\,0)=w_{0}(x),\ w_{\bf
t}(x,\,0)=w_{1}(x),\end{array}\right.$

with $0<x<\ell$ and ${\bf t}>0.$ Here $u=u(x,\,{\bf t})$ and $v=v(x,\,{\bf
t})$ are the longitudinal displacement. $w=w(x,\,{\bf t})$ is the transverse
displacement of the beam. This system is derived from the system (1.6) when we
consider the core material to be linearly elastic, that is,
$\tau=2G_{2}\gamma$ with the shear strain

$\displaystyle\gamma=\frac{1}{2h_{2}}(-u+v+\alpha w_{x})\quad{\rm
and}\quad\alpha=h_{2}=\frac{1}{2}(h_{1}+h_{3}),$

where $k=\frac{G_{2}}{h_{2}},$ the shear modulus
$G_{2}=\frac{E_{2}}{2(1+\nu)},$ and $-1<\nu<\frac{1}{2}$ is the Poisson ratio.

Many authors have been working on the system (1.13) from different points of
view: see [27, 39, 40, 46] and the references therein. In the past three
decades there has been a growing interest by many reserchers for the study of
fractional calculus in different fields of sciences [32, 48, 49]. Several
aspects of engineering, applied sciences, and mathematical physics benefitted
from this ascending wave of applications. Space sciences, fluid mechanics,
porous media flows, viscoelastic and biological processes are but a few areas
in which fractional order differential equations have become a favored tool to
tread new path. To give some examples, fractional derivatives have been used
to model frequency dependent damping behavior of many viscoelastic materials
as well modeling many chemical processes. Many problems in several scientific
applied areas, including the analysis of viscoelastic materials, heat
conduction in materials with memory, electrodynamics with memory, signal
processing, among others, can be modeled by fractional differential calculus.
Indeed, many investigations have shown that models involving fractional
derivatives are more realistic to represent some natural phenomena than models
involving classical derivatives. For more information we refer to [24, 34, 47]
and the references therein. For example, the wave equation with boundary
fractional damping has been treated in [34, 35] where the strong stability and
the lack of uniform stabilization were proved. Applications of this kind are
also found in [33, 36, 42, 44, 50] and the references therein.

In [50] the authors studied the Rao–Nakra system with a boundary dissipation
of fractional derivative type. They established the polynomial stability of
the system. In this article we are interested in studying the polynomial
stabilization to the system (1.22) from a theoretical and numerical point of
view. In order to achieve this, we eliminate the dissipation given on the
border and place the dissipation in the domain. The main idea to reach this
result is to use an idea given in [4]. It is important to note that the
imposed dissipation is of a fractional derivative type. Our method is based on
the fact that the input–output relationship in a certain diffusion equation
introduced by Mbodje in [34] is realized by a fractional derivative operator.

Numerical approximations of derivatives and fractional integrals have been
extensively studied, developed, and refined. Comprehensive reviews on this
topic are available in [11, 26]. Many of these approaches involve direct
approximations of the Riemann–Liouville, Caputo [18, 31, 13, 14, 15], or
Grünwald–Letnikov integral [26], using quadrature methods, finite differences,
and truncated summations [19]. While some progress has been made in enhancing
the convergence rate, such as through trapezoidal approximation [26] or spline
interpolation [17], these methods do not inherently preserve the energy
conservation property of the Rao–Nakra model under study, nor do they maintain
stability in the presence of dissipative terms. In this work, we address this
challenge by employing a suitable combination of the $\beta$-Newmark method
[38], coupled with a Crank–Nicolson method for the Mbodje diffusion equation,
to achieve our objective.

Without lost of generality, we will consider the following system for
$\rho_{1}=\rho_{2}=\rho_{3}=1$:

(1.22) $\displaystyle\left\\{\begin{array}[]{l}u_{\bf tt}-\vartheta
u_{xx}-k(-u+v+\gamma w_{x})+\partial_{\bf
t}^{\alpha,\,\eta}u=0,\,(0,\ell)\times(0,+\infty),\\\ v_{\bf tt}-\chi
v_{xx}+k(-u+v+\gamma w_{x})+\partial_{\bf
t}^{\alpha,\,\eta}v=0,\,(0,\ell)\times(0,+\infty),\\\ w_{\bf tt}+\zeta
w_{xxxx}-k(-u+v+\gamma w_{x})_{x}+\partial_{\bf
t}^{\alpha,\,\eta}w=0,\,(0,\ell)\times(0,+\infty),\\\
u(x,\,0)=u_{0}(x),\;u_{\bf t}(x,\,0)=u_{1}(x),\,(0,\ell),\\\
v(x,\,0)=v_{0}(x),\;v_{\bf t}(x,\,0)=v_{1}(x),\,(0,\ell),\\\
w(x,\,0)=w_{0}(x),\;w_{\bf t}(x,\,0)=w_{1}(x),\,(0,\ell),\\\ u(0,\,{\bf
t})=u(\ell,\,{\bf t})=0,\ v(0,\,{\bf t})=v(\ell,\,{\bf t})=0,\,(0,+\infty),\\\
w(0,\,{\bf t})=w(\ell,\,{\bf t})=w_{x}(0,\,{\bf t})=w_{x}(\ell,\,{\bf
t})=0,\,(0,+\infty),\end{array}\right.$

with $0<x<\ell,$ ${\bf t}>0$, and real-valued functions $u=u(x,\,{\bf t})$ and
$v=v(x,\,{\bf t}),$ $w=w(x,\,{\bf t})$.

The rest of the paper is divided into five sections. In Section 1 we show that
the system (1.22) may be replaced by an augmented system (2.28) obtained by
coupling an equation with a suitable diffusion, and we study the energy
functional associated to system. In Section 3, we establish the existence and
uniqueness of solutions of the system (2.28); for this we use [5, 4]. In
section 4 we prove the strong stability, the lack of uniform stabilization and
the polynomial stability of the system (2.28). In Section 5, we give a
numerical study of the polynomial stability.

Throughout this paper, $C$ is a generic constant, not necessarily the same at
each occasion (it may change from line to line) and depending on the indicated
quantities.

##  2\. Augmented Model

This section is devoted to the reformulation of the model (1.22) with boundary
conditions (1.22)7, 8 as an augmented system. First we provide a brief review
of fractional calculus. There are many definitions of fractional derivatives
[18], among them those of Riemann–Liouville and Caputo are the most widely
used [31]. The latter has the same Laplace transform as the integer order
derivative, so it is widely used in control theory. In this paper the
fractional derivative damping force is regarded as a control force to study
the properties of free damped vibrations of the system, hence the Caputo
definition [13, 14, 15] is used here. Let $0<\alpha<1$. The Caputo fractional
integral of order $\alpha$ is defined by the formula

(2.1) $I^{\alpha}f({\bf t})=\frac{1}{\Gamma(\alpha)}\int_{0}^{\bf t}({\bf
t}-s)^{\alpha-1}f(s)ds,$

where $\Gamma$ is the Gamma function, $f\in L^{1}(0,+\infty)$ and ${\bf t}>0.$

The Caputo fractional derivative operator of order $\alpha$ is defined by the
expression

(2.2) $D^{\alpha}f({\bf t})=I^{1-\alpha}Df({\bf
t}):=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{\bf t}({\bf
t}-s)^{-\alpha}f^{\prime}(s)ds.$

The Caputo definition of fractional derivative possesses a very simple but
interesting interpretation: if the function $f({\bf t})$ represents the strain
history within a viscoelastic material whose relaxation function is
$[\Gamma(1-\alpha){\bf t}^{\alpha}]^{-1}$, then the material will experience
at any time ${\bf t}$ a total stress given by the expression $D^{\alpha}f({\bf
t}).$ Also, it easy to show that $D^{\alpha}$ is a left inverse of
$I^{\alpha},$ but in general it is not a right inverse. More precisely, we
have

$\displaystyle D^{\alpha}I^{\alpha}f=f,\qquad I^{\alpha}D^{\alpha}f({\bf
t})=f({\bf t})-f(0).$

See [47] for the proof of above equalities and for more properties of
fractional calculus.

In this work, we consider a slightly different version of (2.1) and (2.4). In
[16], Choi and MacCamy introduced the following definition of fractional
integro-differential operators with an exponential weight. Given $0<\alpha<1$
and $\eta\geq 0$, the exponential fractional integral of order $\alpha$ is
defined by

(2.3) $I^{\alpha,\ \eta}f({\bf t})=\frac{1}{\Gamma(\alpha)}\int_{0}^{\bf
t}e^{-\eta({\bf t}-s)}({\bf t}-s)^{\alpha-1}f(s)ds,$

and the exponential fractional derivative operator of order $\alpha$ is
defined by

(2.4) $\partial_{\bf t}^{\alpha,\ \eta}f({\bf
t})=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{\bf t}e^{-\eta\,({\bf t}-s)}({\bf
t}-s)^{-\alpha}f^{\prime}(s)ds.$

Note that $\partial_{\bf t}^{\alpha,\ \eta}f({\bf
t})=I^{1-\alpha,\,\eta}f^{\prime}({\bf t}).$

We will need later the following results:

######  Theorem 2.1.

[34] Let $\mu$ be the function

(2.5)
$\displaystyle\mu(\xi)=|\xi|^{(2\alpha-1)/2},\qquad\xi\in\mathbb{R},\qquad
0<\alpha<1.$

Then the relation between the Input $U$ and the Output $O$ is given by the
following system:

(2.6) $\displaystyle\varphi_{\bf t}({\bf t},\,\xi)+|\xi|^{2}\varphi({\bf
t},\,\xi)=\mu(\xi)U({\bf t}),\qquad\xi\in\mathbb{R},\qquad{\bf t}>0,$ (2.7)
$\displaystyle\varphi(0,\,\xi)=0,$ (2.8) $\displaystyle
O=\pi^{-1}\sin(\alpha\pi)\int_{\mathbb{R}}\mu(\xi)\varphi({\bf t},\,\xi)d\xi.$

This implies that

(2.9) $\displaystyle O=I^{1-\alpha}U,$

where $U\in C([0,\,+\infty)).$

######  Lemma 2.2.

[1, Lemma 2.1] If $\lambda\in D=\\{\lambda\in\mathbb{C}:\
Re\,\lambda+\eta>0\text{ or }Im\,\lambda\neq 0\\}$, then

$\displaystyle\int_{\mathbb{R}}\frac{\mu^{2}(\xi)}{\xi^{2}+\eta+\lambda}\,d\xi=\frac{\pi}{\sin(\alpha\pi)}(\eta+\lambda)^{\alpha-1}.$

Our strategy requires the elimination of the fractional derivatives in time
from the domain condition in the system (1.22). To this, setting
$\mu(\xi)=|\xi|^{(2\alpha-1)/2},$ $\xi\in\mathbb{R},$
$\mathfrak{C}=\pi^{-1}\sin(\alpha\,\pi),$ and exploiting the technique from
[21], we reformulate the system (1.22) by using Theorem 2.1. This leads to the
following augmented model, where
$(x,\,t,\,\xi)\in(0,\,\ell)\times(0,\,+\infty)\times\mathbb{R}$ and we
introduce the constant $\mathfrak{C}:=\pi^{-1}\sin(\alpha\,\pi)$:

(2.28) $\displaystyle\left\\{\begin{array}[]{l}u_{\bf tt}(x,\,{\bf
t})-\vartheta u_{xx}(x,\,{\bf t})-k(-u(x,\,{\bf t})+v(x,\,{\bf t})+\gamma
w_{x}(x,\,{\bf t}))\\\ +\
\mathfrak{C}\int_{\mathbb{R}}\mu(\xi)\varphi(x,\,{\bf t},\,\xi)d\xi=0,\\\
\varphi_{\bf t}(x,\,{\bf t},\,\xi)+(\xi^{2}+\eta)\varphi(x,\,{\bf
t},\,\xi)=\mu(\xi)u_{\bf t}(x,\,{\bf t}),\\\ \\\ v_{\bf tt}(x,\,{\bf t})-\chi
v_{xx}(x,\,{\bf t})+k(-u(x,\,{\bf t})+v(x,\,{\bf t})+\gamma w_{x}(x,\,{\bf
t}))\\\ +\ \mathfrak{C}\int_{\mathbb{R}}\mu(\xi)\phi(x,\,{\bf
t},\,\xi)d\xi=0,\\\ \phi_{\bf t}(x,\,{\bf
t},\,\xi)+(\xi^{2}+\eta)\phi(x,\,{\bf t},\,\xi)=\mu(\xi)v_{\bf t}(x,\,{\bf
t}),\\\ \\\ w_{\bf tt}(x,\,{\bf t})+\zeta w_{xxxx}(x,\,{\bf t})-k(-u(x,\,{\bf
t})+v(x,\,{\bf t})+\gamma w_{x}(x,\,{\bf t}))_{x}\\\ +\
\mathfrak{C}\int_{\mathbb{R}}\mu(\xi)\psi(x,\,{\bf t},\,\xi)d\xi=0,\\\
\psi_{\bf t}(x,\,{\bf t},\,\xi)+(\xi^{2}+\eta)\psi(x,\,{\bf
t},\,\xi)=\mu(\xi)w_{\bf t}(x,\,{\bf t}),\\\ \\\ u(x,\,0)=u_{0}(x),\;u_{\bf
t}(x,\,0)=u_{1}(x),\\\ v(x,\,0)=v_{0}(x),\;v_{\bf t}(x,\,0)=v_{1}(x),\\\
w(x,\,0)=w_{0}(x),\;w_{\bf t}(x,\,0)=w_{1}(x),\\\ u(0,\,t)=u(\ell,\,{\bf
t})=0,\ v(0,\,{\bf t})=v(\ell,\,{\bf t})=0,\\\ w(0,\,t)=w(\ell,\,{\bf
t})=w_{x}(0,\,{\bf t})=w_{x}(\ell,\,{\bf t})=0\\\ \varphi(x,0,\xi)=0,\
\phi(x,0,\xi)=0,\ \psi(x,0,\,\xi)=0.\end{array}\right.$

##  3\. Setting of the Semigroup

In this section, we discuss the existence and uniqueness of solutions for the
coupled system (2.28) by the semigroup theory [41, 43].

We will use the standard $L^{2}(0,\,\ell)$ space; the scalar product and the
norm are denoted by

$\langle\varphi,\,\psi\rangle_{L^{2}(0,\,\ell)}=\int_{0}^{\ell}\varphi\,\overline{\psi}\,dx\text{
and }\|\psi\|_{L^{2}(0,\,\ell)}^{2}=\int_{0}^{\ell}|\psi|^{2}dx.$

We introduce the Hilbert space

(3.1) $\displaystyle\mathcal{H}$ $\displaystyle=$ $\displaystyle
H_{0}^{1}(0,\,\ell)\times H_{0}^{1}(0,\,\ell)\times H_{0}^{2}(0,\,\ell)\times
L^{2}(0,\,\ell)\times L^{2}(0,\,\ell)\times L^{2}(0,\,\ell)$
$\displaystyle\qquad\times L^{2}(\mathbb{R},L^{2}(0,\,\ell))\times
L^{2}(\mathbb{R},L^{2}(0,\,\ell))\times L^{2}(\mathbb{R},L^{2}(0,\,\ell))$

equipped with the following inner product:

$\displaystyle\langle{\mathcal{U}},\,\tilde{{\mathcal{U}}}\rangle_{{\mathcal{H}}}=$
$\displaystyle\int_{0}^{\ell}U\overline{\tilde{U}}dx+\int_{0}^{\ell}V\overline{\tilde{V}}dx+\int_{0}^{\ell}W\overline{\tilde{W}}dx$
$\displaystyle+\
\vartheta\int_{0}^{\ell}u_{x}\overline{\tilde{u}}_{x}dx+\chi\int_{0}^{\ell}v_{x}\overline{\tilde{v}}_{x}dx+\zeta\int_{0}^{\ell}w_{xx}\overline{\tilde{w}}_{xx}dx$
$\displaystyle+\
k\int_{0}^{\ell}(-u+v+\gamma\,w_{x})(-\overline{\tilde{u}}+\overline{\tilde{v}}+\gamma\overline{\tilde{w}}_{x})dx$
(3.2) $\displaystyle+\
\mathfrak{C}\int_{\mathbb{R}}\langle\varphi,\,\tilde{\varphi}\rangle_{L^{2}(0,\,\ell)}d\xi+\mathfrak{C}\int_{\mathbb{R}}\langle\phi,\,\tilde{\phi}\rangle_{L^{2}(0,\,\ell)}d\xi+\mathfrak{C}\int_{\mathbb{R}}\langle\psi,\,\tilde{\psi}\rangle_{L^{2}(0,\,\ell)}d\xi,$

where ${\mathcal{U}}=(u,\,v,\,w,\,U,\,V,\,W,\,\varphi,\,\phi,\,\psi)^{T}$ and
$\
\widetilde{{\mathcal{U}}}=(\tilde{u},\,\tilde{v},\,\tilde{w},\,\tilde{U},\,\tilde{V},\,\tilde{W},\,\tilde{\varphi},\,\tilde{\phi},\,\widetilde{\psi})^{T}.$
We wish to transform the initial boundary value problem (2.28) into an
abstract Cauchy problem in the Hilbert space $\mathcal{H}.$ For this we
introduce the functions $u_{\bf t}=U,$ $v_{\bf t}=V,$ $w_{\bf t}=W,$ and we
rewrite the system (2.28) as the following initial value problem:

(3.3) $\displaystyle\frac{d{\mathcal{U}}}{d{\bf t}}({\bf
t})={\mathcal{A}}{\mathcal{U}}({\bf
t}),\quad{\mathcal{U}}(0)={\mathcal{U}}_{0},\quad\forall\;{\bf t}>0.$

Here ${\mathcal{U}}$ is as above,
${\mathcal{U}}_{0}=(u_{0},\,v_{0},\,w_{0},\,u_{1},\,v_{1},\,w_{0},\,0,\,0,\,0)^{T},$
and the operator
$\,\mathcal{A}:\mathcal{D}(\mathcal{A})\subset\mathcal{H}\rightarrow\mathcal{H}$
is given by the formula

(3.4) $\mathcal{A}\left(\begin{array}[]{c}u\\\ \\\ v\\\ \\\ w\\\ \\\ U\\\ \\\
V\\\ \\\ W\\\ \\\ \varphi\\\ \\\ \phi\\\ \\\
\psi\end{array}\right)=\left(\begin{array}[]{c}U\\\ \\\ V\\\ \\\ W\\\ \\\
\vartheta u_{xx}+k(-u+v+\gamma
w_{x})-\mathfrak{C}\int_{\mathbb{R}}\mu(\xi)\varphi(x,\xi)d\xi\\\ \\\ \chi
v_{xx}-k(-u+v+\gamma
w_{x})-\mathfrak{C}\int_{\mathbb{R}}\mu(\xi)\phi(x,\xi)d\xi\\\ \\\ -\zeta
w_{xxxx}+k(-u+v+\gamma
w_{x})_{x}-\mathfrak{C}\int_{\mathbb{R}}\mu(\xi)\psi(x,\xi)d\xi\\\ \\\
-(\xi^{2}+\eta)\varphi(x,\xi)+\mu(\xi)U(x)\\\ \\\
-(\xi^{2}+\eta)\phi(x,\xi)+\mu(\xi)V(x)\\\ \\\
-(\xi^{2}+\eta)\psi(x,\xi)+\mu(\xi)W(x)\end{array}\right).$

with the domain

$\displaystyle\mathcal{D}({\mathcal{A}})=$
$\displaystyle\left\\{{\mathcal{U}}=(u,\,v,\,w,\,U,\,V,\,W,\,\varphi,\,\phi,\,\psi)\in{\mathcal{H}}:\
U,\,V\in H_{0}^{1}(0,\,\ell),W\in H_{0}^{2}(0,\,\ell),\right.$
$\displaystyle\quad\vartheta u_{xx}+k(-u+v+\gamma
w_{x})-\mathfrak{C}\int_{\mathbb{R}}\mu(\xi)\varphi(x,\xi)d\xi\in
L^{2}(0,\,\ell),$ $\displaystyle\quad\chi v_{xx}-k(-u+v+\gamma
w_{x})-\mathfrak{C}\int_{\mathbb{R}}\mu(\xi)\phi(x,\xi)d\xi\in
L^{2}(0,\,\ell),$ $\displaystyle\quad-\zeta w_{xxxx}+k\gamma(-u+v+\gamma
w_{x})_{x}-\mathfrak{C}\int_{\mathbb{R}}\mu(\xi)\psi(x,\xi)d\xi\in
L^{2}(0,\,\ell),$
$\displaystyle\quad-(\xi^{2}+\eta)\varphi(x,\,\xi)+\mu(\xi)U(x)\in
L^{2}({\mathbb{R}},L^{2}(0,\ell))$
$\displaystyle\quad-(\xi^{2}+\eta)\phi(x,\xi)+\mu(\xi)V(x)\in
L^{2}({\mathbb{R}},L^{2}(0,\ell))$
$\displaystyle\quad-(\xi^{2}+\eta)\psi(x,\xi)+\mu(\xi)W(x)\in
L^{2}({\mathbb{R}},L^{2}(0,\ell))$ $\displaystyle\quad\
\quad\xi\varphi(x,\,\xi),\ \xi\phi(x,\,\xi),\ \xi\psi(x,\xi)\in
L^{2}({\mathbb{R}},L^{2}(0,\ell))\\}.$

Let the operator $\,A:\mathcal{D}(A)\subset H\rightarrow H$ which is given by

(3.5) $A\left(\begin{array}[]{c}u\\\ v\\\ w\\\
\end{array}\right)=\left(\begin{array}[]{c}-\vartheta u_{xx}-k(-u+v+\gamma
w_{x})\\\ -\chi v_{xx}+k(-u+v+\gamma w_{x})\\\ \zeta w_{xxxx}-k(-u+v+\gamma
w_{x})_{x}\end{array}\right),$

with the domain $\mathcal{D}(A)=\left[H^{2}(0,\ell)\cap
H^{1}_{0}(0,\ell)\right]\times\left[H^{2}(0,\ell)\cap
H^{1}_{0}(0,\ell)\right]\times\left[H^{4}(0,\ell)\cap
H_{0}^{2}(0,\,\ell)\right],$

$H=(L^{2}(0,\ell))^{3}$.

We note that the operator $A$ is self-adjoint and strictly positive.

So, the problem (1.22) can be rewritten as following, as in [4, 5]:

(3.6) $\left\\{\begin{array}[]{ll}Z_{\bf tt}({\bf t})+AZ({\bf
t})+BB^{*}\partial_{\bf t}^{\alpha,\eta}Z({\bf t})=0,\ {\bf t}>0\\\
Z(0)=Z_{0},Z_{\bf t}(0)=Z_{1},\end{array}\right.$

where $B=B^{*}=I_{H},Z({\bf t})=(u({\bf t}),\,v({\bf t}),\,w({\bf t}))^{T}$
and $Z_{0}=(u_{0},\,v_{0},\,w_{0})^{T},\quad
Z_{1}=(u_{1},\,v_{1},\,w_{1})^{T}.$

We define

(3.7) $\displaystyle\mathcal{H}_{0}$ $\displaystyle=$ $\displaystyle
H_{0}^{1}(0,\,\ell)\times H_{0}^{1}(0,\,\ell)\times H_{0}^{2}(0,\,\ell)\times
L^{2}(0,\,\ell)\times L^{2}(0,\,\ell)\times L^{2}(0,\,\ell)$

equipped with the inner product given by

$\displaystyle\langle{\mathcal{U}},\,\tilde{{\mathcal{U}}}\rangle_{{\mathcal{H}}_{0}}=$
$\displaystyle\int_{0}^{\ell}U\overline{\tilde{U}}dx+\int_{0}^{\ell}V\overline{\tilde{V}}dx+\int_{0}^{\ell}W\overline{\tilde{W}}dx$
$\displaystyle+\
\vartheta\int_{0}^{\ell}u_{x}\overline{\tilde{u}}_{x}dx+\chi\int_{0}^{\ell}v_{x}\overline{\tilde{v}}_{x}dx+\zeta\int_{0}^{\ell}w_{xx}\overline{\tilde{w}}_{xx}dx$
(3.8) $\displaystyle+\
k\int_{0}^{\ell}(-u+v+\gamma\,w_{x})(-\overline{\tilde{u}}+\overline{\tilde{v}}+\gamma\overline{\tilde{w}}_{x})dx,$

where ${\mathcal{U}}=(u,\,v,\,w,\,U,\,V,\,W)^{T}$ and $\
\tilde{{\mathcal{U}}}=(\tilde{u},\,\tilde{v},\,\tilde{w},\,\tilde{U},\,\tilde{V},\,\tilde{W})^{T}.$

Then, the operator
$\,\mathcal{A}_{0}:\mathcal{D}(\mathcal{A}_{0})\subset\mathcal{H}_{0}\rightarrow\mathcal{H}_{0}$
given by

(3.9) $\mathcal{A}_{0}\left(\begin{array}[]{c}u\\\ \\\ v\\\ \\\ w\\\ \\\ U\\\
\\\ V\\\ \\\ W\\\ \end{array}\right)=\left(\begin{array}[]{c}U\\\ \\\ V\\\ \\\
W\\\ \\\ \vartheta u_{xx}+k(-u+v+\gamma w_{x})-U\\\ \\\ \chi
v_{xx}-k(-u+v+\gamma w_{x})-V\\\ \\\ -\zeta w_{xxxx}+k(-u+v+\gamma
w_{x})_{x}-W\end{array}\right).$

with the domain

$\displaystyle\mathcal{D}({\mathcal{A}}_{0})=[H^{2}(0,\,\ell)\times
H_{0}^{1}(0,\,\ell)]^{2}\times[H^{4}(0,\,\ell)\cap H_{0}^{2}(0,\,\ell)]\times
H_{0}^{2}(0,\,\ell).$

We recall two propositions from [4, 5]:

######  Proposition 3.1.

[5, Proposition 2.6.2] The operator ${\mathcal{A}}_{0}$ generates a
C0-semigroup of contractions in the Hilbert space ${\mathcal{H}}_{0}.$
Moreover, the following auxiliary problem:

(3.18) $\displaystyle\left\\{\begin{array}[]{l}u_{\bf tt}({\bf t})-\vartheta
u_{xx}({\bf t})-k(-u({\bf t})+v({\bf t})+\gamma w_{x}({\bf t}))+u_{\bf t}({\bf
t})=0,\\\ v_{\bf tt}({\bf t})-\chi v_{xx}({\bf t})+k(-u({\bf t})+v({\bf
t})+\gamma w_{x}({\bf t}))+v_{\bf t}({\bf t})=0,\\\ w_{\bf tt}({\bf t})+\zeta
w_{xxxx}({\bf t})-k\gamma(-u({\bf t})+v({\bf t})+\gamma w_{x}({\bf
t}))_{x}+w_{\bf t}({\bf t})=0,\\\ u(0,\,{\bf t})=u(\ell,\,{\bf t})=0,\
v(0,\,{\bf t})=v(\ell,\,{\bf t})=0,\\\ w(0,\,{\bf t})=w(\ell,\,{\bf
t})=w_{x}(0,\,{\bf t})=w_{x}(\ell,\,{\bf t})=0,\\\ u(x,\,0)=u_{0}(x),\quad
u_{\bf t}(x,\,0)=u_{1}(x),\\\ v(x,\,0)=v_{0}(x),\quad v_{\bf
t}(x,\,0)=v_{1}(x),\\\ w(x,\,0)=w_{0}(x),\quad w_{\bf
t}(x,\,0)=w_{1}(x),\end{array}\right.$

admits a unique solution $(u(x,\,{\bf t}),\,v(x,\,{\bf t}),\,w(x,\,{\bf t}))$
such that if
$(u_{0},\,v_{0},\,w_{0},\,u_{1},\,v_{1},\,w_{1})\in{\mathcal{D}}({\mathcal{A}}_{0})$,
then the solution $(u(x,\,{\bf t}),\,v(x,\,{\bf t}),\,w(x,\,{\bf t}))$ of
(3.18) verifies the following regularity property:

$\displaystyle{\mathcal{U}}=(u,\,v,\,w,\,u_{\bf t},\,v_{\bf t},\,w_{\bf t})\in
C([0,\,+\infty):\,{\mathcal{D}}({\mathcal{A}}_{0}))\cap
C^{1}([0,\,+\infty),\,{\mathcal{H}}_{0}),$

and when
$(u_{0},\,v_{0},\,w_{0},\,u_{1},\,v_{1},\,w_{1})\in{\mathcal{H}}_{0}$, then

$\displaystyle{\mathcal{U}}=(u,\,v,\,w,\,u_{\bf t},\,v_{\bf t},\,w_{\bf t})\in
C([0,\,+\infty),\,{\mathcal{H}}_{0}).$

The energy of the solution of the system (3.18) is defined as follows:

$\displaystyle E_{0}({\bf t})=$ $\displaystyle\ \frac{1}{2}\left[\|u_{\bf
t}({\bf t})\|_{L^{2}(0,\,\ell)}^{2}+\|v_{\bf t}({\bf
t})\|_{L^{2}(0,\,\ell)}^{2}+\|w_{\bf t}({\bf
t})\|_{L^{2}(0,\,\ell)}^{2}\right.$ $\displaystyle\qquad+\ \vartheta\|u({\bf
t})\|_{H_{0}^{1}(0,\,\ell)}^{2}+\chi\|v({\bf
t})\|_{H_{0}^{1}(0,\,\ell)}^{2}+\zeta\|w({\bf t})\|_{H_{0}^{2}(0,\,\ell)}^{2}$
(3.19) $\displaystyle\qquad\left.+\ k\|-u({\bf t})+v({\bf t})+\gamma
w_{x}({\bf t})\|_{L^{2}(0,\,\ell)}^{2}\right],$

and it is decreasing function of the time ${\bf t}.$ In particular, we have

(3.20) $\displaystyle\frac{dE_{0}({\bf t})}{d{\bf t}}=-\|u_{\bf t}({\bf
t})\|_{L^{2}(0,\,\ell)}^{2}-\|v_{\bf t}({\bf
t})\|_{L^{2}(0,\,\ell)}^{2}-\|w_{\bf t}({\bf t})\|_{L^{2}(0,\,\ell)}^{2}.$

######  Proposition 3.2.

[5, Theorem 2.3.1] The operator ${\mathcal{A}}$ generates a C0-semigroup of
contractions in the Hilbert space ${\mathcal{H}}.$ Hence the problem (3.18)
admits a unique solution

$(u(x,\,{\bf t}),\,v(x,\,{\bf t}),\,w(x,\,{\bf t}),\varphi(x,\,{\bf
t},\,\xi),\,\phi(x,\,{\bf t},\,\xi),\,\psi(x,\,{\bf t},\,\xi))$

such that if
$(u_{0},\,v_{0},\,w_{0},\,u_{1},\,v_{1},\,w_{1},\,0,\,0,\,0)\in{\mathcal{D}}({\mathcal{A}})$,
then the solution $(u,\,v,\,w,\varphi,\,\phi,\,\psi)$ of (3.18) satisfies the
following regularity property:

$\displaystyle{\mathcal{U}}=(u,\,v,\,w,\,u_{\bf t},\,v_{\bf t},\,w_{\bf
t},\,\varphi,\,\phi,\,\psi)\in
C([0,\,+\infty),\,{\mathcal{D}}({\mathcal{A}}))\cap
C^{1}([0,\,+\infty),\,{\mathcal{H}}),$

and when
$(u_{0},\,v_{0},\,w_{0},\,u_{1},\,v_{1},\,w_{1},\,0,\,0,\,0)\in{\mathcal{H}}$,
then

$\displaystyle{\mathcal{U}}=(u,\,v,\,w,\,u_{\bf t},\,v_{\bf t},\,w_{\bf
t},\,\varphi,\,\phi,\,\psi)\in C([0,\,+\infty),\,{\mathcal{H}}).$

The energy of the system (3.18), defined by

(3.21) $\displaystyle E({\bf t}):=\frac{1}{2}\left\|(u,\,v,\,w,\,u_{\bf
t},\,v_{\bf t},\,w_{\bf
t},\,\varphi,\,\phi,\,\psi)\right\|_{\mathcal{H}}^{2},$

is decreasing in time ${\bf t}.$ In particular, we have

$\displaystyle\frac{dE({\bf t})}{d{\bf t}}=$
$\displaystyle-\mathfrak{C}\int_{\mathbb{R}}(|\xi|^{2}+\eta)\|\varphi({\bf
t,\,\xi})\|_{L^{2}(0,\,\ell)}^{2}d\xi-\mathfrak{C}\int_{\mathbb{R}}(|\xi|^{2}+\eta)\|\phi({\bf
t,\,\xi})\|_{L^{2}(0,\,\ell)}^{2}d\xi$ (3.22)
$\displaystyle-\mathfrak{C}\int_{\mathbb{R}}(|\xi|^{2}+\eta)\|\psi({\bf
t,\,\xi})\|_{L^{2}(0,\,\ell)}^{2}d\xi.$

##  4\. Non-uniform stabilization

###  4.1. Polynomial stability for $\eta=0.$

Assume that $\eta=0,$ then the operator ${\mathcal{A}}$ is not onto and
consequently $0\notin\varrho({\mathcal{A}})$ the resolvent set of
${\mathcal{A}}.$

Since the embedding $H_{0}^{1}(0,\,\ell)\times H_{0}^{1}(0,\,\ell)\times
H_{0}^{2}(0,\,\ell)\subset L^{2}(0,\,\ell)\times L^{2}(0,\,\ell)\times
L^{2}(0,\,\ell)$ is compact and the only solution of the following problem:

(4.10) $\displaystyle\left\\{\begin{array}[]{l}u_{\bf tt}-\vartheta
u_{xx}-k(-u+v+\gamma w_{x})=0,\quad(0,\,\ell)\times(0,\,+\infty),\\\ v_{\bf
tt}-\chi v_{xx}+k(-u+v+\gamma w_{x})=0,\quad(0,\,\ell)\times(0,\,+\infty),\\\
w_{\bf tt}+\zeta w_{xxxx}-k(-u+v+\gamma
w_{x})_{x}=0,\quad(0,\,\ell)\times(0,\,+\infty),\\\ u_{\bf t}(x,{\bf
t})=0,v_{\bf t}(x,{\bf t})=0,w_{\bf t}(x,{\bf t})=0,\
(0,\,\ell)\times(0,\,+\infty),\\\ u(0,\,{\bf t})=u(\ell,\,{\bf t})=0,\
v(0,\,{\bf t})=v(\ell,\,{\bf t})=0,\quad(0,\,+\infty),\\\ w(0,\,{\bf
t})=w(\ell,\,{\bf t})=w_{x}(0,\,{\bf t})=w_{x}(\ell,\,{\bf
t})=0,(0,\,+\infty),\\\ u(x,\,0)=u_{0}(x),\quad u_{\bf
t}(x,\,0)=u_{1}(x),\quad(0,\,\ell),\\\ v(x,\,0)=v_{0}(x),\quad v_{\bf
t}(x,\,0)=v_{1}(x),\quad(0,\,\ell),\\\ w(x,\,0)=w_{0}(x),\quad w_{\bf
t}(x,\,0)=w_{1}(x),\quad(0,\,\ell),\end{array}\right.$

is the trivial one.

Then according to [5, Sect. 2.4] the semigroup $(e^{{\bf
t}\mathcal{A}})_{t\geq 0}$ is strongly stable, i.e.,

$\|e^{{\bf t}{\mathcal{A}}}{\mathcal{U}}_{0}\|_{{\mathcal{H}}}\longrightarrow
0\quad{\rm as}\quad{\bf t}\rightarrow+\infty$

for all initial data ${\mathcal{U}}_{0}\in{\mathcal{H}}$.

To determine the rate of stability in this case, we need the following result
[7, Theorem 8.4]:

######  Theorem 4.1 ([7]).

Let $\mathcal{T}(t)$ be a bounded C0-semigroup on a Hilbert space $X$ with
generator $H.$ Assume that $\sigma(H)\cap i\mathbb{R}=\left\\{0\right\\}$ and
that there exist $\beta\geq 1,\gamma>0$ such that

$\left\|(i\lambda
I-H)^{-1}\right\|_{\mathcal{L}(X)}\leq\left\\{\begin{array}[]{ll}O(|\lambda|^{-\beta}),\,\lambda\rightarrow
0,\\\ O(|\lambda|^{\gamma}),\,\lambda\rightarrow\infty.\end{array}\right.$

Then, there exists a constant $C>0$ such that

$\left\|\mathcal{T}(t)z\right\|_{X}\leq\frac{C}{{\bf
t}^{\frac{1}{\max{(\beta,\gamma)}}}}\left\|z\right\|_{{\mathcal{D}}(H)},\;\forall\,t>0,z\in{\mathcal{D}}(H)\cap
R(H),$

where ${\mathcal{D}}(H)$ is the domain of $H$ and $R(H)$ is the range of $H.$

Based on Theorem 4.1, a simple adaptation of the proofs of [8, Theorems 5.8
and 5.9] (see also [4, 5]) lead to the following stability result:

######  Proposition 4.2.

If $\eta=0$, then there exists a $C>0$ such that

(4.11) $\left\|e^{{\bf
t}\mathcal{A}}U_{0}\right\|_{\mathcal{H}}\leq\frac{C}{\sqrt{{\bf
t}}}\left\|U_{0}\right\|_{\mathcal{D}(\mathcal{A})},\,\forall\,{\bf
t}>0,\,U_{0}\in{\mathcal{D}}(\mathcal{A})\cap R(\mathcal{A}),$

where $R(\mathcal{A})$ is the range of $\mathcal{A}.$

###  4.2. Polynomial stability for $\eta>0$

In this case there exists $\delta,\ C>0$, such that the auxiliary dissipative
operator satisfies the following relation:

$\left\|e^{{\bf t}\mathcal{A}_{0}}\right\|_{\mathcal{L}(\mathcal{H}_{0})}\leq
C\,e^{-\delta{\bf t}},\quad\forall\ {\bf t}\geq 0.$

Applying [5, Corollary 2.6.1] hence we obtain the following polynomial decay
result for the systems (1.22) and (2.28):

######  Proposition 4.3.

If $\eta>0$, then the semigroup $(e^{{\bf t}{\mathcal{A}}})_{t\geq 0}$ is
polynomially stable, namely there exists a constant $C>0$ such that

$\displaystyle\|e^{{\bf
t}{\mathcal{A}}}(u_{0},\,v_{0},\,w_{0},\,u_{1},\,v_{1},\,w_{1},\,\varphi_{0},\,\phi_{0},\,\psi_{0})\|_{{\mathcal{H}}}$
$\displaystyle\qquad\qquad\leq\frac{C}{(1+{\bf
t})^{\frac{1}{1-\alpha}}}\|(u_{0},\,v_{0},\,w_{0},\,u_{1},\,v_{1},\,w_{1},\,\varphi_{0},\,\phi_{0},\,\psi_{0})\|_{{\mathcal{D}}({\mathcal{A}})},\quad\forall\;{\bf
t}\geq 0,$

for all initial data
$(u_{0},\,v_{0},\,w_{0},\,u_{1},\,v_{1},\,w_{1},\,\varphi_{0},\,\phi_{0},\,\psi_{0})\in{\mathcal{D}}({\mathcal{A}}).$
In particular, the energy of the strong solution of (1.22) and (2.28)
satisfies the following estimate:

$\displaystyle E({\bf t})\leq\frac{C}{(1+{\bf
t})^{\frac{2}{1-\alpha}}}\|(u_{0},\,v_{0},\,w_{0},\,u_{1},\,v_{1},\,w_{1},\,0,\,0,\,0)\|_{{\mathcal{D}}({\mathcal{A}})}^{2},\quad\forall\,{\bf
t}>0.$

##  5\. Numerical study

In this section we will verify numerically the polynomial rate of decay
obtained in the previous section.

###  5.1. Linear equations of Motion

First, we approximate the longitudinal and transverse displacement vector
$[u,v,w]^{\top}$ in space using an energy-conservative finite difference
method. For $J\in\mathbb{N}$ and $\delta x=\ell/J$, we define $x_{j}$, with
$j=1,\ldots,J$, as a uniform discretization of the interval $(0,\ell)$,
obtaining a vector
$\mathbf{U}=[\mathbf{u}(t),\mathbf{v}(t),\mathbf{w}(t)]^{\top}$ approximation
of $[u,v,w]^{\top}$ in $\mathbb{R}^{3J}$. We obtain the linear equation of
motion

(5.1)
$\mathbf{M}\ddot{\mathbf{U}}(t)+\mathbf{K}{\mathbf{U}}(t)+\mathbf{C}\overset{\alpha,\eta}{\mathbf{U}}(t)=0$

where $\overset{\alpha,\eta}{\mathbf{U}}(t)=D^{\alpha,\eta}{\mathbf{U}}(t)$ is
the generalized Caputo fractional derivative defined in (2.2). In addition,
and according to a discretization consistent with (1.22), we choose
$\mathbf{M}=\mathbf{C}=\mathbf{I}_{3J}$ (identity matrix). The stiffness
matrix is given by
$\mathbf{K}=\mathbf{K}_{\mathrm{stress}}+\mathbf{K}_{\mathrm{coupling}}$,
where

$\mathbf{K}_{\mathrm{stress}}=\begin{pmatrix}-\vartheta\mathbf{D}^{2}&&\\\
&-\chi\mathbf{D}^{2}&\\\
&&\zeta\mathbf{D}^{4}\end{pmatrix},\qquad\mathbf{K}_{\mathrm{coupling}}=k\begin{pmatrix}\mathbf{I}_{J}&-\mathbf{I}_{J}&\gamma\mathbf{R}\\\
-\mathbf{I}_{J}&\mathbf{I}_{J}&-\gamma\mathbf{R}&\\\
\gamma\mathbf{R}^{T}&-\gamma\mathbf{R}^{T}&-\gamma^{2}\mathbf{D}^{2}.\end{pmatrix}$

Here, $\mathbf{D}^{2}$ and $\mathbf{D}^{4}$ are the finite difference matrix
approximation of $\partial_{xx}$ and $\partial_{xxxx}$ respectively with the
boundary conditions (1.22)7,8, and $\mathbf{R}$ is the upper triangular
matrix, obtained from the Cholesky decomposition
$\mathbf{D}^{2}=-\mathbf{R}^{T}\mathbf{R}$. Using uniform finite differences
we have

$\displaystyle\left(\mathbf{D}^{2}\mathbf{u}\right)_{j}$ $\displaystyle=$
$\displaystyle\dfrac{\mathbf{u}_{j-1}-2\mathbf{u}_{j}+\mathbf{u}_{j+1}}{\delta
x^{2}},$ $\displaystyle\left(\mathbf{D}^{4}\mathbf{w}\right)_{j}$
$\displaystyle=$
$\displaystyle\dfrac{\mathbf{w}_{j-2}-4\mathbf{w}_{j-1}+6\mathbf{w}_{j}-4\mathbf{w}_{j+1}+\mathbf{w}_{j+2}}{\delta
x^{4}},\qquad J=1,\ldots,J.$

According to the boundary conditions (1.22)7,8, we seek the solution of Finite
Difference in

$X_{J}:=\left\\{\mathbf{U}=[\mathbf{u}(t),\mathbf{v}(t),\mathbf{w}(t)]^{\top}\in\mathbb{R}^{3J}\right\\},$

considering
$\mathbf{u}_{0}=\mathbf{u}_{J+1}=\mathbf{v}_{0}=\mathbf{v}_{J+1}=\mathbf{w}_{0}=\mathbf{w}_{J+1}=\mathbf{v}_{-1}=\mathbf{v}_{J+2}=0$.

###  5.2. Time discretization

In order to preserve the energy with a second order scheme in time, we choose
a $\beta$-Newmark scheme for $w$. The method consists of updating the
displacement, velocity and acceleration vectors at the current time
$t^{n}=n\delta t$ to the time $t^{n+1}=(n+1)\delta t$, a small time interval
$\delta t$ later. The Newmark algorithm [38] is based on a set of two
relations expressing the forward displacement $\mathbf{U}^{n+1}$ and velocity
$\dot{\mathbf{U}}^{n+1}$ in terms of their current values and the forward and
current values of the acceleration:

(5.2) $\displaystyle\dot{\mathbf{U}}^{n+1}$ $\displaystyle=$
$\displaystyle\dot{\mathbf{U}}^{n}+(1-\widetilde{\gamma})\delta
t\,\ddot{\mathbf{U}}^{n}+\widetilde{\gamma}\delta t\,\ddot{\mathbf{U}}^{n+1},$
(5.3) $\displaystyle{\mathbf{U}}^{n+1}$ $\displaystyle=$
$\displaystyle{\mathbf{U}}^{n}+\delta
t\dot{\mathbf{U}}^{n}+\left(\frac{1}{2}-\widetilde{\beta}\right)\delta
t^{2}\,\ddot{\mathbf{U}}^{n}+\widetilde{\beta}\delta
t^{2}\,\ddot{\mathbf{U}}^{n+1},$

where $\widetilde{\beta}$ and $\widetilde{\gamma}$ are parameters of the
methods that will be fixed later. Replacing (5.2)–(5.3) in the equation of
motion (5.1), we obtain

(5.4) $\left(\mathbf{M}+\widetilde{\beta}\delta
t^{2}\,\mathbf{K}\right)\ddot{\mathbf{U}}^{n+1}+\mathbf{C}\mathord{\hbox{{\vfil#\vfil\cr\hbox{$\overset{\alpha,\eta}{\mathbf{U}}$}\cr\hbox{$\scriptstyle
n+1$}\cr}\kern 0.5pt}}=-\mathbf{K}\left(\mathbf{U}^{n}+\delta
t\dot{\mathbf{U}}^{n}+\left(\dfrac{1}{2}-\widetilde{\beta}\right)\delta
t^{2}\ddot{\mathbf{U}}^{n}\right).$

###  5.3. Approximation of fractional derivatives

At this state we have two possibilities: (1) approximate the fractional
derivative using a classical finite difference method such as the truncation
of the Grünwald–-Letnikov derivative [19], which in turn is an equivalent form
of the Riemann–Liouville derivative; or, (2) propose a conservative numerical
scheme for the Mbodje augmented model (2.28).

####  5.3.1. Approximation using a truncation of the Grünwald-–Letnikov
derivative.

In the case of the classical Grünwald–Letnikov derivative we consider a
fractional trapezoidal formula [26], in order to approximate

(5.5)
$\left[\prescript{}{C}{D}^{\alpha,\eta}_{0,t}f\right]_{t=t_{n}}\approx\sum_{k=0}^{n}a_{k,n}e^{-\eta\delta
t(n-k)}f^{\prime}(t_{k}),$

where

(5.6) $a_{k,n}=\dfrac{\delta
t^{1-\alpha}}{\Gamma(3-\alpha)}\begin{cases}(n-1)^{2-\alpha}-(n-2+\alpha)n^{1-\alpha},&k=0,\\\
(n-k+1)^{2-\alpha}+(n-1-k)^{2-\alpha}-2(n-k)^{2-\alpha},&1\leqslant k\leqslant
n-1,\\\ 1,&k=n.\end{cases}$

Replacing the term
$\mathord{\hbox{{\vfil#\vfil\cr\hbox{$\overset{\alpha,\eta}{\mathbf{U}}$}\cr\hbox{$\scriptstyle
n+1$}\cr}\kern 0.5pt}}$ by the numerical approximation
$\left[\prescript{}{C}{D}^{\alpha,\eta}_{0,t}f\right]_{t=t_{n}}$ given by
(5.5)–(5.6), then the equation (5.4) becomes

(5.7) $\displaystyle\left(\mathbf{M}+\widetilde{\beta}\delta
t^{2}\,\mathbf{K}+\dfrac{\widetilde{\gamma}\delta
t^{2-\alpha}}{\Gamma(3-\alpha)}\,\mathbf{C}\right)\ddot{\mathbf{U}}^{n+1}=-\mathbf{K}\left(\mathbf{U}^{n}+\delta
t\dot{\mathbf{U}}^{n}+\left(\dfrac{1}{2}-\widetilde{\beta}\right)\delta
t^{2}\ddot{\mathbf{U}}^{n}\right)$
$\displaystyle\qquad\qquad\qquad\qquad-\dfrac{\delta
t^{1-\alpha}}{\Gamma(3-\alpha)}\,\mathbf{C}\left(\dot{\mathbf{U}}^{n}+\left(1-\widetilde{\gamma}\right)\delta
t\ddot{\mathbf{U}}^{n}+\sum_{k=0}^{n}a_{k,n+1}\dot{\mathbf{U}}^{k}\right).$

####  5.3.2. Approximation using the Mbodje augmented model

On the other hand, in case of the Mbodje augmented model (2.28) we propose the
following conservative numerical scheme:

(5.8) $\begin{cases}\left(\mathbf{M}+\widetilde{\beta}\delta
t^{2}\,\mathbf{K}\right)\ddot{\mathbf{U}}^{n+1}+\mathfrak{C}\displaystyle\sum_{\ell=1}^{M}\mu_{\ell}\mathbf{\Phi}^{n+1}_{\ell}=-\mathbf{K}\left(\mathbf{U}^{n}+\delta
t\dot{\mathbf{U}}^{n}+\left(\dfrac{1}{2}-\widetilde{\beta}\right)\delta
t^{2}\ddot{\mathbf{U}}^{n}\right),\\\
\mathbf{\Phi}^{n+1}_{\ell}=\mathbf{\Phi}^{n}_{\ell}-\delta
t\left(\xi_{\ell}^{2}+\eta\right)\mathbf{\Phi}^{n+\frac{1}{2}}_{\ell}+\delta
t\mu_{\ell}\dot{\mathbf{U}}^{n+\frac{1}{2}},\end{cases}$

where
$\mathbf{\Phi}^{n+\frac{1}{2}}_{\ell}=\dfrac{\mathbf{\Phi}^{n}_{\ell}+\mathbf{\Phi}^{n+1}_{\ell}}{2}$,
with $\mathbf{\Phi}^{n}_{\ell}$ denoting the approximation of
$[\phi,\varphi,\psi]^{\top}\in\mathbb{R}^{3J}$ evaluated in
$\xi_{\ell}:=\ell\delta\xi$ and $t=n\delta t$, for $\ell=1,\ldots,M$,
$n=1,\ldots,N$ and a fixed $\delta\xi>0$. Using (5.2) again and replacing in
(5.8)2, we can rewrite (5.8) in the following more explicit and computable
way:

(5.9) $\begin{cases}\left(\mathbf{M}+\widetilde{\gamma}\delta
t\,\mathbf{C}_{augm}+\widetilde{\beta}\delta
t^{2}\,\mathbf{K}\right)\ddot{\mathbf{U}}^{n+1}=-\mathfrak{C}\displaystyle\sum_{\ell=1}^{M}\widetilde{\mu}_{\ell}\mathbf{\Phi}^{n}_{\ell}-\mathbf{C}_{augm}\left(\dot{\mathbf{U}}^{n}+\left(1-\widetilde{\gamma}\right)\delta
t\ddot{\mathbf{U}}^{n}\right)\\\
\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\mathbf{K}\left(\mathbf{U}^{n}+\delta
t\dot{\mathbf{U}}^{n}+\left(\dfrac{1}{2}-\widetilde{\beta}\right)\delta
t^{2}\ddot{\mathbf{U}}^{n}\right),\\\
\mathbf{\Phi}^{n+1}_{\ell}=\dfrac{2-\delta
t\left(\xi_{\ell}^{2}+\eta\right)}{2+\delta
t\left(\xi_{\ell}^{2}+\eta\right)}\mathbf{\Phi}^{n}_{\ell}+\dfrac{2\delta
t\mu_{\ell}}{2+\delta
t\left(\xi_{\ell}^{2}+\eta\right)}\dot{\mathbf{U}}^{n+\frac{1}{2}},\end{cases}$

where $\mu_{\ell}=\xi_{\ell}^{(2\alpha-1)/2}$,
$\widetilde{\mu}_{\ell}=\dfrac{2-\delta
t\left(\xi_{\ell}^{2}+\eta\right)}{2+\delta
t\left(\xi_{\ell}^{2}+\eta\right)}\mu_{\ell}$ and $\mathbf{C}_{augm}=\delta
t\mathfrak{C}\left(\displaystyle\sum_{\ell=1}^{M}\dfrac{2\widetilde{\mu}_{\ell}^{2}}{2+\delta
t\left(\xi_{\ell}^{2}+\eta\right)}\right)\textbf{I}_{3J}$. Although the two
approximations are consistently valid from the point of view of finite
differences, we will opt for the second approximation because it conserves the
energy in the non-dissipative case, and its decrease is consistent with the
continuous case. It also involves the calculation of $\mathbf{\Phi}$ which is
inherent to the definition of energy itself.

###  5.4. Decay of the discrete energy using the Mbodje augmented model

Evaluating (5.8)1 in $t=t_{n+\frac{1}{2}}$, multiplying by
$\ddot{\mathbf{U}}^{n+\frac{1}{2}}$, and summing (5.8)2 multiplied by
$\mathfrak{C}\mathbf{\Phi}^{n+\frac{1}{2}}_{\ell}$, we obtain that

(5.10)
$\displaystyle\left[E^{\text{Mbodje}}_{\Delta}\right]_{n}^{n+1}:=\left[\dfrac{1}{2}\dot{\mathbf{U}}^{T}\mathbf{M}\dot{\mathbf{U}}+\dfrac{1}{2}{\mathbf{U}}^{T}\mathbf{K}{\mathbf{U}}+\dfrac{\mathfrak{C}}{2}\displaystyle\sum_{\ell=1}^{M}{\mu}_{\ell}\left|\mathbf{\Phi}_{\ell}\right|^{2}\right]_{n}^{n+1}$
$\displaystyle=$
$\displaystyle\Bigg{[}-\dfrac{\vartheta}{2}\dot{\mathbf{u}}^{T}\mathbf{D}^{2}\dot{\mathbf{u}}-\dfrac{\chi}{2}\dot{\mathbf{v}}^{T}\mathbf{D}^{2}\dot{\mathbf{v}}+\dfrac{\zeta}{2}\dot{\mathbf{w}}^{T}\mathbf{D}^{4}\dot{\mathbf{w}}+\dfrac{k}{2}\|-\mathbf{u}+\mathbf{v}+\gamma\mathbf{R}\mathbf{w}\|^{2}$
$\displaystyle+\
\dfrac{\mathfrak{C}}{2}\displaystyle\sum_{\ell=1}^{M}{\mu}_{\ell}\left(\left|{\phi}_{\ell}\right|^{2}+\left|{\varphi}_{\ell}\right|^{2}+\left|{\psi}_{\ell}\right|^{2}\right)\Bigg{]}_{n}^{n+1}$
$\displaystyle=$ $\displaystyle-\
\mathfrak{C}\displaystyle\sum_{\ell=1}^{M}\left(\xi_{\ell}^{2}+\eta\right)\left(\left|{\phi}^{n+\frac{1}{2}}_{\ell}\right|^{2}+\left|{\varphi}^{n+\frac{1}{2}}_{\ell}\right|^{2}+\left|{\psi}^{n+\frac{1}{2}}_{\ell}\right|^{2}\right)$

which is consistent with the estimates (3.21)–(3), and constitutes a correct
approximation of the energy and its decreasing behavior.

######  Remark 5.1.

In the case of an approximation using a truncation of theGrünwald–-Letnikov
derivative, we have the following estimate of the energy. Evaluating (5.1) in
$t=t_{n+\frac{1}{2}}$ and multiplying by $\ddot{\mathbf{U}}^{n+\frac{1}{2}}$
we obtain

(5.11)
$\displaystyle\left[E^{GL}_{\Delta}\right]_{n}^{n+1}:=\left[\dfrac{1}{2}\dot{\mathbf{U}}^{T}\mathbf{M}\dot{\mathbf{U}}+\dfrac{1}{2}{\mathbf{U}}^{T}\mathbf{K}{\mathbf{U}}\right]_{n}^{n+1}$
$\displaystyle=\left[-\dfrac{\vartheta}{2}\dot{\mathbf{u}}^{T}\mathbf{D}^{2}\dot{\mathbf{u}}-\dfrac{\chi}{2}\dot{\mathbf{v}}^{T}\mathbf{D}^{2}\dot{\mathbf{v}}+\dfrac{\zeta}{2}\dot{\mathbf{w}}^{T}\mathbf{D}^{4}\dot{\mathbf{w}}+\dfrac{k}{2}\|-\mathbf{u}+\mathbf{v}+\gamma\mathbf{R}\mathbf{w}\|^{2}\right]_{n}^{n+1}$
$\displaystyle=-\frac{\delta
t^{2-\alpha}}{\Gamma(3-\alpha)}\left(\dot{\mathbf{U}}^{n+\frac{1}{2}}\right)^{T}\mathbf{C}\dot{\mathbf{U}}^{n+\frac{1}{2}}-(2^{1-\alpha}-1)\frac{\delta
t^{2-\alpha}}{\Gamma(3-\alpha)}\left(\dot{\mathbf{U}}^{n+\frac{1}{2}}\right)^{T}\mathbf{C}\dot{\mathbf{U}}^{n}$
$\displaystyle~{}~{}~{}-\sum_{k=0}^{n-1}a_{k,n+\frac{1}{2}}\left(\dot{\mathbf{U}}^{n+\frac{1}{2}}\right)^{T}\mathbf{C}\dot{\mathbf{U}}^{k}$

where $a_{k,n+\frac{1}{2}}=\dfrac{a_{k,n+1}+a_{k,n}}{2}$. The first term on
the right hand side of (5.11) is obviously strictly negative. However, the
sign of the second and third terms in (5.11) cannot be guaranteed. In fact, it
is verified numerically that the energy discretized in this way is not
necessarily strictly decreasing, unless the solution is sufficiently smooth
and delta t is small enough.

Figure 1. Comparison between the Grünwald–Letnikov derivative approximation
and the discretized Mbodje augmented model. At the top left:
$E^{GL}_{\Delta}(t)=\dfrac{1}{2}\dot{\mathbf{U}}^{T}\mathbf{M}\dot{\mathbf{U}}+\dfrac{1}{2}{\mathbf{U}}^{T}\mathbf{K}{\mathbf{U}}$;
At the top right:
$E^{GL,\Phi}_{\Delta}(t)=E^{GL}_{\Delta}(t)+\dfrac{\mathfrak{C}}{2}\displaystyle\sum_{\ell=1}^{M}{\mu}_{\ell}\left|\mathbf{\Phi}_{\ell}\right|^{2}$;
At the bottom: $E^{Mbodje}_{\Delta}(t)$ defined in (5.10).

Figure 1 shows a comparison of both models. It is proved that the calculation
of energy using the classical Grünwald–-Letnikov derivatives fails, either by
using the derivative defined in (5.11), (in top left graph), as well as by
artificially adding the calculation of $\mathbf{\Phi}$ as seen in the top
right graph. Therefore, the only acceptable method for all our calculations
should be the discretization of the augmented Mbodje model (bottom graph).

###  5.5. Numerical examples

We consider here the value of the parameters
$\vartheta=\chi=\zeta=k=\gamma=1$, a beam of length $\ell=1$, and the initial
conditions

(5.12) $\begin{cases}\displaystyle
u(x,0)=x\left(\frac{\ell^{3}}{8}-\left|x-\frac{\ell}{2}\right|^{3}\right),&u_{t}(x,0)=0,\\\
\displaystyle v(x,0)=0,&v_{t}(x,0)=0,\\\
w(x,0)=x\left(x-\frac{\ell}{2}\right)\left(\displaystyle\frac{\ell^{3}}{8}-\left|x-\frac{\ell}{2}\right|^{3}\right),&w_{t}(x,0)=0.\end{cases}$

We observe that $u(0,t)=u(\ell,t)=0$ and $u(\cdot,0)$ is of class
$\mathcal{C}^{1}(0,\ell)$ but its second derivative has a discontinuity in
$x=\ell/2$, and in turn, $w(0,t)=w(\ell,t)=0$ with $w(\cdot,0)$ of class
$\mathcal{C}^{3}(0,\ell)$ and its fourth derivative has a discontinuity at
$x=\ell/2$. Therefore, we have just $\mathcal{U}_{0}\in\mathcal{D}(A)$, but no
greater regularity than that.

In Figure 2, we consider a fractional derivative with $\alpha=0.5$, and we do
a simulation of the numerical scheme (5.2), (5.3), (5.7) with $J=500$, and
$\delta=T/N_{t}$ with $T=100$ and $N_{t}=1000$. The graphs in (A) of Figure 2
shows the state of the instant of the deformations for $t=10$. The beam is
represented by the thick yellow graph and its deformation corresponds to the
transversal displacements. The longitudinal displacements of $u(x,10)$ and
$v(x,10)$, are represented by the blue and red lines respectively. The graph
(B) of the figure 2 shows the evolution of the transverse displacements $w$ in
term of space and time. Graphs (C) and (D) of the same figure show the
longitudinal displacements $u$ and $v$ (respectively), also as a function of
space and time. In this set of figures, the decay of the quantities, and
therefore of the energy, is visually observed due to the dispative terms of
the fractional derivative type.

(a) State of the deformations at instant $t=10$

(b) Transversal deformations $w$

(c) Longitudinal deformations $u$

(d) Longitudinal deformations $v$

Figure 2. Simulation of the displacements. On the left: the displacements at
$t=2$. On the right: the evolution in time of the transverse displacements.

Figure 3. Behavior of the auxiliary variables $\phi$, $\varphi$ and $\psi$ of
the augmented model (2.28)

On the other hand, we show in figure 3 the graphs of the numerical variables
that approximate the auxiliary functions $\phi(x,t;\xi)$, $\varphi(x,t;\xi)$,
and $\psi(x,t;\xi)$ of the augmented model (2.28), for the instant $t=2$ of
the simulations with the same parameters and data as the example in figure 2.
We numerically observe a decay of these functions with respect to the variable
$\xi$, although not so evident in the case of the variable $\psi$ associated
with the term of the fractional derivative of the transverse displacement
(when truncating at $\xi=10$). In this sense, we must be careful with the
numerical truncation of the variable $\xi\in(0,+\infty)$. For this reason we
will consider a truncation at $\xi=10,000$ for the examples that will follow,
and that will show the polynomial decay of energy for different values of
$\alpha$ and $\eta$.

(a) Energies for different parameters $\alpha$ and $\eta$
regrouped in terms of the values of $\alpha$.

(b) Zoom: Energies for $\alpha=0.25$,
and different values of $\eta$

(c) Comparison between energy curves
with a same $\alpha=0.25$ as a function of time,
and their asymptotes $C(t+1)^{p}$

(d) Comparison between energy curves
with a same $\eta=0.3$ as a function of time,
and their asymptotes $C(t+1)^{p}$

Figure 4.  Energy with different values of $\eta$ and $\alpha$ in loglog scale.  |  | $\alpha$
---|---|---
$C(t+1)^{p}$ | 0.01 | 0.25 | 0.50 | 0.75 | 0.99
$\eta$ | $0$ |  $\begin{array}[]{l}C=\ 0.0300\\\ p=\ -0.6465\end{array}$ |  $\begin{array}[]{l}0.0250\\\ -0.4136\end{array}$ |  $\begin{array}[]{l}0.0600\\\ -0.7497\end{array}$ |  $\begin{array}[]{l}0.0186\\\ -0.9147\end{array}$ |  $\begin{array}[]{l}0.1521\\\ -0.7968\end{array}$
$1.10^{-4}$ |  $\begin{array}[]{l}C=\ 0.0869\\\ p=\ -0.8136\end{array}$ |  $\begin{array}[]{l}0.0408\\\ -0.8630\end{array}$ |  $\begin{array}[]{l}0.0250\\\ -0.9773\end{array}$ |  $\begin{array}[]{l}0.0198\\\ -0.9255\end{array}$ |  $\begin{array}[]{l}0.1521\\\ -0.7968\end{array}$
$3.10^{-4}$ |  $\begin{array}[]{l}C=\ 0.1089\\\ p=\ -0.8600\end{array}$ |  $\begin{array}[]{l}0.4175\\\ -1.2750\end{array}$ |  $\begin{array}[]{l}0.0338\\\ -1.0469\end{array}$ |  $\begin{array}[]{l}0.0198\\\ 0.9266\end{array}$ |  $\begin{array}[]{l}0.1521\\\ -0.7968\end{array}$
$1.10^{-3}$ |  $\begin{array}[]{l}C=\ 0.0894\\\ p=\ -0.8413\end{array}$ |  $\begin{array}[]{l}0.0832\\\ -1.1394\end{array}$ |  $\begin{array}[]{l}0.0338\\\ -0.9985\end{array}$ |  $\begin{array}[]{l}0.0198\\\ -0.9236\end{array}$ |  $\begin{array}[]{l}0.1521\\\ -0.7968\end{array}$
$3.10^{-3}$ |  $\begin{array}[]{l}C=\ 0.0822\\\ p=\ -0.8317\end{array}$ |  $\begin{array}[]{l}0.0258\\\ -1.0075\end{array}$ |  $\begin{array}[]{l}0.0159\\\ -0.9688\end{array}$ |  $\begin{array}[]{l}0.0190\\\ -0.9217\end{array}$ |  $\begin{array}[]{l}0.1522\\\ -0.7968\end{array}$
$1.10^{-2}$ |  $\begin{array}[]{l}C=\ 0.0803\\\ p=\ -0.8290\end{array}$ |  $\begin{array}[]{l}0.0170\\\ -0.9586\end{array}$ |  $\begin{array}[]{l}0.0144\\\ -0.9569\end{array}$ |  $\begin{array}[]{l}0.0188\\\ -0.9207\end{array}$ |  $\begin{array}[]{l}0.1525\\\ -0.7967\end{array}$
$3.10^{-2}$ |  $\begin{array}[]{l}C=\ 0.0801\\\ p=\ -0.8283\end{array}$ |  $\begin{array}[]{l}0.0155\\\ -0.9469\end{array}$ |  $\begin{array}[]{l}0.0141\\\ -0.9534\end{array}$ |  $\begin{array}[]{l}0.0188\\\ -0.9201\end{array}$ |  $\begin{array}[]{l}0.1533\\\ -0.7964\end{array}$
$1.10^{-1}$ |  $\begin{array}[]{l}C=\ 0.0806\\\ p=\ -0.8279\end{array}$ |  $\begin{array}[]{l}0.0152\\\ -0.9431\end{array}$ |  $\begin{array}[]{l}0.0140\\\ -0.9514\end{array}$ |  $\begin{array}[]{l}0.0190\\\ -0.9190\end{array}$ |  $\begin{array}[]{l}0.1562\\\ -0.7955\end{array}$
$3.10^{-1}$ |  $\begin{array}[]{l}C=\ 0.0823\\\ p=\ -0.8269\end{array}$ |  $\begin{array}[]{l}0.0154\\\ -0.9404\end{array}$ |  $\begin{array}[]{l}0.0143\\\ -0.9485\end{array}$ |  $\begin{array}[]{l}0.0196\\\ -0.9161\end{array}$ |  $\begin{array}[]{l}0.1653\\\ -0.7925\end{array}$
Table 1. Constant $C$ and rate $p$ of polynomial decay $C(t+1)^{p}$ for the
energy by a least square method, for different values of the fractional
derivative parameters $\alpha$ and $\eta$.

Finally, in order to numerically appreciate the polynomial decay of the
energy, we do the simulation with the same parameters
$\vartheta=\chi=\zeta=k=\gamma=1$, a beam of length $\ell=1$, and the initial
conditions (5.12), but for a longer time $T=10000$, and a time step $\delta
t=0.1$. Likewise, in order to obtain better precision in the results, we
truncate the variable $\xi$ of the auxiliary functions at $\xi=10000$.

In Figure 4, various plots depict energy on a log-log scale, where curves of
the form $C(t+1)^{p}$ exhibit nearly linear behavior over sufficiently long
times.

The graph in Figure 4(A) illustrates energy curves for 9 different values of
$\eta$ and 5 values of $\alpha$, totaling 45 combinations. Similar curves with
the same $\alpha$ are grouped into 5 distinct colors. While there are
variations in slopes $p$ and constants $C$, graphical observation reveals a
close relationship between these parameters, particularly in their bounded
ranges. The numerical estimation of $C$ and $p$ using a least squares method
is presented in Table 1.

Additionally, from the graph, it is conjectured that energy decays more
rapidly to zero for $\alpha=0.5$ (or similar values) and certain $\eta\neq 0$.
Conversely, $\alpha=0.25$ exhibits intriguing behavior, accentuating
differences across various $\eta$ values. Thus, the curves corresponding to
$\alpha=0.25$ are separately plotted in Figure 4(B). Here, a pattern emerges
with curves lying between two quasi-linear lines with distinct slopes. The
steepest negative slope corresponds to $\eta=0.3$, while the smallest absolute
slope corresponds to $\eta=0$, with intermediate values approaching the former
asymptotically.

Figures 4(C) and 4(D) compare energy curves with their respective asymptotic
forms resembling straight lines in a log-log graph. These asymptotic curves,
of the form $C(1+t)^{p}$, are derived through least squares, effectively
capturing the data’s behavior. In Figure 4(C), distinct $\eta$ values for the
same $\alpha$ exhibit notably different slopes (i.e., decay rates $p$).
Conversely, in Figure 4(D), variations in the multiplicative constant $C$
outweigh changes in the decay rate $p$ across different $\alpha$ values.

## Discussion and Conclusions

We have demonstrated the existence, uniqueness, and stability of a Rao–Nakra
model with dissipative terms of the generalized Caputo fractional derivative
type within the domain, building upon ideas from a proven case for a broader
class of evolution systems [4]. While focusing on the Rao–Nakra beam model,
our treatment of the fractional term through an input–output relationship with
a diffusion equation [34] has illuminated a numerical method of finite-
conservative differences type that ensures stable energy decay, unlike
standard methods for approximating fractional derivatives. This enabled us to
numerically simulate energy decay, confirming its polynomial nature.

One key observation from the data and interpretation of Figure 4 and Table 1
is the clear asymptotic polynomial decay behavior of energy. We note that the
decay rate is notably greater in absolute value when $\eta\neq 0$, and lower
when $\eta=0$, validating the theoretical propositions (Proposition 4.3 and
Theorem 4.2). However, we acknowledge discrepancies between the theoretically
derived rates and the numerically estimated ones. We attribute this to factors
such as the discretization step sizes $\delta t$ and $\delta x$, the order of
convergence of the method introducing additional numerical dissipation, and
the limited simulation times for observing true theoretical decay rates.

Future work should delve into numerical methods enhancing convergence order to
obtain numerical values closer to theoretical ones. Additionally, exploring
inverse or control problems, such as identifying parameter pairs
$(\alpha,\eta)$ for quicker energy decay, merits attention. Notably, from
Table 1, we observe a minimum in the rate $p$ near
$(\alpha,\eta)=(0.5;0.0003)$. An inverse problem could involve minimizing this
slope with respect to $(\alpha,\eta)$. Moreover, considering the importance of
the constant $C$, as evidenced by the curve nearest to zero at $t=10,000$
corresponding to $(\alpha,\eta)=(0.5;0.3)$ in Figure 4(A), remains intriguing.

However, beyond $t=10,000[s]$, maintaining a constant $p$ to ensure optimal
decay is unclear. While minimizing energy decay speed is of interest, defining
a suitable cost function characterizing this decay speed requires further
exploration. Thus, identifying the appropriate cost function to minimize for
solving the inverse problem remains an open task.

## References

  * [1] Z. Achouri, N.E. Amroun and A. Benaissa, The Euler-Bernoulli beam equation with boundary dissipation of fractional derivative type. Methods Appl. Sci. 40 (2017), no.11, 3837–3854.
  * [2] R.A. Adams. Sobolev Spaces. Academic Press,New York. 1975.
  * [3] M.S. Alves, J.E. Muñoz Rivera, M. Sepúlveda, O. Vera. Exponential stability in thermoviscoelastic mixtures of solids. I. J. Solid Struct., 46 (2009), 4151–416.
  * [4] K. Ammari, H. Fathi, L. Robbiano. Fractional-Feedback stabilization for a class of evolution systems. J. Differential Equations. 268 (2020), 5751–5791.
  * [5] K. Ammari, F. Hassine, L. Robbiano, Stabilization for Some Fractional-Evolution Systems. SpringerBriefs in Mathematics. Springer, Cham, 2022.
  * [6] C.J.K. Batty and T. Duyckaerts. Non-Uniform stability for bounded semi-groups on Banach spaces. J. Evol. Equ. (2008), 765–780.
  * [7] C.J.K. Batty, R. Chill, Y. Tomilov, Fine scales of decay of operator semigroups. J. Eur. Math. Soc., 18 (2016), 853–929.
  * [8] A. Benaissa, A. Boudaouad, Some results on the energy decay of solutions for a wave equation with a general internal feedback of diffusive type. Advances in Partial Differential Equations and Control. Trends Math. Birkhäuser/Springer, Cham, 2024.
  * [9] A. Borichev, Y. Tomilov. Optimal polynomial decay of functions and operator semigroups. Math. Ann., Vol. 347. 2 (2010), 455–478.
  * [10] N. Burq. Décroissance de l’energie locale de l’équation des ondes pour le problème extérieur et absence de résonance au voisinage du réel. Acta Math. (1998), 1–29.
  * [11] M. Cai, C. Li. Numerical approaches to fractional integrals and derivatives: A review. Mathematics, Vol. 8(1), 43 (2020).
  * [12] X.-G. Cao, D.-Y. Liu, G.-Q. Xu. Easy test for stability of laminated beams with structural damping and boundary feedback controls. J. Dyn. Control Syst. Vol. 13. 3 (2007), 313–336.
  * [13] M. Caputo. Linear models of dissipation whose Q is almost frequency independent Part II. Geophys. J. R. Astr. Soc.
  * [14] M. Caputo. Elasticità e Dissipazione. Zanichelli, Bologna (1969). (in Italian: Elasticity and Dissipation).
  * [15] M. Caputo, F. Mainardi. Linear models of dissipation in anelastic solids. Riv. Nuovo Cimento. (Ser. II) 1 (1971), 161–198.
  * [16] J. Choi, R. Maccamy. Fractional order Volterra equations with applications to elasticity. J. Math. Anal. Appl., 139 (1989), 448–464.
  * [17] M. Ciesielski. Numerical Algorithms for Approximation of Fractional Integrals and Derivatives Based on Quintic Spline Interpolation. Symmetry, 16(2), 252 (2024).
  * [18] S. Das. Functional Fractional Calculus for System Identification and Control. Springer Science & Business Media, 2011\.
  * [19] M. Duarte Ortigueira, F. Coito From Differences to Derivatives. Fractional Calculus and Applied Analysis. 7 (2004) 459–471.
  * [20] L.M. Gearhart. Spectral theory for contraction semigroups on Hilbert spaces. Trans. Amer. Math. Soc. 236 (1978), 385–394.
  * [21] F.L. Huang. Characteristic condition for exponential stability of linear dynamical systems in Hilbert spaces. Ann. Diff. Eqns., Fuzhou. 1 (1985), 43–56.
  * [22] S.W. Hansen, R. Spies. Structural damping in a laminated beams due to interfacial slip. J. Sound Vib. 204 (1997), 183–202.
  * [23] S.W. Hansen. In Control and Estimation of Distributed Parameter Systems: Non-linear Phenomena. International Series of Numerical Analysis. ISNA 118(1994)143-170. Basel: Birkha$\ddot{u}$ser. A Model for a two-layered plate with interfacial slip.
  * [24] A. Kilbas, H. Srivastava, J. Trujillo. Theory and Applications of Fractional Differential Equations. North-Holland Mathematics Studies, 204. Elsevier Science B.V., Amsterdam, 2006.
  * [25] V. Komornik. Exact controllability and stabilization: The multiplier method. John Wiley & Sons. 1994.
  * [26] C. Li, F. Zeng. Numerical Methods for Fractional Calculus. Boca Raton FL: CRC Press Taylor & Francis Group; 2015.
  * [27] Z. Liu, B. Rao, Q. Zheng. Polynomial stability of the Rao-Nakra beam with a single internal viscous damping. J Differential Equations. Vol. 269 (2020), 6125-6162.
  * [28] Z. Liu, S.A. Trogdon, J. Yong. Modeling and analysis of a laminated beam. Comput Math Model. Vol.30. (1-2)(1999), 149–167.
  * [29] Z. Liu, S. Zheng. Semigroups associated with dissipative systems. New York: Chapman & Hall/CRC. 1999.
  * [30] A. Lo, N-E Tatar. Stabilization of laminated beams with interfacial slip. EJDE. 129 (2015), 1–14.
  * [31] J.A. Machado, I.S. Jesus, R. Barbosa, M. Silva, C. Rei. Application of fractional calculus in engineering, in: Dynamics, Games and Science I. 1 (2011), 619–629.
  * [32] R.L. Magin. Fractional Calculus in Bioengineering. Begell House. Redding. CT. USA. 2006.
  * [33] T. Maryati, J. E. Muñoz Rivera, V. Poblete, O. Vera. Asymptotic behavior in a laminated beams due interfacial slip with a boundary dissipation of fractional derivative. Applied Mathematics and Optimization. 84 (2021), 85–102.
  * [34] B. Mbodje. Wave energy decay under fractional derivative controls. IMA J. Math. Contr. Inform. 23 (2006), 237–257.
  * [35] B. Mbodje, G. Montseny. Boundary fractional derivative control of the wave equation. IEEE Trans. Autom. Control. 40 (1995), 378–382.
  * [36] J. E. Muñoz Rivera, V. Poblete, O. Vera. Stability for an Klein-Gordon equation type with a boundary dissipation of fractional derivative type. Asymptotic Analysis. Vol. 127. 3 (2022), 249–273.
  * [37] M. Nakao. On the decay of solutions of some nonlinear dissipative wave equations. Math. Z. Berlin. 193 (1986), 227–234.
  * [38] N.M. Newmark. A method of computation for structural dynamics. J. Engrg. Mech. Div., ASCE 85 (1959).
  * [39] A. Ozkan ozer, S.W. Hansen. Uniform stabilization of a multilayer Rao-Nakra sandwich beam. Evol. Equ. Control Theory. Vol. 2 (2013), 695–710.
  * [40] A. Ozkan Ozer, S.W. Hansen. Exact boundary controllability results for a multilayer Rao-Nakra sandwich beam. SIAM J Control Optim. Vol. 52 (2014), 1314–1337.
  * [41] A. Pazy. Semigroups of Linear Operators and Applications to Partial Differential Equations. Springer- Verlag, New York, 1983.
  * [42] V. Poblete, F. Toledo, O. Vera. Polynomial stability for a Piezoelectric beam with magnetic effect with magnetic effect and a boundary dissipation of fractional derivative type. Proceedings of the Edinburgh Mathematical Society. In press.
  * [43] J. Prüss. On the spectrum of C0-semigroup. Trans. Am. Math. Soc. 248 (1984), 847–867.
  * [44] J.A. Ramos, C. Nonato, C.A, Raposo, O. Vera. Stability for a weakly coupled equations with a boundary dissipation of fractional derivative type. Rendiconti del Circolo Matematico di Palermo Series 2. (2021) doi.org/10.1007/s12215-021-00703-w .
  * [45] Y.V.K.S Rao, B.C. Nakra. Vibrations of unsymmetrical sanwich beams and plates with viscoelastic cores. J Sound Vibr. Vol. 34. 3 (1974), 309-326.
  * [46] C.A. Raposo, O. Vera, J. Ferreira, E. Piskin. Rao–Nakra sandwich beam with second sound. Partial Differential Equations in Applied Mathematics. 4 (2021), 100053.
  * [47] S. Samko, A. Kilbas, O. Marichev. Integral and Derivatives of Fractional Order. Gordon Breach, New York. 1993.
  * [48] V.E. Tarasov. Fractional Dynamics: Applications of Fractional Calculus to Dynamics of Particles, Fields and Media. Springer. Berlin/Heidelberg. Germany. 2011.
  * [49] D. Valério, J.A.T. Machado, V. Kiryakova. Some pioneers of the applications of fractional calculus. Fract. Calc. Appl. Anal. 17 (2014), 552–578.
  * [50] O. Vera, C.A., Raposo, Carlos A. Nonato, Anderson J.A. Ramos. Stability of solution for Rao-Nakra Sandwich beam with boundary dissipation of fractional derivative type. Journal of Fractional Calculus and Applications. Vol. 13. 2 (2022), 116–143.
  * [51] M. Walter. Dynamical System and evolution Equations, Theory and Applications. Plenum Press. New York. 1980.
  * [52] J.-M. Wang, G.-Q. Xu, S.-P. Yung. Exponential stabilization of laminated beams with structural damping and boundary feedback controls. SIAM J. Control Optim. Vol. 44\. 5 (2005), 1575–1597.