1 Introduction

The modeling of generalized Kirchhoff viscoelastic plate, where a bending moment relation with infinite memory is considered, can be described by the following weakly dissipative viscoelastic equation:

$$\begin{aligned} {\left\{ \begin{array}{ll} u''+\varDelta ^2u+\displaystyle \int _0^\infty g(s)\varDelta u(t-s){\text {d}}s=0 &{}\mathrm {in}\,\,\,\varOmega \times (0,\infty ),\\ u(x,-t)=u_0(x,t),~~ u'(x,0)=u_1(x),&{}\mathrm {in}~~\varOmega ,\\ \end{array}\right. } \end{aligned}$$
(1.1)

subject to the homogeneous conditions \(u=\varDelta u =0\) on \(\partial \varOmega \), where the physical domain \(\varOmega \subset {\mathbb {R}}^n\) is a bounded domain with a smooth (or piecewise smooth) boundary \(\partial \varOmega \). The history condition \(u(x,-t)=u_0(x,t)\) means that we are taking into account all the deformation the material has undergone before the instant \(t=0\). The initial velocity \(u_1\) and the non-negative relaxation function g are given.

The asymptotic or decaying behaviour, of different types of viscoelastic equations, including finite and infinite memories also linear and nonlinear dampings, were the subject of study of many researchers since the pioneer work of Dafermos [4, 5]. The achieved energy decaying in the literature varies between (fractional) polynomial and exponential rates. This depends on the damping term, the memory term, the relaxation function, and the differential operator. In the presence of memory, the main common challenge was showing the energy decaying for most general relaxation function g,  we refer the readers to [1,2,3, 6, 7, 9,10,11, 13].

Assuming that the relaxation function g in (1.1) satisfies \(g'(t)\le -\delta g(t),\) for \(t\ge 0\) and for some positive constant \(\delta \), Revira et al. [15] showed an exponential convergence of the energy decay. The main focus of this work was on investigating the energy decay of problem (1.1) but for a wider class of g, see (2.2) below. Owing to the presence of the weakly dissipative term, a second energy functional is introduced to achieve our goal. To make use of the convolution properties and the history condition, the viscoelastic (infinite memory) term is decomposed accordingly. For numerical illustrations of our theoretical finding, we provide a graphical illustration of the numerical energy decay using an approximate solution of problem (1.1) based on finite differences in time and finite elements in space.

Outline of the paper In the next section, we introduce some necessary notations and assumptions. For our decaying analysis, we state and prove a few technical lemmas. Section 3 is dedicated to show the decaying rates of the energy functional E, see Theorem 3.4. Having a weakly dissipative term in problem (1.1) leads us to introduce a second energy functional \({\mathcal {E}}\) [see (3.4)] to overcome the difficulties in proving the decay of E. For the sake of illustrating the theoretical decaying rate of E numerically, we develop a fully discrete numerical method in Sect. 4. To avoid dealing with \({{\mathcal {C}}}^2\) numerical methods in the spatial variable (which is often not convenient on various physical domains \(\varOmega \)) due to the presence of the biharmonic operator in problem (1.1), we rewrite (1.1) as a coupled system that involves a second-order elliptic operator instead. Then, we apply the \({\mathcal {C}}^0\) Galerkin finite element method to discretize in space. In the time variable, a second central difference is used to handle the second time derivative, while the other terms are approximated appropriately. We show the decaying of both, the numerical solution of problem (1.1) and also the approximation of the energy functional E.

2 Preliminaries

For \(\ell \ge 0,\) \(H^\ell (\varOmega )\) is the standard Sobolev space which reduces to \(L^2(\varOmega )\) space when \(\ell =0.\) On this space, \({\langle \cdot ,\cdot \rangle }\) is the usual inner product and \(\Vert \cdot \Vert \) is the associated \(L^2(\varOmega )\)-norm.

An application of the Poincaré inequality and using the elliptic regularity property,

$$\begin{aligned} \frac{1}{\omega _0} \Vert \nabla w\Vert ^2\le \Vert \varDelta w\Vert ^2\le \omega _1~\Vert \nabla (\varDelta w)\Vert ^2,~\forall w\in \mathcal {H}(\varOmega ), \end{aligned}$$
(2.1)

for some positive constants \(\omega _0\) and \(\omega _1\), where \(\mathcal {H}(\varOmega )=\{u\in H^3(\varOmega ):u=\varDelta u=0~ \mathrm {on}~ \partial \varOmega \}.\)

In the decaying energy analysis (including Lemmas 2.1 and 2.2), we assume that the relaxation function \(g\in \mathcal {C}^1({\mathbb {R}}^+)\) and satisfies

$$\begin{aligned} g(0)>0,~~ 1-\max \{\omega _0,\omega _1\}\mu ^0=:l>0, ~~\mathrm{and} ~~-c_0g(t)<g'(t)\le -\xi (t)g(t), ~~~ \forall ~t\ge 0, \end{aligned}$$
(2.2)

where \(\xi \) is a positive non-increasing \(\mathcal {C}^1\) function, \(\mu ^0=\int _{0}^{\infty }g(s){\text {d}}s\),  and \(c_0\) is a positive constant.

For later use, by (2.1) and the second inequality in (2.2), we have, for \(t\ge 0,\)

$$\begin{aligned} \frac{1}{\mu ^0}\Vert \nabla (\varDelta u(t))\Vert ^2\ge \Vert \varDelta u(t)\Vert ^2\ge \Vert \nabla u(t)\Vert ^2\mu ^0\,. \end{aligned}$$
(2.3)

For convenience, we introduce the following notations: for \(t\ge 0,\)

$$\begin{aligned} (g\circ u)(t):=\int _{0}^{\infty }g(s)\Vert \eta ^t(s)\Vert ^2{\text {d}}s,\quad (g\circ [\nabla u, \varDelta u])(t):=(g\circ \nabla u)(t)+(g\circ \varDelta u)(t),\end{aligned}$$

and for \(0<\varepsilon <1\),

$$\begin{aligned}C_\varepsilon :=\int _0^\infty \frac{g^2(s)}{h_\varepsilon (s)}{\text {d}}s \qquad \mathrm {with}\qquad h_\varepsilon (t):=\varepsilon g(t)-g'(t), \end{aligned}$$

where the function \(\eta ^t\) is the relative history of u [5], defined as \(\eta ^t(s)=u(t)-u(t-s)\).

The next two lemmas will be used in the forthcoming decaying analysis section:

Lemma 2.1

[8] For any \(v \in L^2_{loc} \big ([0,+\infty );L^2(\varOmega )\big )\),

$$\begin{aligned}\Big \Vert \int _0^\infty g(s)(v(t-s)-v(t)){\text {d}}s\Big \Vert ^2\le C_\varepsilon (h_\varepsilon \circ v)(t), \qquad \forall \,t\ge 0.\end{aligned}$$

Lemma 2.2

[14] For any \(\varphi \in H^1 \big ([0,\infty );L^2(\varOmega )\big )\),

$$\begin{aligned}\int _{0}^{\infty }{\langle g(s)\varphi (t-s),\varphi '(t)\rangle }{\text {d}}s =\frac{1}{2}\frac{{\text {d}}}{{\text {d}}t}\Big [\mu ^0\Vert \varphi (t)\Vert ^2-(g\circ \varphi )(t)\Big ]+\frac{1}{2}(g'\circ \varphi )(t).\end{aligned}$$

3 Decay

In this section, we find the energy functional E of problem (1.1). As a starting point, taking the inner product of (1.1) with \(u'\), and then applying Green’s formula (twice for the second term and once for the third term) and using the fact that \(u'=\varDelta u'=0\) on \(\partial \varOmega \), yield the following weak formulation of (1.1):

$$\begin{aligned} {\langle u''(t),u'(t)\rangle }+{\langle \varDelta u(t),\varDelta u'(t)\rangle }-\int _{0}^{\infty }g(s){\langle \nabla u(t-s),\nabla u'(t)\rangle }{\text {d}}s=0\,.\end{aligned}$$

Using Lemma 2.2 with \(\varphi =\nabla u\), this equation can be rewritten as follows:

$$\begin{aligned} E'(t) =\frac{1}{2}(g'\circ \nabla u)(t)\le 0. \end{aligned}$$
(3.1)

where E is the first energy functional given by

$$\begin{aligned} E(t):=\frac{1}{2}\Big [\Vert u'(t)\Vert ^2+\Vert \varDelta u(t)\Vert ^2-\mu ^0\Vert \nabla u(t)\Vert ^2+(g\circ \nabla u)(t)\Big ]\ge 0, \end{aligned}$$
(3.2)

where the non-negative property of E follows from (2.3).

Now, taking the inner product of (1.1) with \(-\varDelta u'\), following the above steps, we get

$$\begin{aligned} {\mathcal {E}}'(t)=\frac{1}{2}(g'\circ \varDelta u)(t)\le 0, \end{aligned}$$
(3.3)

where \({\mathcal {E}}\) is the second energy functional, defined as follows:

$$\begin{aligned} {\mathcal {E}}(t):= \frac{1}{2}\Big [\Vert \nabla u'(t)\Vert ^2+\Vert \nabla (\varDelta u(t))\Vert ^2-\mu ^0\Vert \varDelta u(t)\Vert ^2+(g\circ \varDelta u)(t) \Big ]\ge 0\,,\end{aligned}$$
(3.4)

and its non-negativity follows also from (2.3). In the remaining analysis, c is a generic positive constant. We estimate in the next two lemmas the time derivative of the following functionals:

$$\begin{aligned}I_1(t)={\langle u'(t),u(t)\rangle }\quad \mathrm{and}\quad I_2(t)=-\int _{0}^{\infty }{\langle \eta ^t(s),u'(t)\rangle }g(s){\text {d}}s.\end{aligned}$$

Lemma 3.1

Along the solution of (1.1), we have

$$\begin{aligned} I'_1(t)\le \Vert u'(t)\Vert ^2-\frac{l}{2}\Vert \varDelta u(t)\Vert ^2+c ~C_{\varepsilon }(h_{\varepsilon }\circ \nabla u)(t)\,. \end{aligned}$$
(3.5)

Proof

Since \(I'_1(t)= \Vert u'(t)\Vert ^2+{\langle u''(t),u(t)\rangle },\) using (1.1) and Green’s formula, we get

$$\begin{aligned}I'_1(t)=\Vert u'(t)\Vert ^2-\Vert \varDelta u(t)\Vert ^2 +\mu ^0\Vert \nabla u(t)\Vert ^2 -\int _{0}^{\infty }g(s) {\langle \nabla \eta ^t(s),u(t)\rangle }\,{\text {d}}s. \end{aligned}$$

Young’s inequality, Lemma 2.1, and the inequalities in (2.1) imply that

$$\begin{aligned} \mu ^0\Vert \nabla u(t)\Vert ^2 -\int _{0}^{\infty }g(s) {\langle \nabla \eta ^t(s),u(t)\rangle }\,{\text {d}}s&\le \omega _0\mu ^0\Vert \varDelta u(t)\Vert ^2+ \frac{l}{2\omega _0}\Vert \nabla u(t)\Vert ^2+\frac{\omega _0}{2l}\Big \Vert \int _{0}^{\infty }\nabla \eta ^t(s) g(s){\text {d}}s\Big \Vert ^2\\&\le (1-l)\Vert \varDelta u(t)\Vert ^2+ \frac{l}{2}\Vert \varDelta u(t)\Vert ^2+cC_\varepsilon (h_\varepsilon \circ \nabla u)(t)\,. \end{aligned}$$

Combining the above results completes the proof. \(\square \)

Lemma 3.2

Along the solution of (1.1) and for \(\delta >0,\) we have

$$\begin{aligned} I'_2(t)\le \delta \Vert \varDelta u(t)\Vert ^2-\left( \mu ^0-\delta \right) \Vert u'(t)\Vert ^2+\frac{c}{\delta }(C_{\varepsilon }+1)(h_{\varepsilon }\circ \varDelta u)(t). \end{aligned}$$
(3.6)

Proof

Differentiating \(I_2\) and using the differential equation in (1.1), we get

$$\begin{aligned} I'_2(t)=I_{2,1}(t)+I_{2,2}(t)+I_{2,3}(t)-\mu ^0\Vert u'(t)\Vert ^2, \end{aligned}$$
(3.7)

where

$$\begin{aligned} I_{2,1}(t)&=\int _{0}^{\infty }g(s){\langle \eta ^t(s),\varDelta ^2u(t)\rangle }{\text {d}}s,\quad I_{2,2}(t)={\biggl \langle \int _{0}^{\infty }g(s)\varDelta u(t-s){\text {d}}s,\int _{0}^{\infty }g(s)\eta ^t(s){\text {d}}s\biggr \rangle },\\ I_{2,3}(t)&={\biggl \langle \int _{0}^{\infty }g(s)\eta ^t_s(s)\,{\text {d}}s,u'(t)\biggr \rangle }. \end{aligned}$$

By Green’s formula, Young’s inequality and Lemma 2.1, we have

$$\begin{aligned} I_{2,1}(t) \le \Vert \varDelta u(t)\Vert \int _{0}^{\infty }g(s)\Vert \varDelta \eta ^t(s)\Vert {\text {d}}s \le \frac{\delta }{2}\Vert \varDelta u(t)\Vert ^2+\frac{c}{\delta }C_\varepsilon (h_\varepsilon \circ \varDelta u)(t), \end{aligned}$$

and in addition, using \(u(t-s)=u(t)-\eta ^t(s)\),

$$\begin{aligned} I_{2,2}(t)&=\Big \Vert \int _{0}^{\infty }g(s)\nabla \eta ^t(s){\text {d}}s\Big \Vert ^2 -{\biggl \langle \int _{0}^{\infty }g(s)\nabla u(t){\text {d}}s,\int _{0}^{\infty }g(s)\nabla \eta ^t(s) {\text {d}}s\biggr \rangle }\\&\le C_\varepsilon (h_\varepsilon \circ \nabla u)(t) +c\Vert \nabla u(t)\Vert \int _{0}^{\infty }g(s)\Vert \nabla \eta ^t(s)\Vert {\text {d}}s\\&\le C_\varepsilon (h_\varepsilon \circ \nabla u)(t) +\frac{\delta }{2\omega _0}\Vert \nabla u(t)\Vert ^2+\frac{c}{\delta }C_\varepsilon (h_\varepsilon \circ \nabla u)(t)\,. \end{aligned}$$

To estimate \(I_{2,3}(t)\), we perform integration by parts and then make use of Young’s and Hölder’s inequalities, the fact that \(g'=\varepsilon g-h_\varepsilon \) and Lemma 2.1. So, we obtain

$$\begin{aligned} I_{2,3}(t)&=-{\biggl \langle \int _{0}^{\infty }g'(s)\eta ^t(s){\text {d}}s,u'(t)\biggr \rangle }\\&=-\varepsilon {\biggl \langle \int _{0}^{\infty }g(s)\eta ^t(s){\text {d}}s,u'(t)\biggr \rangle }+ {\biggl \langle \int _{0}^{\infty }h_\varepsilon (s)\eta ^t(s){\text {d}}s,u'(t)\biggr \rangle }\\&\le \frac{\delta }{2}\Vert u'(t)\Vert ^2+\frac{\varepsilon ^2}{2\delta }\Big \Vert \int _{0}^{\infty }g(s)\eta ^t(s){\text {d}}s\Big \Vert ^2+\frac{\delta }{2}\Vert u'(t)\Vert ^2+\frac{1}{2\delta }\Big \Vert \int _{0}^{\infty }h_\varepsilon (s)\eta ^t(s){\text {d}}s\Big \Vert ^2\\&\le \delta \Vert u'(t)\Vert ^2+\frac{c}{\delta }(h_\varepsilon \circ u)(t)+\frac{1}{2\delta }\Big (\int _{0}^{\infty }h^{\frac{1}{2}}_\varepsilon (s)h^{\frac{1}{2}}_\varepsilon (s)\Vert \eta ^t(s)\Vert {\text {d}}s\Big )^2\,. \end{aligned}$$

Inserting the obtained estimates of \(I_{2,1}(t)\), \(I_{2,2}(t),\) and \(I_{2,3}(t)\) in (3.7) and using

$$\begin{aligned}\Big (\int _{0}^{\infty }h^{\frac{1}{2}}_\varepsilon (s)h^{\frac{1}{2}}_\varepsilon (s)\Vert \eta ^t(s)\Vert {\text {d}}s\Big )^2\le \int _{0}^{\infty }h_\varepsilon (s){\text {d}}s\int _{0}^{\infty }h_\varepsilon (s)\Vert \eta ^t(s)\Vert ^2{\text {d}}s\le c(h_\varepsilon \circ u)(t),\end{aligned}$$

(by Young’s and Hölder’s inequalities), in addition to the elliptic regularity property of the operator \(-\varDelta \), the desired bound followed. \(\square \)

Lemma 3.3

For \(\varepsilon _1, \, \varepsilon _2 >0\), the functional \(\mathscr {L}(t):=N(E(t) + {\mathcal {E}}(t)) + \varepsilon _1 I_1(t) + \varepsilon _2 I_2(t)\) satisfies \(\mathscr {L} \sim E + \mathcal {E}\) for a sufficiently large N (the proof is in [12]). Moreover, there exist positive constants \(\alpha _1\) and \(\alpha _1\), such that

$$\begin{aligned} \mathscr {L}'(t) \le -\alpha _1 E(t)+ \alpha _2 \int _{0}^{\infty }g(s)\left( \Vert \nabla \eta ^t(s)\Vert ^2+\Vert \varDelta \eta ^t(s)\Vert ^2\right) {\text {d}}s,\quad \forall t\ge 0. \end{aligned}$$
(3.8)

Proof

To prove the estimate in (3.8), we differentiate \(\mathscr {L}\) and use (3.1) and (3.3), in addition to Lemmas 3.1 and 3.2 . Hence,

$$\begin{aligned}&\mathscr {L}'(t)\le \frac{N}{2}(g'\circ [\nabla u, \varDelta u])(t)+\varepsilon _1\Big [\Vert u'(t)\Vert ^2-\frac{l}{2}\Vert \varDelta u(t)\Vert ^2+c ~C_{\varepsilon }(h_{\varepsilon }\circ \nabla u)(t)\Big ]\\&\quad +\varepsilon _2\Big [\delta \Vert \varDelta u(t)\Vert ^2-\left( \mu ^0-\delta \right) \Vert u'(t)\Vert ^2+\frac{c}{\delta }(C_{\varepsilon }+1)(h_{\varepsilon }\circ \varDelta u)(t)\Big ]\,. \end{aligned}$$

Recalling that \( g'(t)=\varepsilon g(t)-h_\varepsilon (t)\), then rearranging the terms in the above equation, and noting that \(h_\varepsilon >0,\) we arrive at

$$\begin{aligned}&\mathscr {L}'(t)\le -\Big [\left( \mu ^0-\delta \right) \varepsilon _2-\varepsilon _1\Big ]\Vert u'\Vert ^2-\big (\frac{l}{2}\varepsilon _1-\varepsilon _2\delta \big )\Vert \varDelta u\Vert ^2+\frac{N\varepsilon }{2}(g\circ [\nabla u, \varDelta u]) \\&\quad +\Big [\frac{c}{\delta }(1+C_\varepsilon )(\varepsilon _1+\varepsilon _2)-\frac{N}{2}\Big ](h_{\varepsilon }\circ \varDelta u)(t). \end{aligned}$$

Choosing \(\varepsilon =\frac{1}{N}\). Since \(\displaystyle \frac{ g^2(s)}{\varepsilon g(s)-g'(s)}\le \frac{g(s)}{\varepsilon },\) \(C_\varepsilon \le N\mu ^0\). Hence, for \(N\ge 2\frac{c}{\delta }(1+N\mu ^0)(\varepsilon _1+\varepsilon _2)\) (which is valid for \(\varepsilon _1\) and \(\varepsilon _2\) are sufficiently small), we notice from the above equation that

$$\begin{aligned}\mathscr {L}'(t)\le -\Big [\left( \mu ^0-\delta \right) \varepsilon _2-\varepsilon _1\Big ]\Vert u'\Vert ^2-\big (\frac{l}{2}\varepsilon _1-\varepsilon _2\delta \big )\Vert \varDelta u\Vert ^2+\frac{1}{2}(g\circ [\nabla u, \varDelta u]). \end{aligned}$$

Choosing \(\delta<\frac{l}{8}\mu ^0<\frac{1}{8}\mu ^0\) (because \(0<l<1\)) and \(\varepsilon _1=\frac{3}{8}\mu ^0\varepsilon _2\), after some calculations, we have

$$\begin{aligned}(\mu ^0-\delta )\varepsilon _2 - \varepsilon _1> \frac{1}{2}\mu ^0\varepsilon _2,~\frac{l}{2}\varepsilon _1 - \delta \varepsilon _2 > \frac{l}{16}\mu ^0\varepsilon _2,\end{aligned}$$

and consequently,

$$\begin{aligned} \mathscr {L}'(t) \le - \frac{1}{2}\mu ^0\varepsilon _2 \Vert u'\Vert ^2 - \frac{l}{16}\mu ^0\varepsilon _2\Vert \varDelta u\Vert ^2 + \frac{1}{2}(g\circ [\nabla u, \varDelta u])(t). \end{aligned}$$

Finally, recalling the definition of E and \(\eta ^t\), and using the above bound, we obtain the desired estimate. \(\square \)

We are ready now to estimate the energy functional associated with the problem (1.1).

Theorem 3.4

Assume that \(u_0(.,0)\in \mathcal {H}(\varOmega )\) with \(\max \{\Vert \nabla u_0(.,s)\Vert ,\Vert \varDelta u_0(.,s)\Vert \}\le m\) (for some positive constant m), and \(u_1\in H^2(\varOmega )\cap H^1_0(\varOmega )\). Then, for t large enough, there exists a positive constant \(\delta >0\) such that

$$\begin{aligned} E(t)\le \delta \left( \frac{1+\int _{0}^{t}\xi (s)\int _s^\infty g(s){\text {d}}s {\text {d}}s}{1+\int _{0}^{t}\xi (s){\text {d}}s}\right) .\end{aligned}$$

Proof

From the non-increasing property of \(\xi ,\) the inequality \(g'\le -\xi \,g,\) we have for \(t\ge 0;\)

$$\begin{aligned}\xi (t)\int _{0}^{t}g(s)\Vert \nabla \eta ^t(s)\Vert ^2{\text {d}}s \le \int _{0}^{t}\Vert \nabla \eta ^t(s)\Vert ^2\xi (s)g(s){\text {d}}s \le -\int _{0}^{t}\Vert \nabla \eta ^t(s)\Vert ^2g'(s){\text {d}}s \le -c_0E'(t),\end{aligned}$$

for some positive constant \(c_0.\) However, owing to the inequality, \(\left\| \nabla u(t) \right\| ^2\le \frac{2\omega _0}{l}E(t)\le \frac{2\omega _0}{l}E(0)\), we have for \(s>t,\)

$$\begin{aligned} \Vert \nabla \eta ^t(s)\Vert ^2&\le 2\Vert \nabla u(t)\Vert ^2+2\Vert \nabla u(t-s)\Vert ^2 \le \frac{4\omega _0}{l}E(0)+ 2\sup _{\zeta >0}\Vert \nabla u_0(\zeta )\Vert ^2 \le \frac{4\omega _0}{l}E(0)+ 2m^2, \end{aligned}$$

Consequently, we get

$$\begin{aligned} \xi (t)\int _{t}^{\infty }g(s)\Vert \nabla \eta ^t(s)\Vert ^2{\text {d}}s \le \xi (t)\Big (\frac{4\omega _0}{l}E(0)+2m^2\Big ) \int _{t}^{\infty }g(s){\text {d}}s.\end{aligned}$$

Therefore, for \(t\ge 0,\)

$$\begin{aligned} \xi (t)\int _{0}^{\infty }g(s)\Vert \nabla \eta ^t(s)\Vert ^2{\text {d}}s \le -c_0 E'(t)+\xi (t)\Big (\frac{4\omega _0}{l}E(0)+2m^2\Big ) \int _{t}^{\infty }g(s){\text {d}}s. \end{aligned}$$
(3.9)

Following similar arguments, and using the inequality \(\left\| \varDelta u(t) \right\| ^2\le \frac{2}{l}E(t)\le \frac{2}{l}E(0)\), we obtain

$$\begin{aligned} \xi (t)\int _{0}^{\infty }g(s)\Vert \varDelta \eta ^t(s)\Vert ^2{\text {d}}s \le -c_0 \mathcal {E}'(t))+\xi (t)\Big (\frac{4}{l}E(0)+2m^2\Big ) \int _{t}^{\infty }g(s){\text {d}}s. \end{aligned}$$
(3.10)

Multiplying (3.8) by \(\xi (t)\), and using (3.9) and (3.10), we get

$$\begin{aligned} \xi (t)\mathscr {L}'(t) \le -\alpha _1\xi (t) E(t) -\beta _1 (E'(t)+\mathcal {E}'(t))+\beta _2\,\xi (t) \int _{t}^{\infty }g(s){\text {d}}s, \end{aligned}$$
(3.11)

where \(\beta _1=c_0 \alpha _2\) and \(\beta _2=\alpha _2\Big (\frac{4}{l}(1+\omega _0)E(0)+2m^2\Big ).\) Therefore, differentiating the functional \({\mathscr {L}}(t):=\xi \mathscr {L}(t)+\beta _1(E(t)+\mathcal {E}(t))\) and using (3.11), we have

$$\begin{aligned}{\mathscr {L}}'(t) \le \xi (t) \mathscr {L}'(t)+\beta _1(E'(t)+\mathcal {E}'(t)) \le -\alpha _1\xi (t) E(t)+ \beta _2\, \xi (t) \int _{t}^{\infty }g(s){\text {d}}s.\end{aligned}$$

Integrating the above equation over (0, t) yields

$$\begin{aligned} {\mathscr {L}}(t)\le {\mathscr {L}}(0)+\beta _2 \int _{0}^{t}\xi (s) \int _s^{\infty }g(s){\text {d}}s\,{\text {d}}s-\alpha _1\int _{0}^{t}\xi (s)E(s){\text {d}}s. \end{aligned}$$

However, \(E(s)\ge E(t)\) for \(s\le t\) (by the non-increasing property of E), and \(\beta _1E(t)\le {\mathscr {L}}(t)\) (from the definition of \({\mathscr {L}}\)), so

$$\begin{aligned} \left( \beta _1 +\alpha _1\int _{0}^{t}\xi (s){\text {d}}s\right) E(t)\le {\mathscr {L}}(0)+\beta _2\int _{0}^{t}\xi (s) \int _s^{\infty }g(s){\text {d}}s\,{\text {d}}s. \end{aligned}$$

Choosing \(\delta =\frac{\max \left\{ {\mathscr {L}}(0),\beta _2\right\} }{\min \{\beta _1,\alpha _1\}}\), then the desired energy bound is obtained. \(\square \)

4 Numerical study

This section is devoted to illustrate numerically the achieved theoretical decaying results in Theorem 3.4 on a two-dimensional test problem of the form (1.1) with space variables x and y. We develop a numerical scheme for problem (1.1) using finite differences for the time discretization combined with the continuous Galerkin finite element method in space. Applying Galerkin method to problem (1.1) directly forces us to deal with \({\mathcal {C}}^2\) polynomial approximations, which is definitely not convenient owing to the complexity in constructing the basis functions on various physical domains, and also increase the cost of computations. To avoid this, we rewrite (1.1) as a coupled system of lower order elliptic problems:

$$\begin{aligned} {\left\{ \begin{array}{ll} u''-\varDelta w+\displaystyle \int _0^t g(t-s)\varDelta u(s){\text {d}}s=F &{}\mathrm {in}\,\,\,\varOmega \times (0,\infty ),\\ w+\varDelta u=0 &{}\mathrm {in}\,\,\,\varOmega \times (0,\infty ),\\ u(x,y,0)=u_0(x,y,0),~~ u'(x,y,0)=u_1(x,y),&{}\mathrm {in}~~\varOmega ,\\ \end{array}\right. } \end{aligned}$$
(4.1)

with \(u=w =0\) on \(\partial \varOmega \), where \(w=-\varDelta u\), and \(F(x,y,t)=- \int _t^\infty g(s)\varDelta u_0(x,y,s-t){\text {d}}s.\)

Our numerical discretization in space is applicable to any bounded polygonal physical domain \(\varOmega \). However, and for simplicity, we choose \(\varOmega \) to be the unit square \((0,1)\times (0,1).\) Let \(\mathcal {T}_h\) be a family of uniform \(M^2\)-square mesh-cells of \(\varOmega \) with diagonal  \(h=\sqrt{2}/M\) each. Let \(V_h \subset H^1_0(\varOmega )\) denote the usual space of continuous, piecewise-linear polynomials on \({\mathcal {T}}_h\) that vanish on \(\partial \varOmega \). To discretize in time, we truncate the interval \((0,\infty )\) and work instead on the finite interval (0, T] with T that can be sufficiently large. Divide [0, T] uniformly into N subintervals each with size \(\tau :=T/N\) and nodes \(\{t_n\}_{n=0}^N.\) For a given grid function \(v^n\), let

$$\begin{aligned} \delta _{tt} v^n=\frac{v^{n+1}-2v^n+v^{n-1}}{\tau ^2},\quad v^{n+\frac{1}{2}}=\frac{v^n+ v^{n-1}}{2},\quad v^{n+\frac{1}{4}}=\frac{ v^{n+1}+2 v^n+ v^{n-1}}{4}. \end{aligned}$$

Taking the inner product of the first two equations in (4.1) with \(\phi ,\,\psi \in H^1_0(\varOmega )\), respectively, then, applying Green’s formula and using \(u=w =0\) on \(\partial \varOmega \). This implies

$$\begin{aligned}{\left\{ \begin{array}{ll} {\langle u'',\phi \rangle }+{\langle \nabla w, \nabla \phi \rangle }-\displaystyle \int _{0}^t g(t-s){\langle \nabla u(s),\nabla \phi \rangle }\,{\text {d}}s={\langle F,\phi \rangle },\quad &{} \forall ~\phi \in H^1_0(\varOmega ),\\ {\langle w,\psi \rangle }-{\langle \nabla u, \nabla \psi \rangle }=0,\quad &{}\forall ~\psi \in H^1_0(\varOmega )\,. \end{array}\right. }\end{aligned}$$

Motivated by the above formulation, our fully discrete numerical scheme is defined as follows: Find \(U_h^{n+1},\,W_h^{n+1} \in V_h\) such that

$$\begin{aligned} {\left\{ \begin{array}{ll} {\langle \delta _{tt} U^n_h,\phi _h\rangle }+{\langle \nabla W_h^{n+\frac{1}{4}},\nabla \phi _h\rangle }-\displaystyle \int _0^{t_{n+1}} g(t_{n+1}-s){\langle \nabla {\bar{U}}_h(s),\nabla \phi _h\rangle }{\text {d}}s={\langle F^{n+1},\phi _h\rangle }, \\ {\langle W_h^{n+1},\psi _h\rangle }-{\langle \nabla U_h^{n+1}, \nabla \psi _h\rangle }=0, \end{array}\right. } \end{aligned}$$
(4.2)

\(\forall ~\phi _h,\,\psi _h \in V_h\), with \(1\le n\le N-1,\) where the piecewise constant function (in time variable) \({\bar{U}}_h(s)=U_h^{j+\frac{1}{2}}\) for \(t_j<s<t_{j+1}\) with \(0\le j\le N\), and \(F^{n+1}=F(t_{n+1}).\)

Fig. 1
figure 1

The graphical plots of the numerical approximations of the energy E(t) (top) and the weighted energy tE(t) (bottom) the against t

Full size image

Let \(\{\phi _p\}_{p=1}^{d_h}\) be the two-dimensional hat basis functions of \(V_h\). Then, \(U_h^n\) and \(W_h^n\) can be written in term of the basis functions as follows: \(U_h^n =\sum _{i=1}^{d_h} b_i^n {\phi _p}_{i}\) and \(W_h^n =\sum _{i=1}^{d_h} c_i^n {\phi _p}_{i}\,.\) We define \(d_{h}\times d_{h}\) matrices: \(\varvec{A}=\left[ {\langle \nabla \phi _{q},\nabla \phi _p\rangle }\right] \) and \(\varvec{M}=\left[ {\langle \phi _{q},\phi _p\rangle }\right] .\) The \(d_h\)-dimensional constant column vectors \(\mathbf{b}^n\), \(\mathbf{c}^n\), and \({\varvec{F}}^{n+1}\) are, respectively, the transpose of the vectors \([b_1^n,b_2^n,\ldots ,b_{d_h}^n],\) \([c_1^n,c_2^n,\ldots ,c_{d_h}^n]\) and \([{\langle F^{n+1},\phi _1\rangle },{\langle F^{n+1},\phi _2\rangle },\ldots ,{\langle F^{n+1},\phi _{d_h}\rangle }]\,.\)

Therefore, the fully discrete scheme (4.2) has the following matrix representation:

$$\begin{aligned} \varvec{M}\delta _{tt}{} \mathbf{b}^n + \varvec{A}\mathbf{c}^{n+\frac{1}{4}}-g_{n+1,n} \varvec{A}\mathbf{b}^{n+\frac{1}{2}} =\sum _{j=0}^{n-1} g_{n+1,j} \varvec{A}\mathbf{b}^{j+\frac{1}{2}}+{\varvec{F}}^{n+1},\quad \varvec{M}\mathbf{c}^{n+1}-\varvec{A}\mathbf{b}^{n+1}=\mathbf{0},\end{aligned}$$

with \(g_{n+1,j}:=\int _{t_j}^{t_{j+1}}g(t_{n+1}-s)\,{\text {d}}s.\) From the second equation, we notice that \(\varvec{A}\mathbf{c}^{n+1}=\varvec{B}\mathbf{b}^{n+1}\) where \(\varvec{B}=\varvec{A}\,\varvec{M}^{-1} \varvec{A}\). Substitute this in the first equation and rearranging the terms,

$$\begin{aligned}&[4\varvec{M}+ \tau ^2 \varvec{B}-2\tau ^2g_{n+1,n} \varvec{A}] \mathbf{b}^{n+1} =\varvec{M}(8\mathbf{b}^n-4\mathbf{b}^{n-1}) - \tau ^2 \varvec{B}(2\mathbf{b}^n+ \mathbf{b}^{n-1})\nonumber \\&\quad +2\tau ^2 g_{n+1,n} \varvec{A}\mathbf{b}^n+2\tau ^2 \sum _{j=0}^{n-1} g_{n+1,j} \varvec{A}[\mathbf{b}^{j+1}+\mathbf{b}^{j}]+{\varvec{F}}^{n+1}\,. \end{aligned}$$
(4.3)

Therefore, at each time level \(t_{n+1},\) the numerical coupled system (4.2) reduces to a finite square linear system, where the unknown is the column vector \(\mathbf{b}^{n+1}\). Whence the coefficient vector \(\mathbf{b}^{n+1}\) is computed, \(\mathbf{c}^{n+1}\) can be determined by solving \(\varvec{M}\mathbf{c}^{n+1}=\varvec{A}\mathbf{b}^{n+1}\), that is, \(\mathbf{c}^{n+1}=\varvec{M}^{-1} \varvec{A}\mathbf{b}^{n+1}\). This will be used to compute the numerical energy function \(E_h(t_n).\)

To solve linear system in (4.3) for \(\mathbf{b}^{n+1}\), successively, the column vectors \(\mathbf{b}^0\) and \(\mathbf{b}^1\) need to be determined first. In other words, we have to compute the approximate solutions \(U_h^0\) and \(U_h^1\) first. We consider \(U_h^0\) to be the Ritz projection of \(u_0\) on the finite dimensional space \(V_h.\) However, motivated by the Taylor series expansion of u in time about \(t=0,\) we choose \(U_h^1\) to be the Ritz projection \(u_0+t_1 u_1\) on \(V_h.\) For sake of computing efficiently the spatial integrals in the linear system in (4.3), on each cell of our two-dimensional partition, the associated integral is approximated using 4-point (that is, 2-point in each direction) Gauss cubature rule.

In our numerical example, choose \(T=150,\) \(u_0(x,y,t)=t^2\sin (\pi x)\sin (\pi y)\) and \(u_1(x,y)=0\). We run our computer program with \(M=20\) (that is, 400 cells in space) and \(N=120000\). We choose \(g(t)= e^{- t}\); then \(g'(t)=-\xi (t)g(t)\) with \(\xi (t)=1.\) Thus, by Theorem 3.4, we expect our energy to decay monomially, that is, \(tE(t)\le c\) for a sufficiently large t. This is confirmed in Fig. 1. In addition, the plot of the numerical solution \(U_h\) in Fig. 2 shows its convergence to zero as the time t getting far away from 0.

Fig. 2
figure 2

The numerical solution plots for \(t=0\) (top-left), \(t=10\) (top-right), \(t=50\) (bottom-left), and \(t=100\) (bottom-right)

Full size image