# Uniform boundary stabilization of a high-order finite element space discretization of the 1-d wave equation

Tiphaine Delaunay, Sébastien Imperiale, Philippe Moireau

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

## To cite this version:

Tiphaine Delaunay, Sébastien Imperiale, Philippe Moireau. Uniform boundary stabilization of a high-order finite element space discretization of the 1-d wave equation. 2024. hal-04172229v2

HAL Id: hal-04172229

https://inria.hal.science/hal-04172229v2

Preprint submitted on 3 May 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Uniform boundary stabilization of a
high-order finite element space
discretization of the 1-d wave equation

Tiphaine Delaunay¹,², Sébastien Imperiale ◦¹,²✉, Philippe
Moireau ◦¹,²

¹Inria, Team-MΞDISIM, Inria-Saclay Ile de France 91120 Palaiseau; ²LMS,
Ecole Polytechnique, CNRS – Institut Polytechnique de Paris

**Funding:** This work was
supported by ANR
Projet-ANR-19-CE48-0004
ODISSE

Abstract

The aim of this work is to construct and analyze a discretization process that preserves exponential stability at the discrete level for a wave propagation problem with boundary damping when a high-order spectral finite element approximation is used. The main idea is to add a stabilizing term to the wave equation that damps the spurious oscillatory components of the solutions. This term is based on a discrete multiplier analysis that ensures the exponential stability of the semi-discrete problem at any order without compromising the approximation properties. Possible applications of this work concern the solution of boundary control problems or data assimilation problems whose efficiency is based on an exponential stability property of an underlying discrete error system.

1 Introduction

The aim of this paper is to demonstrate uniform exponential stability of a discrete approximation of a one-dimensional wave propagation problem with boundary damping. More precisely, the asymptotic exponential stability is uniform with respect to the time and space discretization. This question typically arises when analyzing boundary control problems [13], data assimilation problems with boundary data [33] or absorbing boundary conditions [23].

Obtaining exponential stability results – by observability condition in some favorable cases – is a rather difficult task at the continuous level and has been extensively studied [16] with various techniques, typically multiplier techniques [1, 6, 18], Carleman estimates [24], spectral approaches [5, 20], or microlocal analyzes [4]. However, discretizing the dynamics so that the stability properties are preserved at the discrete level poses an additional difficulty, as spurious high frequencies are not attenuated when relying on the natural discretization of the stabilization terms [8, 16, 32, 34]. Several strategies – either for finite element or finite difference discretizations – have been proposed in the literature to overcome this problem, namely artificial viscosity [21], Fourier techniques [11, 12, 15], Tikhonov regularization [2], multigrid finite difference methods [14] or mixed finite element methods [17] and their extension for non-uniform meshes [27, 28].

When dealing with the high-performance simulation of wave propagation prob-
lems, high-order finite elements are now considered necessary. Among high-order
finite elements strategies, the high-order spectral finite elements method is widely
used since it provides mass-lumping at any order on quadrangles or hexahedra
meshes, see for instance [10, 26, 29]. However, the discretization strategies that
stabilize spurious high frequencies cannot be easily extended to such high-order dis-
cretizations without losing the advantages of the high-order spectral finite elements

<footnote>✉ For correspondence:
sebastien.imperiale@inria.fr</footnote>
# Uniform boundary stabilization of a high-order finite element space discretization of the 1-d wave equation

Tiphaine Delaunay<sup>1,2</sup>, Sébastien Imperiale<sup>1,2✉</sup>, Philippe Moireau<sup>1,2</sup>

<sup>1</sup>Inria, Team-MΞDISIM, Inria-Saclay Ile de France 91120 Palaiseau; <sup>2</sup>LMS, Ecole Polytechnique, CNRS – Institut Polytechnique de Paris

## Abstract

The aim of this work is to construct and analyze a discretization process that preserves exponential stability at the discrete level for a wave propagation problem with boundary damping when a high-order spectral finite element approximation is used. The main idea is to add a stabilizing term to the wave equation that damps the spurious oscillatory components of the solutions. This term is based on a discrete multiplier analysis that ensures the exponential stability of the semi-discrete problem at any order without compromising the approximation properties. Possible applications of this work concern the solution of boundary control problems or data assimilation problems whose efficiency is based on an exponential stability property of an underlying discrete error system.

## 1 Introduction

The aim of this paper is to demonstrate uniform exponential stability of a discrete approximation of a one-dimensional wave propagation problem with boundary damping. More precisely, the asymptotic exponential stability is uniform with respect to the time and space discretization. This question typically arises when analyzing boundary control problems [13], data assimilation problems with boundary data [33] or absorbing boundary conditions [23].

Obtaining exponential stability results – by observability condition in some favorable cases – is a rather difficult task at the continuous level and has been extensively studied [16] with various techniques, typically multiplier techniques [1, 6, 18], Carleman estimates [24], spectral approaches [5, 20], or microlocal analyzes [4]. However, discretizing the dynamics so that the stability properties are preserved at the discrete level poses an additional difficulty, as spurious high frequencies are not attenuated when relying on the natural discretization of the stabilization terms [8, 16, 32, 34]. Several strategies – either for finite element or finite difference discretizations – have been proposed in the literature to overcome this problem, namely artificial viscosity [21], Fourier techniques [11, 12, 15], Tikhonov regularization [2], multigrid finite difference methods [14] or mixed finite element methods [17] and their extension for non-uniform meshes [27, 28].

When dealing with the high-performance simulation of wave propagation problems, high-order finite elements are now considered necessary. Among high-order finite elements strategies, the high-order spectral finite elements method is widely used since it provides mass-lumping at any order on quadrangles or hexahedra meshes, see for instance [10, 26, 29]. However, the discretization strategies that stabilize spurious high frequencies cannot be easily extended to such high-order discretizations without losing the advantages of the high-order spectral finite elements

<footnote>✉ For correspondence: sebastien.imperiale@inria.fr</footnote>

<footnote>**Funding:** This work was supported by ANR Projet-ANR-19-CE48-0004 ODISSE</footnote>
method, namely convergence order or efficiency.

As a first step towards advanced settings (2D or 3D), in this paper we address the problem of constructing a space-time approximation of a one-dimensional wave-like system with uniform boundary damping using high-order spectral finite elements. The proposed numerical scheme preserves the exponential stability of the problem as well as the convergence order of the method and its efficiency.

We follow a strategy related to the artificial viscosity method, as we incorporate additional stabilization terms – see for example [35]– into the discretization process. These additional terms are constructed in such a way that the classical multiplier analysis – used to prove the exponential stability at the continuous level – can be carried out at the discrete level, hence achieving uniform stabilization with respect to the discretization parameters. Moreover, the stabilization terms are carefully designed to preserve the consistency order. Finally, our result can extend the discretization result for observers proposed in [31] to high-order schemes, which states that using boundary data, one can obtain discretization schemes that do not suffer from the asymptotic time drift of the error.

This work is structured as follows. In Section 2 the motivations for this work are explained in more detail, in particular we focus on spurious waves and their effects on the loss of exponential decay of the energy and on the eigenvalues of the system. In Section 3 we present the continuous setting and the method to prove exponential stability. In Section 4 we prove the uniform exponential stability in h for the semi-discrete problem by first building a general abstract framework and then applying it to our considered case. In Section 5 we describe an application to a data assimilation problem. Then, in Section 6, we present the fully discretized framework using an implicit scheme in time, extending the proof of exponential stability to this particular space-time discretization. The numerical results are presented in Section 7, which is divided into two parts. First, we present an eigenvalue analysis to observe the effects of the proposed stabilization term. In a second step, we introduce an explicit scheme for the time discretization. Using this scheme, we show the property that the additional viscosity terms do not affect the efficiency of the method as well as its high-order accuracy.

# 2 Main contributions

## 2.1 Problem setting

Let $\Omega = (0, 1)$, we study a homogenous wave equation formulated as a first-order system

$$ \begin{cases} \partial_t u(x,t) + \partial_x v(x,t) = 0, & (x,t) \in \Omega \times (0,T), \\ \partial_t v(x,t) + \partial_x u(x,t) = 0, & (x,t) \in \Omega \times (0,T). \end{cases} \qquad (1) $$

The following boundary conditions are used: Neumann atint $x = 0$ and a Robin condition at $x = 1$ with parameter $\gamma > 0$. The system is completed with initial conditions.

$$ \begin{cases} u(1,t) = \gamma v(1,t) & t \in (0,T), \\ v(0,t) = 0, & t \in (0,T), \end{cases} \qquad \begin{cases} u(x,0) = u_0(x), & x \in \Omega, \\ v(x,0) = v_0(x), & x \in \Omega. \end{cases} \qquad (2) $$

Note that the case $\gamma = 1$ corresponds to a transparent boundary condition, it is an absorbing boundary condition otherwise.

## 2.2 Loss of exponential energy decay

The question is whether such a system – which is known to be exponentially stable – retains its exponential stability property after a standard high-order spectral finite element discretization in space and an explicit second-order energy-preserving time discretization [10]. To illustrate, the problem is solved with finite elements (F.E.M.)
method, namely convergence order or efficiency.

As a first step towards advanced settings (2D or 3D), in this paper we address the problem of constructing a space-time approximation of a one-dimensional wave-like system with uniform boundary damping using high-order spectral finite elements. The proposed numerical scheme preserves the exponential stability of the problem as well as the convergence order of the method and its efficiency.

We follow a strategy related to the artificial viscosity method, as we incorporate additional stabilization terms – see for example [35]– into the discretization process. These additional terms are constructed in such a way that the classical multiplier analysis – used to prove the exponential stability at the continuous level – can be carried out at the discrete level, hence achieving uniform stabilization with respect to the discretization parameters. Moreover, the stabilization terms are carefully designed to preserve the consistency order. Finally, our result can extend the discretization result for observers proposed in [31] to high-order schemes, which states that using boundary data, one can obtain discretization schemes that do not suffer from the asymptotic time drift of the error.

This work is structured as follows. In Section 2 the motivations for this work are explained in more detail, in particular we focus on spurious waves and their effects on the loss of exponential decay of the energy and on the eigenvalues of the system. In Section 3 we present the continuous setting and the method to prove exponential stability. In Section 4 we prove the uniform exponential stability in h for the semi-discrete problem by first building a general abstract framework and then applying it to our considered case. In Section 5 we describe an application to a data assimilation problem. Then, in Section 6, we present the fully discretized framework using an implicit scheme in time, extending the proof of exponential stability to this particular space-time discretization. The numerical results are presented in Section 7, which is divided into two parts. First, we present an eigenvalue analysis to observe the effects of the proposed stabilization term. In a second step, we introduce an explicit scheme for the time discretization. Using this scheme, we show the property that the additional viscosity terms do not affect the efficiency of the method as well as its high-order accuracy.

# 2 Main contributions

## 2.1 Problem setting

Let $\Omega = (0, 1)$, we study a homogenous wave equation formulated as a first-order system

$$ \begin{cases} \partial_t u(x,t) + \partial_x v(x,t) = 0, & (x,t) \in \Omega \times (0,T), \\ \partial_t v(x,t) + \partial_x u(x,t) = 0, & (x,t) \in \Omega \times (0,T). \end{cases} \tag{1} $$

The following boundary conditions are used: Neumann atint $x = 0$ and a Robin condition at $x = 1$ with parameter $\gamma > 0$. The system is completed with initial conditions.

$$ \begin{cases} u(1,t) = \gamma v(1,t) & t \in (0,T), \\ v(0,t) = 0, & t \in (0,T), \end{cases} \quad \begin{cases} u(x,0) = u_0(x), & x \in \Omega, \\ v(x,0) = v_0(x), & x \in \Omega. \end{cases} \tag{2} $$

Note that the case $\gamma = 1$ corresponds to a transparent boundary condition, it is an absorbing boundary condition otherwise.

## 2.2 Loss of exponential energy decay

The question is whether such a system – which is known to be exponentially stable – retains its exponential stability property after a standard high-order spectral finite element discretization in space and an explicit second-order energy-preserving time discretization [10]. To illustrate, the problem is solved with finite elements (F.E.M.)
$\mathbb{P}_r$ – the set of polynomials of degree $r$ – and an explicit time scheme and on a uniform mesh. We initialize $u$ as a Gaussian centered in $x = 0.5$ and set $v$ as initial condition to 0. The boundary conditions of the wave propagation have two effects. At $x = 0$ the wave is reflected, while at $x = 1$ part of the wave is absorbed and part is reflected, an exponential convergence of the solution towards 0 should be observed. However, the discretization taints this property. We plot the evolution of the energy of the system over time on a logarithmic scale (Figure 1). We note that the decay rate is no longer exponential after a certain short time.

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Approximation of the logarithm of the square root of the energy obtained with $\gamma = 0.95$, $h = 1/10$ and $\Delta t = 1/1000$ and an explicit leap-frog scheme with different order of F.E.M</figcaption></figure>

## 2.3 Investigation of eigenvalues

As a key indicator of the behavior of the discrete system, we investigate the eigenvalues of the generator of the corresponding semi-group. We formulate (1) as follows,

$$ \dot{z} = Az, \text{ with } z = \begin{pmatrix} u \\ v \end{pmatrix}, \text{ and } A = \begin{pmatrix} 0 & -\partial_x \\ -\partial_x & 0 \end{pmatrix}, \qquad (3) $$

where $A : L^2(\Omega) \times L^2(\Omega) \to L^2(\Omega) \times L^2(\Omega)$ is an unbounded operator with domain $D(A) \subset H^1(\Omega) \times H^1(\Omega)$. The discrete approximation of the above systems is

$$ \dot{z}_h = A_h z_h, $$

where $A_h \in \mathcal{L}(V_h)$ is an approximation of the operator $A$ in a finite-dimensional space $V_h$ that is a subspace of $H^1(\Omega) \times L^2(\Omega)$ obtained using high-order spectral finite elements [10]. Since we look for exponential stability we aim at constructing approximations for which the eigenvalues of $A_h$ lies in the complex plane not close to the imaginary axis.

In Figure 2 (left) we see that for first-order approximations some eigenvalues are close to the imaginary axis, which means that the property of exponential stability is either not uniform in $h$ or even lost. By choosing $h$ smaller, the consistency is improved, i.e. the smallest eigenvalues – in amplitude – are approximated better and better, but the uniform exponential stability property is still not restored. The same behavior can be observed when using high-order approximations, see Figure 2 (right). In fact, it has been shown – at least for $\mathbb{P}_1$ finite elements (see [8, 32]) – that the property of uniform exponential stability at the discrete level cannot be achieved without adding additional dissipative terms. As already mentioned, the main contribution of this work is to construct such terms that stabilize the system – in the sense that the property of uniform exponential stability is guaranteed – and preserve the high order accuracy achieved with $\mathbb{P}_r$ finite elements.
$\mathbb{P}_r$ – the set of polynomials of degree $r$ – and an explicit time scheme and on a uniform mesh. We initialize $u$ as a Gaussian centered in $x = 0.5$ and set $v$ as initial condition to 0. The boundary conditions of the wave propagation have two effects. At $x = 0$ the wave is reflected, while at $x = 1$ part of the wave is absorbed and part is reflected, an exponential convergence of the solution towards 0 should be observed. However, the discretization taints this property. We plot the evolution of the energy of the system over time on a logarithmic scale (Figure 1). We note that the decay rate is no longer exponential after a certain short time.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Approximation of the logarithm of the square root of the energy obtained with $\gamma = 0.95$, $h = 1/10$ and $\Delta t = 1/1000$ and an explicit leap-frog scheme with different order of F.E.M</figcaption></figure>

## 2.3 Investigation of eigenvalues

As a key indicator of the behavior of the discrete system, we investigate the eigenvalues of the generator of the corresponding semi-group. We formulate (1) as follows,

$$ \dot{z} = Az, \text{ with } z = \begin{pmatrix} u \\ v \end{pmatrix}, \text{ and } A = \begin{pmatrix} 0 & -\partial_x \\ -\partial_x & 0 \end{pmatrix}, \qquad (3) $$

where $A: L^2(\Omega) \times L^2(\Omega) \to L^2(\Omega) \times L^2(\Omega)$ is an unbounded operator with domain $D(A) \subset H^1(\Omega) \times H^1(\Omega)$. The discrete approximation of the above systems is

$$ \dot{z}_h = A_h z_h, $$

where $A_h \in \mathcal{L}(V_h)$ is an approximation of the operator $A$ in a finite-dimensional space $V_h$ that is a subspace of $H^1(\Omega) \times L^2(\Omega)$ obtained using high-order spectral finite elements [10]. Since we look for exponential stability we aim at constructing approximations for which the eigenvalues of $A_h$ lies in the complex plane not close to the imaginary axis.

In Figure 2 (left) we see that for first-order approximations some eigenvalues are close to the imaginary axis, which means that the property of exponential stability is either not uniform in $h$ or even lost. By choosing $h$ smaller, the consistency is improved, i.e. the smallest eigenvalues – in amplitude – are approximated better and better, but the uniform exponential stability property is still not restored. The same behavior can be observed when using high-order approximations, see Figure 2 (right). In fact, it has been shown – at least for $\mathbb{P}_1$ finite elements (see [8, 32]) – that the property of uniform exponential stability at the discrete level cannot be achieved without adding additional dissipative terms. As already mentioned, the main contribution of this work is to construct such terms that stabilize the system – in the sense that the property of uniform exponential stability is guaranteed – and preserve the high order accuracy achieved with $\mathbb{P}_r$ finite elements.
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2. Spectrum of $A_h$ for various space discretizations.</figcaption></figure>

# 3 The standard continuous setting

This section presents standard multiplier techniques [6] to show the exponential stability property of System (1). We have chosen to present these classical results, first for the sake of completeness, second because they will be reproduced at the discrete level.

## 3.1 Observability of the undamped system

First, we define the variational formulation of the undamped system ($\gamma = 0$): find $(u(t), v(t)) \in H^1(0, 1) \times H^1(0, 1)$ such that for all $(\tilde{u}, \tilde{v}) \in H^1(0, 1) \times H^1(0, 1)$,

$$ \left\{ \begin{array}{l} \displaystyle \int_0^1 \partial_t u \tilde{u} \, dx - \int_0^1 v \partial_x \tilde{u} \, dx + v(1) \tilde{u}(1) = 0, \\ \displaystyle \int_0^1 \partial_t v \tilde{v} \, dx + \int_0^1 \partial_x u \tilde{v} \, dx - u(1) \tilde{v}(1) = 0, \end{array} \right. \quad \text{and} \quad \left\{ \begin{array}{ll} u(x, 0) = u_0(x), & x \in \Omega, \\ v(x, 0) = v_0(x), & x \in \Omega. \end{array} \right. \qquad (4) $$

For smooth initial data, it is standard to show that (4) has a unique solution,

$$ (u, v) \in C^1([0, T]; L^2(0, 1) \times L^2(0, 1)) \cap C^0([0, T]; H^1(0, 1) \times H^1(0, 1)). \quad (5) $$

Moreover, System (4) is conservative, an energy identity holds as recalled with the following theorem.

**Theorem 3.1.** For all solutions $(u, v)$ to (4), it holds

$$ \frac{d}{dt} \left( \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \right) = 0. \qquad (6) $$

The proof consists in choosing $\tilde{u} = u(t)$ and $\tilde{v} = v(t)$ in (4) and summing the two obtained equations. Note that, as an obvious consequence of (6) we have

$$ \int_0^T \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \, dt = T \left( \|u_0\|_{L^2(0,1)}^2 + \|v_0\|_{L^2(0,1)}^2 \right). \qquad (7) $$

A first step towards our objective is to derive an observability inequality ensuring that the initial condition of the wave propagation problem is controlled by the trace of $v$ at the boundary. This step is not mandatory to prove exponential stability but should rather be seen as a useful intermediate step. The strategy to prove such observability inequality relies here on multipliers [6].

**Theorem 3.2.** For $T$ large enough, there exists $c_{st}(T) > 0$ such that, for all $(u, v)$ solutions to (4),

$$ \frac{1}{2} \left( \|u_0\|_{L^2(0,1)}^2 + \|v_0\|_{L^2(0,1)}^2 \right) \le c_{st}(T) \int_0^T v^2(1, t) \, dt. \qquad (8) $$
<figure><img src="image_4.png" /><figcaption>Figure 2. Spectrum of A<sub>h</sub> for various space discretizations.</figcaption></figure>

3 The standard continuous setting

This section presents standard multiplier techniques [6] to show the exponential
stability property of System (1). We have chosen to present these classical results,
first for the sake of completeness, second because they will be reproduced at the
discrete level.

**3.1 Observability of the undamped system**

First, we define the variational formulation of the undamped system ($\gamma = 0$): find $(u(t), v(t)) \in H^1(0, 1) \times H^1(0, 1)$ such that for all $(\tilde{u}, \tilde{v}) \in H^1(0, 1) \times H^1(0, 1)$,

$$
\left\{
\begin{array}{l}
\displaystyle \int_0^1 \partial_t u \tilde{u} \, dx - \int_0^1 v \partial_x \tilde{u} \, dx + v(1) \tilde{u}(1) = 0, \\
\\
\displaystyle \int_0^1 \partial_t v \tilde{v} \, dx + \int_0^1 \partial_x u \tilde{v} \, dx - u(1) \tilde{v}(1) = 0,
\end{array}
\right.
\quad \text{and} \quad
\left\{
\begin{array}{ll}
u(x, 0) = u_0(x), & x \in \Omega, \\
v(x, 0) = v_0(x), & x \in \Omega.
\end{array}
\right.
\qquad (4)
$$

For smooth initial data, it is standard to show that (4) has a unique solution,

$$
(u, v) \in C^1([0, T]; L^2(0, 1) \times L^2(0, 1)) \cap C^0([0, T]; H^1(0, 1) \times H^1(0, 1)). \quad (5)
$$

Moreover, System (4) is conservative, an energy identity holdsas recalled with
the following theorem.

**Theorem 3.1.** For all solutions (u, v) to (4), it holds

$$
\frac{d}{dt} \left( \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \right) = 0. \quad (6)
$$

The proof consists in choosing $\tilde{u} = u(t)$ and $\tilde{v} = v(t)$ in (4) and summing the
two obtained equations. Note that, as an obvious consequence of (6) we have

$$
\int_0^T \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \, dt = T \left( \|u_0\|_{L^2(0,1)}^2 + \|v_0\|_{L^2(0,1)}^2 \right). \quad (7)
$$

A first step towards our objective is to derive an observability inequality ensuring
that the initial condition of the wave propagation problem is controlled by the trace
of *v* at the boundary. This step is not mandatory to prove exponential stability
but should rather be seen as a useful intermediate step. The strategy to prove such
observability inequality relies here on multipliers [6].

**Theorem 3.2.** For T large enough, there exists $c_{st}(T) > 0$ such that, for all $(u, v)$ solutions to (4),

$$
\frac{1}{2} \left( \|u_0\|_{L^2(0,1)}^2 + \|v_0\|_{L^2(0,1)}^2 \right) \le c_{st}(T) \int_0^T v^2(1,t) \, dt. \quad (8)
$$
Proof. The proof given below is separated in several steps that will all be mimicked at the discrete level.

**Step 1:** We choose $\tilde{u} = -xu(t)$ and $\tilde{v} = -xu(t)$ in (4) and sum the two obtained equations, integrating by parts and neglecting negative terms we obtain,

$$-\frac{d}{dt}\int_{0}^{1}xuv\,dx \leq -\frac{1}{2}(\|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2) + \frac{1}{2}v^2(1,t).$$

**Step 2:** Integrating in time and using the energy conservation property of 3.1 and its direct consequence (7) we get

$$-\left[\int_{0}^{1}xuv\,dx\right]_{0}^{T} \leq -\frac{T}{2}(\|u(0)\|_{L^2(0,1)}^2 + \|v(0)\|_{L^2(0,1)}^2) + \frac{1}{2}\int_{0}^{T}v^2(1,t)\,dt. \quad (9)$$

**Step 3:** We use Cauchy-Schwarz, Young inequalities and Theorem 3.1 to show

$$\left|\int_{0}^{1}xuv\,dx\right| \leq \frac{1}{2}\|v(t)\|_{L^2(0,1)}^2 + \frac{1}{2}\|v(t)\|_{L^2(0,1)}^2 = \frac{1}{2}\|v(0)\|_{L^2(0,1)}^2 + \frac{1}{2}\|v(0)\|_{L^2(0,1)}^2.$$

Finally, using this inequality together with (9) we obtain

$$(T-2) (\|u(0)\|_{L^2(0,1)}^2 + \|v(0)\|_{L^2(0,1)}^2) \leq \int_{0}^{T} v^2(1,t) \, dt. \quad \square$$

The observability inequality (8) does not directly guarantee that the system damped at the boundary is exponentially stable. It is however a good indicator since the damping term corresponds to the right-hand side of (8). This will be used as a clue to define the required damping terms at the discrete level.

## 3.2 Exponential stability for the damped system using Lyapunov functional

We define the variational formulation associated with System (1) with (2): find $(u(t), v(t)) \in H^1(0,1) \times H^1(0,1)$ such that for all $(\tilde{u}, \tilde{v}) \in H^1(0,1) \times H^1(0,1)$,

$$\left\{
\begin{array}{l}
\displaystyle \int_{0}^{1} \partial_{t} u \tilde{u} \, dx - \int_{0}^{1} v \partial_{x} \tilde{u} \, dx + v(1) \tilde{u}(1) = 0, \\
\\
\displaystyle \int_{0}^{1} \partial_{t} v \tilde{v} \, dx + \int_{0}^{1} \partial_{x} u \tilde{v} \, dx - \tilde{v}(1) u(1) + \gamma v(1) \tilde{v}(1) = 0,
\end{array}
\right.
\qquad (10)$$

with the same initial data as (4). For this variational formulation, we have existence and uniqueness of solutions with the regularity (5). We also have an energy relation for this system showing that the system is dissipative.

**Theorem 3.3.** For all solutions $(u, v)$ to (10), it holds

$$\frac{1}{2} \frac{d}{dt} (\|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2) = -\gamma v^2(1,t). \quad (11)$$

Proof. We choose $\tilde{u} = u(t)$ and $\tilde{v} = v(t)$ in (10). Then, adding the two equations in (10) concludes the proof. $\square$

We recall a well-known theorem (see [3] for instance) and prove it for further use.

**Theorem 3.4.** There exist two $c_s > 0$ and $\sigma > 0$ such that, for all $T > 0$ and all $(u, v)$ solutions to (10),

$$\forall t \in [0, T] \quad \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \leq c_s e^{-\sigma t} (\|u_0\|_{L^2(0,1)}^2 + \|v_0\|_{L^2(0,1)}^2). \quad (12)$$
*Proof.* The proof given below is separated in several steps that will all be mimicked at the discrete level.

**Step 1:** We choose $\tilde{u} = -xu(t)$ and $\tilde{v} = -xu(t)$ in (4) and sum the two obtained equations, integrating by parts and neglecting negative terms we obtain,

$$ -\frac{d}{dt} \int_{0}^{1} xuv \, dx \leq -\frac{1}{2} (\|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2) + \frac{1}{2} v^2(1,t). $$

**Step 2:** Integrating in time and using the energy conservation property of 3.1 and its direct consequence (7) we get

$$ - \left[ \int_{0}^{1} xuv \, dx \right]_{0}^{T} \leq -\frac{T}{2} \left( \|u(0)\|_{L^2(0,1)}^2 + \|v(0)\|_{L^2(0,1)}^2 \right) + \frac{1}{2} \int_{0}^{T} v^2(1,t) \, dt. \quad (9) $$

**Step 3:** We use Cauchy-Schwarz, Young inequalities and Theorem 3.1 to show

$$ \left| \int_{0}^{1} xuv \, dx \right| \leq \frac{1}{2} \|v(t)\|_{L^2(0,1)}^2 + \frac{1}{2} \|v(t)\|_{L^2(0,1)}^2 = \frac{1}{2} \|v(0)\|_{L^2(0,1)}^2 + \frac{1}{2} \|v(0)\|_{L^2(0,1)}^2. $$

Finally, using this inequality together with (9) we obtain

$$ (T - 2) (\|u(0)\|_{L^2(0,1)}^2 + \|v(0)\|_{L^2(0,1)}^2) \leq \int_{0}^{T} v^2(1,t) \, dt. \quad \square $$

The observability inequality (8) does not directly guarantee that the system damped at the boundary is exponentially stable. It is however a good indicator since the damping term corresponds to the right-hand side of (8). This will be used as a clue to define the required damping terms at the discrete level.

## 3.2 Exponential stability for the damped system using Lyapunov functional

We define the variational formulation associated with System (1) with (2): find $(u(t), v(t)) \in H^1(0,1) \times H^1(0,1)$ such that for all $(\tilde{u}, \tilde{v}) \in H^1(0,1) \times H^1(0,1)$,

$$ \left\{ \begin{array}{l} \displaystyle \int_{0}^{1} \partial_{t} u \tilde{u} \, dx - \int_{0}^{1} v \partial_{x} \tilde{u} \, dx + v(1) \tilde{u}(1) = 0, \\ \\ \displaystyle \int_{0}^{1} \partial_{t} v \tilde{v} \, dx + \int_{0}^{1} \partial_{x} u \tilde{v} \, dx - \tilde{v}(1) u(1) + \gamma v(1) \tilde{v}(1) = 0, \end{array} \right. \quad (10) $$

with the same initial data as (4). For this variational formulation, we have existence and uniqueness of solutions with the regularity (5). We also have an energy relation for this system showing that the system is dissipative.

**Theorem 3.3.** For all solutions $(u, v)$ to (10), it holds

$$ \frac{1}{2} \frac{d}{dt} (\|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2) = -\gamma v^2(1,t). \quad (11) $$

*Proof.* We choose $\tilde{u} = u(t)$ and $\tilde{v} = v(t)$ in (10). Then, adding the two equations in (10) concludes the proof. $\square$

We recall a well-known theorem (see [3] for instance) and prove it for further use.

**Theorem 3.4.** There exist two $c_s > 0$ and $\sigma > 0$ such that, for all $T > 0$ and all $(u, v)$ solutions to (10),

$$ \forall t \in [0, T] \quad \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \leq c_s e^{-\sigma t} \left( \|u_0\|_{L^2(0,1)}^2 + \|v_0\|_{L^2(0,1)}^2 \right). \quad (12) $$
Proof. The proof given below is separated in several steps that will all be mimicked at the discrete level.

**Step 1:** We choose $\tilde{u} = -xv(t)$ and $\tilde{v} = -xu(t)$ in (10).

Summing the two equations and integrating by parts, we get

$$-\frac{d}{dt} \int_{0}^{1} xuv \, dx = \frac{1}{2} v^{2}(1, t) - \frac{1}{2} u^{2}(1, t) + \gamma v(1, t)u(1, t) - \frac{1}{2} \left( \|u(t)\|_{L^{2}(0,1)}^{2} + \|v(t)\|_{L^{2}(0,1)}^{2} \right).$$

Moreover, thanks to Young's inequality we have $\gamma uv \le \frac{\gamma^2}{2} v^2 + \frac{1}{2} u^2$, finally we obtain

$$-\frac{d}{dt} \int_{0}^{1} xuv \, dx \le \frac{1 + \gamma^2}{2} v^2(1, t) - \frac{1}{2} \left( \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \right). \quad (13)$$

**Step 2:** We add to (13), the energy relation (11) multiplied by $\alpha > 0$,

$$\frac{d}{dt} \left( \frac{\alpha}{2} \left( \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \right) - \int_0^1 xuv \, dx \right) \\ \le \left( \frac{\gamma^2 + 1}{2} - \alpha\gamma \right) v^2(1, t) - \frac{1}{2} \left( \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \right).$$

**Step 3:** Choosing $\alpha$ large enough such that

$$\frac{\gamma^2 + 1}{2} - \alpha\gamma < 0,$$

we get

$$\frac{d}{dt} \left( \frac{\alpha}{2} \left( \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \right) - \int_0^1 xuv \, dx \right) \\ \le -\frac{1}{2} \left( \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \right). \quad (14)$$

**Step 4:** We introduce the Lyapunov functional

$$\mathcal{L}(t) = \frac{\alpha}{2} \left( \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \right) - \int_0^1 xuv \, dx. \quad (15)$$

Using Cauchy-Schwarz and Young inequalities, one can show the equivalence between the energy functional

$$\mathcal{E}(t) := \frac{1}{2} \|u(t)\|_{L^2(0,1)}^2 + \frac{1}{2} \|v(t)\|_{L^2(0,1)}^2$$

and the Lyapunov functional, more precisely,

$$(\alpha - 1)\mathcal{E}(t) \le \mathcal{L}(t) \le (\alpha + 1)\mathcal{E}(t). \quad (16)$$

Using (16) in (14), we get

$$\frac{d}{dt} \mathcal{L}(t) \le -\frac{1}{\alpha + 1} \mathcal{L}(t).$$

**Step 5:** Thanks to the inequation above and using again (16) we finally end up with the estimate

$$\mathcal{E}(t) \le \frac{\alpha + 1}{\alpha - 1} e^{-t(\alpha + 1)^{-1}} \mathcal{E}(0).$$

□
Proof. The proof given below is separated in several steps that will all be mimicked at the discrete level.

**Step 1:** We choose $\tilde{u} = -xv(t)$ and $\tilde{v} = -xu(t)$ in (10).

Summing the two equations and integrating by parts, we get

$$-\frac{d}{dt} \int_{0}^{1} xuv dx = \frac{1}{2} v^{2}(1,t) - \frac{1}{2} u^{2}(1,t) + \gamma v(1,t)u(1,t) - \frac{1}{2} \left( \|u(t)\|_{L^{2}(0,1)}^{2} + \|v(t)\|_{L^{2}(0,1)}^{2} \right).$$

Moreover, thanks to Young's inequality we have $\gamma uv \le \frac{\gamma^2}{2}v^2 + \frac{1}{2}u^2$, finally we obtain

$$-\frac{d}{dt} \int_{0}^{1} xuv \, dx \le \frac{1 + \gamma^2}{2} v^2(1,t) - \frac{1}{2} \left( \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \right). \quad (13)$$

**Step 2:** We add to (13), the energy relation (11) multiplied by $\alpha > 0$,

$$\frac{d}{dt} \left( \frac{\alpha}{2} \left( \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \right) - \int_0^1 xuv \, dx \right) \\ \le \left( \frac{\gamma^2 + 1}{2} - \alpha\gamma \right) v^2(1,t) - \frac{1}{2} \left( \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \right).$$

**Step 3:** Choosing $\alpha$ large enough such that

$$\frac{\gamma^2 + 1}{2} - \alpha\gamma < 0,$$

we get

$$\frac{d}{dt} \left( \frac{\alpha}{2} \left( \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \right) - \int_0^1 xuv \, dx \right) \\ \le -\frac{1}{2} \left( \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \right). \quad (14)$$

**Step 4:** We introduce the Lyapunov functional

$$\mathcal{L}(t) = \frac{\alpha}{2} \left( \|u(t)\|_{L^2(0,1)}^2 + \|v(t)\|_{L^2(0,1)}^2 \right) - \int_0^1 xuv \, dx. \quad (15)$$

Using Cauchy-Schwarz and Young inequalities, one can show the equivalence between the energy functional

$$\mathcal{E}(t) := \frac{1}{2} \|u(t)\|_{L^2(0,1)}^2 + \frac{1}{2} \|v(t)\|_{L^2(0,1)}^2$$

and the Lyapunov functional, more precisely,

$$(\alpha - 1)\mathcal{E}(t) \le \mathcal{L}(t) \le (\alpha + 1)\mathcal{E}(t). \quad (16)$$

Using (16) in (14), we get

$$\frac{d}{dt} \mathcal{L}(t) \le -\frac{1}{\alpha + 1} \mathcal{L}(t).$$

**Step 5:** Thanks to the inequation above and using again (16) we finally end up with the estimate

$$\mathcal{E}(t) \le \frac{\alpha + 1}{\alpha - 1} e^{-t(\alpha+1)^{-1}} \mathcal{E}(0).$$

□
# 4 Uniform stabilization in the semi-discrete setting

As already mentioned, we give a counterpart to Theorem 3.2 and Theorem 3.4 at the semi-discrete level. This is achieved in an abstract framework. Then we show – and this is the most important part of this work – that a well-chosen semi-discretization of the wave propagation problem (10) falls within the described framework.

## 4.1 Abstract framework

Let $\mathcal{Z}_h$ be a family of finite-dimensional Hilbert space parametrized by $h$ and associated with the scalar product $(\cdot, \cdot)_h$. We consider the semi-discretized abstract first-order system: for $z_{h,0}$ given, find $z_h(t) \in \mathcal{Z}_h$, solution to

$$ \begin{cases} \dot{z}_h = A_h z_h, & \text{in } [0, T], \\ z_h(0) = z_{h,0}, \end{cases} \qquad (17) $$

where $A_h \in \mathcal{L}(\mathcal{Z}_h)$ and, when $A_h$ comes from the discretisation of an unbounded operator, $\|A_h\|_h$ goes to infinity when $h$ goes to zero. We assume that the operator $A_h$ satisfies the following property

$$ (A_h z_h, z_h)_h = 0. \qquad (18) $$

Furthermore, we have classically the existence and uniqueness of smooth solutions $(u_h, v_h)$ in the space $C^k([0, T]; \mathcal{Z}_h)$, for all $k \ge 0$. Property (18) implies straightforwardly an energy preservation property.

**Theorem 4.1.** For all $z_h \in \mathcal{Z}_h$ solution to (17), it holds

$$ \frac{1}{2} \frac{d}{dt} \|z\|_h^2 = 0. $$

*Proof.* We multiply (17) by $z_h$, and since $A_h$ satisfies (18), the result follows. □

Let $\mathcal{Y}_h$ be a finite-dimensional space. We have seen – at the continuous level – that a dissipative term supported on the boundary provided exponential stability properties. By analogy, we therefore introduce an operator $C_h$, where $C_h \in \mathcal{L}(\mathcal{Z}_h, \mathcal{Y}_h)$, and consider the system

$$ \begin{cases} \dot{z}_h = (A_h - \gamma C_h^* C_h) z_h, & \text{in } [0, T], \\ z_h(0) = z_{h,0}, \end{cases} \qquad (19) $$

for some $\gamma > 0$ given. By construction, we have an energy dissipation property.

**Theorem 4.2.** Let $z_h \in \mathcal{Z}_h$ the solution of (19), we have the energy dissipation property

$$ \frac{1}{2} \frac{d}{dt} \|z_h\|_h^2 = -\gamma \|C_h z_h\|_h^2. $$

*Proof.* The proof is similar to the proof of Theorem 4.1. □

In order to prove an observability property, we need a fundamental hypothesis.

**Assumption 1.** $\exists \alpha, \beta, c_\pi > 0$, such that, for all $h$, $\exists \Pi_h \in \mathcal{L}(\mathcal{Z}_h)$ symmetric,

$$ \|\Pi_h\|_h \le c_\pi, \quad \text{and} \quad \forall z_h \in \mathcal{Z}_h, \quad (\Pi_h A_h z_h, z_h)_h \le \alpha \|C_h z_h\|_h^2 - \beta \|z_h\|_h^2. $$

From this assumption, we can establish the expected uniform observability property.

**Theorem 4.3.** Let Assumption 1 holds. For $T$ large enough, there exists a constant $c_{st}(T) > 0$ independent of $h$ such that, for any $z_h$ solution to (17),

$$ \|z_{h,0}\|_h^2 \le c_{st}(T) \int_0^T \|C_h z_h(t)\|_h^2 dt. $$
# 4 Uniform stabilization in the semi-discrete setting

As already mentioned, we give a counterpart to Theorem 3.2 and Theorem 3.4 at the semi-discrete level. This is achieved in an abstract framework. Then we show – and this is the most important part of this work – that a well-chosen semi-discretization of the wave propagation problem (10) falls within the described framework.

## 4.1 Abstract framework

Let $\mathcal{Z}_h$ be a family of finite-dimensional Hilbert space parametrized by $h$ and associated with the scalar product $(\cdot, \cdot)_h$. We consider the semi-discretized abstract first-order system: for $z_{h,0}$ given, find $z_h(t) \in \mathcal{Z}_h$, solution to

$$ \begin{cases} \dot{z}_h = A_h z_h, & \text{in } [0, T], \\ z_h(0) = z_{h,0}, \end{cases} \qquad (17) $$

where $A_h \in \mathcal{L}(\mathcal{Z}_h)$ and, when $A_h$ comes from the discretisation of an unbounded operator, $\|A_h\|_h$ goes to infinity when $h$ goes to zero. We assume that the operator $A_h$ satisfies the following property

$$ (A_h z_h, z_h)_h = 0. \qquad (18) $$

Furthermore, we have classically the existence and uniqueness of smooth solutions $(u_h, v_h)$ in the space $C^k([0, T]; \mathcal{Z}_h)$, for all $k \ge 0$. Property (18) implies straightforwardly an energy preservation property.

**Theorem 4.1.** For all $z_h \in \mathcal{Z}_h$ solution to (17), it holds

$$ \frac{1}{2} \frac{d}{dt} \|z\|_h^2 = 0. $$

*Proof.* We multiply (17) by $z_h$, and since $A_h$ satisfies (18), the result follows. □

Let $\mathcal{Y}_h$ be a finite-dimensional space. We have seen – at the continuous level – that a dissipative term supported on the boundary provided exponential stability properties. By analogy, we therefore introduce an operator $C_h$, where $C_h \in \mathcal{L}(\mathcal{Z}_h, \mathcal{Y}_h)$, and consider the system

$$ \begin{cases} \dot{z}_h = (A_h - \gamma C_h^* C_h) z_h, & \text{in } [0, T], \\ z_h(0) = z_{h,0}, \end{cases} \qquad (19) $$

for some $\gamma > 0$ given. By construction, we have an energy dissipation property.

**Theorem 4.2.** Let $z_h \in \mathcal{Z}_h$ the solution of (19), we have the energy dissipation property

$$ \frac{1}{2} \frac{d}{dt} \|z_h\|_h^2 = -\gamma \|C_h z_h\|_h^2. $$

*Proof.* The proof is similar to the proof of Theorem 4.1. □

In order to prove an observability property, we need a fundamental hypothesis.

**Assumption 1.** $\exists \alpha, \beta, c_\pi > 0$, such that, for all $h$, $\exists \Pi_h \in \mathcal{L}(\mathcal{Z}_h)$ symmetric,

$$ \|\Pi_h\|_h \le c_\pi, \quad \text{and} \quad \forall z_h \in \mathcal{Z}_h, \quad (\Pi_h A_h z_h, z_h)_h \le \alpha \|C_h z_h\|_h^2 - \beta \|z_h\|_h^2. $$

From this assumption, we can establish the expected uniform observability property.

**Theorem 4.3.** Let Assumption 1 holds. For $T$ large enough, there exists a constant $c_{st}(T) > 0$ independent of $h$ such that, for any $z_h$ solution to (17),

$$ \|z_{h,0}\|_h^2 \le c_{st}(T) \int_0^T \|C_h z_h(t)\|_h^2 dt. $$
**Proof.** The proof is separated into the same steps as in the proof of Theorem 3.1.

**Step 1:** We multiply the first equation of (17) by $\Pi_h z_h$, using the symmetry of $\Pi_h$, we obtain

$$ (\Pi_h \dot{z}_h, z_h)_h = (\Pi_h A_h z_h, z_h)_h. $$

Using Assumption 1, we get

$$ \frac{1}{2} \frac{d}{dt} (\Pi_h z_h, z_h)_h \leq \alpha \|C_h z_h\|_h^2 - \beta \|z_h\|_h^2. $$

**Step 2:** We integrate in time,

$$ \frac{1}{2} (\Pi_h z_h(T), z_h(T))_h - \frac{1}{2} (\Pi_h z_{h,0}, z_{0,h})_h \leq \alpha \int_0^T \|C_h z_h\|_h^2 dt - \beta \int_0^T \|z_h(t)\|_h^2 dt. $$

and use the energy conservation property of Theorem 4.1,

$$ \frac{1}{2} (\Pi_h z_h(T), z_h(T))_h - \frac{1}{2} (\Pi_h z_{h,0}, z_{0,h})_h \leq \alpha \int_0^T \|C_h z_h\|_h^2 dt - \beta T \|z_h(0)\|_h^2 dt. $$

**Step 3:** Since $\Pi_h$ is bounded, we have

$$ -\frac{1}{2} \|\Pi_h\|_h \|z_h(T)\|_h^2 - \frac{1}{2} \|\Pi_h\| \|z_{h,0}\|_h^2 \leq \alpha \int_0^T \|C_h z_h\|_h^2 dt - \beta T \|z_h(0)\|_h^2 dt. $$

Using again Theorem 4.1, we finally get

$$ (\beta T - \|\Pi_h\|_h) \|z_{0,h}\|_h^2 \leq \alpha \int_0^T \|C_h z_h\|_h^2 dt. $$

The result is obtained choosing $T > \|\Pi_h\|_h \beta^{-1}$. $\square$

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

We need a slightly stronger assumption than Assumption 1 to establish one of the main theorem of this work: a uniform exponential stability property for the damped system (19).

**Assumption 2.** $\exists \chi, \kappa, c_\pi > 0$, such that, for all $h$, $\exists \Pi_h \in \mathcal{L}(\mathcal{Z}_h)$ symmetric,

$$ \|\Pi_h\|_h \leq c_\pi, $$

and

$$ \forall z_h \in \mathcal{Z}_h, \quad (\Pi_h A_h z_h, z_h)_h - \gamma(C_h z_h, C_h \Pi_h z_h)_h \leq \chi \|C_h z_h\|_h^2 - \kappa \|z_h\|_h^2. $$

**Theorem 4.4.** Let Assumption 2 holds. There exist $c_s > 0$ and $\sigma > 0$ independent of $h$ such that, for all $T > 0$ and any $z_h$ solution to (19),

$$ \forall t \in [0, T] \quad \|z_h(t)\|_h \leq c_s \|z_{h,0}\|_h e^{-\sigma t}. $$

**Proof.** The proof follows the same steps as in the proof of Theorem 3.4.

**Step 1:** We multiply the first equation of (19) by $\Pi_h z_h$,

$$ (\Pi_h \dot{z}_h, z_h)_h = (\Pi_h A_h z_h, z_h)_h - \gamma(C_h z_h, C_h \Pi_h z_h)_h. $$

We use Assumption 2 and the symmetry of $\Pi_h$ to obtain

$$ \frac{1}{2} \frac{d}{dt} (\Pi_h z_h, z_h)_h \leq \chi \|C_h z_h\|_h^2 - \kappa \|z_h\|_h^2. $$

**Step 2:** Using the dissipation energy identity of Theorem 4.2 we get, for $\alpha > 0$,

$$ \frac{1}{2} \frac{d}{dt} ((\Pi_h z_h, z_h)_h + \alpha \|z_h\|_h^2) \leq (\chi - \alpha\gamma) \|C_h z_h\|_h^2 - \kappa \|z_h\|_h^2. $$
Proof. The proof is separated into the same steps as in the proof of Theorem 3.1.

**Step 1:** We multiply the first equation of (17) by $\Pi_h z_h$, using the symmetry of $\Pi_h$, we obtain

$$ (\Pi_h \dot{z}_h, z_h)_h = (\Pi_h A_h z_h, z_h)_h. $$

Using Assumption 1, we get

$$ \frac{1}{2} \frac{d}{dt} (\Pi_h z_h, z_h)_h \leq \alpha \|C_h z_h\|_h^2 - \beta \|z_h\|_h^2. $$

**Step 2:** We integrate in time,

$$ \frac{1}{2} (\Pi_h z_h(T), z_h(T))_h - \frac{1}{2} (\Pi_h z_{h,0}, z_{0,h})_h \leq \alpha \int_0^T \|C_h z_h\|_h^2 dt - \beta \int_0^T \|z_h(t)\|_h^2 dt. $$

and use the energy conservation property of Theorem 4.1,

$$ \frac{1}{2} (\Pi_h z_h(T), z_h(T))_h - \frac{1}{2} (\Pi_h z_{h,0}, z_{0,h})_h \leq \alpha \int_0^T \|C_h z_h\|_h^2 dt - \beta T \|z_h(0)\|_h^2 dt. $$

**Step 3:** Since $\Pi_h$ is bounded, we have

$$ -\frac{1}{2} \|\Pi_h\|_h \|z_h(T)\|_h^2 - \frac{1}{2} \|\Pi_h\| \|z_{h,0}\|_h^2 \leq \alpha \int_0^T \|C_h z_h\|_h^2 dt - \beta T \|z_h(0)\|_h^2 dt. $$

Using again Theorem 4.1, we finally get

$$ (\beta T - \|\Pi_h\|_h) \|z_{0,h}\|_h^2 \leq \alpha \int_0^T \|C_h z_h\|_h^2 dt. $$

The result is obtained choosing $T > \|\Pi_h\|_h \beta^{-1}$. $\square$

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

We need a slightly stronger assumption than Assumption 1 to establish one of the main theorem of this work: a uniform exponential stability property for the damped system (19).

**Assumption 2.** $\exists \chi, \kappa, c_\pi > 0$, such that, for all $h$, $\exists \Pi_h \in \mathcal{L}(\mathcal{Z}_h)$ symmetric,

$$ \|\Pi_h\|_h \leq c_\pi, $$

and

$$ \forall z_h \in \mathcal{Z}_h, \quad (\Pi_h A_h z_h, z_h)_h - \gamma(C_h z_h, C_h \Pi_h z_h)_h \leq \chi \|C_h z_h\|_h^2 - \kappa \|z_h\|_h^2. $$

**Theorem 4.4.** Let Assumption 2 holds. There exist $c_s > 0$ and $\sigma > 0$ independent of $h$ such that, for all $T > 0$ and any $z_h$ solution to (19),

$$ \forall t \in [0, T] \quad \|z_h(t)\|_h \leq c_s \|z_{h,0}\|_h e^{-\sigma t}. $$

Proof. The proof follows the same steps as in the proof of Theorem 3.4.

**Step 1:** We multiply the first equation of (19) by $\Pi_h z_h$,

$$ (\Pi_h \dot{z}_h, z_h)_h = (\Pi_h A_h z_h, z_h)_h - \gamma(C_h z_h, C_h \Pi_h z_h)_h. $$

We use Assumption 2 and the symmetry of $\Pi_h$ to obtain

$$ \frac{1}{2} \frac{d}{dt} (\Pi_h z_h, z_h)_h \leq \chi \|C_h z_h\|_h^2 - \kappa \|z_h\|_h^2. $$

**Step 2:** Using the dissipation energy identity of Theorem 4.2 we get, for $\alpha > 0$,

$$ \frac{1}{2} \frac{d}{dt} ((\Pi_h z_h, z_h)_h + \alpha \|z_h\|_h^2) \leq (\chi - \alpha \gamma) \|C_h z_h\|_h^2 - \kappa \|z_h\|_h^2. $$
**Step 3:** We choose $\alpha$ large enough such that $\chi - \alpha\gamma \le 0$. Thus, we have

$$ \frac{1}{2} \frac{d}{dt} \left( (\Pi_h z_h, z_h)_h + \alpha \|z_h\|_h^2 \right) \le -\kappa \|z_h\|_h^2. \quad (20) $$

**Step 4:** We introduce the following discrete Lyapunov functional

$$ \mathcal{L}_h(t) = \frac{1}{2}(\Pi_h z_h, z_h)_h + \frac{\alpha}{2} \|z_h\|_h^2. $$

Since $\Pi_h$ is bounded, we can establish an equivalence between the Lyapunov functional- and the energy functional $\mathcal{E}_h(t) := \|z_h(t)\|_h^2/2$,

$$ (-\|\Pi_h\|_h + \alpha) \mathcal{E}_h(t) \le \mathcal{L}_h(t) \le (\|\Pi_h\|_h + \alpha) \mathcal{E}_h(t). \quad (21) $$

Using (21) with $\alpha > \|\Pi_h\|_h$ and (20), we obtain

$$ \frac{d}{dt} \mathcal{L}_h(t) \le -2\kappa (\|\Pi_h\|_h + \alpha)^{-1} \mathcal{L}_h(t). $$

**Step 5:** The inequality above together with (21) leads to

$$ \|z_h(t)\|_h^2 \le \frac{\alpha + \|\Pi_h\|_h}{\alpha - \|\Pi_h\|_h} \|z_{h,0}\|_h^2 e^{-2\kappa(\|\Pi_h\|_h + \alpha)^{-1}t}. $$

□

## 4.2 Discretization with high-order spectral finite elements

In this section, we show how the semi-discretization in space of the wave propagation problem, using high-order spectral finite elements (1) enters the abstract framework presented above. For any small enough $h$ we assume being given a subdivision of the domain $(0, 1)$ using intervals $(x_i, x_{i+1})$, for $i \in \{0, \dots, N\}$. Moreover, we denote $h_i = |x_{i+1} - x_i|$ for $i \in \{0, \dots, N-1\}$, and assume that

$$ h = \max_{i \in \{0, \dots, N-1\}} h_i. $$

We finally assume that the sequence of meshes is quasi-uniform, i.e., there exists $\rho > 0$ such that

$$ \min_{i \in \{0, N-1\}} h_i \ge \rho h. $$

We then introduce two discrete spaces $U_h^r$ and $V_h^r$ – where $r$ is the order of the used finite elements – defined by

$$ U_h^r = \{u_h \in C^0(0,1) \mid u_h|_{(x_i, x_{i+1})} \in \mathbb{P}_r\}, $$

and

$$ V_h^r = \{v_h \in L^2(0,1) \mid v_h|_{(x_i, x_{i+1})} \in \mathbb{P}_r\}. $$

Note that $U_h^r \subset V_h^r$ and we denote $\mathcal{Z}_h = U_h^r \times V_h^r$. The efficiency of high-order spectral finite elements stems from the fact that mass-lumping can be achieved without loss of consistency using a quadrature rule based on Gauss-Lobatto points. We represent the use of this quadrature by the symbol

$$ \oint_0^1 f(x)dx \sim \int_0^1 f(x)dx \text{ with } \oint_0^1 f(x)dx := \sum_{i=0}^N \sum_{k=0}^r h_i w_k f(x_i + h_i \xi_k), $$

where the $\{w_k\}_{k=0}^r$ are the quadrature weights and $\xi_k \in [0, 1]$ the Gauss-Lobatto quadrature points, that are also interpolation points in the context of high-order spectral finite elements approximations. Note, for latter use, that $\xi_0 = 0$ and $\xi_r = 1$.
**Step 3:** We choose $\alpha$ large enough such that $\chi - \alpha\gamma \le 0$. Thus, we have

$$ \frac{1}{2} \frac{d}{dt} \left( (\Pi_h z_h, z_h)_h + \alpha \|z_h\|_h^2 \right) \le -\kappa \|z_h\|_h^2. \quad (20) $$

**Step 4:** We introduce the following discrete Lyapunov functional

$$ \mathcal{L}_h(t) = \frac{1}{2} (\Pi_h z_h, z_h)_h + \frac{\alpha}{2} \|z_h\|_h^2. $$

Since $\Pi_h$ is bounded, we can establish an equivalence between the Lyapunov functional- and the energy functional $\mathcal{E}_h(t) := \|z_h(t)\|_h^2/2$,

$$ (-\|\Pi_h\|_h + \alpha) \mathcal{E}_h(t) \le \mathcal{L}_h(t) \le (\|\Pi_h\|_h + \alpha) \mathcal{E}_h(t). \quad (21) $$

Using (21) with $\alpha > \|\Pi_h\|_h$ and (20), we obtain

$$ \frac{d}{dt} \mathcal{L}_h(t) \le -2\kappa (\|\Pi_h\|_h + \alpha)^{-1} \mathcal{L}_h(t). $$

**Step 5:** The inequality above together with (21) leads to

$$ \|z_h(t)\|_h^2 \le \frac{\alpha + \|\Pi_h\|_h}{\alpha - \|\Pi_h\|_h} \|z_{h,0}\|_h^2 e^{-2\kappa(\|\Pi_h\|_h + \alpha)^{-1}t}. $$

□

## 4.2 Discretization with high-order spectral finite elements

In this section, we show how the semi-discretization in space of the wave propagation problem, using high-order spectral finite elements (1) enters the abstract framework presented above. For any small enough $h$ we assume being given a subdivision of the domain $(0, 1)$ using intervals $(x_i, x_{i+1})$, for $i \in \{0, \dots, N\}$. Moreover, we denote $h_i = |x_{i+1} - x_i|$ for $i \in \{0, \dots, N-1\}$, and assume that

$$ h = \max_{i \in \{0, \dots, N-1\}} h_i. $$

We finally assume that the sequence of meshes is quasi-uniform, i.e., there exists $\rho > 0$ such that

$$ \min_{i \in \{0, N-1\}} h_i \ge \rho h. $$

We then introduce two discrete spaces $U_h^r$ and $V_h^r$ – where $r$ is the order of the used finite elements – defined by

$$ U_h^r = \{u_h \in C^0(0,1) \mid u_h|_{(x_i, x_{i+1})} \in \mathbb{P}_r\}, $$

and

$$ V_h^r = \{v_h \in L^2(0,1) \mid v_h|_{(x_i, x_{i+1})} \in \mathbb{P}_r\}. $$

Note that $U_h^r \subset V_h^r$ and we denote $\mathcal{Z}_h = U_h^r \times V_h^r$. The efficiency of high-order spectral finite elements stems from the fact that mass-lumping can be achieved without loss of consistency using a quadrature rule based on Gauss-Lobatto points. We represent the use of this quadrature by the symbol

$$ \int_0^1 f(x)dx \sim \int_0^1 f(x)dx \text{ with } \int_0^1 f(x)dx := \sum_{i=0}^N \sum_{k=0}^r h_i w_k f(x_i + h_i \xi_k), $$

where the $\{w_k\}_{k=0}^r$ are the quadrature weights and $\xi_k \in [0, 1]$ the Gauss-Lobatto quadrature points, that are also interpolation points in the context of high-order spectral finite elements approximations. Note, for latter use, that $\xi_0 = 0$ and $\xi_r = 1$.
We have the following bounds (see [19]), there exists $c_{q,\max} > 0$ and $c_{q,\min} > 0$ such that

$$ \forall u \in \mathbb{P}_r, \quad c_{q,\min} \int_0^1 u^2 dx \le \oint_0^1 u^2 dx \le c_{q,\max} \int_0^1 u^2 dx. \quad (22) $$

The semi-discrete formulation of the damped system (10) reads: Find $(u_h, v_h)(t) \in \mathcal{Z}_h$ such that for $(\tilde{u}_h, \tilde{v}_h) \in \mathcal{Z}_h$,

$$ \left\{ \begin{array}{l} \oint_0^1 \partial_t u_h \tilde{u}_h \, dx - \int_0^1 v_h \partial_x \tilde{u}_h \, dx + v_h(1) \tilde{u}_h(1) = 0, \\ \\ \oint_0^1 \partial_t v_h \tilde{v}_h \, dx + \int_0^1 \partial_x u_h \tilde{v}_h \, dx - u_h(1) \tilde{v}_h(1) + \gamma v_h(1) \tilde{v}_h(1) = 0, \end{array} \right. \quad (23) $$

with initial data

$$ u_h(0) = u_{h,0}, \quad v_h(0) = v_{h,0}. $$

Note that when $\gamma = 0$, one recovers the semi-discrete variational formulation of the undamped system (4). In order to relate to the abstract theory presented in Section 4.1, we reformulate the variational formulation (23) as an operator equation. To do so, we introduce the scalar product in $V_h^r$ by

$$ \forall (v_h, \tilde{v}_h) \in V_h^r \times V_h^r, \quad (v_h, \tilde{v}_h)_h := \oint_0^1 v_h \tilde{v}_h \, dx. \quad (24) $$

The scalar product in $U_h^r \subset V_h^r$ is defined similarly, finally we extend this definition to any elements $z_h = (u_h, v_h)^t$ and $\tilde{z}_h = (\tilde{u}_h, \tilde{v}_h)^t$ in $\mathcal{Z}_h = U_h^r \times V_h^r$,

$$ (z_h, \tilde{z}_h)_h = \oint_0^1 u_h \tilde{u}_h \, dx + \oint_0^1 v_h \tilde{v}_h \, dx. $$

The variational formulation (23) is equivalent to

$$ \left\{ \begin{array}{l} \partial_t u_h - R_h^* v_h = 0, \\ \partial_t v_h + R_h u_h + \gamma B_h v_h = 0, \end{array} \right. \quad (u_h, v_h)(0) = (u_{h,0}, v_{h,0}) := z_{h,0}, \quad (25) $$

where, $R_h \in \mathcal{L}(U_h^r, V_h^r)$, $B_h \in \mathcal{L}(V_h^r)$ and for all $(u_h, v_h, \tilde{v}_h) \in U_h^r \times V_h^r \times V_h^r$,

$$ (R_h u_h, \tilde{v}_h)_h = \int_0^1 \tilde{v}_h \frac{d}{dx} u_h \, dx - u_h(1) \tilde{v}_h(1) \quad \text{and} \quad (B_h v_h, \tilde{v}_h)_h = v_h(1) \tilde{v}_h(1). \quad (26) $$

Denoting $z_h = (u_h, v_h)^t$, we can rewrite System (25) as

$$ \left\{ \begin{array}{l} \dot{z}_h = A_h z_h - \gamma \begin{pmatrix} 0 & 0 \\ 0 & B_h \end{pmatrix} z_h, \quad \text{where} \quad A_h = \begin{pmatrix} 0 & R_h^* \\ -R_h & 0 \end{pmatrix}. \quad (27) \\ z_h(0) = z_{h,0}, \end{array} \right. $$

Note that the operator $A_h \in \mathcal{L}(\mathcal{Z}_h)$ satisfies, for all $z_h \in \mathcal{Z}_h$, $(A_h z_h, z_h)_h = 0$, hence is skew-adjoint.

## 4.3 Discrete observability

Our objective is – in a preliminary step – to prove an observability property for the system (27). To do so, Assumption 1 is verified in what follows for a given appropriate operator $\Pi_h \in \mathcal{L}(\mathcal{Z}_h)$.
We have the following bounds (see [19]), there exists $c_{q,\max} > 0$ and $c_{q,\min} > 0$ such that

$$ \forall u \in \mathbb{P}_r, \quad c_{q,\min} \int_0^1 u^2 dx \le \oint_0^1 u^2 dx \le c_{q,\max} \int_0^1 u^2 dx. \quad (22) $$

The semi-discrete formulation of the damped system (10) reads: Find $(u_h, v_h)(t) \in \mathcal{Z}_h$ such that for $(\tilde{u}_h, \tilde{v}_h) \in \mathcal{Z}_h$,

$$ \left\{ \begin{array}{l} \oint_0^1 \partial_t u_h \tilde{u}_h \, dx - \int_0^1 v_h \partial_x \tilde{u}_h \, dx + v_h(1) \tilde{u}_h(1) = 0, \\ \\ \oint_0^1 \partial_t v_h \tilde{v}_h \, dx + \int_0^1 \partial_x u_h \tilde{v}_h \, dx - u_h(1) \tilde{v}_h(1) + \gamma v_h(1) \tilde{v}_h(1) = 0, \end{array} \right. \quad (23) $$

with initial data

$$ u_h(0) = u_{h,0}, \quad v_h(0) = v_{h,0}. $$

Note that when $\gamma = 0$, one recovers the semi-discrete variational formulation of the undamped system (4). In order to relate to the abstract theory presented in Section 4.1, we reformulate the variational formulation (23) as an operator equation. To do so, we introduce the scalar product in $V_h^r$ by

$$ \forall (v_h, \tilde{v}_h) \in V_h^r \times V_h^r, \quad (v_h, \tilde{v}_h)_h := \oint_0^1 v_h \tilde{v}_h \, dx. \quad (24) $$

The scalar product in $U_h^r \subset V_h^r$ is defined similarly, finally we extend this definition to any elements $z_h = (u_h, v_h)^t$ and $\tilde{z}_h = (\tilde{u}_h, \tilde{v}_h)^t$ in $\mathcal{Z}_h = U_h^r \times V_h^r$,

$$ (z_h, \tilde{z}_h)_h = \oint_0^1 u_h \tilde{u}_h \, dx + \oint_0^1 v_h \tilde{v}_h \, dx. $$

The variational formulation (23) is equivalent to

$$ \left\{ \begin{array}{l} \partial_t u_h - R_h^* v_h = 0, \\ \partial_t v_h + R_h u_h + \gamma B_h v_h = 0, \end{array} \right. \quad (u_h, v_h)(0) = (u_{h,0}, v_{h,0}) := z_{h,0}, \quad (25) $$

where, $R_h \in \mathcal{L}(U_h^r, V_h^r)$, $B_h \in \mathcal{L}(V_h^r)$ and for all $(u_h, v_h, \tilde{v}_h) \in U_h^r \times V_h^r \times V_h^r$,

$$ (R_h u_h, \tilde{v}_h)_h = \int_0^1 \tilde{v}_h \frac{d}{dx} u_h \, dx - u_h(1) \tilde{v}_h(1) \quad \text{and} \quad (B_h v_h, \tilde{v}_h)_h = v_h(1) \tilde{v}_h(1). \quad (26) $$

Denoting $z_h = (u_h, v_h)^t$, we can rewrite System (25) as

$$ \left\{ \begin{array}{l} \dot{z}_h = A_h z_h - \gamma \begin{pmatrix} 0 & 0 \\ 0 & B_h \end{pmatrix} z_h, \\ z_h(0) = z_{h,0}, \end{array} \right. \quad \text{where} \quad A_h = \begin{pmatrix} 0 & R_h^* \\ -R_h & 0 \end{pmatrix}. \quad (27) $$

Note that the operator $A_h \in \mathcal{L}(\mathcal{Z}_h)$ satisfies, for all $z_h \in \mathcal{Z}_h$, $(A_h z_h, z_h)_h = 0$, hence is skew-adjoint.

## 4.3 Discrete observability

Our objective is – in a preliminary step – to prove an observability property for the system (27). To do so, Assumption 1 is verified in what follows for a given appropriate operator $\Pi_h \in \mathcal{L}(\mathcal{Z}_h)$.
### 4.3.1 The discrete multiplier

The strategy to define $\Pi_h$ is based upon our analysis presented at the continuous level. More precisely, we construct an operator $\Pi_h$ that is a discrete version of the multiplier used in the continuous case – see the proof of Theorem 3.2 in Section 3 – namely

$$ \Pi_h z_h \approx \begin{pmatrix} -xv_h \\ -xu_h \end{pmatrix}. $$

Note that the previous equation above can not be an equality since in general $xv_h \notin U_h^r$ and $xu_h \notin V_h^r$. A simple remedy is to introduce the interpolation operator $\chi_h \in \mathcal{L}(U_h^r, V_h^r)$ defined by

$$ \forall (u_h, v_h) \in U_h^r \times V_h^r, \quad \oint_0^1 \chi_h(u_h)v_h \, dx = \oint_0^1 xu_h v_h \, dx, \qquad (28) $$

and its adjoint $\chi_h^* \in \mathcal{L}(V_h^r, U_h^r)$ given by

$$ \forall (u_h, v_h) \in U_h^r \times V_h^r, \quad \oint_0^1 \chi_h^*(v_h)u_h \, dx = \oint_0^1 xv_h u_h \, dx. \qquad (29) $$

Due to the specific property of having quadrature points matching interpolation points, it can be shown that the range of $\chi_h$ belongs to $U_h^r$, in particular, for all $u_h \in U_h^r$, $\chi_h(u_h)$ is a continuous function, moreover, $\chi_h(u_h)(x_i) = x_i u_h(x_i)$. We then define the operator $\Pi_h$ as follows

$$ \forall z_h = (u_h, v_h)^t \in \mathcal{Z}_h \quad \Pi_h z_h = - \begin{pmatrix} \chi_h^*(v_h) \\ \chi_h(u_h) \end{pmatrix}. \qquad (30) $$

Note that using (22) one can show that $\Pi_h$ is uniformly bounded with respect to $h$. More precisely, there exists $c_\pi$ independent of $h$ such that

$$ \forall z_h \in \mathcal{Z}_h \quad \|\Pi_h z_h\|_h \le c_\pi \|z_h\|_h. \qquad (31) $$

Note also that $\Pi_h$ is symmetric by construction. The discrepancy between the term $\Pi_h z_h$ and its continuous version is used in what follows. Namely, we need to study

$$ \begin{pmatrix} xu_h \\ xv_h \end{pmatrix} - \begin{pmatrix} \chi_h^*(v_h) \\ \chi_h(u_h) \end{pmatrix}. $$

It turns out that the operator $C_h$ – appearing in Assumption 1 and given later in Eq. (45) and that appears in Assumption 1 – will involve the $H_1$ semi-norm of this discrepancy, on each interval of the subdivision of the propagation domain. That is why, we introduce the following two bilinear forms $d_h$ and $d_{h,*}$, for all $(u_h, \tilde{u}_h) \in U_h^r \times U_h^r$

$$ d_h(u_h, \tilde{u}_h) := \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} \frac{d}{dx} (xu_h - \chi_h(u_h)) \frac{d}{dx} (x\tilde{u}_h - \chi_h(\tilde{u}_h)) \, dx, \qquad (32) $$

and, for all $(v_h, \tilde{v}_h) \in V_h^r \times V_h^r$,

$$ d_{h,*}(v_h, \tilde{v}_h) := \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} \frac{d}{dx} (xv_h - \chi_h^*(v_h)) \frac{d}{dx} (x\tilde{v}_h - \chi_h^*(\tilde{v}_h)) \, dx. \qquad (33) $$

**Proposition 4.5.** Let $d_h$ and $d_{h,*}$ be the positive and symmetric bilinear forms defined by (32) and (33). Then, there exists a constant $c_d > 0$ depending only on $r$ such that

$$ \forall (u_h, v_h) \in U_h^r \times V_h^r, \quad d_h(u_h, u_h) + d_{h,*}(v_h, v_h) \le c_d \int_0^1 u_h^2 + v_h^2 \, dx. $$
### 4.3.1 The discrete multiplier

The strategy to define $\Pi_h$ is based upon our analysis presented at the continuous level. More precisely, we construct an operator $\Pi_h$ that is a discrete version of the multiplier used in the continuous case – see the proof of Theorem 3.2 in Section 3 – namely

$$ \Pi_h z_h \approx \begin{pmatrix} -xv_h \\ -xu_h \end{pmatrix}. $$

Note that the previous equation above can not be an equality since in general $xv_h \notin U_h^r$ and $xu_h \notin V_h^r$. A simple remedy is to introduce the interpolation operator $\chi_h \in \mathcal{L}(U_h^r, V_h^r)$ defined by

$$ \forall (u_h, v_h) \in U_h^r \times V_h^r, \quad \oint_0^1 \chi_h(u_h)v_h \, dx = \oint_0^1 xu_h v_h \, dx, \qquad (28) $$

and its adjoint $\chi_h^* \in \mathcal{L}(V_h^r, U_h^r)$ given by

$$ \forall (u_h, v_h) \in U_h^r \times V_h^r, \quad \oint_0^1 \chi_h^*(v_h)u_h \, dx = \oint_0^1 xv_h u_h \, dx. \qquad (29) $$

Due to the specific property of having quadrature points matching interpolation points, it can be shown that the range of $\chi_h$ belongs to $U_h^r$, in particular, for all $u_h \in U_h^r$, $\chi_h(u_h)$ is a continuous function, moreover, $\chi_h(u_h)(x_i) = x_i u_h(x_i)$. We then define the operator $\Pi_h$ as follows

$$ \forall z_h = (u_h, v_h)^t \in \mathcal{Z}_h \quad \Pi_h z_h = - \begin{pmatrix} \chi_h^*(v_h) \\ \chi_h(u_h) \end{pmatrix}. \qquad (30) $$

Note that using (22) one can show that $\Pi_h$ is uniformly bounded with respect to $h$. More precisely, there exists $c_\pi$ independent of $h$ such that

$$ \forall z_h \in \mathcal{Z}_h \quad \|\Pi_h z_h\|_h \le c_\pi \|z_h\|_h. \qquad (31) $$

Note also that $\Pi_h$ is symmetric by construction. The discrepancy between the term $\Pi_h z_h$ and its continuous version is used in what follows. Namely, we need to study

$$ \begin{pmatrix} xu_h \\ xv_h \end{pmatrix} - \begin{pmatrix} \chi_h^*(v_h) \\ \chi_h(u_h) \end{pmatrix}. $$

It turns out that the operator $C_h$ – appearing in Assumption 1 and given later in Eq. (45) and that appears in Assumption 1 – will involve the $H_1$ semi-norm of this discrepancy, on each interval of the subdivision of the propagation domain. That is why, we introduce the following two bilinear forms $d_h$ and $d_{h,*}$, for all $(u_h, \tilde{u}_h) \in U_h^r \times U_h^r$

$$ d_h(u_h, \tilde{u}_h) := \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} \frac{d}{dx} (xu_h - \chi_h(u_h)) \frac{d}{dx} (x\tilde{u}_h - \chi_h(\tilde{u}_h)) \, dx, \qquad (32) $$

and, for all $(v_h, \tilde{v}_h) \in V_h^r \times V_h^r$,

$$ d_{h,*}(v_h, \tilde{v}_h) := \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} \frac{d}{dx} (xv_h - \chi_h^*(v_h)) \frac{d}{dx} (x\tilde{v}_h - \chi_h^*(\tilde{v}_h)) \, dx. \qquad (33) $$

**Proposition 4.5.** Let $d_h$ and $d_{h,*}$ be the positive and symmetric bilinear forms defined by (32) and (33). Then, there exists a constant $c_d > 0$ depending only on $r$ such that

$$ \forall (u_h, v_h) \in U_h^r \times V_h^r, \quad d_h(u_h, u_h) + d_{h,*}(v_h, v_h) \le c_d \int_0^1 u_h^2 + v_h^2 \, dx. $$
Proof. Let us first estimate $d_h(u_h, u_h)$. We use an interpolation property [22, Chapter 4, Section 4, Theorem 4.4.4]. Since $xu_h|_{K_i}$ belongs to $\mathbb{P}_{r+1}$ with $K_i = (x_i, x_{i+1})$ for any $i \in \{0, N-1\}$, we have

$$ |xu_h - \chi_h(u_h)|_{H^1(K_i)} \leq h_i^r c_r |xu_h|_{H^{r+1}(K_i)}, \quad (34) $$

where $c_r > 0$ depends only on $r$. Now, observe that, as $x \in \mathbb{P}_1$, only its first derivative does not vanish. Moreover as, $u_h \in \mathbb{P}_r$, its $(k+1)$th derivative vanishes. Therefore, the general Leibniz rule leads to

$$ (xu_h)^{(r+1)} = \sum_{j=0}^{r+1} \binom{r+1}{j} x^{(j)} u_h^{(r+1-j)} = (r+1)u_h^{(r)}. $$

Injecting this result in (34), we get

$$ |xu_h - \chi_h(u_h)|_{H^1(K_i)} \leq h_i^r c_r (r+1) |u_h|_{H^r(K_i)}. $$

Thus, using an inverse estimate [22, Chapter 4, Section 5, Lemma 4.5.3], we have that there exists $c_r > 0$, a scalar depending only on $r$, such that

$$ |xu_h - \chi_h(u_h)|_{H^1(K_i)} \leq c_r \|u_h\|_{L^2(K_i)} \quad \Rightarrow \quad d_h(u_h, u_h) \leq c_r^2 \|u_h\|_{L^2(0,1)}^2. $$

The estimate for $d_{h,*}(v_h, v_h)$ can be obtained similarly. □

### 4.3.2 Jump and averages at intervals end-points

Functions in the space $V_h^r$ are discontinuous at each end-points of the intervals $(x_i, x_{i+1})$. Therefore, we introduce the following notations: for all $i \in \{1, N-1\}$, and for any function $v_h \in V_h^r$,

$$ \|v_h(x_i)\| := v_h(x_i^+) - v_h(x_i^-), \quad \|v_h(x_i)\| := \frac{v_h(x_i^+) + v_h(x_i^-)}{2} \quad (35) $$

where,

$$ v_h(x_i^+) = \lim_{\epsilon \to 0} v_h(x_i + \epsilon) \quad \text{and} \quad v_h(x_i^-) = \lim_{\epsilon \to 0} v_h(x_i - \epsilon). $$

The following preliminary estimate is required.

**Lemma 4.6.** There exists $c_a > 0$, such that, for all $h > 0$,

$$ \forall v_h \in V_h^r, \quad \sum_{i=1}^{N-1} \tilde{h}_i \|\{v_h(x_i)\}\|^2 \leq c_a \|v_h\|_{L^2(0,1)}^2 \quad \text{where} \quad \tilde{h}_i = \min(h_{i-1}, h_i) $$

*Proof.* Using the definition of the operator $\|\cdot\|$ and $(a+b)^2 \leq 2(a^2 + b^2)$, we have

$$ \sum_{i=1}^{N-1} \tilde{h}_i \|\{v_h(x_i)\}\|^2 \leq \frac{1}{2} \sum_{i=1}^{N-1} \tilde{h}_i (v_h(x_i^+)^2 + v_h(x_i^-)^2). \quad (36) $$

On any interval $(0, L)$, a trace inequality – [22, Chapter 1, Section 6] – holds. In a one-dimensionnal setting it reads,

$$ \forall v \in H^1(0, L), \quad |v(0)|^2 \leq \|v\|_{L^2(0,L)} \left( \frac{1}{L} \|v\|_{L^2(0,L)} + 2 \left\| \frac{d}{dx} v \right\|_{L^2(0,L)} \right). $$

Using this estimate for a fonction $v_h \in \mathbb{P}_r$ as well as an inverse estimate [22, Chapter 4, Section 5, Lemma 4.5.3], we get, with $c_r > 0$ depending only on $r$,

$$ \forall v_h \in \mathbb{P}_r, \quad |v_h^2(0)|^2 \leq \left( \frac{1 + 2c_r}{L} \right) \|v_h\|_{L^2(0,L)}^2. $$
*Proof.* Let us first estimate $d_h(u_h, u_h)$. We use an interpolation property [22, Chapter 4, Section 4, Theorem 4.4.4]. Since $xu_h|_{K_i}$ belongs to $\mathbb{P}_{r+1}$ with $K_i = (x_i, x_{i+1})$ for any $i \in \{0, N-1\}$, we have

$$ |xu_h - \chi_h(u_h)|_{H^1(K_i)} \leq h_i^r c_r |xu_h|_{H^{r+1}(K_i)}, \quad (34) $$

where $c_r > 0$ depends only on $r$. Now, observe that, as $x \in \mathbb{P}_1$, only its first derivative does not vanish. Moreover as, $u_h \in \mathbb{P}_r$, its $(k+1)$th derivative vanishes. Therefore, the general Leibniz rule leads to

$$ (xu_h)^{(r+1)} = \sum_{j=0}^{r+1} \binom{r+1}{j} x^{(j)} u_h^{(r+1-j)} = (r+1)u_h^{(r)}. $$

Injecting this result in (34), we get

$$ |xu_h - \chi_h(u_h)|_{H^1(K_i)} \leq h_i^r c_r (r+1) |u_h|_{H^r(K_i)}. $$

Thus, using an inverse estimate [22, Chapter 4, Section 5, Lemma 4.5.3], we have that there exists $c_r > 0$, a scalar depending only on $r$, such that

$$ |xu_h - \chi_h(u_h)|_{H^1(K_i)} \leq c_r \|u_h\|_{L^2(K_i)} \quad \Rightarrow \quad d_h(u_h, u_h) \leq c_r^2 \|u_h\|_{L^2(0,1)}^2. $$

The estimate for $d_{h,*}(v_h, v_h)$ can be obtained similarly. □

### 4.3.2 Jump and averages at intervals end-points

Functions in the space $V_h^r$ are discontinuous at each end-points of the intervals $(x_i, x_{i+1})$. Therefore, we introduce the following notations: for all $i \in \{1, N-1\}$, and for any function $v_h \in V_h^r$,

$$ \|v_h(x_i)\| := v_h(x_i^+) - v_h(x_i^-), \quad \|v_h(x_i)\| := \frac{v_h(x_i^+) + v_h(x_i^-)}{2} \quad (35) $$

where,

$$ v_h(x_i^+) = \lim_{\epsilon \to 0} v_h(x_i + \epsilon) \quad \text{and} \quad v_h(x_i^-) = \lim_{\epsilon \to 0} v_h(x_i - \epsilon). $$

The following preliminary estimate is required.

**Lemma 4.6.** There exists $c_a > 0$, such that, for all $h > 0$,

$$ \forall v_h \in V_h^r, \quad \sum_{i=1}^{N-1} \tilde{h}_i \|\{v_h(x_i)\}\|^2 \leq c_a \|v_h\|_{L^2(0,1)}^2 \quad \text{where} \quad \tilde{h}_i = \min(h_{i-1}, h_i) $$

*Proof.* Using the definition of the operator $\|\cdot\|$ and $(a+b)^2 \leq 2(a^2 + b^2)$, we have

$$ \sum_{i=1}^{N-1} \tilde{h}_i \|\{v_h(x_i)\}\|^2 \leq \frac{1}{2} \sum_{i=1}^{N-1} \tilde{h}_i (v_h(x_i^+)^2 + v_h(x_i^-)^2). \quad (36) $$

On any interval $(0, L)$, a trace inequality – [22, Chapter 1, Section 6] – holds. In a one-dimensionnal setting it reads,

$$ \forall v \in H^1(0, L), \quad |v(0)|^2 \leq \|v\|_{L^2(0,L)} \left( \frac{1}{L} \|v\|_{L^2(0,L)} + 2 \left\| \frac{d}{dx} v \right\|_{L^2(0,L)} \right). $$

Using this estimate for a fonction $v_h \in \mathbb{P}_r$ as well as an inverse estimate [22, Chapter 4, Section 5, Lemma 4.5.3], we get, with $c_r > 0$ depending only on $r$,

$$ \forall v_h \in \mathbb{P}_r, \quad |v_h^2(0)|^2 \leq \left( \frac{1 + 2c_r}{L} \right) \|v_h\|_{L^2(0,L)}^2. $$
One can then deduce that for all $i \in \{1, \dots, N-1\}$,

$$|v_h(x_i^+)|^2 \le \left( \frac{1+2c_r}{h_i} \right) \|v_h\|_{L^2(x_i, x_{i+1})}^2 \quad \text{and} \quad |v_h(x_i^-)|^2 \le \left( \frac{1+2c_r}{h_{i-1}} \right) \|v_h\|_{L^2(x_{i-1}, x_i)}^2.$$

Using the two equations above in (36) and using the definition of $\tilde{h}_i$, we get

$$\sum_{i=1}^{N-1} \tilde{h}_i \{\{v(x_i)\}\}^2 \le \frac{1+2c_r}{2} \left( \sum_{i=0}^{N-1} \frac{\tilde{h}_i}{h_i} \|v_h\|_{L^2(x_i, x_{i+1})}^2 + \sum_{i=1}^{N} \frac{\tilde{h}_i}{h_{i-1}} \|v_h\|_{L^2(x_{i-1}, x_i)}^2 \right).$$

The proof is concluded by observing that $\tilde{h}_i/h_i$ and $\tilde{h}_i/h_{i-1}$ are less or equal one. □

### 4.3.3 Main result

The next theorem shows that the multiplier $\Pi_h$ satisfies Assumption 1.

**Theorem 4.7.** $\exists \alpha, \beta > 0$, such that, for all $h$ and all $z_h = (u_h, v_h) \in \mathcal{Z}_h = U_h^r \times V_h^r,$

$$(\Pi_h A_h z_h, z_h)_h \le \alpha \left( d_h(u_h, u_h) + d_{h,*}(v_h, v_h) + v_h^2(1) + \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \left[ \|v_h(x_i)\|^2 \right] \right) - \beta \|z_h\|^2_{\mathcal{Z}}.$$

*Proof.* For all $z_h = (u_h, v_h)^t \in U_h^r \times V_h^r$ we have $(\Pi_h A_h z_h, z_h)_h = (A_h z_h, \Pi_h z_h)_h$ and

$$(A_h z_h, \Pi_h z_h)_h = v_h^2(1) - u_h^2(1) - \int_0^1 v_h \frac{d}{dx} \chi_h^*(v_h) \, dx + \int_0^1 \chi_h(u_h) \frac{d}{dx} u_h \, dx.$$

Inserting $xv_h - xv_h = 0$ in the last terms of this equality, we get

$$(A_h z_h, \Pi_h z_h)_h = v_h^2(1) - u_h^2(1) \\ + \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} v_h \frac{d}{dx} (xv_h - \chi_h^*(v_h)) \, dx - \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} v_h \frac{d}{dx} (xv_h) \, dx \\ - \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} (xu_h - \chi_h(u_h)) \frac{d}{dx} u_h \, dx + \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} xu_h \frac{d}{dx} u_h \, dx. \quad (37)$$

By distributing the derivative with respect to $x$, and then integrating by parts, and finally using the algebraic identity $a^2 - b^2 = (a-b)(a+b)$, we obtain

$$- \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} v_h \frac{d}{dx} (xv_h) \, dx = - \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} v_h^2 \, dx - \frac{1}{2} \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} x \frac{d}{dx} (v_h^2) \, dx \\ = - \frac{1}{2} \int_{0}^{1} v_h^2 \, dx - \frac{1}{2} v_h^2(1) - \sum_{i=1}^{N-1} x_i \{\{v_h(x_i)\}\} \|v_h(x_i)\|. \quad (38)$$

Moreover, using integration by parts and the property that $xu_h - \chi_h(u_h)$ vanishes at $x=0$ and $x=1$ (property of the operator $\chi_h$ mentioned Section 4.3.1), we get

$$- \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} (xu_h - \chi_h(u_h)) \frac{d}{dx} u_h \, dx + \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} xu_h \frac{d}{dx} u_h \, dx \\ = \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} u_h \frac{d}{dx} (xu_h - \chi_h(u_h)) \, dx - \frac{1}{2} \int_{0}^{1} u_h^2 \, dx - \frac{1}{2} u_h^2(1). \quad (39)$$
One can then deduce that for all $i \in \{1, \dots, N-1\}$,

$$|v_h(x_i^+)|^2 \le \left( \frac{1+2c_r}{h_i} \right) \|v_h\|_{L^2(x_i, x_{i+1})}^2 \quad \text{and} \quad |v_h(x_i^-)|^2 \le \left( \frac{1+2c_r}{h_{i-1}} \right) \|v_h\|_{L^2(x_{i-1}, x_i)}^2.$$

Using the two equations above in (36) and using the definition of $\tilde{h}_i$, we get

$$\sum_{i=1}^{N-1} \tilde{h}_i \{\{v(x_i)\}\}^2 \le \frac{1+2c_r}{2} \left( \sum_{i=0}^{N-1} \frac{\tilde{h}_i}{h_i} \|v_h\|_{L^2(x_i, x_{i+1})}^2 + \sum_{i=1}^{N} \frac{\tilde{h}_i}{h_{i-1}} \|v_h\|_{L^2(x_{i-1}, x_i)}^2 \right).$$

The proof is concluded by observing that $\tilde{h}_i/h_i$ and $\tilde{h}_i/h_{i-1}$ are less or equal one. $\square$

### 4.3.3 Main result

The next theorem shows that the multiplier $\Pi_h$ satisfies Assumption 1.

**Theorem 4.7.** $\exists \alpha, \beta > 0$, such that, for all $h$ and all $z_h = (u_h, v_h) \in \mathcal{Z}_h = U_h^r \times V_h^r,$

$$(\Pi_h A_h z_h, z_h)_h \le \alpha \left( d_h(u_h, u_h) + d_{h,*}(v_h, v_h) + v_h^2(1) + \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \left[ \|v_h(x_i)\|^2 \right] \right) - \beta \|z_h\|^2_{\mathcal{Z}}.$$

*Proof.* For all $z_h = (u_h, v_h)^t \in U_h^r \times V_h^r$ we have $(\Pi_h A_h z_h, z_h)_h = (A_h z_h, \Pi_h z_h)_h$ and

$$(A_h z_h, \Pi_h z_h)_h = v_h^2(1) - u_h^2(1) - \int_0^1 v_h \frac{d}{dx} \chi_h^*(v_h) \, dx + \int_0^1 \chi_h(u_h) \frac{d}{dx} u_h \, dx.$$

Inserting $xv_h - xv_h = 0$ in the last terms of this equality, we get

$$(A_h z_h, \Pi_h z_h)_h = v_h^2(1) - u_h^2(1) \\ + \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} v_h \frac{d}{dx} (xv_h - \chi_h^*(v_h)) \, dx - \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} v_h \frac{d}{dx} (xv_h) \, dx \\ - \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} (xu_h - \chi_h(u_h)) \frac{d}{dx} u_h \, dx + \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} xu_h \frac{d}{dx} u_h \, dx. \quad (37)$$

By distributing the derivative with respect to $x$, and then integrating by parts, and finally using the algebraic identity $a^2 - b^2 = (a-b)(a+b)$, we obtain

$$- \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} v_h \frac{d}{dx} (xv_h) \, dx = - \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} v_h^2 \, dx - \frac{1}{2} \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} x \frac{d}{dx} (v_h^2) \, dx \\ = - \frac{1}{2} \int_{0}^{1} v_h^2 \, dx - \frac{1}{2} v_h^2(1) - \sum_{i=1}^{N-1} x_i \{\{v_h(x_i)\}\} \|v_h(x_i)\|. \quad (38)$$

Moreover, using integration by parts and the property that $xu_h - \chi_h(u_h)$ vanishes at $x = 0$ and $x = 1$ (property of the operator $\chi_h$ mentioned Section 4.3.1), we get

$$- \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} (xu_h - \chi_h(u_h)) \frac{d}{dx} u_h \, dx + \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} xu_h \frac{d}{dx} u_h \, dx \\ = \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} u_h \frac{d}{dx} (xu_h - \chi_h(u_h)) \, dx - \frac{1}{2} \int_{0}^{1} u_h^2 \, dx - \frac{1}{2} u_h^2(1). \quad (39)$$
Injecting (38) and (39) into (37) leads to the following equality

$$ (A_h z_h, \Pi_h z_h)_h = \frac{1}{2} v_h^2(1) - \frac{1}{2} u_h^2(1) - \frac{1}{2} \int_0^1 (u_h^2 + v_h^2) dx - \sum_{i=1}^{N-1} x_i \{\{v_h(x_i)\}\} \|v_h(x_i)\| + \sum_{i=0}^N \int_{x_i}^{x_{i+1}} v_h \frac{d}{dx} (xv_h - \chi_h^*(v_h)) dx + \sum_{i=0}^N \int_{x_i}^{x_{i+1}} u_h \frac{d}{dx} (xu_h - \chi_h(u_h)) dx. \quad (40) $$

Using Cauchy-Schwarz and Young inequalities, we obtain for all $\lambda > 0$,

$$ (A_h z_h, \Pi_h z_h)_h \le \frac{1}{2} v_h^2(1) - \frac{1}{2} (1 - \lambda) \int_0^1 (u_h^2 + v_h^2) dx - \sum_{i=1}^{N-1} x_i \{\{v_h(x_i)\}\} \|v_h(x_i)\| + \frac{1}{2\lambda} \sum_{i=0}^N \int_{x_i}^{x_{i+1}} \left[ \frac{d}{dx} (xv_h - \chi_h^*(v_h)) \right]^2 + \left[ \frac{d}{dx} (xu_h - \chi_h(u_h)) \right]^2 dx. \quad (41) $$

Let us focus on the term involving jumps and averages. Using Cauchy-Schwarz and Young inequalities, bounding $x_i$ by 1 as well as Lemma 4.6, we get, for all $\mu > 0$,

$$ - \sum_{i=1}^{N-1} x_i \{\{v_h(x_i)\}\} \|v_h(x_i)\| \le \frac{1}{2\mu} \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 + \frac{\mu}{2} \sum_{i=1}^{N-1} \tilde{h}_i \|v_h(x_i)\|^2 \le \frac{1}{2\mu} \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 + \frac{c_a \mu}{2} \|v_h\|_{L^2(0,1)}^2. $$

Injecting this estimate in (41) and using the definition (32) and (33) of the bilinear forms $d_h$ and $d_{h,*}$ as well as (22), we obtain

$$ (A_h z_h, \Pi_h z_h)_h \le \frac{1}{2} v_h^2(1) - (1 - \lambda) \frac{1}{2c_{q,\max}} \int_0^1 u_h^2 dx - \frac{1}{2c_{q,\max}} (1 - \lambda - \mu c_a) \int_0^1 v_h^2 dx + \frac{1}{2\mu} \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 + \frac{1}{2\lambda} d_{h,*}(v_h, v_h) + \frac{1}{2\lambda} d_h(u_h, u_h). $$

We conclude the proof by choosing $\lambda$ and $\mu$ small enough such that $1 - \lambda$ and $1 - \lambda - \mu c_a$ are strictly positive. $\square$

To formally verify that Assumption 1 holds, the estimation of Theorem 4.7 must be written in the form

$$ (\Pi_h A_h z_h, z_h)_h \le \alpha \|C_h z_h\|_h^2 - \beta \|z_h\|_h^2, \quad (42) $$

for some operator $C_h \in \mathcal{L}(\mathcal{Z}_h, \mathcal{Y}_h)$. In order to define this operator we introduce the operators $D_h \in \mathcal{L}(U_h^r)$ and $D_{h,*} \in \mathcal{L}(V_h^r)$ associated with the bilinear forms $d_h$ and $d_{h,*}$,

$$ \left\{ \begin{array}{ll} \forall (u_h, \tilde{u}_h) \in U_h^r \times U_h^r, & (D_h u_h, \tilde{u}_h)_h = d_h(u_h, \tilde{u}_h), \\ \forall (v_h, \tilde{v}_h) \in V_h^r \times V_h^r, & (D_{h,*} v_h, \tilde{v}_h)_h = d_{h,*}(v_h, \tilde{v}_h). \end{array} \right. \quad (43) $$

We also need to introduce the jump operator involving jumps of discontinuous functions and defined by

$$ \forall (v_h, \tilde{v}_h) \in V_h^r \times V_h^r, \quad (J_h v_h, \tilde{v}_h)_h = \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\| \|\tilde{v}_h(x_i)\|. \quad (44) $$

Note that the operators $J_h, B_h, D_h$ and $D_{h,*}$ are self-adjoint and positive. Therefore, we can introduce their square roots. Setting $\mathcal{Y}_h = U_h^r \times [V_h^r]^3$ (a natural definition
Injecting (38) and (39) into (37) leads to the following equality

$$ (A_h z_h, \Pi_h z_h)_h = \frac{1}{2} v_h^2(1) - \frac{1}{2} u_h^2(1) - \frac{1}{2} \int_0^1 (u_h^2 + v_h^2) \, dx - \sum_{i=1}^{N-1} x_i \{\{v_h(x_i)\}\} \|v_h(x_i)\| \\ + \sum_{i=0}^N \int_{x_i}^{x_{i+1}} v_h \frac{d}{dx} (xv_h - \chi_h^*(v_h)) \, dx + \sum_{i=0}^N \int_{x_i}^{x_{i+1}} u_h \frac{d}{dx} (xu_h - \chi_h(u_h)) \, dx. \quad (40) $$

Using Cauchy-Schwarz and Young inequalities, we obtain for all $\lambda > 0$,

$$ (A_h z_h, \Pi_h z_h)_h \le \frac{1}{2} v_h^2(1) - \frac{1}{2} (1 - \lambda) \int_0^1 (u_h^2 + v_h^2) \, dx - \sum_{i=1}^{N-1} x_i \{\{v_h(x_i)\}\} \|v_h(x_i)\| \\ + \frac{1}{2\lambda} \sum_{i=0}^N \int_{x_i}^{x_{i+1}} \left[ \frac{d}{dx} (xv_h - \chi_h^*(v_h)) \right]^2 + \left[ \frac{d}{dx} (xu_h - \chi_h(u_h)) \right]^2 \, dx. \quad (41) $$

Let us focus on the term involving jumps and averages. Using Cauchy-Schwarz and Young inequalities, bounding $x_i$ by 1 as well as Lemma 4.6, we get, for all $\mu > 0$,

$$ - \sum_{i=1}^{N-1} x_i \{\{v_h(x_i)\}\} \|v_h(x_i)\| \le \frac{1}{2\mu} \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 + \frac{\mu}{2} \sum_{i=1}^{N-1} \tilde{h}_i \|v_h(x_i)\|^2 \\ \le \frac{1}{2\mu} \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 + \frac{c_a \mu}{2} \|v_h\|_{L^2(0,1)}^2. $$

Injecting this estimate in (41) and using the definition (32) and (33) of the bilinear forms $d_h$ and $d_{h,*}$ as well as (22), we obtain

$$ (A_h z_h, \Pi_h z_h)_h \le \frac{1}{2} v_h^2(1) - (1 - \lambda) \frac{1}{2c_{q,\max}} \int_0^1 u_h^2 dx - \frac{1}{2c_{q,\max}} (1 - \lambda - \mu c_a) \int_0^1 v_h^2 dx \\ + \frac{1}{2\mu} \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 + \frac{1}{2\lambda} d_{h,*}(v_h, v_h) + \frac{1}{2\lambda} d_h(u_h, u_h). $$

We conclude the proof by choosing $\lambda$ and $\mu$ small enough such that $1 - \lambda$ and $1 - \lambda - \mu c_a$ are strictly positive. $\square$

To formally verify that Assumption 1 holds, the estimation of Theorem 4.7 must be written in the form

$$ (\Pi_h A_h z_h, z_h)_h \le \alpha \|C_h z_h\|_h^2 - \beta \|z_h\|_h^2, \quad (42) $$

for some operator $C_h \in \mathcal{L}(\mathcal{Z}_h, \mathcal{Y}_h)$. In order to define this operator we introduce the operators $D_h \in \mathcal{L}(U_h^r)$ and $D_{h,*} \in \mathcal{L}(V_h^r)$ associated with the bilinear forms $d_h$ and $d_{h,*}$,

$$ \left\{ \begin{array}{ll} \forall (u_h, \tilde{u}_h) \in U_h^r \times U_h^r, & (D_h u_h, \tilde{u}_h)_h = d_h(u_h, \tilde{u}_h), \\ \forall (v_h, \tilde{v}_h) \in V_h^r \times V_h^r, & (D_{h,*} v_h, \tilde{v}_h)_h = d_{h,*}(v_h, \tilde{v}_h). \end{array} \right. \quad (43) $$

We also need to introduce the jump operator involving jumps of discontinuous functions and defined by

$$ \forall (v_h, \tilde{v}_h) \in V_h^r \times V_h^r, \quad (J_h v_h, \tilde{v}_h)_h = \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\| \|\tilde{v}_h(x_i)\|. \quad (44) $$

Note that the operators $J_h, B_h, D_h$ and $D_{h,*}$ are self-adjoint and positive. Therefore, we can introduce their square roots. Setting $\mathcal{Y}_h = U_h^r \times [V_h^r]^3$ (a natural definition
of its scalar product ensues from its tensorial structure), we define the operator $C_h \in \mathcal{L}(\mathcal{Z}_h, \mathcal{Y}_h)$ by

$$C_h = \begin{pmatrix} D_h^{\frac{1}{2}} & 0 \\ 0 & J_h^{\frac{1}{2}} \\ 0 & B_h^{\frac{1}{2}} \\ 0 & D_{h,*}^{\frac{1}{2}} \end{pmatrix}. \qquad (45)$$

From all these considerations we deduce the following corollary.

**Corollary 4.8.** The operator $\Pi_h$ defined by (30) satisfies Assumption 1 with $C_h$ defined by (45). As a consequence, for $T$ large enough, there exists $c_{st}(T)$ independent of $h$, such that for any $z_h(t) = (u_h(t), v_h(t))^t \in \mathcal{Z}_h$ solution to (27),

$$\|u_{h,0}\|_{h}^{2} + \|v_{h,0}\|_{h}^{2} \le c_{st}(T) \left[ \int_{0}^{T} d_{h}(u_{h}, u_{h}) \, dt + \int_{0}^{T} d_{h,*}(v_{h}, v_{h}) \, dt \right. \\ \left. + \int_{0}^{T} v_{h}^{2}(1) \, dt + \int_{0}^{T} \sum_{i=1}^{N-1} \frac{1}{\bar{h}_{i}} \left[ \|v_{h}(x_{i})\| \right]^{2} \, dt \right] = c_{st}(T) \int_{0}^{T} \|C_{h}z_{h}\|_{h}^{2} \, dt. \quad (46)$$

*Proof.* Observe that Theorem 4.7 indeed implies (42), then Theorem 4.3 can be applied to conclude. □

Following what was done in Section 3, (46) is a guideline to build an exponentially stable system. Using $C_h$ – as defined by (45) – to damp the dynamics (27), we construct an uniformly in $h$ exponentially stable dynamics.

## 4.4 Exponential stability

We incorporate in the system (27) the dissipative terms $C_h$ in the observability inequality (46). In other words, we introduce a perturbation of system (27) by the operator $\gamma C_h^* C_h$, which gives us the following system

$$\begin{cases} \dot{z}_h = A_h z_h - \gamma C_h^* C_h z_h, & \text{in } [0, T], \\ z_h(0) = z_{h,0}, \end{cases} \qquad (47)$$

where $A_h$ and $C_h$ are respectively defined by (27) and (45). We will be able to expose our main result, namely the uniform exponential stability property in the semi-discrete case. First, it is verified that the multiplier $\Pi_h$ satisfies the fundamental Assumption 2.

**Theorem 4.9.** Let $C_h$ be defined by (45). Then $\Pi_h$ defined by (30) satisfies Assumption 2, in particular, $\exists \chi, \kappa > 0$, such that, for all $h$,

$$\forall z_h \in \mathcal{Z}_h, \quad (\Pi_h A_h z_h, z_h)_h - \gamma (C_h z_h, C_h \Pi_h z_h)_h \le \chi \|C_h z_h\|_h^2 - \kappa \|z_h\|_h^2.$$

*Proof.* Since $(\Pi_h A_h z_h, z_h)_h$ has been studied in the proof of Theorem 4.7, we focus on the term $(C_h z_h, C_h \Pi_h z_h)_h$. We have

$$-(C_h z_h, C_h \Pi_h z_h)_h = - \sum_{i=1}^{N-1} \frac{1}{\bar{h}_i} \left[ \|v_h(x_i)\| \left[ \chi_h(u_h)(x_i) \right] + v_h(1)\chi_h(u_h)(1) \right. \\ \left. + d_h(u_h, \chi_h^*(v_h)) + d_{h,*}(v_h, \chi_h(u_h)). \right]$$

Thanks to the properties of the interpolation operator $\chi_h$ discussed in Section 4.3.1, we have $\chi_h(u_h)(1) = u_h(1)$ and $\left[ \chi_h(u_h)(x_i) \right] = 0$. Therefore, we obtain

$$-(C_h z_h, C_h \Pi_h z_h)_h = v_h(1)u_h(1) + d_h(u_h, \chi_h^*(v_h)) + d_{h,*}(v_h, \chi_h(u_h)).$$
of its scalar product ensues from its tensorial structure), we define the operator $C_h \in \mathcal{L}(\mathcal{Z}_h, \mathcal{Y}_h)$ by

$$C_h = \begin{pmatrix} D_h^{\frac{1}{2}} & 0 \\ 0 & J_h^{\frac{1}{2}} \\ 0 & B_h^{\frac{1}{2}} \\ 0 & D_{h,*}^{\frac{1}{2}} \end{pmatrix}. \qquad (45)$$

From all these considerations we deduce the following corollary.

**Corollary 4.8.** The operator $\Pi_h$ defined by (30) satisfies Assumption 1 with $C_h$ defined by (45). As a consequence, for $T$ large enough, there exists $c_{st}(T)$ independent of $h$, such that for any $z_h(t) = (u_h(t), v_h(t))^t \in \mathcal{Z}_h$ solution to (27),

$$\|u_{h,0}\|_{h}^{2} + \|v_{h,0}\|_{h}^{2} \le c_{st}(T) \left[ \int_{0}^{T} d_{h}(u_{h}, u_{h}) dt + \int_{0}^{T} d_{h,*}(v_{h}, v_{h}) dt \right. \\ \left. + \int_{0}^{T} v_{h}^{2}(1) dt + \int_{0}^{T} \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_{i}} \left[ \|v_{h}(x_{i})\|^{2} \right] dt \right] = c_{st}(T) \int_{0}^{T} \|C_{h}z_{h}\|_{h}^{2} dt. \quad (46)$$

*Proof.* Observe that Theorem 4.7 indeed implies (42), then Theorem 4.3 can be applied to conclude. □

Following what was done in Section 3, (46) is a guideline to build an exponentially stable system. Using $C_h$ – as defined by (45) – to damp the dynamics (27), we construct an uniformly in $h$ exponentially stable dynamics.

## 4.4 Exponential stability

We incorporate in the system (27) the dissipative terms $C_h$ in the observability inequality (46). In other words, we introduce a perturbation of system (27) by the operator $\gamma C_h^* C_h$, which gives us the following system

$$\begin{cases} \dot{z}_h = A_h z_h - \gamma C_h^* C_h z_h, & \text{in } [0, T], \\ z_h(0) = z_{h,0}, & \end{cases} \qquad (47)$$

where $A_h$ and $C_h$ are respectively defined by (27) and (45). We will be able to expose our main result, namely the uniform exponential stability property in the semi-discrete case. First, it is verified that the multiplier $\Pi_h$ satisfies the fundamental Assumption 2.

**Theorem 4.9.** Let $C_h$ be defined by (45). Then $\Pi_h$ defined by (30) satisfies Assumption 2, in particular, $\exists \chi, \kappa > 0$, such that, for all $h$,

$$\forall z_h \in \mathcal{Z}_h, \quad (\Pi_h A_h z_h, z_h)_h - \gamma (C_h z_h, C_h \Pi_h z_h)_h \le \chi \|C_h z_h\|_h^2 - \kappa \|z_h\|_h^2.$$

*Proof.* Since $(\Pi_h A_h z_h, z_h)_h$ has been studied in the proof of Theorem 4.7, we focus on the term $(C_h z_h, C_h \Pi_h z_h)_h$. We have

$$-(C_h z_h, C_h \Pi_h z_h)_h = - \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \left[ \|v_h(x_i)\| \left[ \chi_h(u_h)(x_i) \right] + v_h(1) \chi_h(u_h)(1) \right. \\ \left. + d_h(u_h, \chi_h^*(v_h)) + d_{h,*}(v_h, \chi_h(u_h)). \right]$$

Thanks to the properties of the interpolation operator $\chi_h$ discussed in Section 4.3.1, we have $\chi_h(u_h)(1) = u_h(1)$ and $\left[ \chi_h(u_h)(x_i) \right] = 0$. Therefore, we obtain

$$-(C_h z_h, C_h \Pi_h z_h)_h = v_h(1)u_h(1) + d_h(u_h, \chi_h^*(v_h)) + d_{h,*}(v_h, \chi_h(u_h)).$$
Gathering the above equation and (40), we get

$$ (A_h z_h, \Pi_h z_h)_h - \gamma(C_h z_h, C_h \Pi_h z_h)_h = \frac{1}{2} v_h^2(1) - \frac{1}{2} u_h^2(1) + \gamma v_h(1) u_h(1) \\ - \frac{1}{2} \int_0^1 (u_h^2 + v_h^2) dx - \sum_{i=1}^{N-1} x_i \{\{v_h(x_i)\}\} \|v_h(x_i)\| + \gamma d_h(u_h, \chi_h^*(v_h)) + \gamma d_{h,*}(v_h, \chi_h(u_h)) \\ + \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} v_h \frac{d}{dx}(xv_h - \chi_h^*(v_h)) dx + \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} u_h \frac{d}{dx}(xu_h - \chi_h(u_h)) dx. \quad (48) $$

Following the proof of Theorem 4.7, one can deduce from (48), the estimate

$$ (\Pi_h A_h z_h, z_h)_h - \gamma(C_h z_h, C_h \Pi_h z_h)_h \\ \leq -\frac{1}{2c_{q,\max}} (1 - \lambda - \mu c_a) \|z_h\|_h^2 + \frac{1}{2\lambda} d_h(u_h, u_h) + \frac{1}{2\lambda} d_{h,*}(v_h, v_h) \\ + \frac{\gamma^2 + 1}{2} v_h^2(1) + \frac{1}{2\mu} \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 + \gamma d_h(u_h, \chi_h^*(v_h)) + \gamma d_{h,*}(v_h, \chi_h(u_h)), \quad (49) $$

where $\lambda$ and $\mu$ are sufficiently small positive scalars. Thanks to Cauchy-Schwarz and Young inequalities, Proposition 4.5 and (31), one can show that, for any $\delta > 0$,

$$ d_h(u_h, \chi_h^*(v_h)) \leq \frac{1}{2\delta} d_h(u_h, u_h) + \frac{\delta c}{2} \|v_h\|_h^2, \quad (50) $$

with $c = c_d c_\pi^2 / c_{q,\min}$. Similarly, one gets,

$$ d_{h,*}(v_h, \chi_h(u_h)) \leq \frac{1}{2\delta} d_{h,*}(v_h, v_h) + \frac{\delta c}{2} \|u_h\|_h^2. \quad (51) $$

Injecting (50) and (51) in (49) yields

$$ (\Pi_h A_h z_h, z_h)_h - \gamma(C_h z_h, C_h \Pi_h z_h)_h \\ \leq -\chi \|z_h\|_h^2 + \frac{\gamma^2 + 1}{2} v_h^2(1) + \frac{1}{2\mu} \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 \\ + \left( \frac{1}{2\lambda} + \frac{\gamma}{2\delta} \right) (d_h(u_h, u_h) + d_{h,*}(v_h, v_h)), $$

for some scalar $\kappa$, strictly positive for $\lambda, \mu, \delta$ small enough. The proof can be concluded by observing that

$$ \|C_h z_h\|_h^2 = v_h^2(1) + \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 + d_h(u_h, u_h) + d_{h,*}(v_h, v_h). $$

□

We can deduce from Theorem 4.9 the uniform exponential stability property for the damped wave equation (47).

**Corollary 4.10.** There exist $c_s > 0$ and $\sigma > 0$ independent of $h$ such that, for all $T > 0$ and any $z_h(t) = (u_h(t), v_h(t))^t$ solution to (47),

$$ (\|u_h(t)\|_h^2 + \|v_h(t)\|_h^2)^{\frac{1}{2}} \leq c_s e^{-\sigma t} (\|u_{h,0}\|_h^2 + \|v_{h,0}\|_h^2)^{\frac{1}{2}}. $$

*Proof.* Theorem 4.9 shows that the multiplier $\Pi_h$ defined by (30) satisfies Assumption 2. We can then apply Theorem 4.4 with $C_h$ and is defined by (45). □
Gathering the above equation and (40), we get

$$ (A_h z_h, \Pi_h z_h)_h - \gamma(C_h z_h, C_h \Pi_h z_h)_h = \frac{1}{2} v_h^2(1) - \frac{1}{2} u_h^2(1) + \gamma v_h(1) u_h(1) \\ - \frac{1}{2} \int_0^1 (u_h^2 + v_h^2) dx - \sum_{i=1}^{N-1} x_i \{\{v_h(x_i)\}\} \|v_h(x_i)\| + \gamma d_h(u_h, \chi_h^*(v_h)) + \gamma d_{h,*}(v_h, \chi_h(u_h)) \\ + \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} v_h \frac{d}{dx}(xv_h - \chi_h^*(v_h)) dx + \sum_{i=0}^{N} \int_{x_i}^{x_{i+1}} u_h \frac{d}{dx}(xu_h - \chi_h(u_h)) dx. \quad (48) $$

Following the proof of Theorem 4.7, one can deduce from (48), the estimate

$$ (\Pi_h A_h z_h, z_h)_h - \gamma(C_h z_h, C_h \Pi_h z_h)_h \\ \leq -\frac{1}{2c_{q,\max}} (1 - \lambda - \mu c_a) \|z_h\|_h^2 + \frac{1}{2\lambda} d_h(u_h, u_h) + \frac{1}{2\lambda} d_{h,*}(v_h, v_h) \\ + \frac{\gamma^2 + 1}{2} v_h^2(1) + \frac{1}{2\mu} \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 + \gamma d_h(u_h, \chi_h^*(v_h)) + \gamma d_{h,*}(v_h, \chi_h(u_h)), \quad (49) $$

where $\lambda$ and $\mu$ are sufficiently small positive scalars. Thanks to Cauchy-Schwarz and Young inequalities, Proposition 4.5 and (31), one can show that, for any $\delta > 0$,

$$ d_h(u_h, \chi_h^*(v_h)) \leq \frac{1}{2\delta} d_h(u_h, u_h) + \frac{\delta c}{2} \|v_h\|_h^2, \quad (50) $$

with $c = c_d c_\pi^2 / c_{q,\min}$. Similarly, one gets,

$$ d_{h,*}(v_h, \chi_h(u_h)) \leq \frac{1}{2\delta} d_{h,*}(v_h, v_h) + \frac{\delta c}{2} \|u_h\|_h^2. \quad (51) $$

Injecting (50) and (51) in (49) yields

$$ (\Pi_h A_h z_h, z_h)_h - \gamma(C_h z_h, C_h \Pi_h z_h)_h \\ \leq -\chi \|z_h\|_h^2 + \frac{\gamma^2 + 1}{2} v_h^2(1) + \frac{1}{2\mu} \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 \\ + \left( \frac{1}{2\lambda} + \frac{\gamma}{2\delta} \right) (d_h(u_h, u_h) + d_{h,*}(v_h, v_h)), $$

for some scalar $\kappa$, strictly positive for $\lambda, \mu, \delta$ small enough. The proof can be concluded by observing that

$$ \|C_h z_h\|_h^2 = v_h^2(1) + \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 + d_h(u_h, u_h) + d_{h,*}(v_h, v_h). $$

□

We can deduce from Theorem 4.9 the uniform exponential stability property for the damped wave equation (47).

**Corollary 4.10.** There exist $c_s > 0$ and $\sigma > 0$ independent of $h$ such that, for all $T > 0$ and any $z_h(t) = (u_h(t), v_h(t))^t$ solution to (47),

$$ (\|u_h(t)\|_h^2 + \|v_h(t)\|_h^2)^{\frac{1}{2}} \leq c_s e^{-\sigma t} (\|u_{h,0}\|_h^2 + \|v_{h,0}\|_h^2)^{\frac{1}{2}}. $$

*Proof.* Theorem 4.9 shows that the multiplier $\Pi_h$ defined by (30) satisfies Assumption 2. We can then apply Theorem 4.4 with $C_h$ and is defined by (45). □
## 4.5 Additional properties of the stabilization operators

In this section we explain some properties of the operators $D_h$ and $D_{h,*}$ and the corresponding bilinear forms $d_h$ and $d_{h,*}$. In particular, we give another expression for $d_h$ that allows the added stabilization operator $D_h$ to be interpreted as an elementwise $r$-Laplacian operator. The given expression is also more convenient to handle for implementation purposes, but we do not elaborate on this last property. Finally, we point out that in [21] a Laplacian operator is used as a viscous term for the stabilization of finite elements of order 1. Therefore, in a sense, our work generalizes these results by showing that an element-wise $r$-Laplacian operator – scaled by the corresponding power of $h$ – can be used for the stabilization of high-order finite elements

**Proposition 4.11.** For any $u_h \in U_h^r$,

$$d_h(u_h, u_h) = c_r \sum_{i=0}^{N-1} h_i^{2r+1} |u_h^{(r)}(x_{i+\frac{1}{2}})|^2, \quad \text{where} \quad x_{i+\frac{1}{2}} := \frac{x_i + x_{i+1}}{2}.$$

where $c_r$ is a positive scalar depending only on the finite element order $r$.

*Proof.* To show the stated result, it is sufficient to study the following term for $u_h \in \mathbb{P}_r$,

$$\frac{d}{dx}(xu_h - \chi_h(u_h)).$$

Since $u_h \in U_h^r$, we have that $u_i := u_h|_{(x_i, x_{i+1})} \in \mathbb{P}_r$, hence

$$u_i(x) = \sum_{k=0}^{r-1} \frac{x^k}{k!} u_i^{(k)}(x_i) + r_i(x) \quad \text{where} \quad r_i(x) := \frac{x^r}{r!} u_i^{(r)}(x_{i+\frac{1}{2}}).$$

Note that, since $u_i^{(r)}$ is a constant function in space, we simply write $u_i^{(r)}(x_{i+\frac{1}{2}}) = u_i^{(r)}$. Moreover, observe that

$$\chi_h(u_h)|_{(x_i, x_{i+1})} = I_i(xu_h),$$

where $I_i \in \mathcal{L}(C^0([x_i, x_{i+1}]), \mathbb{P}_r)$ is the interpolation operator on the Gauss-Lobatto points $\{x_i + h_i \xi_k\}_{k=0}^r$. The operator $I_i$ is a projection operator on $\mathbb{P}_r$. Therefore, one can deduce that

$$(xu_h - \chi_h(u_h))_{(x_i, x_{i+1})} = xr_i(x) - I_i(xr_i(x)) = \frac{u_i^{(r)}}{r!} (x^{r+1} - I_i(x^{r+1})). \quad (52)$$

and

$$d_h(u_h, u_h) = \sum_{i=0}^{N-1} \frac{(u_i^{(r)})^2}{r!^2} \int_{x_i}^{x_{i+1}} \left| \frac{d}{dx}(x^{r+1} - I_i(x^{r+1})) \right|^2 dx. \quad (53)$$

We now introduce the following change of variable,

$$\phi_i : \begin{cases} [0, 1] \to [x_i, x_{i+1}] \\ \xi \to x = h_i \xi + x_i \end{cases} \Rightarrow (I_{h_i}(x^{r+1}) \circ \phi_i)(\xi) = \sum_{k=0}^{r} \psi_k(\xi) \xi_k^{r+1} h_i^{r+1},$$

where $\{\psi_k(\xi)\}_{k=0}^r$ are a basis for $\mathbb{P}_r$ such that $\psi_k(\xi_\ell) = \delta_{k,\ell}$. The use of the above identity and (52) leads to

$$\left( \frac{d}{dx} (x^{r+1} - I_i(x^{r+1})) \right) \circ \phi_i = \frac{1}{h_i} \frac{d}{d\xi} \left( h_i^{r+1} \xi^{r+1} - \sum_{k=0}^{r} \psi_k(\xi) \xi_k^{r+1} h_i^{r+1} \right).$$
## 4.5 Additional properties of the stabilization operators

In this section we explain some properties of the operators $D_h$ and $D_{h,*}$ and the corresponding bilinear forms $d_h$ and $d_{h,*}$. In particular, we give another expression for $d_h$ that allows the added stabilization operator $D_h$ to be interpreted as an elementwise r-Laplacian operator. The given expression is also more convenient to handle for implementation purposes, but we do not elaborate on this last property. Finally, we point out that in [21] a Laplacian operator is used as a viscous term for the stabilization of finite elements of order 1. Therefore, in a sense, our work generalizes these results by showing that an element-wise r-Laplacian operator – scaled by the corresponding power of h – can be used for the stabilization of high-order finite elements

**Proposition 4.11.** For any $u_h \in U_h^r$,

$$d_h(u_h, u_h) = c_r \sum_{i=0}^{N-1} h_i^{2r+1} |u_h^{(r)}(x_{i+\frac{1}{2}})|^2, \quad \text{where} \quad x_{i+\frac{1}{2}} := \frac{x_i + x_{i+1}}{2}.$$

where $c_r$ is a positive scalar depending only on the finite element order r.

*Proof.* To show the stated result, it is sufficient to study the following term for $u_h \in \mathbb{P}_r$,

$$\frac{d}{dx}(xu_h - \chi_h(u_h)).$$

Since $u_h \in U_h^r$, we have that $u_i := u_h|_{(x_i, x_{i+1})} \in \mathbb{P}_r$, hence

$$u_i(x) = \sum_{k=0}^{r-1} \frac{x^k}{k!} u_i^{(k)}(x_i) + r_i(x) \quad \text{where} \quad r_i(x) := \frac{x^r}{r!} u_i^{(r)}(x_{i+\frac{1}{2}}).$$

Note that, since $u_i^{(r)}$ is a constant function in space, we simply write $u_i^{(r)}(x_{i+\frac{1}{2}}) = u_i^{(r)}$. Moreover, observe that

$$\chi_h(u_h)|_{(x_i, x_{i+1})} = I_i(xu_h),$$

where $I_i \in \mathcal{L}(C^0([x_i, x_{i+1}]), \mathbb{P}_r)$ is the interpolation operator on the Gauss-Lobatto points $\{x_i + h_i \xi_k\}_{k=0}^r$. The operator $I_i$ is a projection operator on $\mathbb{P}_r$. Therefore, one can deduce that

$$(xu_h - \chi_h(u_h))_{(x_i, x_{i+1})} = xr_i(x) - I_i(xr_i(x)) = \frac{u_i^{(r)}}{r!} (x^{r+1} - I_i(x^{r+1})). \quad (52)$$

and

$$d_h(u_h, u_h) = \sum_{i=0}^{N-1} \frac{(u_i^{(r)})^2}{r!^2} \int_{x_i}^{x_{i+1}} \left| \frac{d}{dx}(x^{r+1} - I_i(x^{r+1})) \right|^2 dx. \quad (53)$$

We now introduce the following change of variable,

$$\phi_i : \begin{cases} [0, 1] \to [x_i, x_{i+1}] \\ \xi \to x = h_i \xi + x_i \end{cases} \Rightarrow (I_{h_i}(x^{r+1}) \circ \phi_i)(\xi) = \sum_{k=0}^{r} \psi_k(\xi) \xi_k^{r+1} h_i^{r+1},$$

where $\{\psi_k(\xi)\}_{k=0}^r$ are a basis for $\mathbb{P}_r$ such that $\psi_k(\xi_\ell) = \delta_{k,\ell}$. The use of the above identity and (52) leads to

$$\left( \frac{d}{dx} (x^{r+1} - I_i(x^{r+1})) \right) \circ \phi_i = \frac{1}{h_i} \frac{d}{d\xi} \left( h_i^{r+1} \xi^{r+1} - \sum_{k=0}^{r} \psi_k(\xi) \xi_k^{r+1} h_i^{r+1} \right).$$
Finally, (53) is rewritten

$$d_h(u_h, u_h) = \sum_{i=0}^{N-1} \int_{x_i}^{x_{i+1}} \left| \frac{d}{dx} (xu_h - \chi_h(u_h)) \right|^2 = c_r (u_i^{(r)})^2 \frac{1}{h_i} h_i^{2r+2},$$

$$\text{where } c_r = \frac{1}{r!2} \int_0^1 \left| \frac{d}{d\xi} \left( \xi^{r+1} - \sum_{k=0}^r \psi_k(\xi) \xi_k^{r+1} \right) \right|^2 d\xi. \quad \square$$

The same result cannot be obtained straightforwardly for $d_{h,*}$. In fact, due to the discontinuous nature of the space $V_h^r$, jump terms are involved.

**Proposition 4.12.** For $v_h \in V_h^r$ then

$$d_{h,*}(v_h, v_h) \le c_r^* \left( \sum_{i=0}^{N-1} h_i^{2r+1} \left| v_h^{(r)} \left( x_{i+\frac{1}{2}} \right) \right|^2 + h \sum_{i=1}^{N-1} \left\| v_h(x_i) \right\|^2 \right),$$

where $c_r^*$ is a positive scalar depending only on the finite element order r.

*Proof.* The proof is based on an adequate decomposition of any function $v_h \in V_h^r$. More precisely we introduce $w_h \in V_h^1$ defined by

$$w_h(0) = w_h(1) = 0, \quad w_h(x_i^{\pm}) = \pm \frac{\|v_h(x_i)\|}{2}, \quad i \in \{1, \dots, N-1\} \quad \text{and} \quad w_h|_{(x_i, x_{i+1})} \in \mathbb{P}_1.$$

We then defined $u_h = v_h - w_h$. It can be verified that $u_h \in U_h^r$. Using the Young's inequality, Proposition 4.5 and the property that $u_h \in U_h^r$

$$d_{h,*}(v_h, v_h) \le 2d_h(u_h, u_h) + 2c_d \int_0^1 w_h^2 dx,$$

hence, using Proposition 4.11,

$$d_{h,*}(v_h, v_h) \le 2c_r \sum_{i=0}^{N-1} h_i^{2r+1} \left| u_h^{(r)} \left( x_{i+\frac{1}{2}} \right) \right|^2 + 2c_d \sum_{i=0}^{N-1} \int_{x_i}^{x_{i+1}} w_h^2 dx.$$

It is then possible to show that, there exists $c_j > 0$ independent of $h$, such that

$$\sum_{i=0}^{N-1} \int_{x_i}^{x_{i+1}} w_h^2 dx \le c_j h \sum_{i=1}^{N-1} \|v_h(x_i)\|^2,$$

finally, to conclude, we observe that $u_h^{(r)}(x_{i+\frac{1}{2}}) = (u_h + w_h)^{(r)}(x_{i+\frac{1}{2}}) = v_h^{(r)}(x_{i+\frac{1}{2}})$ by construction. $\square$

**Remark 4.13.** In the numerical results presented below, we used the following simplified expression instead of the term $d_{h,*}$

$$\forall (v_h, \tilde{v}_h) \in V_h^r \times V_h^r, \quad \tilde{d}_h(v_h, \tilde{v}_h) := c_r \sum_{i=0}^{N-1} h_i^{2r+1} v_h^{(r)}(x_{i+\frac{1}{2}}) \tilde{v}_h^{(r)}(x_{i+\frac{1}{2}}), \quad (54)$$

that can be seen as an extension to $V_h^r$ of the bilinear form $d_h$. Numerical evidences show that Theorem 4.9 hold with this substitution.
Finally, (53) is rewritten

$$d_h(u_h, u_h) = \sum_{i=0}^{N-1} \int_{x_i}^{x_{i+1}} \left| \frac{d}{dx} (xu_h - \chi_h(u_h)) \right|^2 = c_r (u_i^{(r)})^2 \frac{1}{h_i} h_i^{2r+2},$$

$$\text{where } c_r = \frac{1}{r!2} \int_0^1 \left| \frac{d}{d\xi} \left( \xi^{r+1} - \sum_{k=0}^r \psi_k(\xi) \xi_k^{r+1} \right) \right|^2 d\xi. \quad \square$$

The same result cannot be obtained straightforwardly for $d_{h,*}$. In fact, due to the discontinuous nature of the space $V_h^r$, jump terms are involved.

**Proposition 4.12.** For $v_h \in V_h^r$ then

$$d_{h,*}(v_h, v_h) \le c_r^* \left( \sum_{i=0}^{N-1} h_i^{2r+1} \left| v_h^{(r)} \left( x_{i+\frac{1}{2}} \right) \right|^2 + h \sum_{i=1}^{N-1} \left| \left[ v_h(x_i) \right] \right|^2 \right),$$

where $c_r^*$ is a positive scalar depending only on the finite element order r.

*Proof.* The proof is based on an adequate decomposition of any function $v_h \in V_h^r$. More precisely we introduce $w_h \in V_h^1$ defined by

$$w_h(0) = w_h(1) = 0, \quad w_h(x_i^{\pm}) = \pm \frac{\left| \left[ v_h(x_i) \right] \right|}{2}, \quad i \in \{1, \dots, N-1\} \quad \text{and} \quad w_h|_{(x_i, x_{i+1})} \in \mathbb{P}_1.$$

We then defined $u_h = v_h - w_h$. It can be verified that $u_h \in U_h^r$. Using the Young's inequality, Proposition 4.5 and the property that $u_h \in U_h^r$

$$d_{h,*}(v_h, v_h) \le 2d_h(u_h, u_h) + 2c_d \int_0^1 w_h^2 dx,$$

hence, using Proposition 4.11,

$$d_{h,*}(v_h, v_h) \le 2c_r \sum_{i=0}^{N-1} h_i^{2r+1} \left| u_h^{(r)} \left( x_{i+\frac{1}{2}} \right) \right|^2 + 2c_d \sum_{i=0}^{N-1} \int_{x_i}^{x_{i+1}} w_h^2 dx.$$

It is then possible to show that, there exists $c_j > 0$ independent of $h$, such that

$$\sum_{i=0}^{N-1} \int_{x_i}^{x_{i+1}} w_h^2 dx \le c_j h \sum_{i=1}^{N-1} \left| \left[ v_h(x_i) \right] \right|^2,$$

finally, to conclude, we observe that $u_h^{(r)}(x_{i+\frac{1}{2}}) = (u_h + w_h)^{(r)}(x_{i+\frac{1}{2}}) = v_h^{(r)}(x_{i+\frac{1}{2}})$ by construction. $\square$

**Remark 4.13.** In the numerical results presented below, we used the following simplified expression instead of the term $d_{h,*}$

$$\forall (v_h, \tilde{v}_h) \in V_h^r \times V_h^r, \quad \tilde{d}_h(v_h, \tilde{v}_h) := c_r \sum_{i=0}^{N-1} h_i^{2r+1} v_h^{(r)}(x_{i+\frac{1}{2}}) \tilde{v}_h^{(r)}(x_{i+\frac{1}{2}}), \quad (54)$$

that can be seen as an extension to $V_h^r$ of the bilinear form $d_h$. Numerical evidences show that Theorem 4.9 hold with this substitution.
# 5 Application to a data assimilation problem

## 5.1 Motivation

One application of the results of exponential stabilization is to construct observers that use available measurements as feedback to improve a physical system approximation. In this case, the goal is to construct a numerical approximation to the solution of a wave propagation problem for which we have boundary measurements. More precisely, we assume that we have access to the boundary data $v(1, t)$ for all times, where $z = (u, v)^t$ is the solution of (1)-(2) with $\gamma = 0$ (hence the problem is conservative). Despite unknown initial data, we want to construct a numerical approximation $z_h = (u_h, v_h)^t$ that converges asymptotically to $z$ in time $t$ up to a discretization error. To do so, one can follow the approach proposed in [7, 30] and compute a so-called observer $\hat{z}_h$ as the solution to

$$ \dot{\hat{z}}_h = A_h \hat{z}_h - \gamma \begin{pmatrix} \tilde{D}_h & 0 \\ 0 & \tilde{D}_{h,*} + J_h \end{pmatrix} \hat{z}_h - \gamma \begin{pmatrix} 0 & 0 \\ 0 & B_h \end{pmatrix} (\hat{z}_h - \bar{z}_h), \quad \hat{z}_h(0) = 0, \quad (55) $$

where $\bar{z}_h$ is the interpolation of the continuous function $z$ solution to (1)-(2). Note that, by assumption, we only have access to $B_h \bar{z}_h(t)$, which corresponds to the measurements at the boundary. These measurements are used in a feedback law for the dynamics of the observer. Note finally that $\bar{z}$ is the solution of a conservative problem (i.e., $\gamma = 0$) and $(u_0, v_0)$ are unknown. The main question to answer concerns the convergence of $\hat{z}_h$ towards $z$.

## 5.2 Convergence result

In this section, we follow the approach of [31] and provide a convergence result of the solution of the exponentially stable discrete problem (55) towards the solution of the continuous problem (1)-(2). The most important feature of the proposed error estimate is that, unlike the standard error estimates, it is uniform in $T$. We introduce the interpolation of $z(t) = (u(t), v(t))^t$ – solution to(1)-(2) – denoted $\bar{z}_h(t) = (I_h u(t), I_h v(t))^t$ where $I_h \in \mathcal{L}(C^0([0, 1]), U_h^r)$, is the interpolation operator. For the sake of simplicity and conciseness, we assume that $\bar{z}_h(t)$ is solution to

$$ \dot{\bar{z}}_h = A_h \bar{z}_h + \varepsilon_h, \quad \bar{z}_h(0) = (I_h u_0, I_h v_0)^t, \quad (56) $$

where $\varepsilon_h(t) \in \mathcal{Z}_h$ is the consistency error, smooth in time, and such that

$$ \sup_{t \in [0,T]} \|\varepsilon_h(t)\|_{L^2(0,1)^2} \le c_{\varepsilon} h^r, \quad (57) $$

with $c_{\varepsilon} > 0$ is independent of $T$. Finally, we assume that the interpolation operator satisfies that there exists $c_I > 0$, such that for $0 \le \ell \le k$ and $0 \le k \le r$, and all $h > 0$,

$$ \forall w \in H^r(0,1)^2, \quad \left( \sum_{i=0}^{N-1} |I_h w - w|^2_{H^\ell(x_i, x_{i+1})} \right)^{\frac{1}{2}} \le c_I h^{k-\ell} \|w\|_{H^k(0,1)}. \quad (58) $$

We refer to [9, 19, 22] for more details on how these properties can be obtained.

**Lemma 5.1.** There exists a constant $c_D > 0$ independent of $h$, such that

$$ \forall w \in H^r(0,1), \quad \|D_h^{\frac{1}{2}} I_h w\|_{L^2(0,1)} = \|D_{h,*}^{\frac{1}{2}} I_h w\|_{L^2(0,1)} \le c_D h^r \|w\|_{H^r(0,1)}, $$

and $J_h I_h w = 0$.
# 5 Application to a data assimilation problem

## 5.1 Motivation

One application of the results of exponential stabilization is to construct observers that use available measurements as feedback to improve a physical system approximation. In this case, the goal is to construct a numerical approximation to the solution of a wave propagation problem for which we have boundary measurements. More precisely, we assume that we have access to the boundary data $v(1, t)$ for all times, where $z = (u, v)^t$ is the solution of (1)-(2) with $\gamma = 0$ (hence the problem is conservative). Despite unknown initial data, we want to construct a numerical approximation $z_h = (u_h, v_h)^t$ that converges asymptotically to $z$ in time $t$ up to a discretization error. To do so, one can follow the approach proposed in [7, 30] and compute a so-called observer $\hat{z}_h$ as the solution to

$$ \dot{\hat{z}}_h = A_h \hat{z}_h - \gamma \begin{pmatrix} \tilde{D}_h & 0 \\ 0 & \tilde{D}_{h,*} + J_h \end{pmatrix} \hat{z}_h - \gamma \begin{pmatrix} 0 & 0 \\ 0 & B_h \end{pmatrix} (\hat{z}_h - \bar{z}_h), \quad \hat{z}_h(0) = 0, \quad (55) $$

where $\bar{z}_h$ is the interpolation of the continuous function $z$ solution to (1)-(2). Note that, by assumption, we only have access to $B_h \bar{z}_h(t)$, which corresponds to the measurements at the boundary. These measurements are used in a feedback law for the dynamics of the observer. Note finally that $\bar{z}$ is the solution of a conservative problem (i.e., $\gamma = 0$) and $(u_0, v_0)$ are unknown. The main question to answer concerns the convergence of $\hat{z}_h$ towards $z$.

## 5.2 Convergence result

In this section, we follow the approach of [31] and provide a convergence result of the solution of the exponentially stable discrete problem (55) towards the solution of the continuous problem (1)-(2). The most important feature of the proposed error estimate is that, unlike the standard error estimates, it is uniform in $T$. We introduce the interpolation of $z(t) = (u(t), v(t))^t$ – solution to(1)-(2) – denoted $\bar{z}_h(t) = (I_h u(t), I_h v(t))^t$ where $I_h \in \mathcal{L}(C^0([0, 1]), U_h^r)$, is the interpolation operator. For the sake of simplicity and conciseness, we assume that $\bar{z}_h(t)$ is solution to

$$ \dot{\bar{z}}_h = A_h \bar{z}_h + \varepsilon_h, \quad \bar{z}_h(0) = (I_h u_0, I_h v_0)^t, \quad (56) $$

where $\varepsilon_h(t) \in \mathcal{Z}_h$ is the consistency error, smooth in time, and such that

$$ \sup_{t \in [0,T]} \|\varepsilon_h(t)\|_{L^2(0,1)^2} \le c_{\varepsilon} h^r, \quad (57) $$

with $c_{\varepsilon} > 0$ is independent of $T$. Finally, we assume that the interpolation operator satisfies that there exists $c_I > 0$, such that for $0 \le \ell \le k$ and $0 \le k \le r$, and all $h > 0$,

$$ \forall w \in H^r(0,1)^2, \quad \left( \sum_{i=0}^{N-1} |I_h w - w|^2_{H^\ell(x_i, x_{i+1})} \right)^{\frac{1}{2}} \le c_I h^{k-\ell} \|w\|_{H^k(0,1)}. \quad (58) $$

We refer to [9, 19, 22] for more details on how these properties can be obtained.

**Lemma 5.1.** There exists a constant $c_D > 0$ independent of $h$, such that

$$ \forall w \in H^r(0,1), \quad \|D_h^{\frac{1}{2}} I_h w\|_{L^2(0,1)} = \|D_{h,*}^{\frac{1}{2}} I_h w\|_{L^2(0,1)} \le c_D h^r \|w\|_{H^r(0,1)}, $$

and $J_h I_h w = 0$.
*Proof.* From the property $\chi_h^*(I_h w) = \chi_h(I_h w)$ and thanks to Proposition 4.11, we have

$$ 
\begin{aligned}
d_{h,*}(I_h w, I_h w) &= d_h(I_h w, I_h w) = c_r \sum_{i=0}^{N-1} h_i^{2r+1} \left| (I_h w)^{(r)}\left(x_{i+\frac{1}{2}}\right) \right|^2 \\
&= c_r \sum_{i=0}^{N-1} \int_{x_i}^{x_{i+1}} h_i^{2r} \left| (I_h w)^{(r)}\left(x_{i+\frac{1}{2}}\right) \right|^2 \\
&\le c_r h^{2r} \sum_{i=0}^{N-1} |I_h w|_{H^r(x_i, x_{i+1})}^2.
\end{aligned}
$$

Using (58) and a triangular inequality, we get the existence of another scalar $c_r$ such that

$$ d_h(I_h w, I_h w) \le c_r h^{2r} \|w\|_{H^r(0,1)}^2. $$

Finally, note that since $I_h w$ is continuous, $J_h I_h w = 0$. □

**Theorem 5.2.** Let $z \in C^1([0, T]; H^r(0, 1)^2)$ be a solution to (10). Then, there exists positive scalars $c_{st} > 0$ and $\sigma > 0$ such that, for any $T > 0$ and $h > 0$, we have

$$ \forall t \in [0, T], \quad \|z(t) - z_h(t)\|_{L^2(0,1)^2} \le c_{st} (e^{-\sigma t} \|\bar{z}_h(0) - z_h(0)\|_h + h^r), $$

where $z_h$ is solution to (55).

*Proof.* Using the triangular inequality we have

$$ \|z(t) - z_h(t)\|_{L^2(0,1)^2} \le \|z(t) - \bar{z}_h\|_{L^2(0,1)^2} + \|\bar{z}_h - z_h(t)\|_{L^2(0,1)^2}. $$

Using the interpolation property (58), we obtain that there exists a positive scalar $c_{st}$ such that

$$ \|z(t) - z_h(t)\|_{L^2(0,1)^2} \le c_{st} h^r + \|e_h(t)\|_{L^2(0,1)^2}, \quad (59) $$

where $e_h = \bar{z}_h - z_h$. The error term $e_h$ satisfies

$$ \dot{e}_h = A_{\gamma,h} e_h + \bar{\varepsilon}_h, \quad e_h(0) = \bar{z}_h(0) - z_h(0), $$

where $A_{\gamma,h} = A_h - \gamma C_h^* C_h$ and

$$ \bar{\varepsilon}_h = \gamma \begin{pmatrix} D_h & 0 \\ 0 & D_{h,*} + J_h \end{pmatrix} \bar{z}_h + \varepsilon_h. $$

Using Lemma 5.1 and estimate (57), we have that, there exists $c_{st}$ independent of $h$ and $T$, such that

$$ \sup_{t \in [0, T]} \|\bar{\varepsilon}_h(t)\|_{L^2(0,1)^2} \le c_{st} h^r. $$

From the Duhamel formula, we recall that

$$ e_h(t) = e^{tA_{\gamma,h}} (\bar{z}_h(0) - z_h(0)) + \int_0^t e^{(t-s)A_{\gamma,h}} \bar{\varepsilon}_h(s) \, ds. $$

Using the exponential stability of the operator $A_{\gamma,h}$ – see Corollary 4.10 – we get that

$$ \|e_h(t)\|_h \le c_{st} \left( e^{-\sigma t} \|\bar{z}_h(0) - z_h(0)\|_h + \sup_{s \in [0,t]} \|\bar{\varepsilon}_h(s)\|_h \int_0^t e^{-\sigma(t-s)} \, ds \right). $$

As a consequence, we finally find that there exist a constant $c_{st} > 0$ independent of $h$ and $T$ such that

$$ \forall t \in [0, T], \quad \|e_h(t)\|_h \le c_{st} \left( e^{-\sigma t} \|\bar{z}_h(0) - z_h(0)\|_h + \frac{1}{\sigma} (1 - e^{-\sigma t}) h^r \right), $$

which combining with (59) and (22) proves the desired result. □
Proof. From the property $\chi_h^*(I_h w) = \chi_h(I_h w)$ and thanks to Proposition 4.11, we have

$$ 
\begin{aligned}
d_{h,*}(I_h w, I_h w) &= d_h(I_h w, I_h w) = c_r \sum_{i=0}^{N-1} h_i^{2r+1} \left| (I_h w)^{(r)}\left(x_{i+\frac{1}{2}}\right) \right|^2 \\
&= c_r \sum_{i=0}^{N-1} \int_{x_i}^{x_{i+1}} h_i^{2r} \left| (I_h w)^{(r)}\left(x_{i+\frac{1}{2}}\right) \right|^2 \\
&\le c_r h^{2r} \sum_{i=0}^{N-1} |I_h w|_{H^r(x_i, x_{i+1})}^2.
\end{aligned}
$$

Using (58) and a triangular inequality, we get the existence of another scalar $c_r$ such that

$$ d_h(I_h w, I_h w) \le c_r h^{2r} \|w\|_{H^r(0,1)}^2. $$

Finally, note that since $I_h w$ is continuous, $J_h I_h w = 0$. □

**Theorem 5.2.** Let $z \in C^1([0, T]; H^r(0,1)^2)$ be a solution to (10). Then, there exists positive scalars $c_{st} > 0$ and $\sigma > 0$ such that, for any $T > 0$ and $h > 0$, we have

$$ \forall t \in [0, T], \quad \|z(t) - z_h(t)\|_{L^2(0,1)^2} \le c_{st} (e^{-\sigma t} \|\bar{z}_h(0) - z_h(0)\|_h + h^r), $$

where $z_h$ is solution to (55).

Proof. Using the triangular inequality we have

$$ \|z(t) - z_h(t)\|_{L^2(0,1)^2} \le \|z(t) - \bar{z}_h\|_{L^2(0,1)^2} + \|\bar{z}_h - z_h(t)\|_{L^2(0,1)^2}. $$

Using the interpolation property (58), we obtain that there exists a positive scalar $c_{st}$ such that

$$ \|z(t) - z_h(t)\|_{L^2(0,1)^2} \le c_{st} h^r + \|e_h(t)\|_{L^2(0,1)^2}, \quad (59) $$

where $e_h = \bar{z}_h - z_h$. The error term $e_h$ satisfies

$$ \dot{e}_h = A_{\gamma,h} e_h + \bar{\varepsilon}_h, \quad e_h(0) = \bar{z}_h(0) - z_h(0), $$

where $A_{\gamma,h} = A_h - \gamma C_h^* C_h$ and

$$ \bar{\varepsilon}_h = \gamma \begin{pmatrix} D_h & 0 \\ 0 & D_{h,*} + J_h \end{pmatrix} \bar{z}_h + \varepsilon_h. $$

Using Lemma 5.1 and estimate (57), we have that, there exists $c_{st}$ independent of $h$ and $T$, such that

$$ \sup_{t \in [0,T]} \|\bar{\varepsilon}_h(t)\|_{L^2(0,1)^2} \le c_{st} h^r. $$

From the Duhamel formula, we recall that

$$ e_h(t) = e^{tA_{\gamma,h}} (\bar{z}_h(0) - z_h(0)) + \int_0^t e^{(t-s)A_{\gamma,h}} \bar{\varepsilon}_h(s) \, ds. $$

Using the exponential stability of the operator $A_{\gamma,h}$ – see Corollary 4.10 – we get that

$$ \|e_h(t)\|_h \le c_{st} \left( e^{-\sigma t} \|\bar{z}_h(0) - z_h(0)\|_h + \sup_{s \in [0,t]} \|\bar{\varepsilon}_h(s)\|_h \int_0^t e^{-\sigma(t-s)} \, ds \right). $$

As a consequence, we finally find that there exist a constant $c_{st} > 0$ independent of $h$ and $T$ such that

$$ \forall t \in [0, T], \quad \|e_h(t)\|_h \le c_{st} \left( e^{-\sigma t} \|\bar{z}_h(0) - z_h(0)\|_h + \frac{1}{\sigma} (1 - e^{-\sigma t}) h^r \right), $$

which combining with (59) and (22) proves the desired result. □
Theorem 5.2 shows that the effect of uncertainty on the initial condition vanishes exponentially fast and $\hat{z}_h$ is asymptotically in time an approximation of $z$ up to a standard discretization error in $O(h^r)$. More precisely, it becomes a better approximation than a standard discretization of the wave equation with known initial condition, since the error does not grow linearly with time.

# 6 The fully discrete setting

## 6.1 Abstract results

In this section, we are interested in the discretization of System (19) with an implicit midpoint scheme. We verify that the use of this time discretization does not alter the stabilization properties. The scheme is given by, for all $n \ge 0$,

$$ \left\{ \begin{array}{l} \frac{z_h^{n+1} - z_h^n}{\Delta t} = A_h z_h^{n+1/2} - \gamma C_h^* C_h z_h^{n+1/2}, \quad z_h^{n+1/2} = \frac{z_h^{n+1} + z_h^n}{2}. \quad (60) \\ z_h^0 = z_{h,0}, \end{array} \right. $$

It is easy to show an energy conservation property for this scheme.

**Theorem 6.1.** Let $\{z_h^n\}_{n \ge 0} \in \mathcal{Z}_h$ be a solution to (60), we have the discrete dissipation energy property

$$ \forall n \ge 0, \quad \frac{1}{2} \frac{\|z_h^{n+1}\|_h^2 - \|z_h^n\|_h^2}{\Delta t} = -\gamma \|C_h z_h^{n+1/2}\|_h^2. \quad (61) $$

*Proof.* We multiply (60) by $z_h^{n+1/2}$. Since $A_h$ satisfies (18) and using $(a-b)(a+b) = a^2 - b^2$ we obtain the result. $\square$

We now demonstrate observability and exponential stability, as in the previous sections, but at the fully discrete level. An important result is that while the scheme (60) is stable for any time step $\Delta t$ (in the sense that the solution remains bounded), exponential stability requires a condition on the time step that is typical of CFL conditions observed when using explicit time discretization.

Following the same reasoning as in the semi-discrete case, we start by demonstrating the observability property for the undamped case $\gamma = 0$.

**Theorem 6.2.** Let Assumption 1 holds. For $N \in \mathbb{N}$ large enough, there exists a constant $c_{st}(N) > 0$ independent of $h$ such that, any solutions $\{z_h^n\}_{n \ge 0} \in \mathcal{Z}_h$ to (60) with $\gamma = 0$ satisfies

$$ \|z_{h,0}\|_h^2 \le c_{st}(N) \Delta t \sum_{n=0}^{N} \|C_h z_h^{n+1/2}\|_h^2. $$

*Proof.* The proof follows the same steps as in the proof of Theorem 4.3

**Step 1:** We multiply the first equation of (60) (with $\gamma = 0$) by $\Pi_h z_h^{n+1/2}$, and we get

$$ \left( \Pi_h \frac{z_h^{n+1} - z_h^n}{\Delta t}, z_h^{n+1/2} \right)_h = \left( \Pi_h A_h z_h^{n+1/2}, z_h^{n+1/2} \right)_h. $$

Using $(a-b)(a+b) = a^2 - b^2$ and Assumption 1 we find,

$$ \frac{1}{2\Delta t} (\Pi_h z_h^{n+1}, z_h^{n+1})_h - \frac{1}{2\Delta t} (\Pi_h z_h^n, z_h^n)_h = \left( \Pi_h A_h z_h^{n+1/2}, z_h^{n+1/2} \right)_h, \quad (62) \\ \le \alpha \|C_h z_h^{n+1/2}\|_h^2 - \beta \|z_h^{n+1/2}\|_h^2. $$

**Step 2:** Summing (62) over $n$ we get,

$$ \frac{1}{2\Delta t} \left( \Pi_h z_h^{N+1}, z_h^{N+1} \right)_h - \frac{1}{2\Delta t} \left( \Pi_h z_h^0, z_h^0 \right)_h \le \alpha \sum_{n=0}^{N} \|C_h z_h^{n+1/2}\|_h^2 - \beta \sum_{n=0}^{N} \|z_h^{n+1/2}\|_h^2. $$
Theorem 5.2 shows that the effect of uncertainty on the initial condition vanishes exponentially fast and $\hat{z}_h$ is asymptotically in time an approximation of $z$ up to a standard discretization error in $O(h^r)$. More precisely, it becomes a better approximation than a standard discretization of the wave equation with known initial condition, since the error does not grow linearly with time.

# 6 The fully discrete setting

## 6.1 Abstract results

In this section, we are interested in the discretization of System (19) with an implicit midpoint scheme. We verify that the use of this time discretization does not alter the stabilization properties. The scheme is given by, for all $n \ge 0$,

$$ \left\{ \begin{array}{l} \frac{z_h^{n+1} - z_h^n}{\Delta t} = A_h z_h^{n+1/2} - \gamma C_h^* C_h z_h^{n+1/2}, \quad z_h^{n+1/2} = \frac{z_h^{n+1} + z_h^n}{2}. \quad (60) \\ z_h^0 = z_{h,0}, \end{array} \right. $$

It is easy to show an energy conservation property for this scheme.

**Theorem 6.1.** Let $\{z_h^n\}_{n \ge 0} \in \mathcal{Z}_h$ be a solution to (60), we have the discrete dissipation energy property

$$ \forall n \ge 0, \quad \frac{1}{2} \frac{\|z_h^{n+1}\|_h^2 - \|z_h^n\|_h^2}{\Delta t} = -\gamma \|C_h z_h^{n+1/2}\|_h^2. \quad (61) $$

*Proof.* We multiply (60) by $z_h^{n+1/2}$. Since $A_h$ satisfies (18) and using $(a-b)(a+b) = a^2 - b^2$ we obtain the result. $\square$

We now demonstrate observability and exponential stability, as in the previous sections, but at the fully discrete level. An important result is that while the scheme (60) is stable for any time step $\Delta t$ (in the sense that the solution remains bounded), exponential stability requires a condition on the time step that is typical of CFL conditions observed when using explicit time discretization.

Following the same reasoning as in the semi-discrete case, we start by demonstrating the observability property for the undamped case $\gamma = 0$.

**Theorem 6.2.** Let Assumption 1 holds. For $N \in \mathbb{N}$ large enough, there exists a constant $c_{st}(N) > 0$ independent of $h$ such that, any solutions $\{z_h^n\}_{n \ge 0} \in \mathcal{Z}_h$ to (60) with $\gamma = 0$ satisfies

$$ \|z_{h,0}\|_h^2 \le c_{st}(N) \Delta t \sum_{n=0}^{N} \|C_h z_h^{n+1/2}\|_h^2. $$

*Proof.* The proof follows the same steps as in the proof of Theorem 4.3

**Step 1:** We multiply the first equation of (60) (with $\gamma = 0$) by $\Pi_h z_h^{n+1/2}$, and we get

$$ \left( \Pi_h \frac{z_h^{n+1} - z_h^n}{\Delta t}, z_h^{n+1/2} \right)_h = \left( \Pi_h A_h z_h^{n+1/2}, z_h^{n+1/2} \right)_h. $$

Using $(a-b)(a+b) = a^2 - b^2$ and Assumption 1 we find,

$$ \frac{1}{2\Delta t} (\Pi_h z_h^{n+1}, z_h^{n+1})_h - \frac{1}{2\Delta t} (\Pi_h z_h^n, z_h^n)_h = \left( \Pi_h A_h z_h^{n+1/2}, z_h^{n+1/2} \right)_h, \quad (62) \\ \le \alpha \|C_h z_h^{n+1/2}\|_h^2 - \beta \|z_h^{n+1/2}\|_h^2. $$

**Step 2:** Summing (62) over $n$ we get,

$$ \frac{1}{2\Delta t} (\Pi_h z_h^{N+1}, z_h^{N+1})_h - \frac{1}{2\Delta t} (\Pi_h z_h^0, z_h^0)_h \le \alpha \sum_{n=0}^{N} \|C_h z_h^{n+1/2}\|_h^2 - \beta \sum_{n=0}^{N} \|z_h^{n+1/2}\|_h^2. $$
The fully discrete energy conservation property of Theorem 6.1 with $\gamma = 0$ gives

$$ \frac{1}{2\Delta t} \left( \Pi_h z_h^{N+1}, z_h^{N+1} \right)_h - \frac{1}{2\Delta t} \left( \Pi_h z_h^0, z_h^0 \right)_h \le \alpha \sum_{n=0}^{N} \|C_h z_h^{n+1/2}\|_h^2 - \beta(N+1) \|z_h^0\|_h^2. $$

**Step 3:** Since $\Pi_h$ is bounded, by the Cauchy-Schwarz inequality and using again the energy conservation theorem 6.1 with $\gamma = 0$, we have

$$ \left( \beta(N+1)\Delta t - \|\Pi_h\|_h \right) \|z_{0,h}\|_h^2 \le \alpha \sum_{n=0}^{N} \Delta t \|C_h z_h^{n+1/2}\|_h^2. $$

Choosing $(N+1)\Delta t > \|\Pi_h\|_h \beta^{-1}$, we obtain the result of the theorem. □

Proving the exponential stability property for the damped system is more involved and is obtained under a condition on the time step.

**Theorem 6.3.** Let Assumption 2 holds with $\gamma > 0$. and assumme that there exists $c_c > 0$ and $\beta > 0$ such that

$$ \Delta t \|C_h\|_h \le c_c, \qquad (63) $$

and

$$ \frac{\Delta t}{2} \|A_h\|_h \le \beta < 1. \qquad (64) $$

Then, there exist $c_s > 0$ and $\sigma > 0$ independent of $h$ such that any solution $\{z_h^n\}_{n \ge 0} \in \mathcal{Z}_h$ to (60) satisfies

$$ \forall n \ge 0, \quad \|z_h^n\|_h^2 \le c_s \|z_{h,0}\|_h^2 e^{-\sigma(1-\beta^2)n\Delta t}. $$

*Proof.* The proof follows the same steps as in the proof of Theorem 3.4.

**Step 1:** We multiply the first equation by (60) and $\Pi_h z_h^{n+1/2}$, and we get

$$ \left( \Pi_h \frac{z_h^{n+1} - z_h^n}{\Delta t}, z_h^{n+1/2} \right)_h = (\Pi_h A_h z_h^{n+1/2}, z_h^{n+1/2})_h - \gamma (C_h z_h^{n+1/2}, C_h \Pi_h z_h^{n+1/2})_h. $$

Using $(a+b)(a-b) = a^2 - b^2$ and Assumption 2, we have

$$ \frac{1}{2\Delta t} (\Pi_h z_h^{n+1}, z_h^{n+1})_h - \frac{1}{2\Delta t} (\Pi_h z_h^n, z_h^n)_h \le \chi \|C_h z_h^{n+1/2}\|_h^2 - \kappa \|z_h^{n+1/2}\|_h^2. $$

**Step 2:** Using the energy identity of Theorem 6.1, we get with $\alpha > 0$,

$$ \frac{1}{2\Delta t} (\Pi_h z_h^{n+1}, z_h^{n+1})_h - \frac{1}{2\Delta t} (\Pi_h z_h^n, z_h^n)_h + \alpha \frac{\|z_h^{n+1}\|_h^2 - \|z_h^n\|_h^2}{2\Delta t} \\ \le (\chi - \alpha\gamma) \|C_h z_h^{n+1/2}\|_h^2 - \kappa \|z_h^{n+1/2}\|_h^2. \quad (65) $$

We intend to bound the last term of this inequality by $-\|z_h^{n+1}\|_h^2 - \|z_h^n\|_h^2$. Note that this intermediate necessary step is specific to the time discretization. Basic algebraic manipulations give

$$ -\|z_h^{n+1/2}\|_h^2 = -\frac{1}{2}\|z_h^{n+1}\|_h^2 - \frac{1}{2}\|z_h^n\|_h^2 + \frac{(\Delta t)^2}{4} \left\| \frac{z_h^{n+1} - z_h^n}{\Delta t} \right\|_h^2. $$

Using the scheme (60) and the Cauchy-Schwarz and Young inequalities, we obtain, for all $\epsilon > 0$,

$$ - \|z_h^{n+1/2}\|_h^2 \le -\frac{1}{2} \|z_h^{n+1}\|_h^2 - \frac{1}{2} \|z_h^n\|_h^2 + \frac{(\Delta t)^2}{4} (1 + \epsilon) \|A_h z_h^{n+1/2}\|_h^2 \\ + \frac{\gamma (\Delta t)^2}{4} (1 + \epsilon^{-1}) \|C_h^*\|_h^2 \|C_h z_h^{n+1/2}\|_h^2. $$
The fully discrete energy conservation property of Theorem 6.1 with $\gamma = 0$ gives

$$ \frac{1}{2\Delta t} \left( \Pi_h z_h^{N+1}, z_h^{N+1} \right)_h - \frac{1}{2\Delta t} \left( \Pi_h z_h^0, z_h^0 \right)_h \le \alpha \sum_{n=0}^{N} \|C_h z_h^{n+1/2}\|_h^2 - \beta(N+1) \|z_h^0\|_h^2. $$

**Step 3:** Since $\Pi_h$ is bounded, by the Cauchy-Schwarz inequality and using again the energy conservation theorem 6.1 with $\gamma = 0$, we have

$$ \left( \beta(N+1)\Delta t - \|\Pi_h\|_h \right) \|z_{0,h}\|_h^2 \le \alpha \sum_{n=0}^{N} \Delta t \|C_h z_h^{n+1/2}\|_h^2. $$

Choosing $(N+1)\Delta t > \|\Pi_h\|_h \beta^{-1}$, we obtain the result of the theorem. □

Proving the exponential stability property for the damped system is more involved and is obtained under a condition on the time step.

**Theorem 6.3.** Let Assumption 2 holds with $\gamma > 0$. and assumme that there exists $c_c > 0$ and $\beta > 0$ such that

$$ \Delta t \|C_h\|_h \le c_c, \qquad (63) $$

and

$$ \frac{\Delta t}{2} \|A_h\|_h \le \beta < 1. \qquad (64) $$

Then, there exist $c_s > 0$ and $\sigma > 0$ independent of $h$ such that any solution $\{z_h^n\}_{n \ge 0} \in \mathcal{Z}_h$ to (60) satisfies

$$ \forall n \ge 0, \quad \|z_h^n\|_h^2 \le c_s \|z_{h,0}\|_h^2 e^{-\sigma(1-\beta^2)n\Delta t}. $$

*Proof.* The proof follows the same steps as in the proof of Theorem 3.4.

**Step 1:** We multiply the first equation by (60) and $\Pi_h z_h^{n+1/2}$, and we get

$$ \left( \Pi_h \frac{z_h^{n+1} - z_h^n}{\Delta t}, z_h^{n+1/2} \right)_h = (\Pi_h A_h z_h^{n+1/2}, z_h^{n+1/2})_h - \gamma (C_h z_h^{n+1/2}, C_h \Pi_h z_h^{n+1/2})_h. $$

Using $(a+b)(a-b) = a^2 - b^2$ and Assumption 2, we have

$$ \frac{1}{2\Delta t} (\Pi_h z_h^{n+1}, z_h^{n+1})_h - \frac{1}{2\Delta t} (\Pi_h z_h^n, z_h^n)_h \le \chi \|C_h z_h^{n+1/2}\|_h^2 - \kappa \|z_h^{n+1/2}\|_h^2. $$

**Step 2:** Using the energy identity of Theorem 6.1, we get with $\alpha > 0$,

$$ \frac{1}{2\Delta t} (\Pi_h z_h^{n+1}, z_h^{n+1})_h - \frac{1}{2\Delta t} (\Pi_h z_h^n, z_h^n)_h + \alpha \frac{\|z_h^{n+1}\|_h^2 - \|z_h^n\|_h^2}{2\Delta t} \\ \le (\chi - \alpha\gamma) \|C_h z_h^{n+1/2}\|_h^2 - \kappa \|z_h^{n+1/2}\|_h^2. \quad (65) $$

We intend to bound the last term of this inequality by $-\|z_h^{n+1}\|_h^2 - \|z_h^n\|_h^2$. Note that this intermediate necessary step is specific to the time discretization. Basic algebraic manipulations give

$$ -\|z_h^{n+1/2}\|_h^2 = -\frac{1}{2}\|z_h^{n+1}\|_h^2 - \frac{1}{2}\|z_h^n\|_h^2 + \frac{(\Delta t)^2}{4} \left\| \frac{z_h^{n+1} - z_h^n}{\Delta t} \right\|_h^2. $$

Using the scheme (60) and the Cauchy-Schwarz and Young inequalities, we obtain, for all $\epsilon > 0$,

$$ - \|z_h^{n+1/2}\|_h^2 \le -\frac{1}{2} \|z_h^{n+1}\|_h^2 - \frac{1}{2} \|z_h^n\|_h^2 + \frac{(\Delta t)^2}{4} (1 + \epsilon) \|A_h z_h^{n+1/2}\|_h^2 \\ + \frac{\gamma (\Delta t)^2}{4} (1 + \epsilon^{-1}) \|C_h^*\|_h^2 \|C_h z_h^{n+1/2}\|_h^2. $$
By triangular inequality and Cauchy-Schwarz inequality and using hypothesis (63), we get

$$ - \|z_h^{n+1/2}\|_h^2 \le -\frac{1}{2} \left( 1 - \frac{(\Delta t)^2}{4} (1+\epsilon) \|A_h\|_h^2 \right) \left( \|z_h^{n+1}\|_h^2 + \|z_h^n\|_h^2 \right) \\ + \frac{\gamma c_c^2}{4} (1+\epsilon^{-1}) \|C_h z_h^{n+1/2}\|_h^2. $$

We choose $\epsilon = \frac{1}{2\beta^2} - \frac{1}{2} > 0$ and use the condition (64), we obtain

$$ -\|z_h^{n+1/2}\|_h^2 \le -\frac{1-\beta^2}{4} (\|z_h^{n+1}\|_h^2 + \|z_h^n\|_h^2) + \frac{\gamma c_c^2}{4} \left(1 + \frac{2\beta^2}{1-\beta^2}\right) \|C_h z_h^{n+1/2}\|_h^2. $$

Using the above equation in (65), we get

$$ \frac{1}{2\Delta t} (\Pi_h z_h^{n+1}, z_h^{n+1})_h - \frac{1}{2\Delta t} (\Pi_h z_h^n, z_h^n)_h + \alpha \frac{\|z_h^{n+1}\|_h^2 - \|z_h^n\|_h^2}{2\Delta t} \\ \le \left( \frac{\kappa \gamma c_c^2}{4} \frac{1+\beta^2}{1-\beta^2} + \chi - \alpha \gamma \right) \|C_h z_h^{n+1/2}\|_h^2 - \frac{\kappa}{4} (1-\beta^2) (\|z_h^{n+1}\|_h^2 + \|z_h^n\|_h^2). $$

**Step 3:** We choose $\alpha$ large enough such that $\alpha \ge \frac{\chi}{\gamma} + \frac{\kappa c_c^2}{4} \frac{1+\beta^2}{1-\beta^2}$. Therefore, we have

$$ \frac{1}{2\Delta t} (\Pi_h z_h^{n+1}, z_h^{n+1})_h - \frac{1}{2\Delta t} (\Pi_h z_h^n, z_h^n)_h + \alpha \frac{\|z_h^{n+1}\|_h^2 - \|z_h^n\|_h^2}{2\Delta t} \\ \le -\frac{\kappa}{4} (1-\beta^2) (\|z_h^{n+1}\|_h^2 + \|z_h^n\|_h^2). \quad (66) $$

**Step 4:** We introduce the fully discrete Lyapunov functional

$$ \mathcal{L}_h^n = \frac{1}{2} (\Pi_h z_h^n, z_h^n)_h + \frac{\alpha}{2} \|z_h^n\|_h^2, $$

and we can also show that it is equivalent to the energy functional $\mathcal{E}_h^n := \frac{1}{2} \|z_h^n\|_h^2$ using the bounded nature of $\Pi_h$ and choosing $\alpha$ sufficiently large,

$$ (\alpha - \|\Pi_h\|_h) \mathcal{E}_h^n \le \mathcal{L}_h^n \le (\alpha + \|\Pi_h\|_h) \mathcal{E}_h^n. \quad (67) $$

Using this equivalence between the Lyapunov functional and the energy, we obtain

$$ \frac{1}{\Delta t} (\mathcal{L}_h^{n+1} - \mathcal{L}_h^n) \le -\frac{\kappa}{2} (1 - \beta^2) (\|\Pi_h\|_h + \alpha)^{-1} (\mathcal{L}_h^{n+1} + \mathcal{L}_h^n). $$

Rearranging the terms, we get

$$ \mathcal{L}_h^{n+1} \le \delta(\Delta t, \alpha) \mathcal{L}_h^n $$

with

$$ \delta(\Delta t, \alpha) = \left( 1 - \frac{\Delta t \frac{\kappa}{2} (1 - \beta^2)}{\|\Pi_h\|_h + \alpha} \right) / \left( 1 + \frac{\Delta t \frac{\kappa}{2} (1 - \beta^2)}{\|\Pi_h\|_h + \alpha} \right). $$

**Step 5:** By induction and using (67), we obtain

$$ \|z_h^n\|_h^2 \le \frac{\alpha + \|\Pi_h\|_h}{\alpha - \|\Pi_h\|_h} \|z_{h,0}\|_h^2 e^{n \ln(\delta(\Delta t, \alpha))}. $$
By triangular inequality and Cauchy-Schwarz inequality and using hypothesis (63), we get

$$ - \|z_h^{n+1/2}\|_h^2 \leq -\frac{1}{2} \left( 1 - \frac{(\Delta t)^2}{4} (1 + \epsilon) \|A_h\|_h^2 \right) \left( \|z_h^{n+1}\|_h^2 + \|z_h^n\|_h^2 \right) \\ + \frac{\gamma c_c^2}{4} (1 + \epsilon^{-1}) \|C_h z_h^{n+1/2}\|_h^2. $$

We choose $\epsilon = \frac{1}{2\beta^2} - \frac{1}{2} > 0$ and use the condition (64), we obtain

$$ -\|z_h^{n+1/2}\|_h^2 \leq -\frac{1-\beta^2}{4} (\|z_h^{n+1}\|_h^2 + \|z_h^n\|_h^2) + \frac{\gamma c_c^2}{4} \left(1 + \frac{2\beta^2}{1-\beta^2}\right) \|C_h z_h^{n+1/2}\|_h^2. $$

Using the above equation in (65), we get

$$ \frac{1}{2\Delta t} (\Pi_h z_h^{n+1}, z_h^{n+1})_h - \frac{1}{2\Delta t} (\Pi_h z_h^n, z_h^n)_h + \alpha \frac{\|z_h^{n+1}\|_h^2 - \|z_h^n\|_h^2}{2\Delta t} \\ \leq \left( \frac{\kappa \gamma c_c^2}{4} \frac{1+\beta^2}{1-\beta^2} + \chi - \alpha \gamma \right) \|C_h z_h^{n+1/2}\|_h^2 - \frac{\kappa}{4} (1-\beta^2) (\|z_h^{n+1}\|_h^2 + \|z_h^n\|_h^2). $$

**Step 3:** We choose $\alpha$ large enough such that $\alpha \geq \frac{\chi}{\gamma} + \frac{\kappa c_c^2}{4} \frac{1+\beta^2}{1-\beta^2}$. Therefore, we have

$$ \frac{1}{2\Delta t} (\Pi_h z_h^{n+1}, z_h^{n+1})_h - \frac{1}{2\Delta t} (\Pi_h z_h^n, z_h^n)_h + \alpha \frac{\|z_h^{n+1}\|_h^2 - \|z_h^n\|_h^2}{2\Delta t} \\ \leq -\frac{\kappa}{4} (1-\beta^2) (\|z_h^{n+1}\|_h^2 + \|z_h^n\|_h^2). \quad (66) $$

**Step 4:** We introduce the fully discrete Lyapunov functional

$$ \mathcal{L}_h^n = \frac{1}{2} (\Pi_h z_h^n, z_h^n)_h + \frac{\alpha}{2} \|z_h^n\|_h^2, $$

and we can also show that it is equivalent to the energy functional $\mathcal{E}_h^n := \frac{1}{2} \|z_h^n\|_h^2$ using the bounded nature of $\Pi_h$ and choosing $\alpha$ sufficiently large,

$$ (\alpha - \|\Pi_h\|_h) \mathcal{E}_h^n \leq \mathcal{L}_h^n \leq (\alpha + \|\Pi_h\|_h) \mathcal{E}_h^n. \quad (67) $$

Using this equivalence between the Lyapunov functional and the energy, we obtain

$$ \frac{1}{\Delta t} (\mathcal{L}_h^{n+1} - \mathcal{L}_h^n) \leq -\frac{\kappa}{2} (1 - \beta^2) (\|\Pi_h\|_h + \alpha)^{-1} (\mathcal{L}_h^{n+1} + \mathcal{L}_h^n). $$

Rearranging the terms, we get

$$ \mathcal{L}_h^{n+1} \leq \delta(\Delta t, \alpha) \mathcal{L}_h^n $$

with

$$ \delta(\Delta t, \alpha) = \left( 1 - \frac{\Delta t \frac{\kappa}{2} (1 - \beta^2)}{\|\Pi_h\|_h + \alpha} \right) / \left( 1 + \frac{\Delta t \frac{\kappa}{2} (1 - \beta^2)}{\|\Pi_h\|_h + \alpha} \right). $$

**Step 5:** By induction and using (67), we obtain

$$ \|z_h^n\|_h^2 \leq \frac{\alpha + \|\Pi_h\|_h}{\alpha - \|\Pi_h\|_h} \|z_{h,0}\|_h^2 e^{n \ln(\delta(\Delta t, \alpha))}. $$
We must now specify the constant ln($\delta$). Since $\Delta t \ge 0$ and $1 - \beta^2 \ge 0$, we get using a series expansion,

$$ \ln(\delta(\Delta t, \alpha)) \le \ln \left( 1 - \Delta t \frac{\kappa(1 - \beta^2)}{\|\Pi_h\|_h + \alpha} \right) = -\frac{\Delta t}{\alpha} \kappa(1 - \beta^2) + o\left(\frac{1}{\alpha}\right). $$

Finally, we can conclude that, for $\alpha$ sufficiently large,

$$ \|z_h^n\|_h^2 \le \left( \frac{\alpha + \|\Pi_h\|_h}{\alpha - \|\Pi_h\|_h} \right) \|z_{h,0}\|_h^2 e^{-\sigma n \Delta t}, $$

with $\sigma > 0$. □

**Remark 6.4.** We show in Section 6.2 that $c_A$ exists independently of $h$ such that,

$$ \|C_h\|_h \le c_A \|A_h\|_h. \qquad (68) $$

In this case, the condition (63) is automatically checked as soon as the condition (64) is fulfilled. Note that, as far as the discretization of hyperbolic PDE problems is concerned, $A_h$ generally satisfies the following conditions

$$ \|A_h\|_h \approx h^{-1}, $$

therefore the condition (64) implies

$$ \Delta t \lesssim h. $$

which is a typical stability condition for explicit schemes. This means that, although the considered implicit scheme (60) is unconditionally stable, we have proved that the exponential stability is obtained under a CFL condition. We do not yet have a proof that this result is optimal. This remains an open question.

## 6.2 Application to the 1D fully discrete wave problem

We can now present the observability theorem in a totally discrete framework.

**Theorem 6.5.** For $N \in \mathbb{N}$ large enough, there exists $c_{st}(N)$ independent of $h$, such that any $\{z_h\}_{n=0}^N = \{(u_h, v_h)^t\}_{n=0}^N \in \mathcal{Z}_h$ solution to (27) with $\gamma = 0$, satisfies

$$ \|u_{h,0}\|_h^2 + \|v_{h,0}\|_h^2 \le c_{st}(N) \sum_{n=1}^N \Delta t \left[ d_h \left( u_h^{n+1/2}, u_h^{n+1/2} \right) + d_{h,*} \left( v_h^{n+1/2}, v_h^{n+1/2} \right) \right. \\ \left. \qquad + \left( v_h^{n+1/2}(1) \right)^2 \right] dt + \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \left[ \|v_i^{n+1/2}\|_h^2 \right]^2, $$

where $d_h$ and $d_h^*$ are defined by (32) and by (33).

*Proof.* Since $\Pi_h$ satisfies Assumption 1, we can apply Theorem 6.2 and the result follows immediately. □

**Theorem 6.6.** Assume that $\exists \beta > 0$ independent of $h$ and $\Delta t$, such that

$$ \frac{\Delta t}{2} \|A_h\|_h \le \beta < 1. \qquad (69) $$

Then, there exist $c_s > 0$ and $\sigma > 0$ independent of $h$ such that any $\{z_h\}_{n \ge 0} = \{(u_h, v_h)^t\}_{n \ge 0} \subset \mathcal{Z}_h$ solution to (27) with $\gamma > 0$ satisfies

$$ \|u_h^n\|_h^2 + \|v_h^n\|_h^2 \le c_s e^{-n \Delta t \sigma(\Delta t)} (\|u_{h,0}\|_h^2 + \|v_{h,0}\|_h^2). $$
We must now specify the constant ln(δ). Since Δt ≥ 0 and 1 − β² ≥ 0 , we get using a series expansion,

$$
\ln(\delta(\Delta t, \alpha)) \le \ln \left( 1 - \Delta t \frac{\kappa(1 - \beta^2)}{\|\Pi_h\|_h + \alpha} \right) = -\frac{\Delta t}{\alpha} \kappa(1 - \beta^2) + o\left(\frac{1}{\alpha}\right).
$$

Finally, we can conclude that, for α sufficiently large,

$$
\|z_h^n\|_h^2 \le \left( \frac{\alpha + \|\Pi_h\|_h}{\alpha - \|\Pi_h\|_h} \right) \|z_{h,0}\|_h^2 e^{-\sigma n \Delta t},
$$

with σ > 0.

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

**Remark 6.4.** We show in Section 6.2 that c_A exists independently of h such that,

$$
\|C_h\|_h \le c_A \|A_h\|_h. \tag{68}
$$

In this case,the condition (63) is automatically checked as soon asthe condition (64)
is fulfilled. Note that, as far as the discretization of hyperbolic PDE problems is
concerned, A_h generally satisfies the following conditions

$$
\|A_h\|_h \approx h^{-1},
$$

therefore the condition (64) implies

$$
\Delta t \lesssim h.
$$

which is a typical stability condition for explicit schemes. This means that, although
the considered implicit scheme (60) is unconditionally stable, we have proved that
the exponential stability is obtained under a CFL condition. We do not yet have a
proof that this result is optimal. This remains an open question.

**6.2 Application to the 1D fully discrete wave problem**

We can now present the observability theorem in a totally discrete framework.

**Theorem 6.5.** For $N \in \mathbb{N}$ large enough, there exists $c_{st}(N)$ independent of $h$, such that any $\{z_h\}_{n=0}^N = \{(u_h, v_h)^t\}_{n=0}^N \in \mathcal{Z}_h$ solution to (27) with $\gamma = 0$, satisfies

$$
\|u_{h,0}\|_h^2 + \|v_{h,0}\|_h^2 \le c_{st}(N) \sum_{n=1}^N \Delta t \left[ d_h \left( u_h^{n+1/2}, u_h^{n+1/2} \right) + d_{h,*} \left( v_h^{n+1/2}, v_h^{n+1/2} \right) \right. \\
\left. \qquad + \left( v_h^{n+1/2}(1) \right)^2 \right] dt + \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \left[ \|v_i^{n+1/2}\|_h^2 \right]^2,
$$

where $d_h$ and $d_h^*$ are defined by (32) and by (33).

Proof. Since Πh satisfies Assumption 1, we can apply Theorem 6.2 and the result
follows immediately.
□

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

**Theorem 6.6.** Assume that $\exists \beta > 0$ independent of $h$ and $\Delta t$, such that

$$
\frac{\Delta t}{2} \|A_h\|_h \le \beta < 1. \tag{69}
$$

Then, there exist $c_s > 0$ and $\sigma > 0$ independent of $h$ such that any $\{z_h\}_{n \ge 0} = \{(u_h, v_h)^t\}_{n \ge 0} \subset \mathcal{Z}_h$ solution to (27) with $\gamma > 0$ satisfies

$$
\|u_h^n\|_h^2 + \|v_h^n\|_h^2 \le c_s e^{-n\Delta t\sigma(\Delta t)} (\|u_{h,0}\|_h^2 + \|v_{h,0}\|_h^2).
$$
Proof. Since $\Pi_h$ satisfies Assumption 2, in order to prove Theorem 6.6 using Theorem 6.3, we need to show that the conditions (63) and (64) hold. Note that we have the following property

$$ \|A_h\|_h \le c_A \|R_h\|_h, c_A > 0, $$

The idea is to use Remark 6.4 – and thus to get rid of the condition (63) –. We need to show that

$$ \|C_h\|_h \le c_c \|R_h\|_h, c_c > 0. $$

The following inequality can be proved

$$ \frac{\alpha}{h} \le \|R_h\|_h, \alpha > 0, $$

so we only need to prove that $\|C_h\|_h \le c_{st}/h$ with $c_{st} > 0$. Let us compute the norm of $C_h$. By definition of the operator $C_h$, we have for $z_h \in \mathcal{Z}_h$,

$$ \|C_h z_h\|_h^2 = d_h(u_h, u_h) + d_{h,*}(v_h, v) + v_h^2(1) + \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2. $$

The two first terms are treated using Proposition 4.5. The third one using a trace inequality and an inverse estimate can be bounded by $O(1/\sqrt{h})$. With computations similar to what was done in the proof of Lemma 4.6 we arrive at,

$$ \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 \le 2 \sum_{i=0}^{N-1} \frac{1+2c_r}{\tilde{h}_i h_i} \|v_h\|_{L^2([x_i, x_{i+1}])}^2 + 2 \sum_{i=1}^{N} \frac{1+2c_r}{\tilde{h}_i h_{i-1}} \|v_h\|_{L^2([x_{i-1}, x_i])}^2. $$

Since the subdivision of the domain is quasi-uniform we have $h_i \ge \rho h$ hence,

$$ \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 \le 4 \frac{1+2c_r}{\rho^2 h^2} \|v_h\|_{L^2(0,1)}^2, $$

this shows that $\|C_h\|_h$ is bounded by $O(1/h)$. $\square$

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

# 7 Numerical results

In a first part, we perform a numerical eigenvalue study where we use our stabilization terms with respect to the ones classically proposed in the literature. In a second part, we introduce a fully explicit leap-frog scheme and provide convergence plots with respect to $h$ as well as energy plots over time. The results obtained illustrate the efficiency of the method. Note that in this section we use the simplified expression for the stabilization term discussed in Remark 4.13 and denote by $\tilde{D}_h$ the operator associated with the bilinear form $\tilde{d}_h$ defined by (54). In the following results, the domain $(0, 1)$ is evenly subdivided, i.e. $h = h_i = \tilde{h}_i$.

## 7.1 Eigenvalues investigation

A first idea is to stabilize the semi-discrete problem (23) is to stabilize the jumps of the discontinuous functions, as seen in [25] for Maxwell's equations. We define

$$ A_{h,\beta} = A_h - \begin{pmatrix} 0 & 0 \\ 0 & B_h \end{pmatrix} $$

and define the operator

$$ \forall z_h = \begin{pmatrix} u_h \\ v_h \end{pmatrix} \in \mathcal{Z}_h, \quad \tilde{z}_h = \begin{pmatrix} \tilde{u}_h \\ \tilde{v}_h \end{pmatrix} \in \mathcal{Z}_h, \quad (J_h z_h, \tilde{z}_h)_h = \frac{1}{h} \sum_{i=1}^{N} \|v_h(x_i)\| \|\tilde{v}_h(x_i)\| $$
*Proof.* Since $\Pi_h$ satisfies Assumption 2, in order to prove Theorem 6.6 using Theorem 6.3, we need to show that the conditions (63) and (64) hold. Note that we have the following property

$$ \|A_h\|_h \le c_A \|R_h\|_h, c_A > 0, $$

The idea is to use Remark 6.4 – and thus to get rid of the condition (63) –. We need to show that

$$ \|C_h\|_h \le c_c \|R_h\|_h, c_c > 0. $$

The following inequality can be proved

$$ \frac{\alpha}{h} \le \|R_h\|_h, \alpha > 0, $$

so we only need to prove that $\|C_h\|_h \le c_{st}/h$ with $c_{st} > 0$. Let us compute the norm of $C_h$. By definition of the operator $C_h$, we have for $z_h \in \mathcal{Z}_h$,

$$ \|C_h z_h\|_h^2 = d_h(u_h, u_h) + d_{h,*}(v_h, v) + v_h^2(1) + \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2. $$

The two first terms are treated using Proposition 4.5. The third one using a trace inequality and an inverse estimate can be bounded by $O(1/\sqrt{h})$. With computations similar to what was done in the proof of Lemma 4.6 we arrive at,

$$ \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 \le 2 \sum_{i=0}^{N-1} \frac{1+2c_r}{\tilde{h}_i h_i} \|v_h\|_{L^2([x_i, x_{i+1}])}^2 + 2 \sum_{i=1}^{N} \frac{1+2c_r}{\tilde{h}_i h_{i-1}} \|v_h\|_{L^2([x_{i-1}, x_i])}^2. $$

Since the subdivision of the domain is quasi-uniform we have $h_i \ge \rho h$ hence,

$$ \sum_{i=1}^{N-1} \frac{1}{\tilde{h}_i} \|v_h(x_i)\|^2 \le 4 \frac{1+2c_r}{\rho^2 h^2} \|v_h\|_{L^2(0,1)}^2, $$

this shows that $\|C_h\|_h$ is bounded by $O(1/h)$. $\square$

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

# 7 Numerical results

In a first part, we perform a numerical eigenvalue study where we use our stabilization terms with respect to the ones classically proposed in the literature. In a second part, we introduce a fully explicit leap-frog scheme and provide convergence plots with respect to $h$ as well as energy plots over time. The results obtained illustrate the efficiency of the method. Note that in this section we use the simplified expression for the stabilization term discussed in Remark 4.13 and denote by $\tilde{D}_h$ the operator associated with the bilinear form $\tilde{d}_h$ defined by (54). In the following results, the domain $(0, 1)$ is evenly subdivided, i.e. $h = h_i = \tilde{h}_i$.

## 7.1 Eigenvalues investigation

A first idea is to stabilize the semi-discrete problem (23) is to stabilize the jumps of the discontinuous functions, as seen in [25] for Maxwell's equations. We define

$$ A_{h,\beta} = A_h - \begin{pmatrix} 0 & 0 \\ 0 & B_h \end{pmatrix} $$

and define the operator

$$ \forall z_h = \begin{pmatrix} u_h \\ v_h \end{pmatrix} \in \mathcal{Z}_h, \quad \tilde{z}_h = \begin{pmatrix} \tilde{u}_h \\ \tilde{v}_h \end{pmatrix} \in \mathcal{Z}_h, \quad (J_h z_h, \tilde{z}_h)_h = \frac{1}{h} \sum_{i=1}^{N} \|v_h(x_i)\| \|\tilde{v}_h(x_i)\| $$
in the dynamics (23) which corresponds to a stabilization of the jumps of $v_h$ that are discontinuous. Namely, we represent in Figure 3 the spectrum of the operator associated with the following equation ($\gamma = 1$),

$$ \left\{ \begin{array}{l} \oint_0^1 \partial_t u_h \tilde{u}_h \, dx - \int_0^1 v_h \partial_x \tilde{u}_h \, dx + v_h(1) \tilde{u}_h(1) = 0, \\ \\ \oint_0^1 \partial_t v_h \tilde{v}_h \, dx + \int_0^1 \partial_x u_h \tilde{v}_h \, dx - u_h(1) \tilde{v}_h(1) + v_h(1) \tilde{v}_h(1) \\ \qquad + \frac{1}{h} \sum_{i=1}^N \llbracket v_h(x_i) \rrbracket \llbracket \tilde{v}_h(x_i) \rrbracket = 0. \end{array} \right. $$

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

<figure><img src="image_7.png" /><figcaption>Figure 3. Spectrum of $A_{h,b} - J_h$.</figcaption></figure>

We have shown in Figure 3 (left) some eigenvalues of the continuous operator (they all have a real part equal to -0.549) and the spectrum of the operator $A_{h,b} - J_h$ for two $\mathbb{P}_1$ discretizations. The eigenvalues of the operator $A_h - J_h$ all have a strictly negative real part, so that the dynamics is dissipative but suboptimal, since the real parts of many eigenvalues are too closed to the imaginary axis. Figure 3 (right) shows the same kind of results, but for higher order finite elements. The situation is even worse as some eigenvalues are very close to the imaginary axis. We also represent in Figure 4 the spectrum obtained when viscosity is added. Such a strategy has been studied in [21]. We introduce the operator

$$ \forall z_h = \begin{pmatrix} u_h \\ v_h \end{pmatrix} \in \mathcal{Z}_h, \quad \tilde{z}_h = \begin{pmatrix} \tilde{u}_h \\ \tilde{v}_h \end{pmatrix} \in \mathcal{Z}_h, \quad (V_h z_h, \tilde{z}_h)_h = \int_0^1 \partial_x u_h \partial_x \tilde{u}_h dx, $$

and represents the spectrum of the operator $A_{h,b} - V_h$ associated with the following variational formulations:

$$ \left\{ \begin{array}{l} \oint_0^1 \partial_t u_h \tilde{u}_h \, dx - \int_0^1 v_h \partial_x \tilde{u}_h \, dx + v_h(1) \tilde{u}_h(1) + h^k \int_0^1 \partial_x u_h \partial_x \tilde{u}_h \, dx = 0, \\ \\ \oint_0^1 \partial_t v_h \tilde{v}_h \, dx + \int_0^1 \partial_x u_h \tilde{v}_h \, dx - u_h(1) \tilde{v}_h(1) + v_h(1) \tilde{v}_h(1) \\ \qquad + \frac{1}{h} \sum_{i=1}^N \llbracket v_h(x_i) \rrbracket \llbracket \tilde{v}_h(x_i) \rrbracket = 0, \end{array} \right. \quad (70) $$

where $k \in \mathbb{N}^*$ should be chosen large to preserve consistency.

Thus, as we can see in Figure 4 (left), a second order perturbation gives the exponential stability for $\mathbb{P}_1$ finite elements. For high-order finite elements the consistency is deteriorated. Furthermore, If we soften the impact of the viscous term – Figure 4 (right) – for instance by choosing $k = r + 1$ where r is the order of the finite element method then some eigenvalues are close to the imaginary axis (however consistency is recovered). Finally we introduce the operator

$$ D_h = \begin{pmatrix} D_h & 0 \\ 0 & \tilde{D}_h \end{pmatrix}, $$
in the dynamics (23) which corresponds to a stabilization of the jumps of $v_h$ that are discontinuous. Namely, we represent in Figure 3 the spectrum of the operator associated with the following equation ($\gamma = 1$),

$$ \left\{ \begin{array}{l} \oint_0^1 \partial_t u_h \tilde{u}_h \, dx - \int_0^1 v_h \partial_x \tilde{u}_h \, dx + v_h(1) \tilde{u}_h(1) = 0, \\ \oint_0^1 \partial_t v_h \tilde{v}_h \, dx + \int_0^1 \partial_x u_h \tilde{v}_h \, dx - u_h(1) \tilde{v}_h(1) + v_h(1) \tilde{v}_h(1) \\ \qquad + \frac{1}{h} \sum_{i=1}^N \llbracket v_h(x_i) \rrbracket \llbracket \tilde{v}_h(x_i) \rrbracket = 0. \end{array} \right. $$

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

<figure><img src="image_10.png" /><figcaption>Figure 3. Spectrum of $A_{h,b} - J_h$.</figcaption></figure>

We have shown in Figure 3 (left) some eigenvalues of the continuous operator (they all have a real part equal to -0.549) and the spectrum of the operator $A_{h,b} - J_h$ for two $\mathbb{P}_1$ discretizations. The eigenvalues of the operator $A_h - J_h$ all have a strictly negative real part, so that the dynamics is dissipative but suboptimal, since the real parts of many eigenvalues are too closed to the imaginary axis. Figure 3 (right) shows the same kind of results, but for higher order finite elements. The situation is even worse as some eigenvalues are very close to the imaginary axis. We also represent in Figure 4 the spectrum obtained when viscosity is added. Such a strategy has been studied in [21]. We introduce the operator

$$ \forall z_h = \begin{pmatrix} u_h \\ v_h \end{pmatrix} \in \mathcal{Z}_h, \quad \tilde{z}_h = \begin{pmatrix} \tilde{u}_h \\ \tilde{v}_h \end{pmatrix} \in \mathcal{Z}_h, \quad (V_h z_h, \tilde{z}_h)_h = \int_0^1 \partial_x u_h \partial_x \tilde{u}_h dx, $$

and represents the spectrum of the operator $A_{h,b} - V_h$ associated with the following variational formulations:

$$ \left\{ \begin{array}{l} \oint_0^1 \partial_t u_h \tilde{u}_h \, dx - \int_0^1 v_h \partial_x \tilde{u}_h \, dx + v_h(1) \tilde{u}_h(1) + h^k \int_0^1 \partial_x u_h \partial_x \tilde{u}_h \, dx = 0, \\ \oint_0^1 \partial_t v_h \tilde{v}_h \, dx + \int_0^1 \partial_x u_h \tilde{v}_h \, dx - u_h(1) \tilde{v}_h(1) + v_h(1) \tilde{v}_h(1) \\ \qquad + \frac{1}{h} \sum_{i=1}^N \llbracket v_h(x_i) \rrbracket \llbracket \tilde{v}_h(x_i) \rrbracket = 0, \end{array} \right. \quad (70) $$

where $k \in \mathbb{N}^*$ should be chosen large to preserve consistency.

Thus, as we can see in Figure 4 (left), a second order perturbation gives the exponential stability for $\mathbb{P}_1$ finite elements. For high-order finite elements the consistency is deteriorated. Furthermore, If we soften the impact of the viscous term – Figure 4 (right) – for instance by choosing $k = r + 1$ where r is the order of the finite element method then some eigenvalues are close to the imaginary axis (however consistency is recovered). Finally we introduce the operator

$$ D_h = \begin{pmatrix} D_h & 0 \\ 0 & \tilde{D}_h \end{pmatrix}, $$
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Figure 4. Spectrum of $A_{h,\blacksquare} - J_h - h^2V_h$ (left) and of $A_h - J_h - h^{r+1}V_h$ (right), r is the order of the finite element method.</figcaption></figure>

and represents Figure 5 the eigenvalues of the operator $A_h - J_h - D_h$. One can notice that consistency and exponential stability are preserved.

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

<figure><img src="image_11.png" /><figcaption>Figure 5. Spectrum of $A_h - J_h - D_h$</figcaption></figure>

## 7.2 A numerical example with an explicit leap-frog scheme

One of our objectives being to propse an effective numerical strategy, we now present the leap-frog scheme used in the numerical simulations. We consider the damped wave equation (47) discretized with an explicit scheme where the dissipation terms are off-centered.

$$ \left\{ \begin{array}{l} \displaystyle \frac{u_h^{n+1} - u_h^n}{\Delta t} - R_h^* v_h^{n+1/2} + \gamma D_h u_h^n = 0 \\[2ex] \displaystyle \frac{v_h^{n+3/2} - v_h^{n+1/2}}{\Delta t} - R_h u_h^{n+1} + \gamma B_h v_h^{n+1} \\ \qquad + \gamma J_h \frac{v_h^{n+3/2} + v_h^{n+1/2}}{2} + \gamma \tilde{D}_h v_h^{n+1/2} = 0. \end{array} \right. \qquad (71) $$

Note that when no-stabilization is used $J_h$, $D_h$ and $\tilde{D}_h$ are neglected. In general, the scheme is “almost” fully explicit thanks to the mass-lumping strategy and off-centering of the newly added dissipative terms $D_h$ and $\tilde{D}_h$. Note that:

* centering of the terms involving $J_h$, implies at the algebraic level, to solve $(N-1)$ linear systems of size $2 \times 2$ at each time iteration (N is the number of sub-intervals);

* off-centering does not imply here formally a loss of consistency in time since it is used on terms that are small when apply to interpolation of smooth solutions of the problem.
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 4. Spectrum of $A_{h,\blacksquare} - J_h - h^2V_h$ (left) and of $A_h - J_h - h^{r+1}V_h$ (right), r is the order of the finite element method.</figcaption></figure>

and represents Figure 5 the eigenvalues of the operator $A_h - J_h - D_h$. One can notice that consistency and exponential stability are preserved.

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

<figure><img src="image_14.png" /><figcaption>Figure 5. Spectrum of $A_h - J_h - D_h$</figcaption></figure>

## 7.2 A numerical example with an explicit leap-frog scheme

One of our objectives being to propse an effective numerical strategy, we now present the leap-frog scheme used in the numerical simulations. We consider the damped wave equation (47) discretized with an explicit scheme where the dissipation terms are off-centered.

$$ \left\{ \begin{array}{l} \displaystyle \frac{u_h^{n+1} - u_h^n}{\Delta t} - R_h^* v_h^{n+1/2} + \gamma D_h u_h^n = 0 \\[1.5em] \displaystyle \frac{v_h^{n+3/2} - v_h^{n+1/2}}{\Delta t} - R_h u_h^{n+1} + \gamma B_h v_h^{n+1} \\ \qquad + \gamma J_h \frac{v_h^{n+3/2} + v_h^{n+1/2}}{2} + \gamma \tilde{D}_h v_h^{n+1/2} = 0. \end{array} \right. \qquad (71) $$

Note that when no-stabilization is used $J_h$, $D_h$ and $\tilde{D}_h$ are neglected. In general, the scheme is “almost” fully explicit thanks to the mass-lumping strategy and off-centering of the newly added dissipative terms $D_h$ and $\tilde{D}_h$. Note that:

* centering of the terms involving $J_h$, implies at the algebraic level, to solve $(N-1)$ linear systems of size $2 \times 2$ at each time iteration (N is the number of sub-intervals);

* off-centering does not imply here formally a loss of consistency in time since it is used on terms that are small when apply to interpolation of smooth solutions of the problem.
The scheme preserves an energy identity. To show this energy identity, we rewrite the scheme using the following trick

$$ \begin{cases} D_h u_h^n = D_h u_h^{n+1/2} - \frac{\Delta t}{2} D_h \frac{u_h^{n+1} - u_h^n}{\Delta t}, \\ \tilde{D}_h v_h^{n+1/2} = \tilde{D}_h v_h^{n+1} - \frac{\Delta t}{2} \tilde{D}_h \frac{v_h^{n+3/2} - v_h^{n+1/2}}{\Delta t}, \end{cases} $$

thus, we obtain

$$ \left\{ \begin{aligned} & \widehat{M}_h \frac{u_h^{n+1} - u_h^n}{\Delta t} - R_h^* v_h^{n+1/2} + \gamma D_h u_h^{n+1/2} = 0, \\ & \widetilde{M}_h \frac{v_h^{n+3/2} - v_h^{n+1/2}}{\Delta t} + R_h u_h^{n+1} + \gamma B_h v_h^{n+1} + \gamma J_h v_h^{n+1} + \gamma \tilde{D}_h v_h^{n+1} = 0, \end{aligned} \right. \quad (72) $$

where

$$ \widehat{M}_h = I_h - \gamma \frac{\Delta t}{2} D_h \quad \text{and} \quad \widetilde{M}_h = I_h - \gamma \frac{\Delta t}{2} \tilde{D}_h, $$

where $I_h$ are (different) identity operators in $U_h^r$ and $V_h^r$. The energy identity for System (71) involves a new energy functional that is given for $u_h^n \in U_h^r$ and $v_h^{n+1/2} \in V_h^r$ by

$$ \mathcal{E}_{h,ex}^n = \frac{1}{2} \left( \widehat{M}_h u_h^n, u_h^n \right)_h + \frac{1}{2} \left( \widetilde{M}_h v_h^{n+1/2}, v_h^{n+1/2} \right)_h - \frac{\Delta t}{2} \left( R_h u_h^n, v_h^{n+1/2} \right)_h. \quad (73) $$

**Theorem 7.1.** Let $\{u_h^n, v_h^n\}_{n \ge 0}$ belonging to $U_h^r \times V_h^r$ be the solutions to (72). The following energy identity holds,

$$ \frac{\mathcal{E}_{h,ex}^{n+1} - \mathcal{E}_{h,ex}^n}{\Delta t} = -\gamma (B_h v_h^{n+1}, v_h^{n+1})_h - \gamma (J_h v_h^{n+1}, v_h^{n+1})_h \\ - \gamma (\tilde{D}_h v_h^{n+1}, v_h^{n+1})_h - \gamma (D_h u_h^{n+1/2}, u_h^{n+1/2})_h. $$

Assume that there exists $\eta > 0$ such that

$$ \gamma \frac{\Delta t}{2} \|D_h\|_h \le \eta < 1 \quad \text{and} \quad \gamma \frac{\Delta t}{2} \|\tilde{D}_h\|_h \le \eta < 1 - \eta, \quad (74) $$

then, when

$$ \frac{\Delta t}{2} \|R_h\|_h < 1 - \eta, \quad (75) $$

then, for all $n \ge 0$ we have $\mathcal{E}_{h,ex}^n \ge 0$ and the scheme is stable.

**Proof.** First we prove the energy identity. We multiply the two equations of (72) by $u_h^{n+1/2}$ and $v_h^{n+1}$ respectively, and we obtain,

$$ \left( \widehat{M}_h \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} + u_h^n}{2} \right)_h - \left( R_h^* v_h^{n+1/2}, \frac{u_h^{n+1} + u_h^n}{2} \right)_h + \gamma \left( D_h u_h^{n+1/2}, u_h^{n+1/2} \right)_h \\ + \left( \widetilde{M}_h \frac{v_h^{n+3/2} - v_h^{n+1/2}}{\Delta t}, \frac{v_h^{n+3/2} + v_h^{n+1/2}}{2} \right)_h + \left( R_h u_h^{n+1}, \frac{v_h^{n+3/2} + v_h^{n+1/2}}{2} \right)_h \\ + \gamma (B_h v_h^{n+1}, v_h^{n+1})_h + \gamma (J_h v_h^{n+1}, v_h^{n+1})_h + \gamma (\tilde{D}_h v_h^{n+1}, v_h^{n+1})_h = 0. $$

Then we use the following identity

$$ \begin{aligned} & - \left( R_h^* v_h^{n+1/2}, \frac{u_h^{n+1} + u_h^n}{2} \right)_h + \left( R_h u_h^{n+1}, \frac{v_h^{n+3/2} + v_h^{n+1/2}}{2} \right)_h \\ & \qquad = -\frac{1}{2} \left( R_h^* v_h^{n+1/2}, u_h^n \right)_h + \frac{1}{2} \left( R_h^* v_h^{n+3/2}, u_h^{n+1} \right)_h, \end{aligned} $$
The scheme preserves an energy identity. To show this energy identity, we rewrite the scheme using the following trick

$$ egin{cases} D_h u_h^n = D_h u_h^{n+1/2} - \frac{\Delta t}{2} D_h \frac{u_h^{n+1} - u_h^n}{\Delta t}, \\ \tilde{D}_h v_h^{n+1/2} = \tilde{D}_h v_h^{n+1} - \frac{\Delta t}{2} \tilde{D}_h \frac{v_h^{n+3/2} - v_h^{n+1/2}}{\Delta t}, \end{cases} $$

thus, we obtain

$$ egin{cases} \widehat{M}_h \frac{u_h^{n+1} - u_h^n}{\Delta t} - R_h^* v_h^{n+1/2} + \gamma D_h u_h^{n+1/2} = 0, \\ \widetilde{M}_h \frac{v_h^{n+3/2} - v_h^{n+1/2}}{\Delta t} + R_h u_h^{n+1} + \gamma B_h v_h^{n+1} + \gamma J_h v_h^{n+1} + \gamma \tilde{D}_h v_h^{n+1} = 0, \end{cases} \quad (72) $$

where

$$ \widehat{M}_h = I_h - \gamma \frac{\Delta t}{2} D_h \quad \text{and} \quad \widetilde{M}_h = I_h - \gamma \frac{\Delta t}{2} \tilde{D}_h, $$

where $I_h$ are (different) identity operators in $U_h^r$ and $V_h^r$. The energy identity for System (71) involves a new energy functional that is given for $u_h^n \in U_h^r$ and $v_h^{n+1/2} \in V_h^r$ by

$$ \mathcal{E}_{h,ex}^n = \frac{1}{2} (\widehat{M}_h u_h^n, u_h^n)_h + \frac{1}{2} (\widetilde{M}_h v_h^{n+1/2}, v_h^{n+1/2})_h - \frac{\Delta t}{2} (R_h u_h^n, v_h^{n+1/2})_h. \quad (73) $$

**Theorem 7.1.** Let $\{u_h^n, v_h^n\}_{n \ge 0}$ belonging to $U_h^r \times V_h^r$ be the solutions to (72). The following energy identity holds,

$$ \frac{\mathcal{E}_{h,ex}^{n+1} - \mathcal{E}_{h,ex}^n}{\Delta t} = -\gamma (B_h v_h^{n+1}, v_h^{n+1})_h - \gamma (J_h v_h^{n+1}, v_h^{n+1})_h \\ - \gamma (\tilde{D}_h v_h^{n+1}, v_h^{n+1})_h - \gamma (D_h u_h^{n+1/2}, u_h^{n+1/2})_h. $$

Assume that there exists $\eta > 0$ such that

$$ \gamma \frac{\Delta t}{2} \|D_h\|_h \le \eta < 1 \quad \text{and} \quad \gamma \frac{\Delta t}{2} \|\tilde{D}_h\|_h \le \eta < 1 - \eta, \quad (74) $$

then, when

$$ \frac{\Delta t}{2} \|R_h\|_h < 1 - \eta, \quad (75) $$

then, for all $n \ge 0$ we have $\mathcal{E}_{h,ex}^n \ge 0$ and the scheme is stable.

**Proof.** First we prove the energy identity. We multiply the two equations of (72) by $u_h^{n+1/2}$ and $v_h^{n+1}$ respectively, and we obtain,

$$ \left( \widehat{M}_h \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} + u_h^n}{2} \right)_h - \left( R_h^* v_h^{n+1/2}, \frac{u_h^{n+1} + u_h^n}{2} \right)_h + \gamma \left( D_h u_h^{n+1/2}, u_h^{n+1/2} \right)_h \\ + \left( \widetilde{M}_h \frac{v_h^{n+3/2} - v_h^{n+1/2}}{\Delta t}, \frac{v_h^{n+3/2} + v_h^{n+1/2}}{2} \right)_h + \left( R_h u_h^{n+1}, \frac{v_h^{n+3/2} + v_h^{n+1/2}}{2} \right)_h \\ + \gamma (B_h v_h^{n+1}, v_h^{n+1})_h + \gamma (J_h v_h^{n+1}, v_h^{n+1})_h + \gamma (\tilde{D}_h v_h^{n+1}, v_h^{n+1})_h = 0. $$

Then we use the following identity

$$ - \left( R_h^* v_h^{n+1/2}, \frac{u_h^{n+1} + u_h^n}{2} \right)_h + \left( R_h u_h^{n+1}, \frac{v_h^{n+3/2} + v_h^{n+1/2}}{2} \right)_h \\ = -\frac{1}{2} \left( R_h^* v_h^{n+1/2}, u_h^n \right)_h + \frac{1}{2} \left( R_h^* v_h^{n+3/2}, u_h^{n+1} \right)_h, $$
to obtain

$$ 
\begin{aligned}
& \frac{1}{2\Delta t} (\widehat{M}_h u_h^{n+1}, u_h^{n+1})_h - \frac{1}{2\Delta t} (\widehat{M}_h u_h^n, u_h^n)_h + \frac{1}{2\Delta t} (\widetilde{M}_h v_h^{n+3/2}, v_h^{n+3/2})_h \\
& - \frac{1}{2\Delta t} (\widetilde{M}_h v_h^{n+1/2}, v_h^{n+1/2})_h + \frac{\Delta t}{2} (R_h u_h^{n+1}, v_h^{n+3/2})_h - \frac{\Delta t}{2} (R_h u_h^n, v_h^{n+1/2})_h \\
& = -\gamma (B_h v_h^{n+1}, v_h^{n+1})_h + \gamma (J_h v_h^{n+1}, v_h^{n+1})_h \\
& \quad - \gamma (\widetilde{D}_h v_h^{n+1}, v_h^{n+1})_h - \gamma (D_h u_h^{n+1/2}, u_h^{n+1/2})_h.
\end{aligned}
$$

We recognize the expression of the explicit energy (73). To prove the positivity of the energy, we use that

$$ |(R_h u_h^n, v_h^{n+1/2})_h| \le \|R_h\|_h \|u_h^n\|_h \|v_h^{n+1/2}\|_h \le \frac{\|R_h\|_h}{2} (\|u_h^n\|_h^2 + \|v_h^{n+1/2}\|_h^2) $$

and

$$ (\widehat{M}_h u_h^n, u_h^n)_h \ge (1-\eta) \|u_h^n\|_h^2, \quad (\widetilde{M}_h v_h^{n+1/2}, v_h^{n+1/2})_h \ge (1-\eta) \|v_h^{n+1/2}\|_h^2 $$

Injecting these inequalities in (73) leads to

$$ \mathcal{E}_{h,ex}^n \ge \frac{1}{2} (1 - \eta - \Delta t \frac{\|R_h\|_h}{2}) (\|u_h^n\|^2 + \|v_h^{n+1/2}\|^2) $$

the positiveness follows from the conditions (74) and (75). The stability of the scheme is a consequence of the property that if (75) holds, the energy is in fact a norm (uniform in h) for $(u_h^n, v_h^{n+1/2})^t$. □

**Remark 7.2.** Note that since $D_h$ and $\widetilde{D}_h$ are bounded with bounds independent of h, the condition (74) is not restrictive. Condition (75) corresponds to a classical CFL condition for explicit schemes.

As in Section 2, we now solve the transient problem with the explicit scheme given above (with and without stabilization) and look at the evolution of the solution over time and its energy. The numerical pollution due to spurious waves that can be observed for the non-stabilized numerical experiment is not observed in the stabilized version, see Figure 6.

<figure><img src="image_12.png" /><figcaption>**Figure 6.** Snapshot of the solution $u_h^n$ discretized with $\mathbb{P}_4$ F.E.M. with $h = 1/10$ and $dt = 1/1000$ around time $t = 2$ (left) and at time $t = 4$ (right) and computed with and without stabilization</figcaption></figure>

While the system without the stabilization terms had an energy that reaches a plateau asymptotically in time (Figure 7), we notice that with the addition of the
to obtain

$$ 
\begin{aligned}
& \frac{1}{2\Delta t} (\widehat{M}_h u_h^{n+1}, u_h^{n+1})_h - \frac{1}{2\Delta t} (\widehat{M}_h u_h^n, u_h^n)_h + \frac{1}{2\Delta t} (\widetilde{M}_h v_h^{n+3/2}, v_h^{n+3/2})_h \\
& \quad - \frac{1}{2\Delta t} (\widetilde{M}_h v_h^{n+1/2}, v_h^{n+1/2})_h + \frac{\Delta t}{2} (R_h u_h^{n+1}, v_h^{n+3/2})_h - \frac{\Delta t}{2} (R_h u_h^n, v_h^{n+1/2})_h \\
& \quad = -\gamma (B_h v_h^{n+1}, v_h^{n+1})_h + \gamma (J_h v_h^{n+1}, v_h^{n+1})_h \\
& \quad \quad - \gamma (\widetilde{D}_h v_h^{n+1}, v_h^{n+1})_h - \gamma (D_h u_h^{n+1/2}, u_h^{n+1/2})_h.
\end{aligned}
$$

We recognize the expression of the explicit energy (73). To prove the positivity of the energy, we use that

$$ |(R_h u_h^n, v_h^{n+1/2})_h| \le \|R_h\|_h \|u_h^n\|_h \|v_h^{n+1/2}\|_h \le \frac{\|R_h\|_h}{2} (\|u_h^n\|_h^2 + \|v_h^{n+1/2}\|_h^2) $$

and

$$ (\widehat{M}_h u_h^n, u_h^n)_h \ge (1-\eta) \|u_h^n\|_h^2, \quad (\widetilde{M}_h v_h^{n+1/2}, v_h^{n+1/2})_h \ge (1-\eta) \|v_h^{n+1/2}\|_h^2 $$

Injecting these inequalities in (73) leads to

$$ \mathcal{E}_{h,ex}^n \ge \frac{1}{2} (1 - \eta - \Delta t \frac{\|R_h\|_h}{2}) (\|u_h^n\|^2 + \|v_h^{n+1/2}\|^2) $$

the positiveness follows from the conditions (74) and (75). The stability of the scheme is a consequence of the property that if (75) holds, the energy is in fact a norm (uniform in h) for $(u_h^n, v_h^{n+1/2})^t$. □

**Remark 7.2.** Note that since $D_h$ and $\widetilde{D}_h$ are bounded with bounds independent of h, the condition (74) is not restrictive. Condition (75) corresponds to a classical CFL condition for explicit schemes.

As in Section 2, we now solve the transient problem with the explicit scheme given above (with and without stabilization) and look at the evolution of the solution over time and its energy. The numerical pollution due to spurious waves that can be observed for the non-stabilized numerical experiment is not observed in the stabilized version, see Figure 6.

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

<figure><img src="image_16.png" /><figcaption>**Figure 6.** Snapshot of the solution $u_h^n$ discretized with $\mathbb{P}_4$ F.E.M. with $h = 1/10$ and $dt = 1/1000$ around time $t = 2$ (left) and at time $t = 4$ (right) and computed with and without stabilization</figcaption></figure>

While the system without the stabilization terms had an energy that reaches a plateau asymptotically in time (Figure 7), we notice that with the addition of the
<figure><img src="image_13.png" /><figcaption>**Figure 7.** Plot of $\ln (\|u_h^n\|_h^2 + \|v_h^n\|_h^2)^{1/2}$ with respect to time $t^n = n\Delta t$.</figcaption></figure>

stabilization terms, we obtain an exponentially decreasing energy as can be seen in Figure 7. This confirms numerically illustrate the result of Theorem 6.6.

Furthermore, Figure 8 shows numerically that the order of the scheme is not affected by the addition of the stabilization terms. We plot the relative difference in the norm $L^\infty((0,T), L^2(0,1))$ with $T = 1$ between the fully discrete solution $u_h^n$ and the interpolation of the exact solution $u(t)$ as a function of different space steps $h$. We observe a convergence of order $r+1$ – where $r$ is the order of the finite elements used – until a plateau is reached. This is explained by the fact that the explicit scheme is of lower order in time. Nevertheless, we can conclude that the additional stabilization terms do not affect the accuracy of the high-order scheme.

<figure><img src="image_14.png" /><figcaption>**Figure 8.** Evolution of the error according of the space steps $h$ with finite elements of order 3 (left) and order 4 (right)</figcaption></figure>

## References

[1] J.-L. Lions. *Contrôlabilité Exacte, Perturbations et Stabilisation de Systèmes Distribués*. Tome 1. Vol. 8. Recherches En Mathématiques Appliquées [Research in Applied Mathematics]. Masson, Paris, 1988, pp. x+541.

[2] R. Glowinski, C.-H. Li, and J.-L. Lions. "A Numerical Approach to the Exact Boundary Controllability of the Wave Equation (I) Dirichlet Controls: Description of the Numerical Methods". In: *Japan Journal of Applied Mathematics* 7 (1990), pp. 1–76.

[3] Enrike Zuazua. "Uniform stabilization of the wave equation by nonlinear boundary feedback". In: *SIAM Journal on Control and Optimization* 28.2 (1990), pp. 466–477.

[4] G. Lebeau, C. Bardos, and J. Rauch. "Sharp Sufficient Conditions for the Observation, Control, and Stabilization of Waves from the Boundary". In: *SIAM Journal on Control and Optimization* 30.5 (1992), pp. 1024–1065.
<figure><img src="image_17.png" /><figcaption>**Figure 7.** Plot of $\ln (\|u_h^n\|_h^2 + \|v_h^n\|_h^2)^{1/2}$ with respect to time $t^n = n\Delta t$.</figcaption></figure>

stabilization terms, we obtain an exponentially decreasing energy as can be seen in Figure 7. This confirms numerically illustrate the result of Theorem 6.6.

Furthermore, Figure 8 shows numerically that the order of the scheme is not affected by the addition of the stabilization terms. We plot the relative difference in the norm $L^\infty((0,T), L^2(0,1))$ with $T = 1$ between the fully discrete solution $u_h^n$ and the interpolation of the exact solution $u(t)$ as a function of different space steps $h$. We observe a convergence of order $r+1$ – where $r$ is the order of the finite elements used – until a plateau is reached. This is explained by the fact that the explicit scheme is of lower order in time. Nevertheless, we can conclude that the additional stabilization terms do not affect the accuracy of the high-order scheme.

<figure><img src="image_18.png" /><figcaption>**Figure 8.** Evolution of the error according of the space steps $h$ with finite elements of order 3 (left) and order 4 (right)</figcaption></figure>

## References

[1] J.-L. Lions. *Contrôlabilité Exacte, Perturbations et Stabilisation de Systèmes Distribués*. Tome 1. Vol. 8. Recherches En Mathématiques Appliquées [Research in Applied Mathematics]. Masson, Paris, 1988, pp. x+541.

[2] R. Glowinski, C.-H. Li, and J.-L. Lions. "A Numerical Approach to the Exact Boundary Controllability of the Wave Equation (I) Dirichlet Controls: Description of the Numerical Methods". In: *Japan Journal of Applied Mathematics* 7 (1990), pp. 1-76.

[3] Enrike Zuazua. "Uniform stabilization of the wave equation by nonlinear boundary feedback". In: *SIAM Journal on Control and Optimization* 28.2 (1990), pp. 466-477.

[4] G. Lebeau, C. Bardos, and J. Rauch. "Sharp Sufficient Conditions for the Observation, Control, and Stabilization of Waves from the Boundary". In: *SIAM Journal on Control and Optimization* 30.5 (1992), pp. 1024-1065.
[5] S. Cox and E. Zuazua. “The Rate At Which Energy Decays In A String Damped At One End”. In: *Communications in Partial Differential Equations* 19 (1994), pp. 213–243.

[6] V. Komornik. “Exact Controllability and Stabilization. The Multiplier Method”. In: (1997), p. 161.

[7] A. J. Krener. “A Lyapunov of nonlinear observers”. In: (1998).

[8] J. A. Infante and E. Zuazua. “Boundary Observability for the Space Semi-Discretizations of the 1-d Wave Equation”. In: *ESAIM: Mathematical Modelling and Numerical Analysis* 33.2 (1999), pp. 407–438.

[9] P. G. Ciarlet. *The finite element method for elliptic problems*. SIAM, 2002.

[10] G. C. Cohen. *Higher-order numerical methods for transient wave equations*. Vol. 5. Springer, 2002.

[11] S. Micu. “Uniform Boundary Controllability of a Semi-Discrete 1-D Wave Equation”. In: *Numerische Mathematik* 91 (2002), pp. 723–768.

[12] M. Negreanu and E. Zuazua. “Uniform Boundary Controllability of a Discrete 1-D Wave Equation”. In: *Systems & Control Letters* 48.3-4 (2003), pp. 261–279.

[13] I. Lasiecka and R. Triggiani. “The operator B*L for the wave equation with Dirichlet control”. In: *Abstract and Applied Analysis* 2004.7 (2004), pp. 625–634.

[14] M. Negreanu and E. Zuazua. “Convergence of a multigrid method for the controllability of a 1-d wave equation”. In: *Comptes Rendus Mathematique* 338.5 (2004), pp. 413–418.

[15] A. Münch. “A uniformly controllable and implicit scheme for the 1-D wave equation”. In: *ESAIM: Mathematical Modelling and Numerical Analysis* 39.2 (2005), pp. 377–418.

[16] E. Zuazua. “Propagation, observation, control and numerical approximation of waves”. In: *SIAM Review* 47.2 (2005), pp. 197–243.

[17] C. Castro and S. Micu. “Boundary Controllability of a Linear Semi-Discrete 1-D Wave Equation Derived from a Mixed Finite Element Method”. In: *Numerische Mathematik* 102.3 (2006), pp. 413–462.

[18] R. Joly. “Convergence of the wave equation damped on the interior to the one damped on the boundary”. In: *Journal of Differential Equations* 229.2 (2006), pp. 588–653.

[19] P. Joly. “Numerical Methods for Elastic Wave Propagation”. In: *Waves in nonlinear pre-stressed materials* (2007), pp. 181–281.

[20] K. Ramdani, T Takahashi, and M. Tucsnak. “Uniformly exponentially stable approximations for a class of second order evolution equations—application to LQR problems”. In: *ESAIM. Control, Optimisation and Calculus of Variations. European Series in Applied and Industrial Mathematics* 13.3 (2007), pp. 503–527.

[21] L. Tebou and E. Zuazua. “Uniform Boundary Stabilization of the Finite Difference Space Discretization of the 1-d wave equation”. In: 26 (2007), pp. 337–365.

[22] S. C. Brenner and L. R. Scott. *The Mathematical Theory of Finite Element Methods*. Red. by J. E. Marsden, L. Sirovich, and S. S. Antman. Vol. 15. Texts in Applied Mathematics. New York, NY: Springer, 2008.

[23] S. Ervedoza and E. Zuazua. “Perfectly matched layers in 1-d: energy decay for continuous and semi-discrete waves”. In: *Numerische Mathematik* 109.4 (2008), pp. 597–634.
[5] S. Cox and E. Zuazua. "The Rate At Which Energy Decays In A String Damped At One End". In: *Communications in Partial Differential Equations* 19 (1994), pp. 213-243.

[6] V. Komornik. "Exact Controllability and Stabilization. The Multiplier Method". In: (1997), p. 161.

[7] A. J. Krener. "A Lyapunov of nonlinear observers". In: (1998).

[8] J. A. Infante and E. Zuazua. "Boundary Observability for the Space Semi-Discretizations of the 1-d Wave Equation". In: *ESAIM: Mathematical Modelling and Numerical Analysis* 33.2 (1999), pp. 407-438.

[9] P. G. Ciarlet. *The finite element method for elliptic problems*. SIAM, 2002.

[10] G. C. Cohen. *Higher-order numerical methods for transient wave equations*. Vol. 5. Springer, 2002.

[11] S. Micu. "Uniform Boundary Controllability of a Semi-Discrete 1-D Wave Equation". In: *Numerische Mathematik* 91 (2002), pp. 723-768.

[12] M. Negreanu and E. Zuazua. "Uniform Boundary Controllability of a Discrete 1-D Wave Equation". In: *Systems & Control Letters* 48.3-4 (2003), pp. 261-279.

[13] I. Lasiecka and R. Triggiani. "The operator B*L for the wave equation with Dirichlet control". In: *Abstract and Applied Analysis* 2004.7 (2004), pp. 625-634.

[14] M. Negreanu and E. Zuazua. "Convergence of a multigrid method for the controllability of a 1-d wave equation". In: *Comptes Rendus Mathematique* 338.5 (2004), pp. 413-418.

[15] A. Münch. "A uniformly controllable and implicit scheme for the 1-D wave equation". In: *ESAIM: Mathematical Modelling and Numerical Analysis* 39.2 (2005), pp. 377-418.

[16] E. Zuazua. "Propagation, observation, control and numerical approximation of waves". In: *SIAM Review* 47.2 (2005), pp. 197-243.

[17] C. Castro and S. Micu. "Boundary Controllability of a Linear Semi-Discrete 1-D Wave Equation Derived from a Mixed Finite Element Method". In: *Numerische Mathematik* 102.3 (2006), pp. 413-462.

[18] R. Joly. "Convergence of the wave equation damped on the interior to the one damped on the boundary". In: *Journal of Differential Equations* 229.2 (2006), pp. 588-653.

[19] P. Joly. "Numerical Methods for Elastic Wave Propagation". In: *Waves in nonlinear pre-stressed materials* (2007), pp. 181-281.

[20] K. Ramdani, T Takahashi, and M. Tucsnak. "Uniformly exponentially stable approximations for a class of second order evolution equations—application to LQR problems". In: *ESAIM. Control, Optimisation and Calculus of Variations. European Series in Applied and Industrial Mathematics* 13.3 (2007), pp. 503-527.

[21] L. Tebou and E. Zuazua. "Uniform Boundary Stabilization of the Finite Difference Space Discretization of the 1-d wave equation". In: 26 (2007), pp. 337-365.

[22] S. C. Brenner and L. R. Scott. *The Mathematical Theory of Finite Element Methods*. Red. by J. E. Marsden, L. Sirovich, and S. S. Antman. Vol. 15. Texts in Applied Mathematics. New York, NY: Springer, 2008.

[23] S. Ervedoza and E. Zuazua. "Perfectly matched layers in 1-d: energy decay for continuous and semi-discrete waves". In: *Numerische Mathematik* 109.4 (2008), pp. 597-634.
[24] L. Tebou. "A Carleman estimates based approach for the stabilization of some locally damped semilinear hyperbolic equations". In: *ESAIM: COCV* 14.3 (2008), pp. 561-574.

[25] G. Cohen and A. Sinding. "Non Spurious Mixed Spectral Element Methods for Maxwell's Equations". In: *The 9th International Conference on Mathematical and Numerical Aspects of Waves Propagation*. Pau, France, 2009.

[26] D. Komatitsch, D. Michéa, and G. Erlebacher. "Porting a High-Order Finite-Element Earthquake Modeling Application to NVIDIA Graphics Cards Using CUDA". In: *Journal of Parallel and Distributed Computing* 69.5 (2009), pp. 451-460.

[27] S. Ervedoza. "Observability Properties of a Semi-Discrete 1D Wave Equation Derived from a Mixed Finite Element Method on Nonuniform Meshes". In: *ESAIM: Control, Optimisation and Calculus of Variations* 16.2 (2010), pp. 298-326.

[28] S. Ervedoza. "Observability properties of a semi-discrete 1D wave equation derived from a mixed finite element method on nonuniform meshes". In: *ESAIM: Control, Optimisation and Calculus of Variations* 16.2 (2010), pp. 298-326.

[29] D. Komatitsch, G. Erlebacher, D. Göddeke, and D. Michéa. "High-Order Finite-Element Seismic Wave Propagation Modeling with MPI on a Large GPU Cluster". In: *Journal of computational physics* 229.20 (2010), pp. 7692-7714.

[30] K. Ramdani, M. Tucsnak, and G. Weiss. "Recovering the Initial State of an Infinite-Dimensional System Using Observers". In: *Automatica* 46.10 (2010), pp. 1616-1625.

[31] D. Chapelle, N. Cîndea, and P. Moireau. "Improving convergence in numerical analysis using observers—the wave-like equation case". In: *Mathematical Models and Methods in Applied Sciences* 22.12 (2012), p. 1250040.

[32] S. Ervedoza and E. Zuazua. *The Wave Equation: Control and Numerics*. Vol. 2048. Lecture Notes in Mathematics. Berlin, Heidelberg: Springer Berlin Heidelberg, 2012, pp. 245-339.

[33] G. Haine. "Recovering the observable part of the initial data of an infinite-dimensional linear system with skew-adjoint generator". In: *Mathematics of Control Signals and Systems* 26.3 (2014), pp. 435-462.

[34] A. Marica and E. Zuazua. "Propagation of 1D Waves in Regular Discrete Heterogeneous Media: A Wigner Measure Approach". In: *Foundations of Computational Mathematics* 15.6 (2015), pp. 1571-1636.

[35] E. Burman, A. Feizmohammadi, and L. Oksanen. "A finite element data assimilation method for the wave equation". In: *Mathematics of Computation* 89.324 (2020), pp. 1681-1709.
[24] L. Tebou. "A Carleman estimates based approach for the stabilization of some locally damped semilinear hyperbolic equations". In: *ESAIM: COCV* 14.3 (2008), pp. 561-574.

[25] G. Cohen and A. Sinding. "Non Spurious Mixed Spectral Element Methods for Maxwell's Equations". In: *The 9th International Conference on Mathematical and Numerical Aspects of Waves Propagation*. Pau, France, 2009.

[26] D. Komatitsch, D. Michéa, and G. Erlebacher. "Porting a High-Order Finite-Element Earthquake Modeling Application to NVIDIA Graphics Cards Using CUDA". In: *Journal of Parallel and Distributed Computing* 69.5 (2009), pp. 451-460.

[27] S. Ervedoza. "Observability Properties of a Semi-Discrete 1D Wave Equation Derived from a Mixed Finite Element Method on Nonuniform Meshes". In: *ESAIM: Control, Optimisation and Calculus of Variations* 16.2 (2010), pp. 298-326.

[28] S. Ervedoza. "Observability properties of a semi-discrete 1D wave equation derived from a mixed finite element method on nonuniform meshes". In: *ESAIM: Control, Optimisation and Calculus of Variations* 16.2 (2010), pp. 298-326.

[29] D. Komatitsch, G. Erlebacher, D. Göddeke, and D. Michéa. "High-Order Finite-Element Seismic Wave Propagation Modeling with MPI on a Large GPU Cluster". In: *Journal of computational physics* 229.20 (2010), pp. 7692-7714.

[30] K. Ramdani, M. Tucsnak, and G. Weiss. "Recovering the Initial State of an Infinite-Dimensional System Using Observers". In: *Automatica* 46.10 (2010), pp. 1616-1625.

[31] D. Chapelle, N. Cîndea, and P. Moireau. "Improving convergence in numerical analysis using observers—the wave-like equation case". In: *Mathematical Models and Methods in Applied Sciences* 22.12 (2012), p. 1250040.

[32] S. Ervedoza and E. Zuazua. *The Wave Equation: Control and Numerics*. Vol. 2048. Lecture Notes in Mathematics. Berlin, Heidelberg: Springer Berlin Heidelberg, 2012, pp. 245-339.

[33] G. Haine. "Recovering the observable part of the initial data of an infinite-dimensional linear system with skew-adjoint generator". In: *Mathematics of Control Signals and Systems* 26.3 (2014), pp. 435-462.

[34] A. Marica and E. Zuazua. "Propagation of 1D Waves in Regular Discrete Heterogeneous Media: A Wigner Measure Approach". In: *Foundations of Computational Mathematics* 15.6 (2015), pp. 1571-1636.

[35] E. Burman, A. Feizmohammadi, and L. Oksanen. "A finite element data assimilation method for the wave equation". In: *Mathematics of Computation* 89.324 (2020), pp. 1681-1709.