# Analysis and numerical simulation of the nonlinear beam equation with moving
ends

N. Quintino1, M. A. Rincon1

1Instituto de Matemática, Universidade Federal do Rio de Janeiro, Rio de
Janeiro, Brazil E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The numerical analysis for the small amplitude motion of an elastic beam with
internal damping is investigated in domain with moving ends. An efficient
numerical method is constructed to solve this moving boundary problem. The
stability and convergence of the method is studied, and the errors of both the
semi-discrete and fully-discrete schemes are derived, using Hermite’s
polynomials as a base function, proving that the method has order of quadratic
convergence in space and time.

Numerical simulations using the finite element method associated with the
finite difference method (Newmark’s method) are employed, for one-dimensional
and two-dimensional cases. To validate the theoretical results, tables and
figures are shown comparing approximate and exact solutions. In addition,
numerically the uniform decay rate for energy and the order of convergence of
the approximate solution are also shown.

Keywords: Beam equation; Hermite’s polynomials, Non-cylindrical domain Error
estimates, Order of convergence, Numerical simulation, Newton’s method,
Newmark’s method.

AMS Subject Classification: 35L20; 35R37; 65M60; 65M06.

##  1 Introduction

The equation of motion of a thin beam with weak-internal damping undergoing
cylindrical bending can be written as

$u^{\prime\prime}(x,t)-\left(\zeta_{0}+\zeta_{1}\int_{\Omega_{t}}\big{|}\nabla
u(x,t)\big{|}^{2}dx\right)\Delta u(x,t)+\Delta^{2}u(x,t)+\nu
u^{\prime}(x,t)=0,\mbox{ in }Q_{t},$ (1)

where $\,u\,$ is the transverse displacement, $x\in\mathbb{R}^{n}$ is the
vector of spatial coordinates, $\,t\,$ is the time. The aerodynamic damping
term is denoted by $\,\nu,\,$ $\,\zeta_{1}\,$ is the nonlinear stiffness of
beam, $\,\zeta_{0}\,$ is an in-plane tensile load, and $\,(x,t)\,$ belongs to
the noncylindrical domain, defined by

$Q_{t}=\big{\\{}(x,t)\in\mathbb{R}^{n+1};x=K(t)y,\leavevmode\nobreak\
y\in\Omega,\leavevmode\nobreak\ 0<t<T\big{\\}},$ (2)

with $\Omega=(0,1)^{n}\subset\mathbb{R}^{n}$ is an open, limited and regular
and $K:[0,\infty[\longrightarrow\mathbb{R}$ is a $C^{2}$\- real function. All
quantities are physically nondimensionalized, $\,\nu,\,$ $\,\zeta_{1}\,$ are
fixed positive and $\,\zeta_{0}$ is not necessarily positive.

We will study the equation (1) with zero boundary conditions $\forall
t\in[0,T[$, as follows,

$u(x,t)=0\quad\mbox{and}\quad\nabla u\equiv 0,\leavevmode\nobreak\
\mbox{on}\leavevmode\nobreak\ \partial\Omega_{t}\times[0,T[,$ (3)

where $\partial\Omega_{t}=\\{x=K(t)y;\leavevmode\nobreak\
y\in\partial\Omega\\}$.

In the absence of the damping term $\nu\leavevmode\nobreak\ \\!u_{t}(x,t)$,
many authors have already studied Cauchy’s problems and the mixed problems
associated with the equation (1) on bounded and unlimited cylinders. The
existence and uniqueness of solution was demonstrated in [9].

The Cauchy problem associated with the (1) equation in abstract terms on a
Hilbert space is studied, among other authors, by [4, 18]. In addition, in the
last were established results on existence, unicity and asymptotic stability
of the solution.

For the one-dimensional case, omitting the term $u_{xxxx}$ and considering the
cylinder $Q=\Omega\times[0,\infty[$, with $\Omega\subset\mathbb{R}$ open and
limited, we obtain the Kirchhoff equation with internal damping, namely,

$u^{\prime\prime}(x,t)-\left(\zeta_{0}+\zeta_{1}\int_{\Omega}\big{|}u_{x}(t)\big{|}^{2}dx\right)u_{xx}(x,t)+\nu
u^{\prime}(x,t)=0,\mbox{ in }Q.$ (4)

The equation (4) has been extensively studied by a wide variety of authors on
$n$-dimensional cases and generic mathematical models defined in Hilbert
spaces. The existence of a local and global solution can be found in several
physical-mathematical contexts, such as [14, 15].

In the article [15] a new model associated with the Kirchhoff equation (4) was
introduced, called the Medeiros-Kirchhoff equation, given by

$\displaystyle
u^{\prime\prime}(x,t)-\left(a(t)+b(t)\int_{\Omega}\big{|}u_{x}(t)\big{|}^{2}dx\right)u_{xx}(x,t)=0,\mbox{
in }Q_{t}.$

There are several areas that work with the application of the Beam equation,
for example, in civil and naval engineering, as well as in the aerospace area
as the reaction of rocket and satellite metal structures to various space
situations, as can be observed in the articles [12, 13]. In addition to these
areas, it is possible to associate Maxwell’s equations, which govern
electromagnetic theory, and the Euler-Bernoulli Beam Equation, in the study
presented by article [16]

The equation (1) follows the Euler-Bernoulli Beam model. In addition to this
model, there are other models, extensions of the Euler-Bernoulli model, which
consider rotational movements besides transverse ones, being possible to
create models closer to reality. The article [11] presents a more detailed
study on three of these models.

The existence and uniqueness of the solution, the asymptotic decay and the
numerical simulations of the beam equation with movinf boundary were studied
in [Raquel] for the one-dimensional case.

In this paper, we present a family of numerical method for the one-dimensional
and two-dimensional cases, based on the finite element method and the finite
difference method. In addition, we are doing numerical analysis for both the
semi-discrete problem and the fully discrete problem, showing that the order
of convergence is quadratic in space and time. For this we will use the Faedo-
Galerkin method, using the Hermite polynomials as base functions. In addition,
the nonlinear system associated with the system of ordinary differential
equations is being solved by Newton’s method and for temporal discretization
is applying the Neumark’s method, where we will show the convergence of the
approximate solution to a family of numerical methods depending on the choice
of $\theta\in]1/4;1]$. Numerical simulations are presented as examples, with
different types of boundary and tables showing the order of numerical
convergence. In addition, we present the asymptotic decay, and the graphs of
the solutions for the one-dimensional and two-dimensional cases.

##  2 Analytical results

Our section will present, without demonstrating, the existence and uniqueness
results of the mixed problem solution for the equation (1), considering the
boundary conditions (3), given by:

$\left\\{\begin{aligned} &\begin{aligned}
u^{\prime\prime}(x,t)-\left(\zeta_{0}+\zeta_{1}\int_{\Omega_{t}}\big{|}\nabla
u(x,t)\big{|}^{2}dx\right)\Delta u(x,t)+\Delta^{2}u(x,t)+\nu
u^{\prime}(x,t)=0,\mbox{in}\leavevmode\nobreak\ Q_{t}\end{aligned}\\\
&\leavevmode\nobreak\ \\!u(x,0)=u_{0}(x),\leavevmode\nobreak\
\leavevmode\nobreak\ u^{\prime}(x,0)=u_{1}(x),\leavevmode\nobreak\
\mbox{in}\leavevmode\nobreak\ \Omega_{0},\\\ &\leavevmode\nobreak\
\\!u(x,t)=0\leavevmode\nobreak\ \mbox{and}\leavevmode\nobreak\ \nabla u\equiv
0,\leavevmode\nobreak\ \mbox{in}\leavevmode\nobreak\
\partial\Omega_{t}\times[0,T[,\end{aligned}\right.$ (5)

where $Q_{t}$ is defined in 2 and satisfying the hypotheses

$\left\\{\begin{aligned}
&Q_{t}\subset\Omega_{t}\times]0,T[,\leavevmode\nobreak\ \mbox{with
}\Omega_{t}\subset\mathbb{R}^{n}\leavevmode\nobreak\ \mbox{open and limited
with }0\in{\Omega_{t}},\\\ &K\in
C^{2}\big{(}[0,+\infty[\big{)},\leavevmode\nobreak\ \mbox{with }0<K_{0}\leq
K(t)\leq K_{1},\quad 0<K^{\prime}(t)\leq K_{2},\\\ \end{aligned}\right.$ (H1)

where $K_{0},K_{1}$ and $K_{2}$ are positive real constants.

Let us consider the diffeomorphism $\tau:Q_{t}\rightarrow
Q=(0,1)^{n}\times(0,T)$, defined by $\tau(x,t)=(y,t)$, with $x=K(t)y$ .
Therefore, the change of variable $u(x,t)=v(y,t)$ transforms the problem(5)
into the equivalent problem with cylindrical boundary

$\left\\{\begin{aligned} &\begin{aligned}
&v^{\prime\prime}(y,t)-b_{1}(t)\big{|}\nabla v\big{|}_{0}^{2}\Delta
v(y,t)+b_{2}(t)\Delta^{2}v(y,t)+\nu\leavevmode\nobreak\
\\!v^{\prime}(y,t)+a^{(2)}_{ij}(y,t)\nabla_{y_{i},y_{j}}v(y,t)\\\ &\quad-
a_{i}^{(1)}(y,t)\Delta_{y_{i}}v(y,t)-a^{(3)}_{i}(y,t)\nabla_{y_{i}}v(y,t)-a^{(4)}_{i}(y,t)\nabla_{y_{i}}v^{\prime}(y,t)=0,\mbox{
in }Q\\\ \end{aligned}\\\ &\leavevmode\nobreak\
v(y,0)=v_{0}(y),\leavevmode\nobreak\ \leavevmode\nobreak\
v^{\prime}(y,0)=v_{1}(y)\quad\mbox{in }\leavevmode\nobreak\ \Omega,\\\
&\leavevmode\nobreak\ v(y,t)=0\leavevmode\nobreak\
\mbox{and}\leavevmode\nobreak\ \nabla v\equiv 0,\leavevmode\nobreak\
\mbox{in}\leavevmode\nobreak\ \partial\Omega\times[0,T[,\end{aligned}\right.$
(6)

where, for simplicity, we denote $K=K(t)$ and

$\displaystyle b_{1}(t)=\frac{\zeta_{1}}{K^{4}},\leavevmode\nobreak\
b_{2}(t)=\frac{1}{K^{4}},\leavevmode\nobreak\
a^{(1)}_{i}(y,t)=\frac{1}{K^{2}}\Big{[}\zeta_{0}-4(y_{i}K^{\prime})^{2}\Big{]},\leavevmode\nobreak\
a^{(2)}_{ij}(y,t)=4y_{i}y_{j}{(K^{\prime}/K)^{2}},$ (7) $\displaystyle
a^{(3)}_{i}(y,t)=\frac{1}{K^{2}}\Big{[}2y_{i}(K^{\prime})^{2}-y_{i}K(\nu
K^{\prime}+K^{\prime\prime})\Big{]},\leavevmode\nobreak\
a^{(4)}_{i}(y,t)=-2y_{i}({K^{\prime}}/{K}).$

In the following, we present the results on the existence and uniqueness of
the weak solution of problems (5) and (6) and asymptotic decay, whose
demonstrations can be found in [8] for the one-dimensional case and [5] for
the two-dimensional case.

######  Theorem 1.

(Existence and Uniqueness of solution in nocylindrical domain.)

Let $u_{0}\in H_{0}^{2}(\Omega_{0})$, $u_{1}\in L^{2}(\Omega_{0})$ and the
hypothesis (H1). Then there is a unique weak solution
$u:\widehat{Q}\longrightarrow\mathbb{R}$ for the problem (5) satisfying the
following conditions:

$\displaystyle 1.\quad u\in L^{2}(0,T;H_{0}^{2}(\Omega_{t})),$ $\displaystyle
2.\quad u^{\prime}\in L^{2}(0,T;L^{2}(\Omega_{t}))$

and the weak formulation of (5)1 is verified in the sense of
$L^{2}(0,T;H^{-2}(\Omega_{t}))$.

######  Theorem 2.

(Existence and Uniqueness of solution in cylindrical domain)

Let $v_{0}\in H_{0}^{2}(\Omega_{0})$, $v_{1}\in L^{2}(\Omega_{0})$ and the
hypothesis (H1). Then , Then there is a unique weak solution
$v:Q\longrightarrow\mathbb{R}$ for the problem (5) satisfying the following
conditions:

$\displaystyle 1.\quad v\in L^{2}(0,T;H_{0}^{2}(\Omega)),$ $\displaystyle
2.\quad v^{\prime}\in L^{2}(0,T;L^{2}(\Omega))$

and the weak formulation of (6)1 is verified in the sense of
$L^{2}(0,T;H^{-2}(\Omega))$.

######  Theorem 3.

(Asymptotic Decay: Two-dimensional case)

Consider the hypothesis H1, then we have the energy (5), given by

$\displaystyle
E(t)=\frac{1}{2}\int_{\Omega}\Big{[}\big{|}u^{\prime}(x,t)\big{|}^{2}+\big{|}\Delta
u(x,t)\big{|}^{2}+\zeta_{0}\big{|}\nabla
u(t)\big{|}_{0}^{2}+\frac{\zeta_{1}}{2}\big{|}\nabla
u(x,t)\big{|}^{4}\Big{]}dx,$

satisfie

$\displaystyle E(t)\leq A_{0}\leavevmode\nobreak\
e^{-A_{1}t},\leavevmode\nobreak\ \quad\forall t\geq 0,$

where $A_{0}$ and $A_{1}$ are positive constants.

From the important theoretical results, we can now develop numerical methods
to determine an approximate numerical solution and verify the consistency of
the numerical methods used.

Due to the equivalence of the problems, the methods used will be for the
problem (6) that is defined in a cylindrical domain.

##  3 Numerical Method

The purpose of this section is to develop a numerical method, based on the
finite element method and the finite difference method to make the numerical
simulations in order to obtain a approximate solution of the problem (6).

Denoting the Hilbert space $H_{0}^{2}(\Omega)=V$, then weak formulation of the
problem (6) can be written in the follows form, when we integrate by parts the
second, third, fifth and sixth terms:

$\left\\{\begin{aligned}
&\big{(}v^{\prime\prime}(t),w\big{)}+b_{1}(t)\big{|}\nabla
v(t)\big{|}_{0}^{2}\big{(}\nabla v(t),\nabla w\big{)}+b_{2}(t)\big{(}\Delta
v(t),\Delta w\big{)}+\nu\leavevmode\nobreak\
\\!\big{(}v^{\prime}(t),w\big{)}\\\
&\quad+\big{(}a^{(1)}_{i}(t)\nabla_{y_{i}}v(t),\nabla_{y_{i}}w\big{)}+\big{(}a^{(2)}_{ij}(t)\nabla_{y_{i},y_{j}}v(t),w\big{)}\\\
&\quad+\big{(}a^{(4)}_{i}(t)\nabla_{y_{i}}v^{\prime}(t),w\big{)}+\big{(}a^{(5)}_{i}(t)\nabla_{y_{i}}v(t),w\big{)}=0,\mbox{
in }Q\\\ &\big{(}v(y,0),w\big{)}=\big{(}v_{0}(y),w\big{)},\leavevmode\nobreak\
\leavevmode\nobreak\
\big{(}v^{\prime}(y,0),w\big{)}=\big{(}v_{1}(y),w\big{)},\leavevmode\nobreak\
\mbox{in }\leavevmode\nobreak\ \Omega,\\\
&\big{(}v,w\big{)}=0\leavevmode\nobreak\ \mbox{and}\leavevmode\nobreak\
\big{(}\nabla v,w\big{)}\equiv 0,\leavevmode\nobreak\ \mbox{in
}\partial\Omega\times[0,T[,\leavevmode\nobreak\ \forall w\in
V,\end{aligned}\right.$ (8)

where $\displaystyle
a^{(5)}_{i}(t)=a^{(3)}_{i}(t)+2\frac{K^{\prime}}{K}a^{(4)}(t)$.

###  3.1 Approximated Problem

Let $V_{h}=[\varphi_{1},\varphi_{1},\cdots,\varphi_{m}]$ be a subspace
generated by the $m$ first vectors of the basis of Hilbert space
$V=H_{0}^{2}(\Omega)$. Then we want to determine an approximate solution
$v_{h}(t)\in V_{h}$ of the following approximate system, given by:

$\displaystyle\big{(}v_{h}^{\prime\prime}(t),w_{h}\big{)}+b_{1}(t)\big{|}\nabla
v_{h}(t)\big{|}_{0}^{2}\big{(}\nabla v_{h}(t),\nabla
w_{h}\big{)}+b_{2}(t)\big{(}\Delta v_{h}(t),\Delta
w_{h}\big{)}+\nu\leavevmode\nobreak\ \\!\big{(}v_{h}^{\prime}(t),w_{h}\big{)}$
(9) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+\big{(}a^{(1)}_{i}(t)\nabla_{y_{i}}v_{h}(t),\nabla_{y_{i}}w_{h}\big{)}-\big{(}a^{(2)}_{ij}(t)\nabla_{y_{i}}v_{h}(t),\nabla_{y_{j}}w_{h}\big{)}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+\big{(}a^{(4)}_{i}(t)\nabla_{y_{i}}v_{h}^{\prime}(t),w_{h}\big{)}+\big{(}a^{(5)}_{i}(t)\nabla_{y_{i}}v_{h}(t),w_{h}\big{)}=0,\leavevmode\nobreak\
\forall w_{h}\in V_{h},$

where we integrate by parts the sixth term.

Since that, $v_{h}(t)\in V_{h}$, then the function can be represented by

$v_{h}(y,t)=\sum_{i=1}^{m}d_{i}(t)\varphi_{i}(y),\quad\varphi_{i}(y)\in
V_{h}.$ (10)

Taking $w_{h}=\varphi_{l}(y)\in V_{h}$, for some $l\in\\{1,\cdots,m\\}$, and
substituting (10) in (9), we obtain

$\displaystyle\sum_{k=1}^{m}\Big{\\{}d^{\prime\prime}_{k}(t)\big{(}\varphi_{k}(y),\varphi_{l}(y)\big{)}+b_{1}(t)d_{k}(t)\Big{|}\sum_{k=1}^{m}d_{k}(t)\nabla\varphi_{k}(y)\Big{|}_{0}^{2}\big{(}\nabla\varphi_{k}(y),\nabla\varphi_{l}(y)\big{)}$
$\displaystyle\quad+b_{2}(t)d_{k}(t)\big{(}\Delta\varphi_{k}(y),\Delta\varphi_{l}(y)\big{)}+d^{\prime}_{k}(t)\Big{[}\nu\leavevmode\nobreak\
\\!\big{(}\varphi_{k}(y),\varphi_{l}(y)\big{)}+\Big{(}a^{(4)}_{i}(t)\nabla_{y_{i}}\varphi_{k}(y),\varphi_{l}(y)\Big{)}\Big{]}$
(11)
$\displaystyle\quad+d_{k}(t)\Big{[}\big{(}a^{(1)}_{i}(t)\nabla_{y_{i}}\varphi_{k},\nabla_{y_{i}}\varphi_{l}\big{)}-\big{(}a^{(2)}_{ij}(t)\nabla_{y_{i}}\varphi_{k},\nabla_{y_{j}}\varphi_{l}\big{)}+\big{(}a^{(5)}_{i}(t)\nabla_{y_{i}}\varphi_{k},\varphi_{l}\big{)}\Big{]}\Big{\\}}=0.$

In order to validate the method (11), we add the term source
$\big{(}f(y,t),w\big{)}$ on the right side of the equation (9) and, taking
$w=\varphi_{l}$, we have the following definitions for the matrices and
vectors:

$\displaystyle A_{kl}=\Big{(}\varphi_{k}(y),\varphi_{l}(y)\Big{)},\quad
K_{1(kl)}=\Big{(}\nabla\varphi_{k}(y),\nabla\varphi_{l}(y)\Big{)},\quad
K_{2(kl)}=\Big{(}\Delta\varphi_{k}(y),\Delta\varphi_{l}(y)\Big{)},$ (12)
$\displaystyle
B_{1(ikl)}(t)=\Big{(}a_{i}^{(1)}(t)\nabla_{y_{i}}\varphi_{k}(y),\nabla_{y_{i}}\varphi_{l}(y)\Big{)},\quad
B_{2(ijkl)}(t)=\Big{(}a_{ij}^{(2)}(t)\nabla_{y_{i}}\varphi_{k}(y),\nabla_{y_{j}}\varphi_{l}(y)\Big{)},$
$\displaystyle
B_{3(ikl)}(t)=\Big{(}a_{i}^{(4)}(t)\nabla_{y_{i}}\varphi_{k}(y),\varphi_{l}(y)\Big{)},\quad
B_{4(ikl)}(t)=\Big{(}a_{i}^{(5)}(t)\nabla_{y_{i}}\varphi_{k}(y),\varphi_{l}(y)\Big{)},$
$\displaystyle F_{l}(t)=\big{(}f(t),\varphi_{l}(y)\big{)},\quad
d_{0(l)}=\big{(}v_{0}(y),\varphi_{l}(y)\big{)},\quad
d_{1(l)}=\big{(}v_{1}(y),\varphi_{l}(y)\big{)}.$

Considering the definitions (12) and varying $l=1,2,\cdots,m$, we obtain the
following nonlinear system of second order ordinary differential equations:

$\left\\{\begin{aligned}
&Ad^{\prime\prime}(t)+G\big{(}t,d\big{)}\leavevmode\nobreak\
\\!K_{1}\leavevmode\nobreak\ \\!d(t)+L_{1}(t)\leavevmode\nobreak\
\\!d^{\prime}(t)+L_{2}(t)\leavevmode\nobreak\ \\!d(t)=F(t),\\\
&d(0)=d_{0}\quad\mbox{and}\quad d^{\prime}(0)=d_{1},\end{aligned}\right.$ (13)

with

$\displaystyle
G\big{(}t,d\big{)}=b_{1}(t)\Big{|}\sum_{k=1}^{m}d_{k}(t)\varphi_{k}\Big{|}_{0}^{2},\quad
L_{1}(t)=\nu A+\widehat{B}_{3(i)}(t),$ (14) $\displaystyle
L_{2}(t)=b_{2}(t)K_{2}+\Big{(}\widehat{B}_{1(i)}(t)+\widehat{B}_{4(i)}(t)-\widehat{B}_{2(ij)}(t)\Big{)},$

where we denote by $\widehat{B}$ the transpose of the matrix $B$.

Now, to solve the system of ordinary nonlinear equations (13) in discrete
time, we will use the finite difference method.

###  3.2 The Newmark’s method

In order, consider the discretize the fixed time interval $[0,T]$, using a
uniform time step of size $\Delta t=T/N$, then $t_{\eta}=\eta\Delta t$, for
$\eta=0,1,\cdots,N$ and $\chi^{\eta}=\chi(\cdot,t_{\eta})$. Consider the
Neumark’s approximations:

$\chi^{\eta+\theta}=\theta\chi^{\eta-1}+(1-2\theta)\chi^{\eta}+\theta\chi^{\eta+1}.$

Substituting the approximations in the system (13) and taking the equation in
the discrete time $t=t_{\eta}=\eta\Delta t$ , we obtain

$\displaystyle\Big{(}M_{1}^{\eta+1}+\theta\leavevmode\nobreak\ \\!(\Delta
t)^{2}G^{\eta+1}(d)K_{1}\Big{)}d^{\eta+1}+M_{2}^{\eta}d^{\eta}$ (15)
$\displaystyle\qquad+\Big{(}M_{3}^{\eta-1}+\theta\leavevmode\nobreak\
\\!(\Delta t)^{2}G^{\eta-1}(d)K_{1}\Big{)}d^{\eta-1}=(\Delta
t)^{2}F^{\eta+\theta},$

or equivalently

$\displaystyle\Big{(}M_{1}^{\eta+1}+\theta\leavevmode\nobreak\ \\!(\Delta
t)^{2}G^{\eta+1}(d)K_{1}\Big{)}d^{\eta+1}=-M_{2}^{\eta}d^{\eta}-\Big{(}M_{3}^{\eta-1}+\theta\leavevmode\nobreak\
\\!(\Delta t)^{2}G^{\eta-1}(d)K_{1}\Big{)}d^{\eta-1}+(\Delta
t)^{2}F^{\eta+\theta},$ (16)

where

$\displaystyle M_{1}^{\eta+1}$ $\displaystyle=A+\frac{\Delta
t}{2}L_{1}^{\eta+1}+\theta(\Delta t)^{2}\leavevmode\nobreak\
\\!L_{2}^{\eta+1},\quad M_{2}^{\eta}=(\Delta
t)^{2}(1-2\theta)\big{(}G^{\eta}(d)K_{1}+L_{2}^{\eta}\big{)}-2A,$
$\displaystyle M_{3}^{\eta-1}$ $\displaystyle=A-\frac{\Delta
t}{2}L_{1}^{\eta-1}+\theta(\Delta t)^{2}\leavevmode\nobreak\
\\!L_{2}^{\eta-1},\quad F^{\eta+\theta}=\theta
F^{\eta-1}+(1-2\theta)F^{\eta}+\theta F^{\eta+1}.$

#### Remark.

For the initial step, that is, for $\eta=0$, appears in the equation the term
$d^{-1}$, that does not exist. Proceeding formally, we use a quadratic order
approximation for the initial data $d^{\prime}(0)$. Therefore, we consider
$d^{-1}=d^{1}-2\Delta td^{\prime}(0)$, so we have $d^{\eta+\theta}=2\theta
d^{\eta+1}+(1-2\theta)d^{\eta}-2\Delta t\theta d^{\prime}(0)$. So, for
$\eta=0$ in (16), we have

$\Big{[}M_{1}^{1}+M_{3}^{0}+\theta\leavevmode\nobreak\ \\!(\Delta
t)^{2}\Big{(}G^{1}(d)+G^{-1}(d)\Big{)}K_{1}\Big{]}d^{1}=-M_{2}^{0}d^{0}+2\Delta
t\leavevmode\nobreak\ \\!M_{3}^{0}d_{1}+(\Delta t)^{2}F^{0+\theta},$ (17)

where we assume the approximation $M_{3}^{-1}=M_{3}^{0},\quad
F^{0+\theta}=\theta F^{1}+(1-\theta)F^{0}\leavevmode\nobreak\ $ and

$G^{-1}(X)=b_{1}^{0}\Big{|}\sum_{k=1}^{m}\big{(}d_{k}^{1}-2\Delta
t\leavevmode\nobreak\ d_{1(k)}\big{)}\nabla\varphi_{k}(y)\Big{|}_{0}^{2},$

considering the approximation $b_{1}^{-1}\approx b_{1}^{0}$.

To determine the solution vector
$d^{\eta}=(d_{1}^{\eta},d_{2}^{\eta},\cdots,d_{m}^{\eta})$, we must now solve
the algebraic system (17) and (16), for $\eta=0,1,\cdots,N$. As the system is
not linear, we will employ the Newton Method, which we briefly describe below,
aiming at an optimization of the computational cost in the calculations.

###  3.3 Newton’s method

For each $\eta=1,2,\cdots,N$, consider $\mathcal{F}(X)$ the left side of the
equation in (17). In this way, solve (17) is equivalent to finding the root of
$\mathcal{F}(X)$. By Newton’s method, given the initial approximation and
denoting $X=d^{n+1}$,, we obtain $X_{k+1}$ a new approximate solution of
$\mathcal{F}(X)=0$, given by $X_{k+1}=X_{k}+s_{k}$, where $s_{k}$ is the
linear system solution $J\mathcal{F}(X_{k})s_{k}=-\mathcal{F}(X_{k})$. The key
point of the method is the calculation of the Jacobian matrix
$J\mathcal{F}(X)$. This is because, if it has a high computational cost, it
can make the method unfeasible.

Define

$\mathcal{F}^{\eta}(X)=\Big{(}M_{1}^{\eta+1}+\theta(\Delta
t)^{2}G^{\eta+1}(X)K_{1}\Big{)}X+\Gamma^{\eta},$ (18)

where $\leavevmode\nobreak\
\Gamma^{\eta}=M_{2}^{\eta}d^{\eta}+\Big{(}M_{3}^{\eta-1}+\theta\leavevmode\nobreak\
\\!(\Delta t)^{2}G^{\eta-1}(d)K_{1}\Big{)}d^{\eta-1}-(\Delta
t)^{2}F^{\eta+\theta}$.

For $\eta=0$, follows

$\mathcal{F}^{0}(X)=\Big{[}M_{1}^{1}+M_{3}^{0}+\theta(\Delta
t)^{2}\big{(}G^{1}(X)+G^{-1}(X)\big{)}K_{1}\Big{]}X+2\theta(\Delta
t)^{3}G^{-1}(X)K_{1}d_{1}+\Gamma^{0},$ (19)

where $\leavevmode\nobreak\ \Gamma^{0}=M_{2}^{0}d^{0}-2\Delta
t\leavevmode\nobreak\ \\!M_{3}^{0}d_{1}-(\Delta t)^{2}F^{0+\theta}$.

For $\eta>0$:

$\displaystyle J{\mathcal{F}}^{\eta}(X)=\frac{d\leavevmode\nobreak\
\\!\mathcal{F}^{\eta}\\!\\!}{dX}(X)$
$\displaystyle=M_{1}^{\eta+1}+\theta(\Delta
t)^{2}\frac{d}{dX}\Big{[}G^{\eta+1}(X)K_{1}X\Big{]}$ (20)
$\displaystyle=M_{1}^{\eta+1}+\theta(\Delta
t)^{2}\Big{(}K_{1}X\leavevmode\nobreak\
\\!\frac{dG^{\eta+1}}{dX}(X)+G^{\eta+1}(X)K_{1}\Big{)},$

and for $\eta=0$:

$\displaystyle J{\mathcal{F}}^{0}(X)$
$\displaystyle=M_{1}^{1}+M_{3}^{0}+\theta(\Delta
t)^{2}\Big{[}K_{1}X\leavevmode\nobreak\
\\!\Big{(}\frac{dG^{1}}{dX}(X)+\frac{dG^{-1}}{dX}(X)\Big{)}+\Big{(}G^{1}(X)+G^{-1}(X)\Big{)}K_{1}\Big{]}$
(21) $\displaystyle\quad+2\theta(\Delta t)^{3}\leavevmode\nobreak\
\\!K_{1}d_{1}\leavevmode\nobreak\ \\!\frac{dG^{-1}}{dX}(X),$

where

$\displaystyle\frac{dG^{\eta}}{dX}(X)=\Big{[}\frac{\partial G^{\eta}}{\partial
X_{1}}(X)\leavevmode\nobreak\ \leavevmode\nobreak\ \cdots\leavevmode\nobreak\
\leavevmode\nobreak\ \frac{\partial G^{\eta}}{\partial
X_{m}}(X)\Big{]}_{m\times
1},\quad\forall\eta\in\\{-1\\}\cup\\{1,\cdots,N-1\\}.$

Using (14), we obtain

$\displaystyle\frac{\partial G^{\eta}}{\partial
X_{k}}(X)=2\leavevmode\nobreak\ \\!b_{1}^{\eta}\leavevmode\nobreak\
X_{k}\leavevmode\nobreak\
\Big{[}\int_{\Omega}\big{(}\nabla_{y_{i}}\varphi_{k}(y)\big{)}^{2}\leavevmode\nobreak\
dx\Big{]},$ $\displaystyle\frac{\partial G^{-1}}{\partial
X_{k}}(X)=2\leavevmode\nobreak\ \\!b_{1}^{0}\leavevmode\nobreak\
\big{(}X_{k}-2\Delta t\leavevmode\nobreak\
\\!d_{1k}\big{)}\leavevmode\nobreak\
\Big{[}\int_{\Omega}\big{(}\nabla_{y_{i}}\varphi_{k}(y)\big{)}^{2}\leavevmode\nobreak\
dx\Big{]},$

where $d_{1k}$ is the $k$-th coordinate of the initial velocity vector
$d_{1}$.

From there, we apply Newton’s refinement algorithm, for $\eta=1,2,\ldots,N-1$,
until one of the stopping criteria is satisfied:
$\|X_{k+1}-X_{k}\|_{\infty}<10^{-14}$ or $|\mathcal{F}(X_{k})|<10^{-14}$,
where $X_{k}$ is the solution obtained in the $i$-th iteration of Newton. We
also define a maximum limit for iterations in order to avoid infinite
refinements.

###  3.4 Finite Element Method- Hermite Base

Let $[T_{h}]$ be a polygonalization family of $\Omega$, where
$T_{h}=\\{K;\leavevmode\nobreak\ K\subset\Omega\\}$, satisfying the minimum
angle condition (see [7]) and indexed by the parameter $h$, representing the
maximum diameter of the elements $K\in T_{h}$. Given an integer $l\geq 1$, we
introduce the following finite element space $V_{h}$ spanned by the Hermite
interpolation polynomials basis functions. More specifically, we let

$\displaystyle V_{h}=\Big{\\{}\psi_{h}\in V_{h}^{3}\cap
C^{1}(\Omega);\leavevmode\nobreak\ \psi_{h}\big{|}_{K}\in
P_{l}(K),\leavevmode\nobreak\ \forall\leavevmode\nobreak\ K\in
T_{h}\Big{\\}}\subset V,$

where $V^{3}_{h}$ represents a finite-dimensional function space generated by
piecewise polynomial $P_{l}(K)$ functions of degree less or equal $3$
associated with a discretization by finite elements of size $h$. In
particular, in this paper we will use the Hermite polynomials of degree equal
$l=3$.

From the Lemma of Douglas-Dupont (see in [7]), it follows that, given a
function $w:]0,T[\longrightarrow H^{l+1}(\Omega)$, there is an interpolator
$\widehat{w}_{h}:]0,T[\longrightarrow V_{h}^{l}$ such that

$\displaystyle\|w(t)-\widehat{w}_{h}(t)\|_{m}\leq C_{1}\leavevmode\nobreak\
\\!h^{l+1-m}\|w(t)\|_{l+1},$

where $\|\cdot\|_{m}$ denotes the seminorm over Hilbert space $H^{m}(\Omega)$
and $0\leq m\leq l$.

Knowing the Hermite polynomials, then all the matrices of the system (16) can
be calculated explicitly and solving the nonlinear system by Newton’s method
gives the approximate solution, given by

$v_{h}(y,t)=\sum_{i=1}^{m}d_{i}(t)\varphi_{i}(y),\quad\varphi_{i}(y)\in
V_{h}.$

We note that since the Hermite polynomials are orthogonal then
$v^{n}_{h}(y_{i})=d_{i}^{n}$

##  4 Numerical analysis

In this section we present error estimates for the finite element semi-
discrete and fully discrete approximations of the nonlinear beam equation (8)
or equivalently for the weak formulation of the nonlinear beam equation with
moving ends.

Remarks on regularity: In order to get an error estimate, we require more
regularity for the function $v=v(y,t)$ and its derivatives in time:

$v^{\prime},v^{\prime\prime}\in L^{2}(0,T;H^{4}(\Omega)),\leavevmode\nobreak\
\leavevmode\nobreak\ v^{\prime\prime\prime}\in
L^{\infty}(0,T;H^{1}(\Omega))\leavevmode\nobreak\ \leavevmode\nobreak\
\mbox{and}\leavevmode\nobreak\ \leavevmode\nobreak\ v^{(iv)}\in
L^{\infty}(0,T;L^{2}(\Omega)).$ (H2)

###  4.1 Error estimates for the semi-discrete problem

The semi-discrete problem consists of finding $v_{h}:\leavevmode\nobreak\
[0,T[\longrightarrow V_{h}$ such that satisfies the problem (8) restricted to
$w=w_{h}\in V_{h}$, that is, to find $v_{h}(t)$ for all
$t\in\leavevmode\nobreak\ \\!]0,T[$ such that

$\left\\{\begin{aligned} &\begin{aligned}
&\big{(}v_{h}^{\prime\prime}(t),w_{h}\big{)}+b_{1}(t)\big{|}\nabla
v_{h}(t)\big{|}_{0}^{2}\big{(}\nabla v_{h}(t),\nabla
w_{h}\big{)}+b_{2}(t)\big{(}\Delta v_{h}(t),\Delta
w_{h}\big{)}+\nu\leavevmode\nobreak\
\\!\big{(}v_{h}^{\prime}(t),w_{h}\big{)}\\\
&\quad+\big{(}a^{(1)}_{i}(t)\nabla_{y_{i}}v_{h}(t),\nabla_{y_{i}}w_{h}\big{)}+\big{(}a^{(2)}_{ij}(t)\nabla_{y_{i},y_{j}}v_{h}(t),w_{h}\big{)}\\\
&\quad+\big{(}a^{(4)}_{i}(t)\nabla_{y_{i}}v_{h}^{\prime}(t),w_{h}\big{)}+\big{(}a^{(5)}_{i}(t)\nabla_{y_{i}}v_{h}(t),w_{h}\big{)}=0,\mbox{
in }Q,\end{aligned}\\\ &\leavevmode\nobreak\
\big{(}v_{h}(y,0),w_{h}\big{)}=\big{(}v_{0h}(y),w_{h}\big{)},\leavevmode\nobreak\
\leavevmode\nobreak\
\big{(}v_{h}^{\prime}(y,0),w_{h}\big{)}=\big{(}v_{1h}(y),w_{h}\big{)}\quad\mbox{in
}\leavevmode\nobreak\ \Omega,\\\ &\leavevmode\nobreak\
\big{(}v_{h},w_{h}\big{)}=0\leavevmode\nobreak\ \mbox{and}\leavevmode\nobreak\
\big{(}\nabla v_{h},w_{h}\big{)}\equiv 0,\leavevmode\nobreak\
\mbox{on}\leavevmode\nobreak\ \partial\Omega\times[0,T[.\end{aligned}\right.$
(22)

For the estimates, we will need the following hypothesis about the projections
of $V$ in $V_{h}$, denoted by $R_{h}$ and defined in (27), of the initial data

$|v_{1h}-R_{h}v_{1}|_{0}\leavevmode\nobreak\ \leq\leavevmode\nobreak\
c_{1}h^{4},\quad|\nabla v_{0h}-R_{h}\nabla v_{0}|_{0}\leavevmode\nobreak\
\leq\leavevmode\nobreak\ c_{2}h^{3}\quad\mbox{and}\quad|\Delta
v_{0h}-R_{h}\Delta v_{0}|_{0}\leavevmode\nobreak\ \leq\leavevmode\nobreak\
c_{3}h^{2},$ (H3)

where $R_{h}\Delta v_{0},\leavevmode\nobreak\ R_{h}\nabla v_{0}$ and
$R_{h}v_{1}$ are the projections of $\Delta v_{0},\leavevmode\nobreak\ \nabla
v_{0}$ and $v_{1}$ in $V_{h}$, respectively, and $c_{1},c_{2},c_{3}\geq 0$ are
constants independent of $h$.

Let us assume that,

$\displaystyle\max_{t\in[0,T]}(K^{\prime}(t))^{2}<\frac{\zeta_{0}}{4}$ (H4)

which requires that the speed of the end points be smaller than the
characteristic speed of the equation.

Let us denote by $e(t)=v(t)-v_{h}(t),\forall t\in]0,T[$, the error between the
exact and approximate solution. Then, we have the following error estimate:

######  Theorem 4.

(Semi-discrete problem - Estimate)

Let $v=v(y,t)$ the solution of the problem (8). Under hypotheses H1– (H4) and
given the initial data $v_{0},v_{1}\in H_{0}^{2}(\Omega)\cap H^{4}(\Omega)$,
there exists a constant $C>0$, such that

$\|e^{\prime}\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\|e\|_{L^{\infty}(0,T;H^{2}_{0}(\Omega))}\leavevmode\nobreak\
\leq\leavevmode\nobreak\ Ch^{2}.$ (23)

where the constant $C$ is dependent on $v$ and independent of $h$.

###### Proof.

We note that, as $V_{h}\subset V=H_{0}^{2}(\Omega)$, taking in particular
$w=w_{h}$ in (8)1 and doing the difference with (22)1, denoting $e=e(t)$,
$v=v(t)$ and $v_{h}=v_{h}(t)$ , we obtain:

$\displaystyle\big{(}e^{\prime\prime},w_{h}\big{)}-b_{1}(t)\big{(}\big{[}\big{|}\nabla
v\big{|}_{0}^{2}\nabla v-\big{|}\nabla v_{h}\big{|}_{0}^{2}\nabla
v_{h}\big{]},\nabla w_{h}\big{)}+b_{2}(t)\big{(}\Delta{e},\Delta w_{h}\big{)}$
(24) $\displaystyle\quad+\nu\leavevmode\nobreak\
\\!\big{(}e^{\prime},w_{h}\big{)}+\big{(}a^{(1)}_{i}(t)\nabla_{y_{i}}e\leavevmode\nobreak\
,\nabla_{y_{i}}w_{h}\big{)}+\big{(}a^{(2)}_{ij}(t)\nabla_{y_{i},y_{j}}e\leavevmode\nobreak\
,w_{h}\big{)}$
$\displaystyle\quad+\big{(}a^{(4)}_{i}(t)\nabla_{y_{i}}e^{\prime},w_{h}\big{)}+\big{(}a^{(5)}_{i}(t)\nabla_{y_{i}}e\leavevmode\nobreak\
,w_{h}\big{)}=0,\leavevmode\nobreak\ \forall\leavevmode\nobreak\ w_{h}\in
V_{h}.$

We use the standard decomposition of the error

$e(t)=v(t)-v_{h}(t)=(v(t)-R_{h}v(t))+(R_{h}v(t)-v_{h}(t))=\rho(t)+\xi(t)=\rho+\xi$
(25)

where $\rho=v(t)-R_{h}v(t)$ and $\xi=R_{h}v(t)-v_{h}(t)$, with $R_{h}v(t)$
being the Ritz projection of $v(t)\in V_{h}$, to be defined as follows:

To define the Ritz projection $R_{h}$ we will follow the analogue of the ideas
of the elliptic approximation theory. Consider the following definitions.

######  Definition 1.

(Bilinear form)

Let the functional $a^{(\vartheta,t)}:V\times V\rightarrow\mathbb{R}$, for all
$\vartheta\in\\{v,v_{h}\\}$ and $t\in[0,T[$ both fixed, such that

$a^{(\vartheta,t)}(v,w)=(a^{(1)}_{i}(t)\nabla_{y_{i}}v\leavevmode\nobreak\
,\nabla_{y_{i}}w)+b_{1}(t)|\nabla\vartheta(t)|_{0}^{2}(\nabla v,\nabla
w)+b_{2}(t)(\Delta{v},\Delta{w}).$ (26)

######  Definition 2.

(Ritz Projection on $a^{(\vartheta,t)}(\cdot,\cdot)$) We define the Ritz
projection on the functional $a^{(\vartheta,t)}(\cdot,\cdot)$, defined in
(26), as follows:

$\displaystyle R_{h}:\leavevmode\nobreak\ V$ $\displaystyle\longrightarrow
V_{h}$ $\displaystyle v$ $\displaystyle\longmapsto R_{h}v$

satisfying the orthogonality condition

$a^{(\vartheta,t)}\big{(}v(t)-R_{h}v(t),w_{h}\big{)}=0,\leavevmode\nobreak\
\forall w_{h}\in V_{h}.$ (27)

For the Ritz projection to satisfy the equation (27), it is necessary that the
functional (26) satisfies the following Lemma:

######  Lemma 1.

(Functional Bilinear, Continuous and Coercive)

The functional $a^{(\vartheta,t)}(\cdot,\cdot)$, defined by (26), is bilinear,
continuous and coercive, for all $t\in[0,T[$ and $\vartheta\in\\{v,v_{h}\\}$.

###### Proof.

The proof is given in three steps:

  1. 1.

The bilinearity of the functional $a^{(\vartheta,t)}(\cdot,\cdot)$ is
immediate.

  2. 2.

The coercivity of the functional $a^{(\vartheta,t)}(\cdot,\cdot)$:

$\displaystyle a^{(\vartheta,t)}(w,w)$
$\displaystyle=\big{(}a^{(1)}_{i}(t),(\nabla_{y_{i}}w)^{2}\big{)}+b_{1}(t)|\vartheta(t)|_{0}^{2}\leavevmode\nobreak\
\\!|\nabla w|_{0}^{2}+b_{2}(t)|\Delta{w}|_{0}^{2}$ (28)
$\displaystyle\geq\kappa_{0}\leavevmode\nobreak\
\\!\big{[}|\nabla_{y_{i}}w|_{0}^{2}+|\nabla
w|_{0}^{2}+|\Delta{w}|_{0}^{2}\big{]}\geq c\leavevmode\nobreak\
\\!\|w\|_{V}^{2}.$

The constants are given by

$\displaystyle\kappa_{0}=\min\big{\\{}\min_{t\in[0,T]}b_{1}(t)|\vartheta(t)|_{0}^{2},\leavevmode\nobreak\
\min_{t\in[0,T]}b_{2}(t),\leavevmode\nobreak\
\min_{i\in\\{1,\cdots,n\\}}\big{[}\min_{t\in[0,T]}\big{(}\inf_{y\in\Omega}a^{(1)}_{i}(y,t)\big{)}\big{]}\big{\\}}.$

We note that, from (7) and $|\vartheta|_{0}^{2}\not\equiv 0$, once
$\vartheta\in\\{v,v_{h}\\}$ and $v,v_{h}$ are non-trivial solutions,

$\displaystyle\min_{t\in[0,T]}b_{1}(t)|\vartheta|_{0}^{2}$
$\displaystyle=\min_{t\in[0,T]}\frac{\zeta_{1}|\vartheta|_{0}^{2}}{K^{4}}>\frac{\zeta_{1}|\vartheta|_{0}^{2}}{K_{1}^{4}}>0$
(29) $\displaystyle\min_{t\in[0,T]}b_{2}(t)$
$\displaystyle=\min_{t\in[0,T]}\leavevmode\nobreak\ \leavevmode\nobreak\
\frac{1}{K^{4}}\leavevmode\nobreak\ \leavevmode\nobreak\
\\!>\leavevmode\nobreak\ \leavevmode\nobreak\
\\!\frac{1}{K_{1}^{4}}\leavevmode\nobreak\ \leavevmode\nobreak\
\\!>0\leavevmode\nobreak\ \\!.$

And from hypothesis (H4),

$\displaystyle\min_{t\in[0,T]}\big{(}\inf_{y\in\Omega}a^{(1)}_{i}(y,t)\big{)}$
$\displaystyle=\min_{t\in[0,T]}\Big{\\{}\frac{1}{K^{2}}\big{[}\zeta_{0}-4(K^{\prime})^{2}\big{]}\Big{\\}}>0.$
(30)

Therefore, by (29) and (30), $\kappa_{0}>0$, concluding the coercivity.

  3. 3.

The continuity of the functional $a^{(\vartheta,t)}(\cdot,\cdot)$:

$\displaystyle a^{(\vartheta,t)}(v,v)$
$\displaystyle\leq|(a^{(1)}_{i}(t)\nabla_{y_{i}}v\leavevmode\nobreak\
,\nabla_{y_{i}}w)|+b_{1}(t)|\nabla\vartheta(t)|_{0}^{2}|(\nabla v,\nabla
w)|+b_{2}(t)|(\Delta{v},\Delta{w})|$
$\displaystyle\leq\kappa_{1}\big{[}|\nabla_{y_{i}}v|_{0}\leavevmode\nobreak\
|\nabla_{y_{i}}w|_{0}+|\nabla v|_{0}|\nabla w|_{0}^{2}+|\Delta
v|_{0}\leavevmode\nobreak\ |\Delta w|_{0}\big{]}$ (31)
$\displaystyle\leq\leavevmode\nobreak\ c\leavevmode\nobreak\
\\!||v||_{V}||w||_{V},$

where we consider $|\nabla_{y_{i}}v|_{0}\leq|\nabla v|_{0}$ and apply the
inequalities of Cauchy-Schwarz and Poincare-Friedrich and the equivalence of
norms $|\Delta v|_{0}$ and $\|v\|_{2}$ in $H_{0}^{2}(\Omega)$. The constant
are given by

$\displaystyle\kappa_{1}=\max_{i\in\\{1,\cdots,n\\}}\big{\\{}\max_{t\in[0,T]}\big{[}\sup_{y\leavevmode\nobreak\
\\!\in\leavevmode\nobreak\
\\!\Omega}|a_{i}^{(1)}(y,t)|+b_{1}(t)(|v(t)|_{0}^{2}+|v_{h}(t)|_{0}^{2})+b_{2}(t)\big{]}\big{\\}}.$

Once $\\{v(t),v_{h}(t)\\}\subset V=H_{0}^{2}(\Omega)$, we have
$|v|_{0}+|v_{h}|_{0}^{2}<+\infty$.

Thus, we conclude that the functional $a^{(\vartheta,t)}(\cdot,\cdot)$, is
bilinear, continuous and coercive, $\forall t\in[0,T[$ and
$\vartheta\in\\{v,v_{h}\\}$. ∎

We return the proof of the theorem: Summing and subtracting the Ritz
projection of $v$ and its derivatives conveniently on each term of the
equation (24), using the decomposition of the error (25), and taking in
particular $w_{h}=\xi^{\prime}\in V_{h}$ and using the (27), we get:

$\Gamma_{1}(\xi,v_{h},t)+\Gamma_{2}(\xi,t)=-\big{[}\Gamma_{1}(\rho,v,t)+\Gamma_{2}(\rho,t)+\mu(t)\big{]},$
(32)

where $\Gamma_{1}:V\times\\{v,v_{h}\\}\times[0,T[\leavevmode\nobreak\
\longrightarrow\leavevmode\nobreak\ \mathbb{R}$, $\leavevmode\nobreak\
\Gamma_{2}:V\times[0,T[\leavevmode\nobreak\
\longrightarrow\leavevmode\nobreak\ \mathbb{R}\leavevmode\nobreak\ $ and
$\leavevmode\nobreak\ \mu:[0,T[\leavevmode\nobreak\
\longrightarrow\leavevmode\nobreak\ \mathbb{R}\leavevmode\nobreak\ $ such that

$\displaystyle\Gamma_{1}(\chi,\vartheta,t)$
$\displaystyle=\big{(}\chi^{\prime\prime},\xi^{\prime}\big{)}+a^{(\vartheta,t)}(\chi,\xi^{\prime})+\nu\leavevmode\nobreak\
\\!\big{(}\chi^{\prime},\xi^{\prime}\big{)}+\big{(}a^{(4)}_{i}(t)\nabla_{y_{i}}\chi^{\prime},\xi^{\prime}\big{)},$
$\displaystyle\Gamma_{2}(\chi,t)$
$\displaystyle=\big{(}a^{(2)}_{ij}(t)\nabla_{y_{i},y_{j}}\chi,\xi^{\prime}\big{)}+\big{(}a^{(5)}_{i}(t)\nabla_{y_{i}}\chi,\xi^{\prime}\big{)},$
$\displaystyle\mu(t)$ $\displaystyle=b_{1}(t)\big{[}|\nabla
v_{h}|_{0}^{2}-|\nabla v|_{0}^{2}\big{]}\big{(}R_{h}\nabla
v,\nabla\xi^{\prime}\big{)}.$

Next, we will analyze each term of the equation (32)

Analisys: Terms dependent of $\xi$

By definition, we have

$\displaystyle\Gamma_{1}(\xi,v_{h},t)=\big{(}\xi^{\prime\prime},\xi^{\prime}\big{)}+a^{(v_{h},t)}(\xi,\xi^{\prime})+\nu\leavevmode\nobreak\
\\!\big{(}\xi^{\prime},\xi^{\prime}\big{)}+\big{(}a^{(4)}_{i}(t)\nabla_{y_{i}}\xi^{\prime},\xi^{\prime}\big{)}.$

Analysing each term, we have

$\big{(}\xi^{\prime\prime},\xi^{\prime}\big{)}=\frac{1}{2}\frac{d}{dt}|\xi^{\prime}|_{0}^{2},\qquad\nu\big{(}\xi^{\prime},\xi^{\prime}\big{)}=\nu\leavevmode\nobreak\
\\!|\xi^{\prime}|_{0}^{2}\geq 0.$ (33)

Integrating by parts and using the definition (7), follows that

$\big{(}a^{(4)}_{i}(t)\nabla_{y_{i}}\xi^{\prime},\xi^{\prime}\big{)}=-\frac{1}{2}\int_{\Omega}\big{(}\nabla_{y_{i}}a^{(4)}_{i}(y,t)\big{)}\big{(}\xi^{\prime}\big{)}^{2}dy=\frac{K^{\prime}(t)}{K(t)}|\xi^{\prime}|_{0}^{2}\geq
0,$ (34)

since the
$\displaystyle\nabla_{y_{i}}a^{(4)}_{i}(y,t)=-\frac{2K^{\prime}(t)}{K(t)}$ and
$K^{\prime}(t)>0$, for the hypotesis (H1).

Now using the definition (7) and the mean value theorem for integrals, we
obtain

$\displaystyle a^{(v_{h},t)}(\xi,\xi^{\prime})$
$\displaystyle=\frac{1}{2}\leavevmode\nobreak\
\Big{[}\int_{\Omega}a^{(1)}_{i}(t)\frac{d}{dt}(\nabla_{y_{i}}\xi)^{2}\leavevmode\nobreak\
dy\Big{]}+\frac{b_{1}(t)}{2}|v_{h}|_{0}^{2}\leavevmode\nobreak\
\\!\frac{d}{dt}|\nabla\xi|_{0}^{2}+\frac{b_{2}(t)}{2}\frac{d}{dt}|\Delta\xi|_{0}^{2}$
(35)
$\displaystyle=\frac{1}{2}\frac{d}{dt}\big{[}a^{(1)}_{i}(\sigma_{1},t)\leavevmode\nobreak\
\\!|\nabla_{y_{i}}\xi|_{0}^{2}+b_{1}(t)|v_{h}|_{0}^{2}|\nabla\xi|_{0}^{2}+b_{2}(t)|\Delta\xi|_{0}^{2}\big{]}$
$\displaystyle\qquad-\frac{1}{2}\leavevmode\nobreak\
\big{[}a_{i}^{\prime}\leavevmode\nobreak\
\\!{}^{(1)}(\sigma_{2},t)\leavevmode\nobreak\
\\!\big{|}\nabla_{y_{i}}\xi\big{|}_{0}^{2}+\frac{d}{dt}\big{(}b_{1}(t)|v_{h}|_{0}^{2}\big{)}|\nabla\xi|_{0}^{2}+b^{\prime}_{2}(t)|\Delta\xi|_{0}^{2}\big{]},$

where $\sigma_{1},\sigma_{2}\in\Omega$. So, by (33) – (35) and decreasing the
term $\displaystyle\frac{\kappa_{2}}{2}|\xi^{\prime}|_{0}^{2}$ :

$\displaystyle\Gamma_{1}(\xi,v_{h},t)$
$\displaystyle\geq\frac{1}{2}\frac{d}{dt}\big{[}|\xi^{\prime}|_{0}^{2}+a^{(1)}_{i}(\sigma_{1},t)\leavevmode\nobreak\
\\!|\nabla_{y_{i}}\xi|_{0}^{2}+b_{1}(t)|v_{h}|_{0}^{2}|\nabla\xi|_{0}^{2}+b_{2}(t)|\Delta\xi|_{0}^{2}\big{]}-\frac{\kappa_{2}}{2}E(t),$
(36)

where

$\displaystyle
E(t)=|\nabla_{y_{i}}\xi|_{0}^{2}+|\nabla\xi|_{0}^{2}+|\Delta\xi|_{0}^{2}+|\xi^{\prime}|_{0}^{2},$

and
$\displaystyle\kappa_{2}=\max_{i,j\in\\{1,\cdots,n\\}}\big{\\{}\max_{t\in[0,T]}\big{[}\sup_{y\leavevmode\nobreak\
\\!\in\leavevmode\nobreak\
\\!\Omega}\big{(}|a_{i}^{\prime}\leavevmode\nobreak\
\\!{}^{(1)}(y,t)|+\Big{|}\frac{d}{dt}\big{(}b_{1}(t)|v_{h}(t)|_{0}^{2}\big{)}\Big{|}+|b^{\prime}_{2}(t)|\big{)}\big{]}\big{\\}}$.

We note that, by definition (7) and applying the inequality of Cauchy-Schwarz,

$\displaystyle\Big{|}\frac{d}{dt}\big{(}b_{1}(t)|v_{h}|_{0}^{2}\big{)}\Big{|}$
$\displaystyle=b^{\prime}_{1}(t)|v_{h}|_{0}^{2}+b_{1}(t)\frac{d}{dt}\big{(}v_{h},v_{h}\big{)}\leq
c\big{[}|v_{h}|_{0}^{2}+|v_{h}|_{0}|v^{\prime}_{h}|_{0}\big{]}<+\infty$ (37)

where $\displaystyle c=\frac{4\leavevmode\nobreak\
\\!\zeta_{1}}{K_{0}^{5}}(K_{2}+K_{0})$. By Theorem (2), we have
$v^{\prime}_{h}(t)\in L^{2}(\Omega)$. Soon, by (37), we ensure $\kappa_{2}$
limited.

Now, analyzing each term of $\Gamma_{2}(\xi,t)$, give by

$\displaystyle\Gamma_{2}(\xi,t)=\big{(}a^{(2)}_{ij}(t)\nabla_{y_{i},y_{j}}\xi,\xi^{\prime}\big{)}+\big{(}a^{(5)}_{i}(t)\nabla_{y_{i}}\xi,\xi^{\prime}\big{)},$

we have, using the inequality of Cauchy-Schwarz and Young, the definition of
$|\Delta\xi|_{0}$ and increasing the norm $|\nabla\xi|_{0}^{2}$, follows that

$\displaystyle\Gamma_{2}(\xi,t)\leq\frac{\kappa_{3}}{2}E(t),$ (38)

where $\kappa_{3}=\displaystyle
2\\!\\!\\!\max_{i,j\in\\{1,\cdots,n\\}}\big{\\{}\max_{t\in[0,T]}\big{[}\sup_{y\leavevmode\nobreak\
\\!\in\leavevmode\nobreak\
\\!\Omega}\big{(}|a_{ij}^{(2)}(y,t)|+|a_{i}^{(4)}(y,t)|+|a_{i}^{(5)}(y,t)|\big{)}\big{]}\big{\\}}$.

Analisys: Terms dependent on $\rho$

By definition, we have

$\displaystyle\Gamma_{1}(\rho,v,t)+\Gamma_{2}(\rho,t)$
$\displaystyle=\big{(}\rho^{\prime\prime},\xi^{\prime}\big{)}+a^{(v,t)}(\rho,\xi^{\prime})+\nu\leavevmode\nobreak\
\\!\big{(}\rho^{\prime},\xi^{\prime}\big{)}+\big{(}a^{(4)}_{i}(t)\nabla_{y_{i}}\rho^{\prime},\xi^{\prime}\big{)}$
$\displaystyle\qquad+\big{(}a^{(2)}_{ij}(t)\nabla_{y_{i},y_{j}}\rho,\xi^{\prime}\big{)}+\big{(}a^{(5)}_{i}(t)\nabla_{y_{i}}\rho,\xi^{\prime}\big{)}.$

Analysing the terms of $\Gamma_{1}(\rho,v,t)$ and $\Gamma_{2}(\rho,t)$
applying the inequality of Cauchy-Schwarz and Young, follows that

$\displaystyle\Gamma_{1}(\rho,v,t)+\Gamma_{2}(\rho,t)$
$\displaystyle\leq\frac{\kappa_{4}}{2}\big{[}|\rho^{\prime\prime}|_{0}^{2}+|\rho^{\prime}|_{0}^{2}+|\nabla\rho^{\prime}|_{0}^{2}+|\nabla\rho|_{0}^{2}+|\Delta\rho|_{0}^{2}+|\xi^{\prime}|_{0}^{2}\big{]}$
(39)
$\displaystyle\leq\frac{\kappa_{4}}{2}\big{[}|\rho^{\prime\prime}|_{0}^{2}+\|\rho^{\prime}\|_{1}^{2}+\|\rho\|_{2}^{2}+|\xi^{\prime}|_{0}^{2}\big{]},$

once
$|\rho^{\prime}|_{0}^{2}+|\nabla\rho^{\prime}|_{0}^{2}\leq\|\rho^{\prime}\|_{1}^{2}$
and $|\nabla\rho|_{0}^{2}+|\Delta\rho|_{0}^{2}\leq\|\rho\|_{2}^{2}$, and where

$\kappa_{4}=5(\kappa_{1}+\kappa_{2}+\kappa_{3}+1+\nu).$

We note that by definition of Ritz projection in (27), we get
$a^{(v,t)}(\rho,\xi^{\prime})=a^{(v,t)}(v-R_{h}v,\xi^{\prime})=0$.

Analisys: Nonlinear term

For the term $\mu(t)$, after integrating by parts, applying the module and the
inequalities of Cauchy-Schwarz and Young, we get

$\displaystyle\mu(t)\leq b_{1}(t)\leavevmode\nobreak\ \\!|\nabla v-\nabla
v_{h}|_{0}\big{[}|\nabla v|_{0}+|\nabla v_{h}|_{0}\big{]}\big{|}\big{(}\nabla
R_{h}\nabla v,\xi^{\prime}\big{)}\big{|}$ (40)
$\displaystyle\qquad\\!\leq\leavevmode\nobreak\
\kappa_{1}\big{[}|\nabla\rho|_{0}+|\nabla\xi|_{0}\big{]}\leavevmode\nobreak\
\\!|\nabla R_{h}\nabla v|_{0}|\xi^{\prime}|_{0}\leavevmode\nobreak\
\leq\frac{\kappa_{5}}{2}\big{[}|\nabla\rho|_{0}^{2}+|\nabla\xi|_{0}^{2}+|\xi^{\prime}|_{0}^{2}\big{]},$

where $\kappa_{5}=\kappa_{4}+2\leavevmode\nobreak\
\\!\kappa_{1}\leavevmode\nobreak\ \\!|\nabla R_{h}\nabla v|_{0}<+\infty$, once
$R_{h}\nabla v\in V_{h}\subset V=H_{0}^{2}(\Omega)$ implies $\nabla
R_{h}\nabla v\in L^{2}(\Omega)$.

Final estimate

Now, passing the term $\Gamma_{2}(\xi,t)$ to the right side of (32) and using
the estimates (36) – (40), we obtain

$\frac{1}{2}\Big{[}\frac{d}{dt}\Theta_{1}(t)-\kappa_{2}\leavevmode\nobreak\
\\!E(t)\Big{]}\leq\Gamma_{1}(\xi,v_{h},t)=-\big{[}\Gamma_{1}(\rho,v,t)+\Gamma_{2}(\rho+\xi,t)+\mu(t)\big{]}\leq\frac{\kappa_{5}}{2}\leavevmode\nobreak\
\\!\Theta_{3}(t)$ (41)

where

$\displaystyle\Theta_{1}(t)=|\xi^{\prime}|_{0}^{2}+a^{(1)}_{i}(\sigma_{1},t)\leavevmode\nobreak\
\\!|\nabla_{y_{i}}\xi|_{0}^{2}+b_{1}(t)|v_{h}|_{0}^{2}|\nabla\xi|_{0}^{2}+b_{2}(t)|\Delta\xi|_{0}^{2},$
$\displaystyle\Theta_{2}(t)=E(t)+|\rho^{\prime\prime}|_{0}^{2}+\|\rho^{\prime}\|_{1}^{2}+\|\rho\|_{2}^{2}.$

Passing the term $\kappa_{2}\leavevmode\nobreak\ \\!E(t)$ of left to the right
side of (41), integrating from $0$ to $t<T$ and multiplying by both side by
$2$, we get

$\kappa_{0}\leavevmode\nobreak\ \\!E(t)\leq\kappa_{1}\leavevmode\nobreak\
\\!E(0)+\kappa_{5}\Big{[}\int_{0}^{t}|\rho^{\prime\prime}|_{0}^{2}+\|\rho^{\prime}\|_{1}^{2}+\|\rho\|_{2}^{2}\leavevmode\nobreak\
dt+\int_{0}^{t}E(s)\leavevmode\nobreak\ dt\Big{]}.$ (42)

From hypothesis (H3), follows that

$E(0)\leq C_{1}h^{4},\quad\mbox{for some}\quad C_{1}>0.$ (43)

On the other hand, by the Céa Lemma and Douglas-Dupont Theorem, follows that

$\displaystyle\int_{0}^{t}|\rho^{\prime\prime}(s)|_{0}^{2}+\|\rho^{\prime}(s)\|_{1}^{2}+\|\rho(s)\|_{2}^{2}\leavevmode\nobreak\
ds\leq\leavevmode\nobreak\ C_{2}\leavevmode\nobreak\
\\!h^{4}\sum_{k=0}^{2}\Big{\|}\frac{d^{k}v}{dt^{k}}\Big{\|}^{2}_{L^{2}(0,T;H^{4}(\Omega))},\leavevmode\nobreak\
\mbox{for some}\leavevmode\nobreak\ C_{2}>0.$ (44)

Substituting the estimations (43) and (44) in (42) and dividing both sides by
$\kappa_{0}$, since $\kappa_{0}>0$ by (29) and (30), we obtain

$E(t)\leq C_{3}\leavevmode\nobreak\
\\!h^{4}+C_{3}\int_{0}^{t}E(s)\leavevmode\nobreak\ ds,$ (45)

where $\leavevmode\nobreak\
C_{3}=(\kappa_{5}/\kappa_{0})(C_{1}+C_{2})\leavevmode\nobreak\ \\!$. As the
estimate (45) holds for all $t\in[0,T[$, for Gronwall-Bellman’s Lemma, we have
to

$\frac{1}{2}\big{(}|\xi^{\prime}(t)|_{0}+|\Delta\xi(t)|_{0}\big{)}^{2}\leq|\xi^{\prime}(t)|_{0}^{2}+|\Delta\xi(t)|_{0}^{2}\leq
E(t)\leq C_{5}\leavevmode\nobreak\ \\!h^{4},\leavevmode\nobreak\
\mbox{where}\leavevmode\nobreak\ C_{5}=C_{3}e^{C_{3}T}.$ (46)

Therefore, from the decomposition of the error (25) and the estimates (44) and
(46), we obtain

$|e^{\prime}(t)|_{0}+|\Delta
e(t)|_{0}\leq|\rho^{\prime}(t)|_{0}+|\Delta\rho(t)|_{0}+|\xi^{\prime}(t)|_{0}+|\Delta\xi(t)|_{0}\leq
C\leavevmode\nobreak\ \\!h^{2},\leavevmode\nobreak\ \forall t\in[0,T[,$ (47)

for some $C>0$ independent of $t$ and $h$.

Finally, we just take (47) the essential supreme in time and using the
equivalence of norms $|\Delta e(t)|_{0}$ and $\|e(t)\|_{2}$ in
$H_{0}^{2}(\Omega)$ to obtain the estimate (23). ∎

###  4.2 Fully-Discrete Problem

In order to obtain an error estimate in the norm
$L^{\infty}(0,T);L^{2}(\Omega)$ for the fully-discrete problem, we discretize
the fixed time interval $[0,T]$, using the discretization defined in (3.2).

Let $\theta\in[0,1]$ and consider the notation $v^{\eta}=v(\cdot,t_{\eta})$,
for $\eta\in\\{1,\cdots,N-1\\}$ and the operators differences

$\displaystyle\delta^{2}v^{\eta}=\frac{v^{\eta+1}-2\leavevmode\nobreak\
\\!v^{\eta}+v^{\eta-1}}{(\Delta t)^{2}},\quad\delta
v^{\eta}=\frac{v^{\eta+1}-v^{\eta-1}}{2\leavevmode\nobreak\ \\!\Delta
t},\quad\delta v^{\eta+1/2}=\frac{v^{\eta+1}-v^{\eta}}{\Delta t},$ (48)
$\displaystyle\qquad\delta v^{\eta-1/2}=\frac{v^{\eta}-v^{\eta-1}}{\Delta
t},\quad v^{\eta+\theta}=\theta(v^{\eta+1}+v^{\eta-1})+(1-2\theta)v^{\eta}.$

Using the operators (48) , we introduce the following fully discrete method
based on the semi-discrete problem (8), given by the following equation:

$\displaystyle\big{(}\delta^{2}v^{\eta},w\big{)}+a^{\eta+\theta}(v,w)+\nu\leavevmode\nobreak\
\\!\big{(}\delta
v^{\eta},w\big{)}+\big{(}\big{[}a^{(2)}_{ij}\nabla_{y_{i},y_{j}}\leavevmode\nobreak\
\\!v\big{]}^{\eta+\theta},w\big{)}$ (49)
$\displaystyle\quad+\big{(}a^{(4,\eta)}_{i}\nabla_{y_{i}}\delta
v^{\eta},w\big{)}+\big{(}\big{[}a^{(5)}_{i}\nabla_{y_{i}}v\big{]}^{\eta+\theta},w\big{)}=\big{(}\psi^{\eta},w\big{)},\leavevmode\nobreak\
\forall\leavevmode\nobreak\ w\in V,$

where $a_{i}^{(4,\eta)}=a_{i}^{(4)}(y,t_{\eta})\leavevmode\nobreak\ $ and

$\displaystyle a^{(\vartheta,\leavevmode\nobreak\
\\!\eta+\theta)}(v,w)=\sum_{k=1}^{3}\big{(}\big{[}f_{k}D_{k}v\big{]}^{\eta+\theta},D_{k}w\big{)},$
(50)
$\displaystyle\psi^{\eta}=\delta^{2}v^{\eta}-v_{tt}^{\eta+\theta}+\nu\leavevmode\nobreak\
\\!\big{(}\delta
v^{\eta}-v_{t}^{\eta+\theta}\big{)}+\big{(}a^{(4,\eta)}_{i}\nabla_{y_{i}}\delta
v^{\eta}-\big{[}a^{(4)}_{i}\nabla_{y_{i}}v_{t}\big{]}^{\eta+\theta}\big{)},$

were $f_{1}:\Omega\times\\{1,\cdots,N\\}\longrightarrow\mathbb{R}$,
$f_{2}:\\{v,\leavevmode\nobreak\
\\!v_{h}\\}\times\\{1,\cdots,N\\}\longrightarrow\mathbb{R}$ and
$f_{3}:\\{1,\cdots,N\\}\longrightarrow\mathbb{R}$, such that
$f_{1}^{\eta}(y)=a_{i}^{(1,\eta)}(y),\leavevmode\nobreak\
f_{2}^{\eta}(\vartheta)=b_{1}^{\eta}|\vartheta^{\eta}|_{0}^{2}$ and
$f_{3}^{\eta}=b_{2}^{\eta}$, and $D_{k}$ representing the operators
$\nabla_{y_{i}},\leavevmode\nobreak\ \nabla$ and $\Delta$, for the $k=1,2,3$,
respectively.

In the following theorem, we will present an error estimate for the discrete
time associated with the problem (49):

######  Theorem 5.

(Error estimates for the fully-discrete problem)

Let $v$ be the solution to the problem (49). Assuming $v_{0},v_{1}\in
H_{0}^{2}(\Omega)\cap H^{4}(\Omega)$, and considering the hypotheses of the
theorem 4, we have the following estimate for the approximation error
$\forall\leavevmode\nobreak\ \\!\theta\in\leavevmode\nobreak\ \\!]1/4,1]$:

$\displaystyle\|\delta
e\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\|e\|_{L^{\infty}(0,T;H_{0}^{2}(\Omega))}\leq
C\leavevmode\nobreak\ \\!\big{[}h^{2}+(\Delta t)^{2}\big{]},$

where $C>0$ is constant and independent of $h$ and $\Delta t$.

###### Proof.

The decomposition of the error (25), for the discrete case, is given by
$e^{\eta}=\rho^{\eta}+\xi^{\eta}=(v^{\eta}-R_{h}v^{\eta})+(R_{h}v^{\eta}-v_{h}^{\eta})$.
Making the difference between the problem (49) restricted to $V_{h}$ and its
approximate problem, taking in particular $w_{h}=\delta\xi^{\eta}$ and using
the decomposition of the error and the Ritz projection, we obtain

$\Gamma_{1}^{\eta}(\xi,v_{h})+\Gamma_{2}^{\eta}(\xi)=(\psi^{\eta},\delta\xi^{\eta})-\big{[}\Gamma_{1}^{\eta}(\rho,v)+\Gamma_{2}^{\eta}(\rho)+\mu^{\eta}\big{]},$
(51)

where
$\Gamma_{1}:V\times\\{v,v_{h}\\}\times\\{0,\cdots,N\\}\longrightarrow\mathbb{R}$,
$\leavevmode\nobreak\
\Gamma_{2}:V\times\\{0,\cdots,N\\}\longrightarrow\mathbb{R}\leavevmode\nobreak\
$ and $\leavevmode\nobreak\
\mu:\\{0,\cdots,N\\}\longrightarrow\mathbb{R}\leavevmode\nobreak\ $ such that

$\displaystyle\Gamma_{1}^{\eta}(\chi,\vartheta)$
$\displaystyle=\big{(}\delta^{2}\chi^{\eta},\delta\xi^{\eta}\big{)}+a^{(\vartheta,\leavevmode\nobreak\
\\!\eta+\theta)}(\chi,\delta\xi^{\eta})+\nu\leavevmode\nobreak\
\\!\big{(}\delta\chi^{\eta},\delta\xi^{\eta}\big{)}+\big{(}a^{(4,\eta)}_{i}(t)\nabla_{y_{i}}\delta\chi^{\eta},\delta\xi^{\eta}\big{)},$
$\displaystyle\Gamma_{2}^{\eta}(\chi)$
$\displaystyle=\big{(}\big{[}a^{(2)}_{ij}\nabla_{y_{i},y_{j}}\chi\big{]}^{\eta+\theta},\delta\xi^{\eta}\big{)}+\big{(}\big{[}a^{(5)}_{i}\nabla_{y_{i}}\chi\big{]}^{\eta+\theta},\delta\xi^{\eta}\big{)},$
$\displaystyle\mu^{\eta}$ $\displaystyle=\big{(}\big{[}b_{1}\big{(}|\nabla
v_{h}|_{0}^{2}-|\nabla v|_{0}^{2}\big{)}R_{h}\nabla
v\big{]}^{\eta+\theta},\nabla\delta\xi^{\eta}\big{)}.$

Analysis: Terms dependent on $\xi$

By definition, we have

$\displaystyle\Gamma_{1}^{\eta}(\xi,v_{h})$
$\displaystyle=\big{(}\delta^{2}\xi^{\eta},\delta\xi^{\eta}\big{)}+a^{(v_{h},\leavevmode\nobreak\
\\!\eta+\theta)}(\xi,\delta\xi^{\eta})+\nu\leavevmode\nobreak\
\\!\big{(}\delta\xi^{\eta},\delta\xi^{\eta}\big{)}+\big{(}a^{(4,\eta)}_{i}(t)\nabla_{y_{i}}\delta\xi^{\eta},\delta\xi^{\eta}\big{)}.$

Analysing each term, we get

$\displaystyle(\delta^{2}\xi^{\eta},\delta\xi^{\eta})=\frac{1}{2\Delta
t}\big{(}|\delta\xi^{\eta+1/2}|^{2}_{0}-|\delta\xi^{\eta-1/2}|^{2}_{0}\big{)}.$

Summing of $\eta=1$ at $N_{0}$, for some $1<N_{0}\leq N$, we have

$\sum_{\eta=1}^{N_{0}-1}(\delta^{2}\xi^{\eta},\delta\xi^{\eta})=\frac{1}{2\Delta
t}\big{(}|\delta\xi^{N_{0}-1/2}|^{2}_{0}-|\delta\xi^{1/2}|^{2}_{0}\big{)}.$
(52)

We also have to

$\nu\leavevmode\nobreak\
\\!\big{(}\delta\xi^{\eta},\delta\xi^{\eta}\big{)}=\nu\leavevmode\nobreak\
\\!|\delta\xi^{\eta}|_{0}^{2}\geq 0.$ (53)

Integrating by parts and using the definitions (7), we have the following
equality:

$\big{(}a^{(4,\eta)}_{i}\nabla_{y_{i}}\delta\xi^{\eta},\delta\xi^{\eta}\big{)}=-\frac{1}{2}\int_{\Omega}\Big{(}\nabla_{y_{i}}a^{(4,\eta)}_{i}(y)\Big{)}\big{(}\delta\xi^{\eta}\big{)}^{2}dy=\frac{K^{\prime}\leavevmode\nobreak\
\\!{}^{\eta}}{K^{\eta}}|\delta\xi^{\eta}|_{0}^{2}\geq 0,$ (54)

where, we use that
$\displaystyle\nabla_{y_{i}}a^{(4,\eta)}_{i}(y)=-2(K^{\prime}\leavevmode\nobreak\
\\!{}^{\eta}/K^{\eta})$.

For to analyse the second term of $\Gamma_{1}(\xi,v_{h})$, we note that

$\displaystyle\big{(}[f\leavevmode\nobreak\
\\!D\xi]^{\eta+\theta},D\delta\xi^{\eta}\big{)}$
$\displaystyle=\frac{\theta}{2\Delta
t}\int_{\Omega}f^{\eta+1}\big{(}D\xi^{\eta+1}\big{)}^{2}-f^{\eta-1}\big{(}D\xi^{\eta-1}\big{)}^{2}\leavevmode\nobreak\
dy-\theta\big{(}f^{\prime}(\sigma_{k1})D\xi^{\eta+1},D\xi^{\eta-1}\big{)}$
(55) $\displaystyle\qquad+\frac{(1-2\theta)}{2\Delta
t}\big{[}\big{(}f^{\eta}D\xi^{\eta+1},D\xi^{\eta}\big{)}-\big{(}f^{\eta}D\xi^{\eta},D\xi^{\eta-1}\big{)}\big{]},$

for some $\sigma_{k1}\in(t_{\eta-1},t_{\eta+1})$, varying
$k\in\\{1,\cdots,3\\}$. The terms with the derivative in time were obtained
from the application of the Theorem Fundamental of Calculus and the Theorem of
the Mean Value as follow

$f^{\eta-1}-f^{\eta+1}=-\int_{t_{\eta-1}}^{t_{\eta+1}}f^{\prime}(t)\leavevmode\nobreak\
dt=-f^{\prime}(\sigma_{k1})\int_{t_{\eta-1}}^{t_{\eta+1}}\\!dt=-2\leavevmode\nobreak\
\\!\Delta t\leavevmode\nobreak\ \\!f^{\prime}(\sigma_{k1}).$ (56)

Therefore by (55), we get that

$\displaystyle\begin{aligned} &a^{(v_{h},\leavevmode\nobreak\
\\!\eta+\theta)}(\xi,\delta\xi^{\eta})=\sum_{k=1}^{3}\bigg{\\{}\frac{\theta}{2\Delta
t}\int_{\Omega}f_{k}^{\eta+1}\big{(}D_{k}\xi^{\eta+1}\big{)}^{2}-f_{k}^{\eta-1}\big{(}D_{k}\xi^{\eta-1}\big{)}^{2}\leavevmode\nobreak\
dy\\\
&\quad-\theta\big{(}f_{k}^{\prime}(\sigma_{k1})D_{k}\xi^{\eta+1},D_{k}\xi^{\eta-1}\big{)}+\frac{(1-2\theta)}{2\Delta
t}\big{[}\big{(}f_{k}^{\eta}D_{k}\xi^{\eta+1},D_{k}\xi^{\eta}\big{)}-\big{(}f_{k}^{\eta}D_{k}\xi^{\eta},D_{k}\xi^{\eta-1}\big{)}\big{]}\bigg{\\}}.\end{aligned}$

Summing of $\eta=1$ at $N_{0}$, we have

$\displaystyle\sum_{\eta=1}^{N_{0}-1}a^{(v_{h},\leavevmode\nobreak\
\\!\eta+\theta)}(\xi,\delta\xi^{\eta})=\frac{1}{2\Delta
t}\sum_{k=1}^{3}\Big{\\{}\theta\int_{\Omega}f_{k}^{N_{0}}\big{(}D_{k}\xi^{N_{0}}\big{)}^{2}+f_{k}^{N_{0}-1}\big{(}D_{k}\xi^{N_{0}-1}\big{)}^{2}\leavevmode\nobreak\
dy$ (57)
$\displaystyle\quad+(1-2\theta)\big{(}f_{k}^{N_{0}-1}D_{k}\xi^{N_{0}},D_{k}\xi^{N_{0}-1}\big{)}-\theta\int_{\Omega}f_{k}^{1}\big{(}D_{k}\xi^{1}\big{)}^{2}+f_{k}^{0}\big{(}D_{k}\xi^{0}\big{)}^{2}\leavevmode\nobreak\
dy$
$\displaystyle\quad-(1-2\theta)\big{(}f_{k}^{1}D_{k}\xi^{1},D_{k}\xi^{0}\big{)}-\Delta
t\sum_{\eta=1}^{N_{0}-1}\big{[}2\leavevmode\nobreak\
\\!\theta\big{(}f_{k}^{\prime}(\sigma_{k1})D_{k}\xi^{\eta+1},D_{k}\xi^{\eta-1}\big{)}$
$\displaystyle\quad-\Delta t\leavevmode\nobreak\ \\!(1-2\leavevmode\nobreak\
\\!\theta)\big{(}f_{k}^{\prime}(\sigma_{k2})D_{k}\xi^{\eta},D_{k}\xi^{\eta-1}\big{)}\big{]}\Big{\\}},$

for some $\sigma_{k2}\in(t_{\eta-1},t_{\eta})$, varying
$k\in\\{1,\cdots,3\\}$.

So, by analysis (52), (54) and (57) we obtain

$\displaystyle\sum_{\eta=1}^{N_{0}-1}\Gamma_{1}^{\eta}(\xi,v_{h})\geq\frac{1}{2\Delta
t}\Big{[}\Theta_{1}^{N_{0},N_{0}-1}-\leavevmode\nobreak\
\\!\Theta_{1}^{1,1}-\Delta t\sum_{\eta=1}^{N_{0}-1}\Theta_{2}^{\eta}\Big{]},$
(58)

where

$\displaystyle\Theta_{1}^{\eta,\lambda}$
$\displaystyle=|\delta\xi^{\eta-1/2}|^{2}_{0}+\sum_{k=1}^{3}\Big{[}\theta\\!\int_{\Omega}f_{k}^{\eta}\big{(}D_{k}\xi^{\eta}\big{)}^{2}+f_{k}^{\eta-1}\big{(}D_{k}\xi^{\eta-1}\big{)}^{2}\leavevmode\nobreak\
dy+(1-2\theta)\big{(}f_{k}^{\lambda}D_{k}\xi^{\eta},D_{k}\xi^{\eta-1}\big{)}\Big{]},$
$\displaystyle\Theta_{2}^{\eta}$
$\displaystyle=\sum_{k=1}^{3}\big{[}2\leavevmode\nobreak\
\\!\theta\big{(}f_{k}^{\prime}(\sigma_{k1})D_{k}\xi^{\eta+1},D_{k}\xi^{\eta-1}\big{)}+(1-2\leavevmode\nobreak\
\\!\theta)\big{(}f_{k}^{\prime}(\sigma_{k2})D_{k}\xi^{\eta},D_{k}\xi^{\eta-1}\big{)}\big{]}.$

Now, analysing each term of $\Gamma_{2}^{\eta}(\xi)$, give by

$\displaystyle\Gamma_{2}^{\eta}(\xi)=\big{(}\big{[}a^{(2)}_{ij}\nabla_{y_{i},y_{j}}\xi\big{]}^{\eta+\theta},\delta\xi^{\eta}\big{)}+\big{(}\big{[}a^{(5)}_{i}\nabla_{y_{i}}\xi\big{]}^{\eta+\theta},\delta\xi^{\eta}\big{)},$

we have, using the inequality of Cauchy-Schwarz and Young, the definition of
$|\nabla\xi^{\eta}|_{0}$, increasing the norm
$\big{[}|\nabla\xi|_{0}^{2}\big{]}^{\eta+\theta}$ and consider $\theta\leq 1$,
follows that

$\displaystyle\Gamma_{2}^{\eta}(\xi)\leq\frac{\kappa_{3}}{2}\big{[}\varepsilon\big{(}\bar{E}^{\eta+1}+\bar{E}^{\eta}+\bar{E}^{\eta-1}\big{)}+\frac{1}{\varepsilon}|\delta\xi^{\eta}|_{0}^{2}\big{]},$

for some $\varepsilon>0$ and for

$\displaystyle\bar{E}\leavevmode\nobreak\
\\!^{\eta}=|\nabla_{y_{i}}\xi^{\eta}|_{0}^{2}+|\nabla\xi^{\eta}|_{0}^{2}+|\Delta\xi^{\eta}|_{0}^{2}\leavevmode\nobreak\
.$

Summing of $\eta=1$ at $N_{0}-1$, we get

$\sum_{\eta=1}^{N_{0}-1}\Gamma_{2}^{\eta}(\xi)\leq\frac{\kappa_{3}}{2}\leavevmode\nobreak\
\\!\varepsilon\leavevmode\nobreak\
\\!\bar{E}^{N_{0}}+\frac{\kappa_{3}}{2}\sum_{\eta=1}^{N_{0}-1}\big{[}\varepsilon\leavevmode\nobreak\
\\!\bar{E}^{\eta}+\frac{1}{\varepsilon}|\delta\xi^{\eta}|_{0}^{2}\big{]}$ (59)

Analisys: Terms dependent on $\rho$

By definition, we have

$\displaystyle\Gamma_{1}^{\eta}(\rho,v)+\Gamma_{2}^{\eta}(\rho)=\big{(}\delta^{2}\rho^{\eta},\delta\xi^{\eta}\big{)}+a^{(v,\leavevmode\nobreak\
\\!\eta+\theta)}(\rho,\delta\xi^{\eta})+\nu\leavevmode\nobreak\
\\!\big{(}\delta\rho^{\eta},\delta\xi^{\eta}\big{)}+\big{(}a^{(4,\eta)}_{i}(t)\nabla_{y_{i}}\delta\rho^{\eta},\delta\xi^{\eta}\big{)}$
$\displaystyle\qquad+\big{(}\big{[}a^{(2)}_{ij}\nabla_{y_{i},y_{j}}\rho\big{]}^{\eta+\theta},\delta\xi^{\eta}\big{)}+\big{(}\big{[}a^{(5)}_{i}\nabla_{y_{i}}\rho\big{]}^{\eta+\theta},\delta\xi^{\eta}\big{)}.$

Analysing the terms of $\Gamma_{1}^{\eta}(\rho,v)$ and
$\Gamma_{2}^{\eta}(\rho)$ applying the inequality of Cauchy-Schwarz and Young,
and considering $\theta\leq 1$, follows that

$\displaystyle\Gamma_{1}^{\eta}(\rho,v)+\Gamma_{2}^{\eta}(\rho)$
$\displaystyle\leq\frac{\kappa_{4}}{2}\Big{\\{}\frac{1}{\varepsilon}\big{[}|\delta^{2}\rho^{\eta}|_{0}^{2}+|\delta\rho^{\eta}|_{0}^{2}+|\nabla\delta\rho^{\eta}|_{0}^{2}+\big{(}|\nabla\rho|_{0}^{2}+|\Delta\rho|_{0}^{2}\big{)}^{\eta+1}$
$\displaystyle\qquad+\big{(}|\nabla\rho|_{0}^{2}+|\Delta\rho|_{0}^{2}\big{)}^{\eta}+\big{(}|\nabla\rho|_{0}^{2}+|\Delta\rho|_{0}^{2}\big{)}^{\eta-1}\big{]}+\varepsilon\leavevmode\nobreak\
\\!|\delta\xi^{\eta}|_{0}^{2}\Big{\\}}$
$\displaystyle\leq\frac{\kappa_{4}}{2}\Big{[}\frac{1}{\varepsilon}\big{(}|\delta^{2}\rho^{\eta}|_{0}^{2}+\|\delta\rho^{\eta}\|_{1}^{2}+\|\rho^{\eta+1}\|_{2}^{2}+\|\rho^{\eta}\|_{2}^{2}+\|\rho^{\eta-1}\|_{2}^{2}\big{)}+\varepsilon\leavevmode\nobreak\
\\!|\delta\xi^{\eta}|_{0}^{2}\Big{]}.$

This is similar to (39). Summing of $\eta=1$ at $N_{0}$ we get

$\sum_{\eta=1}^{N_{0}-1}\big{[}\Gamma_{1}^{\eta}(\rho,v)+\Gamma_{2}^{\eta}(\rho)\big{]}\leq\frac{\kappa_{4}}{2}\Big{\\{}\frac{1}{\varepsilon}\sum_{\eta=0}^{N_{0}}\big{[}|\delta^{2}\rho^{\eta}|_{0}^{2}+\|\delta\rho^{\eta}\|_{1}^{2}+\|\rho^{\eta}\|_{2}^{2}\big{]}+\varepsilon\sum_{\eta=1}^{N_{0}-1}|\delta\xi^{\eta}|_{0}^{2}\Big{\\}}.$
(60)

We note that by definition of Ritz projection in (27), for each discrete time
$t=t_{\eta-1},\leavevmode\nobreak\ t_{\eta},\leavevmode\nobreak\ t_{\eta+1}$,
we get
$a^{(v,\eta)}(\rho,\delta\xi^{\eta})=a^{(v,\eta)}(v-R_{h}v,\delta\xi^{\eta})=0$
.

Analisys: Nonlinear term

The calculation procedure for the nonlinear term is analogous to the estimates
(40), that is, for each of the discrete times
$t=t_{\eta-1},\leavevmode\nobreak\ t_{\eta},\leavevmode\nobreak\ t_{\eta+1}$,
we have

$\displaystyle\mu^{\eta}\leq\frac{\kappa_{5}}{2}\Big{[}\frac{1}{\varepsilon}\big{(}|\nabla\rho^{\eta+1}|_{0}^{2}+|\nabla\rho^{\eta}|_{0}^{2}+|\nabla\rho^{\eta-1}|_{0}^{2}\big{)}+\varepsilon\big{(}|\nabla\xi^{\eta+1}|_{0}^{2}+|\nabla\xi^{\eta}|_{0}^{2}+|\nabla\xi^{\eta-1}|_{0}^{2}+|\delta\xi^{\eta}|_{0}^{2}\big{)}\Big{]}.$

Summing of $\eta=1$ at $N_{0}-1$, we get

$\sum_{\eta=1}^{N_{0}-1}\mu^{\eta}\leq\frac{\kappa_{5}}{2}\Big{\\{}\varepsilon\leavevmode\nobreak\
\\!|\nabla\xi^{N_{0}}|_{0}^{2}+\frac{1}{\varepsilon}\sum_{\eta=0}^{N_{0}}\big{[}\|\rho^{\eta}\|_{1}^{2}\big{]}+\varepsilon\sum_{\eta=1}^{N_{0}-1}\big{[}|\nabla\xi^{\eta}|_{0}^{2}+|\delta\xi^{\eta}|_{0}^{2}\big{]}\Big{\\}}.$
(61)

Analysis: Term of finite differences

By definition, from (50)2, we have

$\psi^{\eta}=\delta^{2}v^{\eta}-v_{tt}^{\eta+\theta}+\nu\leavevmode\nobreak\
\\!\big{(}\delta
v^{\eta}-v_{t}^{\eta+\theta}\big{)}+\Big{(}a^{(4,\eta)}_{i}\nabla_{y_{i}}\delta
v^{\eta}-\big{[}a^{(4)}_{i}\nabla_{y_{i}}v_{t}\big{]}^{\eta+\theta}\Big{)}.$

Using integration by parts we can obtain

$\displaystyle\psi^{\eta}$ $\displaystyle=\frac{1}{(\Delta
t)^{2}}\bigg{[}\int_{t_{\eta-1}}^{t_{\eta}}v^{(iv)}(\cdot,s)\frac{(s-t_{\eta-1})^{3}}{3!}\leavevmode\nobreak\
ds+\int_{t_{\eta}}^{t_{\eta+1}}v^{(iv)}(\cdot,s)\frac{(t_{\eta+1}-s)^{3}}{3!}\leavevmode\nobreak\
ds\bigg{]}$ $\displaystyle\quad+\frac{1}{2\Delta
t}\bigg{[}\int_{t_{\eta-1}}^{t_{\eta}}v^{\prime\prime\prime}(\cdot,s)\frac{(s-t_{\eta-1})^{2}}{2}\leavevmode\nobreak\
ds+\int_{t_{\eta}}^{t_{\eta+1}}v^{\prime\prime\prime}(\cdot,s)\frac{(t_{\eta+1}-s)^{2}}{2}\leavevmode\nobreak\
ds\bigg{]}$ $\displaystyle\quad+\frac{a_{i}^{(4,\eta)}(y)}{2\Delta
t}\bigg{[}\int_{t_{\eta-1}}^{t_{\eta}}\nabla_{y_{i}}v^{\prime\prime\prime}(\cdot,s)\frac{(s-t_{\eta-1})^{2}}{2}\leavevmode\nobreak\
ds+\int_{t_{\eta}}^{t_{\eta+1}}\nabla_{y_{i}}v^{\prime\prime\prime}(\cdot,s)\frac{(t_{\eta+1}-s)^{2}}{2}\leavevmode\nobreak\
ds\bigg{]}$
$\displaystyle\quad+\theta\bigg{[}\int_{t_{\eta-1}}^{t_{\eta}}\big{(}v^{(iv)}-v^{\prime\prime\prime}+(a_{i}^{(4)}\nabla_{y_{i}}v^{\prime})^{\prime\prime}\big{)}(\cdot,s)(s-t_{\eta-1})\leavevmode\nobreak\
ds$
$\displaystyle\quad\qquad+\int_{t_{\eta}}^{t_{\eta+1}}\big{(}v^{(iv)}-v^{\prime\prime\prime}+(a_{i}^{(4)}\nabla_{y_{i}}v^{\prime})^{\prime\prime}\big{)}(\cdot,s)(t_{\eta+1}-s)\leavevmode\nobreak\
ds\bigg{]}.$

Considering that $0\leq(s-t_{\eta-1})\leq\Delta t$, for all
$s\in[t_{\eta-1},t_{\eta}]$ and that $0\leq(t_{\eta+1}-s)\leq\Delta t$, for
all $s\in[t_{\eta},t_{\eta+1}]$, applying the supreme essential in the
derivatives of $v$ and then resolving the integrals, considering $\theta\leq
1$, we get

$\psi^{\eta}\leq\frac{\kappa_{6}}{2}\leavevmode\nobreak\ \\!(\Delta
t)^{2}\Big{\\{}\sum_{k=1}^{3}\Big{\|}\frac{d^{k}\nabla
v}{dt^{k}}\Big{\|}_{L^{\infty}(t_{\eta-1},t_{\eta+1})}+\sum_{k=3}^{4}\Big{\|}\frac{d^{k}v}{dt^{k}}\Big{\|}_{L^{\infty}(t_{\eta-1},t_{\eta+1})}\Big{\\}},$
(62)

where $\kappa_{6}=2\leavevmode\nobreak\
\\!\kappa_{5}+\displaystyle\max_{i\in\\{1,\cdots,n\\}}\Big{\\{}\max_{t\in[0,T]}\Big{[}2\leavevmode\nobreak\
\\!\sup_{y\in\Omega}\big{(}|a_{i}^{\prime\prime}\leavevmode\nobreak\
\\!{}^{(1)}(y,t)|+|a_{i}^{\prime}\leavevmode\nobreak\
\\!{}^{(1)}(y,t)|\big{)}\Big{]}\Big{\\}}$. Once that

$\displaystyle(a_{i}^{(4)}\nabla_{y_{i}}v^{\prime})^{\prime\prime}$
$\displaystyle=a_{i}^{\prime\prime}\leavevmode\nobreak\
\\!{}^{(4)}\nabla_{y_{i}}v^{\prime}+2\leavevmode\nobreak\
\\!a_{i}^{\prime}\leavevmode\nobreak\
\\!{}^{(4)}\nabla_{y_{i}}v^{\prime\prime}+a_{i}^{(4)}\nabla_{y_{i}}v^{\prime\prime\prime}\leq\frac{\kappa_{6}}{2}\sum_{k=1}^{3}\Big{\|}\frac{d^{k}\nabla
v}{dt^{k}}\Big{\|}_{L^{\infty}(t_{\eta-1},t_{\eta+1})}.$

Thus, by Cauchy-Schwarz, Young and (62), we have that

$\displaystyle\big{(}\psi^{\eta},\delta\xi^{\eta}\big{)}$
$\displaystyle\leq\frac{1}{2}\Big{[}\frac{1}{\varepsilon}|\psi^{\eta}|_{0}^{2}+\varepsilon\leavevmode\nobreak\
\\!|\delta\xi^{\eta}|_{0}^{2}\Big{]}$
$\displaystyle\leq\frac{(\kappa_{6})^{2}\\!\\!}{4\leavevmode\nobreak\
\\!\varepsilon}\leavevmode\nobreak\ (\Delta
t)^{4}\Big{[}\sum_{k=1}^{3}\Big{\|}\frac{d^{k}v}{dt^{k}}\Big{\|}^{2}_{L^{\infty}(t_{\eta-1},t_{\eta+1};H^{1}(\Omega))}+\|v^{(iv)}\|^{2}_{L^{\infty}(t_{\eta-1},t_{\eta+1};L^{2}(\Omega))}\Big{]}+\frac{\varepsilon}{2}|\delta\xi^{\eta}|_{0}^{2}.$

Summing of $\eta=1$ at $N_{0}-1$, we get

$\sum_{\eta=1}^{N_{0}-1}\big{(}\psi^{\eta},\delta\xi^{\eta}\big{)}\leq\frac{(\kappa_{6})^{2}\\!\\!}{4\leavevmode\nobreak\
\\!\varepsilon}\leavevmode\nobreak\ \\!(\Delta
t)^{4}\Big{[}\sum_{k=1}^{3}\Big{\|}\frac{d^{k}v}{dt^{k}}\Big{\|}^{2}_{L^{\infty}(0,T;H^{1}(\Omega))}+\|v^{(iv)}\|^{2}_{L^{\infty}(0,T;L^{2}(\Omega))}\Big{]}+\frac{\varepsilon}{2}\sum_{\eta=1}^{N_{0}-1}|\delta\xi^{\eta}|_{0}^{2}$
(63)

Final estimate

Now, passing the term $\Gamma_{2}^{\eta}(\xi)$ to the right side of (51),
summing of $\eta=1$ at $N_{0}-1$ and using the estimates (58) – (63), we
obtain

$\frac{1}{2\Delta t}\Big{[}\Theta_{1}^{N_{0},N_{0}-1}-\Theta_{1}^{1,1}-\Delta
t\sum_{\eta=1}^{N_{0}-1}\Theta_{2}^{\eta}\Big{]}\leq\frac{\kappa_{7}}{2}\leavevmode\nobreak\
\\!\Theta_{3}^{\eta}$ (64)

where $\kappa_{7}=\kappa_{6}(1+\kappa_{6}/2)\leavevmode\nobreak\ $ and

$\displaystyle\Theta_{3}^{\eta}=\varepsilon\leavevmode\nobreak\
\\!\bar{E}^{N_{0}}+\frac{1}{\varepsilon}\sum_{\eta=0}^{N_{0}}\big{[}|\delta^{2}\rho^{\eta}|_{0}^{2}+\|\delta\rho^{\eta}\|_{1}^{2}+\|\rho^{\eta}\|_{2}^{2}\big{]}+\Big{(}\varepsilon+\frac{1}{\varepsilon}\Big{)}\sum_{\eta=1}^{N_{0}-1}\big{[}\bar{E}^{\eta}+|\delta\xi^{\eta}|_{0}^{2}\big{]}$
(65) $\displaystyle\qquad+\frac{(\Delta
t)^{4}}{\varepsilon}\Big{[}\sum_{k=1}^{3}\Big{\|}\frac{d^{k}v}{dt^{k}}\Big{\|}^{2}_{L^{\infty}(0,T;H^{1}(\Omega))}+\|v^{(iv)}\|^{2}_{L^{\infty}(0,T;L^{2}(\Omega))}\Big{]}.$

We need to analyze the term $\Theta_{2}^{\eta}$ before proceeding. We note
that, by Cauchy-Schwarz and Young, considering $\theta\leq 1$:

$\displaystyle\Theta_{2}^{\eta}$
$\displaystyle=\sum_{k=1}^{3}\big{[}2\leavevmode\nobreak\
\\!\theta\big{(}f_{k}^{\prime}(\sigma_{k1})D_{k}\xi^{\eta+1},D_{k}\xi^{\eta-1}\big{)}+(1-2\leavevmode\nobreak\
\\!\theta)\big{(}f_{k}^{\prime}(\sigma_{k2})D_{k}\xi^{\eta},D_{k}\xi^{\eta-1}\big{)}\big{]}$
$\displaystyle\leq 2\leavevmode\nobreak\
\\!\kappa_{2}\sum_{k=1}^{3}\big{[}\varepsilon\big{(}|D_{k}\xi^{\eta+1}|_{0}^{2}+|D_{k}\xi^{\eta}|_{0}^{2}\big{)}+\frac{1}{\varepsilon}|D_{k}\xi^{\eta-1}|_{0}^{2}\big{]}=2\leavevmode\nobreak\
\\!\kappa_{2}\big{[}\varepsilon\big{(}\bar{E}^{\eta+1}+\bar{E}^{\eta}\big{)}+\frac{1}{\varepsilon}\bar{E}^{\eta-1}\big{]}.$

Therefore,

$\Delta
t\sum_{\eta=1}^{N_{0}-1}\Theta_{2}^{\eta}\leq\kappa_{7}\leavevmode\nobreak\
\\!\varepsilon\leavevmode\nobreak\ \\!\Delta t\leavevmode\nobreak\
\\!\bar{E}^{N_{0}}+\kappa_{8}(\varepsilon)\leavevmode\nobreak\ \\!\Delta
t\sum_{\eta=0}^{N_{0}-1}\bar{E}^{\eta}\leq\kappa_{7}\leavevmode\nobreak\
\\!\varepsilon\leavevmode\nobreak\
\\!\bar{E}^{N_{0}}+\kappa_{8}(\varepsilon)\leavevmode\nobreak\ \\!\Delta
t\sum_{\eta=0}^{N_{0}-1}\bar{E}^{\eta},$ (66)

where $\kappa_{8}(\varepsilon)=\kappa_{7}(2+\varepsilon+1/\varepsilon)$ and in
the last inequality we consider $\Delta t\leq 1$ for the first parcel.

On the other hand, for the term $\Theta_{1}^{1,1}$, considering $\theta\leq 1$
and applying Cauchy-Schwarz and Young, we get

$\displaystyle\Theta_{1}^{1,1}$
$\displaystyle=|\delta\xi^{1/2}|^{2}_{0}+\sum_{k=1}^{3}\Big{[}\theta\\!\int_{\Omega}f_{k}^{1}\big{(}D_{k}\xi^{1}\big{)}^{2}+f_{k}^{0}\big{(}D_{k}\xi^{0}\big{)}^{2}\leavevmode\nobreak\
dy+(1-2\theta)\big{(}f_{k}^{1}D_{k}\xi^{1},D_{k}\xi^{0}\big{)}\Big{]}$ (67)
$\displaystyle\leq\kappa_{1}\leavevmode\nobreak\
\\!\Big{[}|\delta\xi^{1/2}|_{0}^{2}+\sum_{k=1}^{3}\big{[}|D_{k}\xi^{1}|_{0}^{2}+|D_{k}\xi^{0}|_{0}^{2}\big{]}\Big{]}\leq
2\leavevmode\nobreak\ \\!\kappa_{1}\leavevmode\nobreak\
\\!\big{[}|\delta\xi^{1/2}|_{0}^{2}+\|\xi^{1}\|_{2}^{2}+\|\xi^{0}\|_{0}^{2}\big{]}.$

In the last inequality we use that

$\displaystyle\sum_{k=1}^{3}\big{[}|D_{k}\xi^{\eta}|_{0}^{2}\big{]}=|\nabla_{y_{i}}\xi^{\eta}|_{0}^{2}+|\nabla\xi^{\eta}|_{0}^{2}+|\Delta\xi^{\eta}|_{0}^{2}\leq
2\big{[}|\xi^{\eta}|_{0}^{2}+|\nabla\xi^{\eta}|_{0}^{2}+|\Delta\xi^{\eta}|_{0}^{2}\big{]}=2\leavevmode\nobreak\
\\!\|\xi^{\eta}\|_{2}^{2},\leavevmode\nobreak\ \forall\leavevmode\nobreak\
\\!\eta.$

Following the procedures of article [6] and considering (67) and the
hypothesis (H3) we get

$\Theta_{1}^{1,1}\leq 2\leavevmode\nobreak\
\\!\kappa_{1}\big{[}|\delta\xi^{1/2}|_{0}^{2}+\|\xi^{1}\|_{2}^{2}+\|\xi^{0}\|_{2}^{2}\big{]}\leq
2\leavevmode\nobreak\ \\!\kappa_{1}\big{[}(\Delta t)^{4}+h^{4}\big{]}.$ (68)

Now, analysing the $\rho$ norms of term $\Theta_{3}^{\eta}$, we note that, by
the Lemma of Céa and Douglas-Dupont,

$\displaystyle\|\rho^{\eta}\|_{m}\leq\|v^{\eta}-\hat{v}^{\eta}\|_{m}\leq
C\leavevmode\nobreak\ \\!h^{4-m}\|v\|_{4},$

were $\hat{v}$ is the interpolator of $v$ in $V_{h}$. So we have, considering
$h^{4}\leq h^{3}\leq h^{2}$ and multiply by $\Delta t$:

$\displaystyle\Delta
t\sum_{\eta=0}^{N^{*}}\Big{[}|\delta^{2}\rho^{\eta}|_{0}^{2}+\|\delta\rho^{\eta}\|_{1}^{2}+\|\rho^{\eta}\|_{2}^{2}\Big{]}$
$\displaystyle\leq\int_{0}^{T}|\rho^{\prime\prime}(s)|_{0}^{2}+\|\rho^{\prime}(s)\|_{1}^{2}+\|\rho(s)\|_{2}^{2}\leavevmode\nobreak\
ds$ (69) $\displaystyle\leq
h^{4}\sum_{k=0}^{2}\Big{\|}\frac{d^{k}\leavevmode\nobreak\
\\!v}{dt^{k}}\Big{\|}_{L^{2}(0,T;H^{4}(\Omega))}^{2}.$

The hypothesis (H2) ensures that the norms of $v$ and its derivatives in the
time of (69), and the norms of last parcels of $\Theta_{3}^{\eta}$ in (65) are
constants.

Multiply both side of (64) by $2\leavevmode\nobreak\ \\!\Delta t$ and passing
the terms $\Theta_{1}^{1,1}$ and the sum of $\Theta_{2}^{\eta}$ to the right
side, and applying the estimates (66) – (69) we have

$\Theta_{1}^{N_{0},N_{0}-1}\leq\kappa_{7}\leavevmode\nobreak\
\\!\varepsilon\leavevmode\nobreak\
\\!\bar{E}^{N_{0}}+\kappa_{9}(\varepsilon)\Big{\\{}\big{[}h^{4}+(\Delta
t)^{4}\big{]}+\sum_{\eta=1}^{N_{0}-1}\bar{E}^{\eta}+|\delta\xi^{\eta}|_{0}^{2}\Big{\\}}.$
(70)

were $\kappa_{9}(\varepsilon)$ is the product between
$\kappa_{8}(\varepsilon)$ and maximum of norms of $v$ and its derivatives in
time.

By definition, we have

$\displaystyle|\delta\xi^{\eta}|_{0}=\frac{1}{2}|\delta\xi^{\eta+1/2}+\delta\xi^{\eta-1/2}|_{0}\leq\frac{1}{2}\big{[}|\delta\xi^{\eta+1/2}|_{0}+|\delta\xi^{\eta-1/2}|_{0}\big{]}.$

So, summing of $\eta=1$ at $N_{0}-1$, we get

$\sum_{\eta=1}^{N_{0}-1}|\delta\xi^{\eta}|_{0}\leq\frac{1}{2}|\delta\xi^{N_{0}-1/2}|_{0}^{2}+\sum_{\eta=1}^{N_{0}-1}|\delta\xi^{\eta-1/2}|_{0}.$
(71)

Therefore substituing (71) in (70), we have

$\Theta_{1}^{N_{0},N_{0}-1}\leq\kappa_{7}\leavevmode\nobreak\
\\!\varepsilon\leavevmode\nobreak\
\\!E^{N_{0}}+\kappa_{9}(\varepsilon)\Big{\\{}\big{[}h^{4}+(\Delta
t)^{4}\big{]}+\sum_{\eta=1}^{N_{0}-1}E^{\eta}\Big{\\}}.$ (72)

were

$\displaystyle
E^{\eta}=|\delta\xi^{\eta-1/2}|_{0}+|\nabla_{y_{i}}\xi^{\eta}|_{0}^{2}+|\nabla\xi^{\eta}|_{0}^{2}+|\Delta\xi^{\eta}|_{0}^{2}.$

Finally, we will analyze the term $\Theta_{1}^{N_{0},N_{0}-1}$. By definition
(58) we have

$\displaystyle\Theta_{1}^{N_{0},N_{0}-1}=|\delta\xi^{N_{0}-1/2}|^{2}_{0}+\sum_{k=1}^{3}\Big{[}\theta\int_{\Omega}f_{k}^{N_{0}}\big{(}D_{k}\xi^{N_{0}}\big{)}^{2}+f_{k}^{N_{0}-1}\big{(}D_{k}\xi^{N_{0}-1}\big{)}^{2}\leavevmode\nobreak\
dy$
$\displaystyle\quad+(1-2\theta)\big{(}f_{k}^{N_{0}-1}D_{k}\xi^{N_{0}},D_{k}\xi^{N_{0}-1}\big{)}\Big{]}.$

Adding and subtracting the terms $\theta\leavevmode\nobreak\
\\!f_{k}^{N_{0}}\big{(}D_{k}\xi^{N_{0}-1}\big{)}^{2}$ and
$(1-2\theta)\big{(}f_{k}^{N_{0}}D_{k}\xi^{N_{0}},D_{k}\xi^{N_{0}-1}\big{)}$,
we get

$\displaystyle\Theta_{1}^{N_{0},N_{0}-1}$
$\displaystyle=|\delta\xi^{N_{0}-1/2}|^{2}_{0}+\sum_{k=1}^{3}\Big{[}\theta\int_{\Omega}f_{k}^{N_{0}}\big{[}\big{(}D_{k}\xi^{N_{0}}\big{)}^{2}+\big{(}D_{k}\xi^{N_{0}-1}\big{)}^{2}\big{]}\leavevmode\nobreak\
dy$ (73)
$\displaystyle\quad+(1-2\theta)\big{(}f_{k}^{N_{0}}D_{k}\xi^{N_{0}},D_{k}\xi^{N_{0}-1}\big{)}-\Theta_{4}^{N_{0},k}\Big{]}$
$\displaystyle=|\delta\xi^{N_{0}-1/2}|^{2}_{0}+\sum_{k=1}^{3}\Big{[}\theta\int_{\Omega}f_{k}^{N_{0}}\big{(}D_{k}\xi^{N_{0}}-D_{k}\xi^{N_{0}-1}\big{)}^{2}\leavevmode\nobreak\
dy$
$\displaystyle\quad+\big{(}f_{k}^{N_{0}}D_{k}\xi^{N_{0}},D_{k}\xi^{N_{0}-1}\big{)}-\Theta_{4}^{N_{0},k}\Big{]},$

were, by application of procedures in (56), we obtain to follow the equality
and the estimate, applying Cauchy-Schwarz and Young, considering
$\theta,\Delta t\leq 1$ :

$\displaystyle\Theta_{4}^{N_{0},k}$ $\displaystyle=\Delta
t\Big{[}\theta\int_{\Omega}f^{\prime}_{k}(\sigma_{k3})\big{(}D_{k}\xi^{N_{0}-1}\big{)}^{2}\leavevmode\nobreak\
dy+(1-2\leavevmode\nobreak\
\\!\theta)\big{(}f^{\prime}_{k}(\sigma_{k3})D_{k}\xi^{N_{0}},D_{k}\xi^{N_{0}-1}\big{)}\Big{]}$
$\displaystyle\leq\kappa_{2}\leavevmode\nobreak\
\\!|D_{k}\xi^{N_{0}}|_{0}^{2}+2\leavevmode\nobreak\
\\!\kappa_{2}\leavevmode\nobreak\ \\!\Delta t|D_{k}\xi^{N_{0}-1}|_{0}^{2}.$

for some $\sigma_{k3}\in\leavevmode\nobreak\ \\!(t_{N_{0}-1},t_{N_{0}})$.

Now, by the polar identity, we have to

$\big{(}f_{k}^{N_{0}}D_{k}\xi^{N_{0}},D_{k}\xi^{N_{0}-1}\big{)}=\frac{1}{4}\int_{\Omega}f_{k}^{N_{0}}\big{[}(D_{k}\xi^{N_{0}}+D_{k}\xi^{N_{0}-1})^{2}-(D_{k}\xi^{N_{0}}-D_{k}\xi^{N_{0}-1})^{2}\big{]}dy.$
(74)

So, replacing (74) in (73) we have

$\displaystyle\Theta_{1}^{N_{0},N_{0}-1}$
$\displaystyle\geq|\delta\xi^{N_{0}-1/2}|^{2}_{0}+\sum_{k=1}^{3}\Big{[}\int_{\Omega}f_{k}^{N_{0}}\Big{[}\Big{(}\theta-\frac{1}{4}\Big{)}\big{(}D_{k}\xi^{N_{0}}-D_{k}\xi^{N_{0}-1}\big{)}^{2}$
$\displaystyle\quad+\frac{1}{4}\big{(}D_{k}\xi^{N_{0}}+D_{k}\xi^{N_{0}-1}\big{)}^{2}dy-\Theta_{4}^{N_{0},k}\Big{]}.$

Taking $\theta\in\leavevmode\nobreak\ \\!]1/4,1]$, we get

$\displaystyle\Theta_{1}^{N_{0},N_{0}-1}$
$\displaystyle\geq|\delta\xi^{N_{0}-1/2}|^{2}_{0}+\sum_{k=1}^{3}\big{[}\leavevmode\nobreak\
\\!\kappa_{10}\big{(}\big{|}D_{k}\xi^{N_{0}}-D_{k}\xi^{N_{0}-1}\big{|}_{0}^{2}+\big{|}D_{k}\xi^{N_{0}}+D_{k}\xi^{N_{0}-1}\big{|}_{0}^{2}\big{)}-\Theta_{4}^{N_{0},k}\leavevmode\nobreak\
\\!\big{]}\Big{\\}}$ $\displaystyle>\kappa_{10}\leavevmode\nobreak\
\\!E^{N_{0}}-\sum_{k=1}^{3}\big{[}\Theta_{4}^{N_{0},k}\leavevmode\nobreak\
\\!\big{]},$ (75)

were $\kappa_{10}=\min\\{1,\leavevmode\nobreak\
\kappa_{0}(\theta-1/4),\leavevmode\nobreak\ \kappa_{0}/4\\}>0$ and once that

$\displaystyle|D_{k}\xi^{N_{0}}-D_{k}\xi^{N_{0}-1}\big{|}_{0}^{2}+\big{|}D_{k}\xi^{N_{0}}+D_{k}\xi^{N_{0}-1}\big{|}_{0}^{2}=2\big{(}|D_{k}\xi^{N_{0}}|_{0}^{2}+|D_{k}\xi^{N_{0}-1}|_{0}^{2}\big{)}>|D_{k}\xi^{N_{0}}|_{0}^{2}.$

Substituing (75) in (72) and passing the term $\Theta_{4}^{N_{0},k}$ to the
right side, we get

$(\kappa_{10}-\kappa_{7}\leavevmode\nobreak\
\\!\varepsilon)\leavevmode\nobreak\
\\!E^{N_{0}}\leq\kappa_{9}(\varepsilon)\Big{\\{}\big{[}h^{4}+(\Delta
t)^{4}\big{]}+\sum_{\eta=1}^{N_{0}-1}E^{\eta}\Big{\\}},\leavevmode\nobreak\
\forall\leavevmode\nobreak\ \\!N_{0}\in\\{1,\cdots,N\\}.$ (76)

Taking $\varepsilon=(\kappa_{10}/2\leavevmode\nobreak\ \\!\kappa_{7})>0$,
applying the discrete Gronwall inequality and the equivalence of norms, we
have

$|\delta\xi^{\eta-1/2}|_{0}^{2}+\|\xi^{\eta}\|_{2}^{2}\leq E^{\eta}\leq
C\big{[}(\Delta t)^{4}+h^{4}\big{]},\leavevmode\nobreak\
\forall\leavevmode\nobreak\ \eta=1,\cdots,N.$ (77)

where $\displaystyle
C=({\kappa_{9}}/{\kappa_{10}})\exp\big{(}({\kappa_{9}}/{\kappa_{10}})T\big{)}$.

From (77), using the decomposition of error, the Douglas-Dupont Lemma and
taking the essential supreme in time, we conclude the theorem, i.e.

$\displaystyle\|\delta
e\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\|e\|_{L^{\infty}(0,T;H_{0}^{2}(\Omega))}\leq
C\leavevmode\nobreak\ \\!\big{[}h^{2}+(\Delta t)^{2}\big{]},$

Then, we obtain that the order of convergence of the numerical method is
quadratic in time and space. ∎

##  5 Numerical Simulations

In this section, the accuracy of the discrete methods (18) is tested by
comparing the approximate numerical solutions with the exact solutions and
computing the corresponding approximation error in the norm
$L^{\infty}(0,T;L^{2}(0,L))$, for one-dimensional and two-dimensional cases.
The results of two numerical examples are presented below.

In order to compare the approximate and exact solutions, we consider the
inhomogeneous versions of equations (16) and (17) by introducing the functions
$f(y,t)$, see (12), in their right-hand side, respectively. For these new
problems, we can construct exact solutions by appropriately choosing the
right-hand sides functions.

In our analysis, we use two solutions for one-dimensional and two-dimensional
cases, and two examples of moving boundarys, as shown in Table 1 and Table 2,
respectively.

$v(y,t)$
---
| $\mathbb{R}$ | $\mathbb{R}^{2}$
S1 | $10^{-1}(y^{2}-1)^{2}\cos(2\pi t)$ | $10^{-1}\big{[}(y_{1}^{2}-1)(y_{2}^{2}-1)\big{]}^{2}\cos(2\pi t)$
S2 | $10^{-3}(y^{2}-1)^{2}\sin(2\pi t)$ | $10^{-7}\big{[}(y_{1}^{2}-1)(y_{2}^{2}-1)\big{]}^{2}\sin(2\pi t)$
Table 1: Examples of exact solutions $v(y,t)$ $K(t)$
---
| $\mathbb{R}$ | $\mathbb{R}^{2}$
B${}_{1}\leavevmode\nobreak\ $ | $64+t/2^{7}$ | $64+t/2^{17}$
B${}_{2}\leavevmode\nobreak\ $ | $64+2(1-e^{-t})$ | $64+(1-e^{-t})/2^{17}$
Table 2: Examples of Moving boundaries $K(t)$

We emphasize that the boundaries of the Table 2, satisfy the hypothesis (H1).
The initial conditions immediately follow the exact solutions given.

The Theorem (5), shows that the solution is unconditionally convergent with
order of quadratic convergence in space and time for
$\forall\leavevmode\nobreak\ \\!\theta\in\leavevmode\nobreak\ \\!]1/4,1]$.
However, we know from the literature that for $\forall\leavevmode\nobreak\
\\!\theta\in[0,1/4[$, the system is conditionally convergent and we can see
this fact in numerical simulations in one-dimensional case for the example S1
and moving boundary B1, see Table 3. In addition, we can see that the error is
inversely proportional to the parameter $\forall\leavevmode\nobreak\
\\!\theta\in\leavevmode\nobreak\ \\!]1/4,1]$, that is, the smaller $\theta$
the minor error, as can be observed in the Table 3.

$E_{L^{\infty}(0,T;L^{2}(\Omega))}$
---
$h\diagdown\theta$ | $0$ | $0.25$ | $0.5$ | $0.75$ | $1$
$2^{-1}$ | $5.495e-3$ | $5.463e-3$ | $5.434e-3$ | $5.406e-3$ | $5.379e-3$
$2^{-2}$ | $2.759e-3$ | $2.750e-3$ | $2.745e-3$ | $2.743e-3$ | $2.745e-3$
$2^{-3}$ | $6.404e-4$ | $6.417e-4$ | $6.572e-4$ | $6.859e-4$ | $7.264e-4$
$2^{-4}$ | $2.016e-4$ | $2.076e-4$ | $2.532e-4$ | $3.220e-4$ | $4.022e-4$
$2^{-5}$ | $4.266e-5$ | $5.507e-5$ | $1.510e-4$ | $2.472e-4$ | $3.434e-4$
$2^{-6}$ | diverge | $5.380e-5$ | $1.498e-4$ | $2.460e-4$ | $3.422e-4$
Table 3: Error with $\Delta t=2^{-7}$ fixed and varying $\theta$, for
$\Omega\in\mathbb{R}$, example S1 and boundary B1

Let’s set, from now on, $\theta=1/4$, for generating a minor error. In the
Table 4, the parameter $h=2^{-6}$ is fixed and the $\Delta t$ value are
varying in $\Delta t=2^{-(i+1)},\leavevmode\nobreak\ i=1,\cdots,6$.

$E_{L^{\infty}(0,T;L^{2}(\Omega))}$
---
| $\mathbb{R}$ |  | $\mathbb{R}^{2}$
| B1 | B2 |  | B1 | B2
| $9.651e-2$ | $9.647e-2$ |  | $8.711e-3$ | $8.711e-3$
| $1.750e-2$ | $1.757e-2$ |  | $1.577e-3$ | $1.577e-3$
S1 | $3.173e-3$ | $3.200e-3$ |  | $2.853e-4$ | $2.853e-4$
| $7.713e-4$ | $7.670e-4$ |  | $6.954e-5$ | $6.954e-5$
| $2.063e-4$ | $2.051e-4$ |  | $1.857e-5$ | $1.857e-5$
| $5.380e-5$ | $5.350e-5$ |  | $4.837e-6$ | $4.837e-6$
| $3.181e-4$ | $3.169e-4$ |  | $9.562e-2$ | $9.570e-2$
| $7.296e-5$ | $7.273e-5$ |  | $2.385e-2$ | $2.386e-2$
S2 | $1.700e-5$ | $1.695e-5$ |  | $5.978e-3$ | $5.980e-3$
| $4.072e-6$ | $4.060e-6$ |  | $1.515e-3$ | $1.514e-3$
| $9.946e-7$ | $9.917e-7$ |  | $3.805e-4$ | $3.802e-4$
| $2.457e-7$ | $2.450e-7$ |  | $3.805e-4$ | $3.802e-4$
Table 4: Error with $h=2^{-6}$ fixed and $\Delta
t=2^{-(i+1)},\leavevmode\nobreak\ i=1,\cdots,6$

In the Table 5 is the opposite, the parameter $\Delta t=2^{-7}$ is fixed and
the $h$ value are varying $h=2^{-i},\leavevmode\nobreak\ i=1,\cdots,6$.

$E_{L^{\infty}(0,T;L^{2}(\Omega))}$
---
| $\mathbb{R}$ |  | $\mathbb{R}^{2}$
| B1 | B2 |  | B1 | B2
| $5.463e-3$ | $5.462e-3$ |  | $9.895e-3$ | $9.895e-3$
| $2.750e-3$ | $2.751e-3$ |  | $2.105e-3$ | $2.105e-3$
S1 | $6.417e-4$ | $6.419e-4$ |  | $6.283e-4$ | $6.283e-4$
| $2.076e-4$ | $2.077e-4$ |  | $1.889e-4$ | $1.889e-4$
| $5.507e-5$ | $5.477e-5$ |  | $5.837e-5$ | $5.837e-5$
| $5.380e-5$ | $5.350e-5$ |  | $4.837e-5$ | $4.837e-5$
| $5.254e-5$ | $5.254e-5$ |  | $1.000e-9$ | $1.000e-9$
| $2.745e-5$ | $2.745e-5$ |  | $2.200e-10$ | $2.200e-10$
S2 | $6.383e-6$ | $6.383e-6$ |  | $7.000e-11$ | $7.000e-11$
| $1.984e-6$ | $1.984e-6$ |  | $2.335e-11$ | $2.335e-11$
| $2.460e-7$ | $2.453e-7$ |  | $2.205e-11$ | $2.205e-11$
| $2.457e-7$ | $2.450e-7$ |  | $2.001e-11$ | $2.001e-11$
Table 5: Error with $\Delta t=2^{-7}$ fixed and $h=2^{-i},\leavevmode\nobreak\
i=1,\cdots,6$

Finally, we present the order of numerical convergence in the Table 6,
$h=2\Delta t$ and varying $\Delta t=2^{-(i+1)},\leavevmode\nobreak\
i=1,\cdots,6$, which coincides with the expected theoretical results obtained
in Theorem 4 and Theorem 5 demonstrated in this paper, what suggests to us,
that these results are valid also for $\theta=1/4$.

Convergence Rate
---
| $\mathbb{R}$ |  | $\mathbb{R}^{2}$
| B1 | B2 |  | B1 | B2
| $-$ | $-$ |  | $-$ | $-$
| $2.466$ | $2.460$ |  | $2.208$ | $2.208$
S1 | $2.434$ | $2.429$ |  | $2.475$ | $2.475$
| $2.014$ | $2.031$ |  | $2.095$ | $2.095$
| $1.953$ | $1.954$ |  | $1.884$ | $1.884$
| $1.947$ | $1.947$ |  | $1.955$ | $1.955$
| $-$ | $-$ |  | $-$ | $-$
| $2.087$ | $2.085$ |  | $1.552$ | $1.552$
S2 | $2.093$ | $2.093$ |  | $1.922$ | $1.922$
| $2.061$ | $2.061$ |  | $2.002$ | $2.002$
| $2.037$ | $2.037$ |  | $1.982$ | $1.982$
| $2.018$ | $2.018$ |  | $2.009$ | $2.009$
Table 6: Convergence rate with $h=2\Delta t$ and $\Delta
t=2^{-(i+1)},\leavevmode\nobreak\ i=1,\cdots,6$

###  5.1 Asymptotic Behavior

In Figure 1 we present the decay of the energy of the homogeneous solution of
the mobile problem (5) for all $t\in[0,T^{*}]$, were $T^{*}$ is the time when
the energy is of the order of $10^{-10}$. We consider the initial conditions
taking $t=0$ in the exact functions of the Table 1 and the moving borders of
the Table 2.

(a) Example S1

(b) Example S2

Figure 1: Asymptotic behavior of energy

In the simulations we obtained an energy of the order of $10^{-10}$ at the
approximate times $T^{*}=770$, in example S1, and $T^{*}=560$, in example S2 ,
for both moving borders. In Figure 2 we present the evolution of the moving
boundaries.

Figure 2: Moving boundaries

###  5.2 Homogeneos Solution

Now we will present the evolution of the numerical solution of the problem
(6), obtained from the (18), by considering $f(y,t)=0$ and the initial
conditions from taking $t=0$ in exact functions of Table 1.

In the simulations were considered $h=2^{-6}$ and $\Delta t=2^{-7}$. The
following parameters were prefixed $\zeta_{0}=128$, $\zeta_{1}=2$ and $\nu=1$.

Figure 3: Homogeneous solutions for example S1

(a) Under the effect of the initial velocity

(b) Post-effect of initial velocity

Figure 4: Homogeneous solutions for example S2

##  6 Conclusions

In this paper, we focused our study on the numerical analysis of Problem (1),
introducing a numerical method and demonstrating the order of convergence for
both the semi-discrete problem and the totally discrete problem. Initially, we
announced the known theoretical results in the presence of moving boundary and
damping. In addition we announce the asymptotic decay. We introduce a family
of numeric methods, based on a theta parameter, and show numerically that for
$\theta\in[0,1/4[$, the numerical method is conditionally convergent (as
expected) and for $\theta\in[1/4,1]$ the convergence is unconditionally
convergent (see the Tables when $h$ or $\Delta t$ are varying). In addition,
we can verify that the numerical error is smaller the smaller the $\theta$
limited to $[1/4,1]$. Finally, let us show the decay of the system’s energy.

It is worth emphasizing that the presented numerical error results were in
complete agreement with the results presented in the theoretical part, and
they also indicate the reliability of this numerical method for obtaining good
approximate solutions of these nonlinear problems.

## References

  * [1] S. Bagheri, A. Nikkar e H. Ghaffarzadeh, “Study of nonlinear vibration of Euler-Bernoulli beams by using analytical approximate techniques”, Latin American Journal of Solids and Structures, 11 (2014) 157–168.
  * [2] P. Z. Bar-Yoseph, D. Fisher e O. Gottlieb, “Spectral element methods for nonlinear spatio-temporal dynamics of an Euler-Bernoulli beam”, Computational Mechanics, 19 (1996) 136–15.
  * [3] A. Barari, H. D. Kaliji, M.Ghadimi e G. Domairry, “Non-linear vibration of Euler-Bernoulli beams”, Latin American Journal of Solids and Structures, 8 (2011) [S.p.].
  * [4] P. Biler, “Remark on the decay for damped string and beam equations”, Nonlinear Analisys, 10 (1986) 839–842.
  * [5] C. S. Q. Caldas, J. Limaco e R. K. Barreto, “Remarks on the beam evolution equations in noncylindrical domains”, Nonlinear Analysis, 70 (2009) 693–710.
  * [6] S.-I. Chou e C.-C.Wang, “Error estimates of finite element approximations for problems in linear elasticity - Part 3. Problems in elastodynamics discrete time approximations”, Archive for Rational Mechanics and Analysis, 73 (1980) 159–182.
  * [7] P.G. Ciarlet, “The Finite Element Method for Elliptic Problems”, North-Holland, Amsterdam, 1978.
  * [8] H. R. Clark, M. A. Rincon e R. D. Rodrigues, “Beam Equation with weak-internal damping in domain with moving boudary”, Applied Numerical Mathematics, 47 (2003) 139–157.
  * [9] R. W. Dickey, “The initial value problem for a nonlinear semi-infinite string”, Proceedings of the Royal Society Edinburgh Section A, 82 (1978) 19–26.
  * [10] S. R. Gunakata, D. M. G. Comissiong, K. Jordan e A. Sankar, “A Finite Element Solution of the Beam Equation via Matlab”, International Journal of Applied Science and Technology, 2 (2012) 80–88.
  * [11] S. M. Han, H. Benaroya e T. Wei, “Dynamics of Transversely Vibrating Beams Using Four Engineering Theories”, Journal of Sound and Vibration, 225 (1999) 935–988.
  * [12] D. H. Hodgese e E. H. Dowell, “Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades”, [S.l.]; National Aeronautics and Space Adiministration (NASA), 1974 (NASA TN D-7818).
  * [13] D. H. Hodgese e R. A. Ormiston, “On the Nonlinear Deformation Geometry of Euler-Bernoulli”, [S.l.]; National Aeronautics and Space Adiministration (NASA), 1960 (NASA A-7985).
  * [14] G. Kirchhoff, “Vorlesungen über Mechanik”, Tauber, Leipzig, 1883.
  * [15] L. A. Medeiros, J. L. Ferrel e S. B. Menezes, “Vibrations of elastic strings: Mathematical aspects I and II”, Journal of Computational Analysis and Applications, 4(2) (2002) 91–127; 4(3) (2002) 211-263.
  * [16] S. Narendar, S. S. Gupta e S. Gopalakrishnan, “Wave Propagation in Single-Walled Carbon Nanotube Under Longitudinal Magnetic Field Using Nonlocal Euler-Bernoulli Beam Theory”, Applied Mathematical Modelling, 36 (2012) 4529–4538.
  * [17] G. S. Payette e J. N. Reddy, “Nonlinear quasi-static finite element formulations for viscoelastic Euler-Bernoulli and Timoshenko beams”, International Journal for Numerical Methods in Biomedical Engineering, 26 (2010) 1736–1755.
  * [18] D. C. Pereira, “Existence, uniqueness and asymptotic behavior for solutions of the nonlinear bema equation”, Nonlinear Analisys, 8 (1990) 613–623.
  * [19] M. Tucsnak, “Semi-internal Stabilization for a Nonlinear Bernoulli-Euler Equation”, Mathematical Methods in the Applied Sciences, 19 (1996) 897-907.
  * [20] O. Weeger, U. Wever and B. Simeon, “Isogeometric analysis of nonlinear Euler-Bernoulli beam vibrations”, Nonlinear Dynamics, 72 (2013) 813–835.