hypothesisHypothesis claimClaim DG methods for the acoustic vibration
problemFelipe Lepe, David Mora, Jesus Vellojin

# Discontinuous Galerkin methods for the acoustic vibration problem††thanks:
Submitted to the editors DATE. The first author was partially supported by
DIUBB through project 2120173 GI/C Universidad del Bío-Bío and by ANID-Chile
through FONDECYT project 11200529 (Chile). The second author was partially
supported by DIUBB through project 2120173 GI/C Universidad del Bío-Bío and
project Centro de Modelamiento Matemático (CMM), FB210005, BASAL funds for
centers of excellence. The second and third author were partially supported by
ANID-Chile through project Anillo of Computational Mathematics for
Desalination Processes ACT210087.

Felipe Lepe GIMNAP-Departamento de Matemática, Universidad del Bío - Bío,
Casilla 5-C, Concepción, Chile<EMAIL_ADDRESS>David Mora GIMNAP-
Departamento de Matemática, Universidad del Bío - Bío, Casilla 5-C,
Concepción, Chile and C$\textrm{I}^{2}$MA, Universidad de Concepción, Chile.
<EMAIL_ADDRESS>Jesus Vellojin GIMNAP-Departamento de Matemática, Universidad
del Bío - Bío, Casilla 5-C, Concepción, Chile<EMAIL_ADDRESS>
###### Abstract

In two and three dimension we analyze discontinuous Galerkin methods for the
acoustic problem. The acoustic fluid that we consider on this paper is
inviscid, leading to a linear eigenvalue problem. The acoustic problem is
written, in first place, in terms of the displacement. Under the approach of
the non-compact operators theory, we prove convergence and error estimates for
the method when the displacement formulation is considered. We analyze the
influence of the stabilization parameter on the computation of the spectrum,
where spurious eigenmodes arise when this parameter is not correctly chosen.
Alternatively we present the formulation depending only on the pressure,
comparing the performance of the DG methods with the pure displacement
formulation. Computationally, we study the influence of the stabilization
parameter on the arising of spurious eigenvalues when the spectrum is
computed. Also, we report tests in two and three dimensions where convergence
rates are reported, together with a comparison between the displacement and
pressure formulations for the proposed DG methods.

##  1 Introduction

In this paper we are interested in the analysis of a discontinuous Galerkin
method (DG) to approximate the natural frequencies of the acoustic problem.
The method is based in a interior penalization strategy, in which we focus our
attention since the aforementioned parameter may introduce spurious
eigenvalues on the computed spectrum.

The DG discretization for eigenvalues is a ongoing research topic where the
pioneer work related to this subject is [1], where under suitable norms that
depend on jumps and averages, the authors proved spectral correctness and
convergence of the DG method for the Laplace eigenvalue problem for both, the
primal and mixed formulations. Also, this paper consider symmetric and
nonsymmetric DG methods. This study inspired a number of works where DG
methods have been considered for other spectral problems [8, 9, 10, 19, 20].

The contribution of the present work is the development of DG methods for the
acoustic problem. This problem is one of the most important on engineering and
physics, since the correct knowledge of some acoustic fluids, allows to design
different structures as containers, motors, vehicles, aircrafts, etc., or
noise reduction devices, necessary on the industry. These applications have
motivated the development of different numerical methods and approaches such
as those of [4, 7, 6, 14, 4, 26, 25, 17]. Some specific applications that
deserve to be mentioned are those related to space programs, where mechanical
effects of shock, acceleration, and acoustic vibrations are taken into account
on the design of durable structures [24]. Also, there is the ultrasonic
welding, where ultrasonic acoustic vibrations are used to generate heat and
pressure to weld two work pieces together [15]. Another use is found in the
acoustic emission detectors, which deals with leaking on pipelines by using
the energy signal produced by acoustic vibrations when the fluid escapes [3].

In particular, we are interested in the interior penalization discontinuous
Galerkin method (IPDG) in order to obtain a suitable alternative to
approximate accurately the eigenfunctions and eigenvalues of the acoustic
vibration problem. In recent years, the IPDG has been applied to different
eigenvalue problems as for instance, the Maxwell spectral problem [9, 8], the
Stokes eigenvalue problem [20], or the elasticity eigenvalue problem [19],
where the accuracy on the approximation of the spectrum is remarkable. The
IPDG method results to be easy to handle from the computational point of view.
This allows the consideration of arbitrary polynomial degrees in two
dimensional (resp. three dimensional) convex (resp. non-convex) domain,
together with single or mixed boundary conditions. However, the main drawback
of the IPDG methods for eigenvalue problems is the correct choice of the
penalization parameter, since an incorrect selection of the parameter may
introduce spurious eigenvalues. This is a inherent concern that comes with
numerical methods that depend on some stabilizations, as for instance the
virtual element method (VEM), where in [4, 23, 21], the authors show in their
numerical tests that spurious eigenvalues arise when the stabilization
parameter changes. The same is observed with IPDG methods such as those of
[18, 20], where an analysis of spurious eigenvalues is performed. This is
important to take into account, since in real applications, the correct
frequencies must be well approximated to avoid undesirable physical
phenomenons, as for instance, the resonance of certain structures, or the
incorrect development of noise insulators.

The available literature for IPDG methods for eigenvalue problems indicates
that the correct choice of the aforementioned stabilization parameters depend
strongly on the configuration of the problem of interest. This implies that,
as the development of the theoretical well posedness of the problem is
analyzed, the domain, physical parameters and boundary conditions, need to be
taken into account. Also, the polynomial degree for the approximation plays an
important role when the stabilization is chosen, since when the polynomial
approximation increases, the pollution of the spectrum vanishes, as the
computational results available on, for instance, [18, 20] suggests.

With these considerations at hand, our aim is to prove that a IPDG approach is
a suitable technique to approximate the natural frequencies of the acoustic
eigenvalue problem, where the numerical method, whether symmetric or
nonsymmetric, approximates accurately the spectrum, and the stabilization
parameter can be chosen with no major problems, according to the needs of the
problem. In this sense, we will apply our IPDG methods for the two classic
formulations of the acoustic problem: the displacement and pressure
formulations, in order to compare the performance of the discontinuous
numerical schemes when different unknowns are considered in the mathematical
formulation.

The paper is organized as follows: in Section 2, we present the acoustic
eigenvalue problem, the displacement formulation and summarize some important
results such as the solution operator, regularity properties and the
corresponding spectral characterization. Section 3 contains the development of
the DG methods, presenting definitions associated to the mesh elements, norms
and technical results that are needed for the analysis. Also, we introduce the
DG discretization for the acoustic problem and the discrete solution operator
that will approximate the continuous one. Also, the spectral analysis is
presented, where approximation results and error estimates for the eigenvalues
and eigenfunctions are presented. In Section 4 we introduce the pressure
formulation and its DG discretizations, which also are symmetric and
nonsymmetric. Finally, in Section 5, we perform a number of numerical tests in
two and three dimensions where we analyze, in one hand, the influence of the
stabilization parameters on the computation of the spectrum, performing a
spurious analysis and, on the other hand, computation of orders of convergence
for the eigenvalues. In this section we also compare the performance of the
displacement and pressure DG formulations with a benchmark test, where the
computational cost is mesaured trough the cpu-time, sparsity and problem size.

###  1.1 Notations

Throughout this work, $\Omega$ is a generic Lipschitz bounded domain of
$\mathbb{R}^{\texttt{d}}$, where $\texttt{d}\in\\{2,3\\}$. For $s\geq 0$,
$\lVert\cdot\rVert_{s,\Omega}$ stands indistinctly for the norm of the
Hilbertian Sobolev spaces ${\mathrm{H}^{s}(\Omega)}$ or
${\mathrm{H}^{s}(\Omega)}^{\texttt{d}}$ with the convention
$\mathrm{H}^{0}(\Omega):=\mathrm{L}^{2}(\Omega)$. If ${X}$ and ${Y}$ are
normed vector spaces, we write ${X}\hookrightarrow{Y}$ to denote that ${X}$ is
continuously embedded in ${Y}$. We denote by ${X}^{\prime}$ and
$\|\cdot\|_{{X}}$ the dual and the norm of ${X}$, respectively. Finally, we
employ $\boldsymbol{0}$ to denote a generic null vector and the relation
$\texttt{a}\lesssim\texttt{b}$ indicates that $\texttt{a}\leq C\texttt{b}$,
with a positive constant $C$ which is independent of a, b, and the size of the
elements in the mesh. The value of $C$ might change at each occurrence. We
remark that we will write the constant $C$ only when is needed.

##  2 The model problem

The acoustic problem reads as follows: Find $\omega\in\mathbb{R}$, the
displacement $\boldsymbol{u}$ and the pressure $p$ on a domain
$\Omega\subset\mathbb{R}^{\texttt{d}}$, such that

(2.1) $\left\\{\begin{array}[]{rcll}\nabla
p-\omega^{2}\rho(x)\boldsymbol{u}&=&\boldsymbol{0}&\text{in}\,\Omega\\\
p+\rho(x)c^{2}\mathop{\mathrm{div}}\nolimits\boldsymbol{u}&=&0&\text{in}\,\Omega\\\
\boldsymbol{u}\cdot\boldsymbol{n}&=&0&\text{on}\,\partial\Omega,\end{array}\right.$

where $\rho(x)$ is the density defined in $x\in\Omega$, $c$ is the sound speed
and $\boldsymbol{n}$ is the outward unitary vector. As an assumption on the
density, we consider two constants $\overline{\rho}$ and $\underline{\rho}$,
both positives, such that, $\underline{\rho}\leq\rho(x)\leq\overline{\rho}$,
for $x\in\Omega$.

Let us introduce the displacement formulation for (2.1). Let us define the
space

$\mathcal{V}:=\\{\boldsymbol{v}\in\mathrm{H}(\mathop{\mathrm{div}}\nolimits,\Omega)\,:\,\boldsymbol{v}\cdot\boldsymbol{n}=0\,\,\text{on}\,\,\partial\Omega\\}.$

With this space at hand, a variational formulation for the acoustic problem is
the following: Find $(\omega,\boldsymbol{u})\in\mathbb{R}\times\mathcal{V}$,
with $\boldsymbol{u}\neq\boldsymbol{0}$, such that

(2.2)
$\displaystyle\int_{\Omega}\rho(x)c^{2}\mathop{\mathrm{div}}\nolimits\boldsymbol{u}\mathop{\mathrm{div}}\nolimits\boldsymbol{v}=\omega^{2}\int_{\Omega}\rho(x)\boldsymbol{u}\cdot\boldsymbol{v}\quad\forall\boldsymbol{v}\in\mathcal{V}.$

Defining $\widehat{\lambda}:=\omega^{2}$, we immediately deduce that the
eigenspace associated to $\widehat{\lambda}=0$ is

$\mathcal{K}:=\\{\boldsymbol{v}\in\mathcal{V}\,:\,\mathop{\mathrm{div}}\nolimits\boldsymbol{v}=0\,\,\text{in}\,\Omega\\}.$

We observe that (2.2) is not an elliptic problem. Hence, with a shift argument
we obtain the following problem: Find
$(\lambda,\boldsymbol{u})\in\mathbb{R}\times\mathcal{V}$, with
$\boldsymbol{u}\neq\boldsymbol{0}$, such that

$\displaystyle\int_{\Omega}\rho(x)c^{2}\mathop{\mathrm{div}}\nolimits\boldsymbol{u}\mathop{\mathrm{div}}\nolimits\boldsymbol{v}+\int_{\Omega}\rho(x)\boldsymbol{u}\cdot\boldsymbol{v}=\lambda\int_{\Omega}\rho(x)\boldsymbol{u}\cdot\boldsymbol{v}\quad\forall\boldsymbol{v}\in\mathcal{V},$

where $\lambda:=1+\widehat{\lambda}$.

Let us introduce the bilinear forms
$a:\mathcal{V}\times\mathcal{V}\rightarrow\mathbb{R}$ and
$b:\mathcal{V}\times\mathcal{V}\rightarrow\mathbb{R}$ defined by

$\displaystyle
a(\boldsymbol{v},\boldsymbol{w}):=\int_{\Omega}\rho(x)c^{2}\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\mathop{\mathrm{div}}\nolimits\boldsymbol{w}+\int_{\Omega}\rho(x)\boldsymbol{v}\cdot\boldsymbol{w}\quad\forall\boldsymbol{v},\boldsymbol{w}\in\mathcal{V}.$

and

$b(\boldsymbol{v},\boldsymbol{w}):=\int_{\Omega}\rho(x)\boldsymbol{v}\cdot\boldsymbol{w},\quad\forall\boldsymbol{v},\boldsymbol{w}\in\mathcal{V}.$

With these bilinear forms at hand, problem (2.2) is rewritten as follows: Find
$(\lambda,\boldsymbol{u})\in\mathbb{R}\times\mathcal{V}$, with
$\boldsymbol{u}\neq\boldsymbol{0}$, such that

(2.3) $a(\boldsymbol{u},\boldsymbol{v})=\lambda
b(\boldsymbol{u},\boldsymbol{v})\quad\forall\boldsymbol{v}\in\mathcal{V}.$

Clearly $a(\cdot,\cdot)$ is continuous, in the sense that there exists a
positive constant $M:=\|\rho\|_{\infty,\Omega}\min\\{c^{2},1\\}$ such that
$|a(\boldsymbol{u},\boldsymbol{v})|\leq
M\|\boldsymbol{u}\|_{\mathop{\mathrm{div}}\nolimits,\Omega}\|\boldsymbol{v}\|_{\mathop{\mathrm{div}}\nolimits,\Omega}$.
Observe that for any $\boldsymbol{v}\in\mathcal{V}$ there holds

(2.4)
$a(\boldsymbol{v},\boldsymbol{v})\geq\underline{\rho}c^{2}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}+\underline{\rho}\|\boldsymbol{v}\|_{0,\Omega}^{2}\geq\underbrace{\min\\{\underline{\rho}c^{2},\underline{\rho}\\}}_{\alpha^{c}}\|\boldsymbol{v}\|_{\mathop{\mathrm{div}}\nolimits,\Omega}^{2},$

implying the coercivity of $a(\cdot,\cdot)$ with a constant
$\alpha^{c}:=\min\\{\underline{\rho}c^{2},\underline{\rho}\\}$. This allows us
to introduce the so called solution operator associated to (2.3), defined by
$\boldsymbol{T}:\mathcal{V}\rightarrow\mathcal{V}$ where for a given
$\boldsymbol{f}\in\mathcal{V}$, we have
$\boldsymbol{T}\boldsymbol{f}=\widehat{\boldsymbol{u}}$, where
$\widehat{\boldsymbol{u}}$ is the solution of the source problem

$a(\widehat{\boldsymbol{u}},\boldsymbol{v})=b(\boldsymbol{f},\boldsymbol{v})\quad\forall\boldsymbol{v}\in\mathcal{V}.$

Note that operator $\boldsymbol{T}$ is well defined and bounded thanks to Lax-
Milgram’s lemma, non-compact, and selfadjoint respect to $a(\cdot,\cdot)$ and
$b(\cdot,\cdot)$. Moreover, $(\lambda,\boldsymbol{u})$ is a solution of (2.2)
if and only if $(1/\lambda,\boldsymbol{u})$ is an eigenpair of
$\boldsymbol{T}$.

Let us remark that since we have the space $\mathcal{K}$ at hand, there exists
an operator $\boldsymbol{P}:\mathcal{V}\rightarrow\mathcal{V}$ such that the
space $\mathcal{V}$ can be decomposed as a direct sum given by
$\mathcal{V}=\mathcal{K}\oplus\boldsymbol{P}(\mathcal{V})$, where
$\boldsymbol{P}(\mathcal{V})$ is the orthogonal complement of $\mathcal{K}$
with respect to the inner product defined in $\mathcal{V}$.

The following technical results is available in [4, 5].

######  Lemma 2.1 (Orthogonal complement).

The space $\boldsymbol{P}(\mathcal{V})$ is defined by
$\boldsymbol{P}(\mathcal{V}):=\\{\nabla\,q\,:\,q\in\mathrm{H}^{1}(\Omega)\\}$
and satisfies $\mathcal{V}=\mathcal{K}\oplus\boldsymbol{P}(\mathcal{V})$,
which is a orthogonal decomposition in $\mathrm{L}^{2}(\Omega)^{\texttt{d}}$
and $\mathrm{H}(\mathop{\mathrm{div}}\nolimits,\Omega)$. Moreover, there
exists $s\in(1/2,1]$ such that for all $\boldsymbol{v}\in\mathcal{V}$, if
$\boldsymbol{v}=\boldsymbol{\varphi}+\nabla q$, with
$\boldsymbol{\varphi}\in\mathcal{K}$ and $\nabla
q\in\boldsymbol{P}(\mathcal{V})$, then $\nabla
q\in\mathrm{H}^{s}(\Omega)^{\texttt{d}}$ and $\|\nabla
q\|_{s,\Omega}\lesssim\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|_{0,\Omega}.$

Now we recall some important properties that $\boldsymbol{T}$ satisfies.

######  Lemma 2.2 (Invariance of $\boldsymbol{T}$ on
$\boldsymbol{P}(\mathcal{V})$).

There holds
$\boldsymbol{T}(\boldsymbol{P}(\mathcal{V}))\subset\boldsymbol{P}(\mathcal{V}).$

Also, there holds

$\boldsymbol{T}(P(\mathcal{V}))\subset\left\\{\boldsymbol{v}\in\mathrm{H}^{s}(\Omega)^{\texttt{d}}\,:\,\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\in\mathrm{H}^{1}(\Omega)\right\\},$

and for all $\boldsymbol{f}\in\mathcal{G}\cap\mathcal{V}$, if
$\boldsymbol{u}=\boldsymbol{T}\boldsymbol{f}$, the following estimate holds

$\|\boldsymbol{u}\|_{s,\Omega}+\|\mathop{\mathrm{div}}\nolimits\boldsymbol{u}\|_{1+s,\Omega}\lesssim\|\boldsymbol{f}\|_{\mathop{\mathrm{div}}\nolimits,\Omega},$

where this estimate is obtained with a bootstrap argument. Moreover, we have
that
$\boldsymbol{T}|_{\boldsymbol{P}(\mathcal{V})}:\boldsymbol{P}(\mathcal{V})\rightarrow\boldsymbol{P}(\mathcal{V})$
is compact.

Finally, we present the following spectral characterization of
$\boldsymbol{T}$.

######  Theorem 2.3.

The spectrum of $\boldsymbol{T}$ decomposes as
$\mathop{\mathrm{sp}}\nolimits(\boldsymbol{T})=\\{0,1\\}\cup\\{\mu_{k}\\}_{k\in\mathbb{N}}$
where

  1. 1.

$\mu=1$ is an infinite-multiplicity eigenvalue of $\boldsymbol{T}$ and its
associated eigenspace is $\mathcal{K}$,

  2. 2.

$\\{\mu_{k}\\}_{k\in\mathbb{N}}\subset(0,1)$ is a sequence of finite-
multiplicity eigenvalues if $\boldsymbol{T}$ which converge to zero and if
$\boldsymbol{u}$ is an eigenfunction of $\boldsymbol{T}$ associated with such
an eigenvalue, then there exists $\widehat{s}>1/2$, depending on the
eigenvalue, such that

$\|\boldsymbol{u}\|_{\widehat{s},\Omega}+\|\mathop{\mathrm{div}}\nolimits\boldsymbol{u}\|_{1+\widehat{s},\Omega}\lesssim\|\boldsymbol{u}\|_{\mathop{\mathrm{div}}\nolimits,\Omega},$

where the hidden constant also depends on the eigenvalue.

  3. 3.

$\mu=0$ is not an eigenvalue of $\boldsymbol{T}$.

With this result at hand, we are in a position to introduce the numerical
method of our interest.

##  3 The DG method

To begin, we need to set some definitions and preliminaries. For
$\Omega\subset\mathbb{R}^{\texttt{d}}$, for $\texttt{d}\in\\{2,3\\}$, let
$\mathcal{T}_{h}$ be a shape regular family of meshes which subdivide the
domain $\bar{\Omega}$ into triangles/tetrahedra $K$. We denote by $h_{K}$ the
diameter of the element $K\in\mathcal{T}_{h}$ and $h$ the maximum of the
diameters of all the elements of the mesh, i.e.
$h:=\max_{K\in\mathcal{T}_{h}}\\{h_{K}\\}$. Let $F$ be a closed set such that
$F\subset\overline{\Omega}$ is an interior edge/face if $F$ has a positive
$(n-1)$-dimensional measure and if there are distinct elements $K$ and
$K^{\prime}$ such that $F=\bar{K}\cap\bar{K}^{\prime}$. A closed subset
$F\subset\overline{\Omega}$ is a boundary edge/face if there exists
$K\in\mathcal{T}_{h}$ such that $F$ is an edge/face of $K$ and
$F=\bar{K}\cap\partial\Omega$. Let $\mathcal{F}_{h}^{*}$ and
$\mathcal{F}_{h}^{\partial}$ be the sets of interior edges/faces and boundary
edges/face, respectively.

For any $t\geq 0$, we define the following broken Sobolev space

$\mathrm{H}^{t}(\mathcal{T}_{h})^{\texttt{d}}:=\\{\boldsymbol{v}\in\mathrm{L}^{2}(\Omega)^{\texttt{d}}:\quad\boldsymbol{v}|_{K}\in\mathrm{H}^{t}(K)^{\texttt{d}}\quad\forall
K\in\mathcal{T}_{h}\\}.$

Also, the space of the skeletons of the triangulations $\mathcal{T}_{h}$ is
defined by
$\mathrm{L}^{2}(\mathcal{F}_{h}):=\prod_{F\in\mathcal{F}_{h}}\mathrm{L}^{2}(F),$
where $\mathcal{F}_{h}:=\mathcal{F}_{h}^{*}\cup\mathcal{F}_{h}^{\partial}$.

In the forthcoming analysis,
$h_{\mathcal{F}}\in\mathrm{L}^{2}(\mathcal{F}_{h})$ will represent the
piecewise constant function defined by $h_{\mathcal{F}}|_{F}:=h_{F}$ for all
$F\in\mathcal{F}_{h}$, where $h_{F}$ denotes the diameter of edge/face $F$.

Let $\mathcal{P}_{m}(\mathcal{T}_{h})$ be the space of piecewise polynomials
respect with to $\mathcal{T}_{h}$ of degree at most $m\geq 0$; namely,

$\mathcal{P}_{m}(\mathcal{T}_{h}):=\left\\{v\in\mathrm{L}^{2}(\Omega):\,v|_{K}\in\mathcal{P}_{m}(K)\,\,\,\,\,\forall
K\in\mathcal{T}_{h}\right\\}.$

For any $k\geq 1$, we define the finite element spaces
$\mathcal{V}_{h}:=\mathcal{P}_{k}(\mathcal{T}_{h})^{\texttt{d}}$ and
$\mathcal{V}_{h}^{c}:=\mathcal{V}_{h}\cap\mathcal{V}$. We observe that the
space $\mathcal{V}_{h}^{c}$ is the Brezzi-Douglas-Marini (BDM) finite element
space. Now, we recall some well-known properties of the space
$\mathcal{V}_{h}^{c}$.

We define averages
$\\{\boldsymbol{v}\\}\in\mathrm{L}^{2}(\mathcal{F}_{h})^{\texttt{d}}$ and
jumps $\llbracket\boldsymbol{v}\rrbracket\in\mathrm{L}^{2}(\mathcal{F}_{h})$
as follows

$\\{\boldsymbol{v}\\}_{F}:=(\boldsymbol{v}_{K}+\boldsymbol{v}_{K^{\prime}})/2\quad\text{and}\quad\llbracket\boldsymbol{v}\rrbracket_{F}:=\boldsymbol{v}_{K}\cdot\boldsymbol{n}_{K}+\boldsymbol{v}_{K^{\prime}}\cdot\boldsymbol{n}_{K^{\prime}}\quad\forall
F\in\mathcal{F}(K)\cap\mathcal{F}(K^{\prime}),$

where $\boldsymbol{n}_{K}$ is the outward unit normal vector to $\partial K$.
Also, on the boundary $\partial\Omega$ and for all
$F\in\mathcal{F}(K)\cap\partial\Omega$, the averages and jumps are defined by
$\\{\boldsymbol{v}\\}_{F}:=\boldsymbol{v}_{K}$ and
$\llbracket\boldsymbol{v}\rrbracket_{F}:=\boldsymbol{v}_{K}\cdot\boldsymbol{n}$,
respectively.

Before concluding this section, we shall note some results required for the
analysis that follows. For a vector field $\boldsymbol{v}\in\mathcal{V}_{h}$
we define
$\mathop{\mathrm{div}}\nolimits_{h}\boldsymbol{v}\in\mathrm{L}^{2}(\Omega)$ by
$\mathop{\mathrm{div}}\nolimits_{h}\boldsymbol{v}|_{K}:=\mathop{\mathrm{div}}\nolimits(\boldsymbol{v}|_{K})$
for all $K\in\mathcal{T}_{h}$ and we endow
$\mathcal{V}(h):=\mathcal{V}+\mathcal{V}_{h}$ with the following seminorm

$|\boldsymbol{v}|^{2}_{\mathcal{V}(h)}:=\|\mathop{\mathrm{div}}\nolimits_{h}\boldsymbol{v}\|^{2}_{0,\Omega}+\|h_{\mathcal{F}}^{-1/2}\llbracket\boldsymbol{v}\rrbracket\|^{2}_{0,\mathcal{F}^{*}_{h}},$

which is well defined in $\mathcal{V}(h)$ and the norm
$\|\boldsymbol{v}\|^{2}_{\textrm{DG}}:=|\boldsymbol{v}|^{2}_{\mathcal{V}(h)}+\|\boldsymbol{v}\|^{2}_{0,\Omega}.$
Another important ingredient for the analysis is the following discrete trace
inequality [13]:

(3.5) $\|h^{1/2}\\{v\\}\|_{0,\mathcal{F}}\lesssim\|v\|_{0,\Omega}\qquad\forall
v\in\mathcal{P}_{k}(\mathcal{T}_{h}).$

###  3.1 The DG acoustic problem

Now we introduce the DG discretization for the acoustic problem (2.2), which
reads as follows: Find $\lambda_{h}:=1+\widehat{\lambda}_{h}\in\mathbb{C}$ and
$\boldsymbol{0}\neq\boldsymbol{u}_{h}\in\mathcal{V}_{h}$ such that

(3.6)
$a_{h}(\boldsymbol{u}_{h},\boldsymbol{\tau}_{h})=\lambda_{h}b(\boldsymbol{u}_{h},\boldsymbol{\tau}_{h})\quad\forall\boldsymbol{\tau}_{h}\in\mathcal{V}_{h},$

where the bilinear form
$a_{h}:\mathcal{V}_{h}\times\mathcal{V}_{h}\rightarrow\mathbb{C}$ is defined
by

(3.7)
$a_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h}):=\int_{\Omega}\rho(x)c^{2}\mathop{\mathrm{div}}\nolimits_{h}\boldsymbol{u}_{h}\mathop{\mathrm{div}}\nolimits_{h}\boldsymbol{v}_{h}+\int_{\Omega}\rho(x)\boldsymbol{u}_{h}\cdot\boldsymbol{v}_{h}\\\
+\int_{\mathcal{F}^{*}_{h}}\frac{\texttt{a}_{S}}{h_{\mathcal{F}}}\,\llbracket\boldsymbol{u}_{h}\rrbracket\cdot\llbracket\boldsymbol{v}_{h}\rrbracket-\int_{\mathcal{F}^{*}_{h}}\\{\rho(x)c^{2}\mathop{\mathrm{div}}\nolimits_{h}\boldsymbol{u}_{h}\\}\cdot\llbracket\boldsymbol{v}_{h}\rrbracket-\varepsilon\int_{\mathcal{F}^{*}_{h}}\\{\rho(x)c^{2}\mathop{\mathrm{div}}\nolimits_{h}\boldsymbol{v}_{h}\\}\cdot\llbracket\boldsymbol{u}_{h}\rrbracket,$

where $\texttt{a}_{S}>0$ is the so called stabilization parameter and
$\varepsilon\in\\{-1,0,1\\}$. These parameters are relevant for our DG
methods. In one hand, the stabilization $\texttt{a}_{S}$ is a parameter that,
if is not correctly chosen, could yield to instabilities for the method. In
general, this parameter depends on theoretical and physical constants and, in
the numerical performance, also on the polynomial degrees. (cf. Section 5). On
the other hand, $\varepsilon$ establishes if the DG method is symmetric or
nonsymmetric. More precisely, if $\varepsilon=1$ we obtain the classic
symmetric interior penalty method (SIP), if $\varepsilon=-1$ we obtain the
non-symmetric interior penalty method (NIP) and if $\varepsilon=0$ the
incomplete interior penalty method (IIP) (see [1, 20, 19] for instance).

For the analysis, we introduce the following norm

$\|\boldsymbol{u}\|^{*}_{\textrm{DG}}:=\big{(}\|\boldsymbol{u}\|^{2}_{\textrm{DG}}+\|h_{\mathcal{F}_{h}}^{1/2}\\{\mathop{\mathrm{div}}\nolimits_{h}\boldsymbol{u}\\}\|^{2}_{0,\mathcal{F}_{h}^{*}}\big{)}^{1/2}.$

It is easy to check that for all
$\boldsymbol{u},\boldsymbol{v}\in\mathcal{V}(h)$, and
$\mathop{\mathrm{div}}\nolimits\boldsymbol{u},\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\in\mathrm{H}^{t}(\Omega)$
with $t>1/2$, the bilinear form $a_{h}(\cdot,\cdot)$ is bounded. In fact,
there exists a constant $C>0$, independent of $h$, such that
$\big{|}a_{h}(\boldsymbol{u},\boldsymbol{v})\big{|}\leq
C\|\boldsymbol{u}\|_{\textrm{DG}}^{*}\|\boldsymbol{v}\|_{\textrm{DG}}^{*}$.
Moreover, by means of (3.5) it is possible to prove that for all
$\boldsymbol{v}_{h}\in\mathcal{V}_{h}$ there exists a positive constant
$M_{\textrm{DG}}$ such that $|a_{h}(\boldsymbol{u},\boldsymbol{v}_{h})|\leq
M_{\textrm{DG}}\|\boldsymbol{u}\|_{\textrm{DG}}^{*}\|\boldsymbol{v}_{h}\|_{\textrm{DG}}$.

In order to analyze the discrete eigenvalue problem (3.6), we need to
decompose the space $\mathcal{V}^{c}_{h}$. With this aim, we consider the
following subspace of $\mathcal{K}$

$\mathcal{K}_{h}:=\left\\{\boldsymbol{v}_{h}\in\mathcal{V}^{c}_{h}:\,\,\mathop{\mathrm{div}}\nolimits\boldsymbol{v}_{h}=0\right\\}.$

With $\mathcal{K}_{h}$ at hand, we are interested in its orthogonal complement
at the discrete level. With this in mind, is necessary to introduce the
discrete counterpart of the operator $\boldsymbol{P}$, which we denote by
$\boldsymbol{P}_{h}$ and is defined by
$\boldsymbol{P}_{h}:\mathcal{V}_{h}^{c}\rightarrow\mathcal{V}_{h}^{c}$. With
this operator at hand, the following approximation result holds, an its proof
is identical to [20, Lemma 3.1].

######  Lemma 3.1.

The operator
$\boldsymbol{P}_{h}:\mathcal{V}_{h}^{c}\rightarrow\mathcal{V}_{h}^{c}$ is a
projection with kernel $\mathcal{K}_{h}$ and satisfies the following estimate

$\|(\boldsymbol{P}-\boldsymbol{P}_{h})\boldsymbol{v}_{h}\|_{\mathop{\mathrm{div}}\nolimits,\Omega}\lesssim
h^{s}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}_{h}\|_{0,\Omega}\quad\forall\boldsymbol{v}_{h}\in\mathcal{V}_{h}^{c},$

where the hidden constant is independent of $h$ and $s$ is as in Lemma 2.1.

The following technical result proved in [22, Proposition 5.2] for tensorial
fields, and that is easily adapted for vectorial fields, will be useful in the
forthcoming analysis.

######  Proposition 3.2.

There exist a projection
$\mathcal{I}_{h}:\,\mathcal{V}_{h}\to\mathcal{V}_{h}^{c}$ and two constants
$\underline{C},\bar{C}>0$, independent of $h$, such that

(3.8)
$\underline{C}\,\,\|\boldsymbol{v}_{h}\|_{\mathrm{DG}}\leq\Big{(}\|\mathcal{I}_{h}\boldsymbol{v}_{h}\|^{2}_{\mathop{\mathrm{div}}\nolimits,\Omega}+\|h_{\mathcal{F}}^{-1/2}\llbracket\boldsymbol{v}_{h}\rrbracket\|^{2}_{0,\mathcal{F}^{*}_{h}}\Big{)}^{1/2}\leq\bar{C}\,\,\|\boldsymbol{v}_{h}\|_{\mathrm{DG}}\quad\forall\boldsymbol{v}_{h}\in\mathcal{V}_{h}.$

Moreover, we have that

(3.9)
$\|\mathop{\mathrm{div}}\nolimits_{h}(\boldsymbol{v}_{h}-\mathcal{I}_{h}\boldsymbol{v}_{h})\|^{2}_{0,\Omega}+\sum_{K\in\mathcal{T}_{h}}h_{K}^{-2}\|\boldsymbol{v}_{h}-\mathcal{I}_{h}\boldsymbol{v}_{h}\|^{2}_{0,K}\leq
C\,\,\|h_{\mathcal{F}}^{-1/2}\llbracket\boldsymbol{v}_{h}\rrbracket\|^{2}_{0,\mathcal{F}^{*}_{h}},$

with $C>0$ independent of $h$.

Now, with this result at hand, we will prove that $a_{h}(\cdot,\cdot)$ is
coercive in $\mathcal{V}_{h}$, independent of the DG method that we are
considering.

######  Lemma 3.3.

For any $\varepsilon\in\\{-1,0,1\\}$, there exists a parameter
$\textup{{a}}^{*}>0$ such that for all $\textup{{a}}_{S}\geq\textup{{a}}^{*}$
there holds

$a_{h}(\boldsymbol{v}_{h},\boldsymbol{v}_{h})\geq\alpha_{\mathrm{DG}}\|\boldsymbol{v}_{h}\|_{\mathrm{DG}}^{2}\qquad\forall\boldsymbol{v}_{h}\in\mathcal{V}_{h},$

where $\alpha_{\mathrm{DG}}>0$ is independent of $h$.

######  Proof 3.4.

The proof follows the same arguments of those in [20, Lemma 3.2], taking into
account for our case, the presence of density and sound speed of the fluid.

######  Remark 3.5.

Notice that the coercivity constant of $a_{h}(\cdot,\cdot)$ depends strongly
on the physical parameters, namely the density and sound speed. Hence, the
stability of the method will depend on the configuration of the problem where
these parameters play an important role, together with the geometry of domain,
since the constants that appear on the discrete trace inequality (3.5) and
Proposition 3.2 depend on $\Omega$. This is important to have in mind, since
as we will analyze in the numerical section (cf. Section 5), the appearance of
spurious modes depend on the choice of a and hence, on $\texttt{a}^{*}$.

We notice that the discrete coercivity of $a_{h}(\cdot,\cdot)$ holds true for
any $\varepsilon\in\\{-1,0,1\\}$. Implying that symmetric and nonsymmetric
methods are well posed with our DG approach. Moreover, we introduce the
discrete solution operator associated to (3.6):

$\boldsymbol{T}^{\varepsilon}_{h}:\mathcal{V}\rightarrow\mathcal{V}_{h},\qquad\boldsymbol{f}\mapsto\boldsymbol{T}^{\varepsilon}_{h}\boldsymbol{f}:=\widetilde{\boldsymbol{u}}^{\varepsilon}_{h},$

where $\widetilde{\boldsymbol{u}}^{\varepsilon}_{h}\in\mathcal{V}_{h}$ is the
unique solution of the following discrete source problem:

$a_{h}(\widetilde{\boldsymbol{u}}^{\varepsilon}_{h},\boldsymbol{v}_{h})=b(\boldsymbol{f},\boldsymbol{v}_{h})\qquad\forall\boldsymbol{v}_{h}\in\mathcal{V}_{h}.$

Observe that the superindex $\varepsilon$ denotes if the operator is
associated to a symmetric or nonsymmetric method. This is important to have in
mind since the choice of $\varepsilon$ leads to different solutions that lie
on real or complex spaces.

Clearly $\boldsymbol{T}^{\varepsilon}_{h}$ is well defined thanks to Lax-
Milgram’s lemma. Moreover, there exists a constant $C>0$ independent of $h$
such that
$\lVert\boldsymbol{T}^{\varepsilon}_{h}\boldsymbol{f}\rVert_{\mathrm{DG}}\leq
C\lVert\boldsymbol{f}\rVert_{\mathop{\mathrm{div}}\nolimits,\Omega}$, for all
$\boldsymbol{f}\in\mathcal{V}$. Also, It is easy to check that
$(\lambda_{h},\boldsymbol{u}_{h})\in\mathbb{C}\times\mathcal{V}_{h}$ is a
solution of problem (3.6) if and only if
$(\mu_{h},\boldsymbol{u}_{h})\in\mathbb{C}\times\mathcal{V}_{h}$ with
$\mu_{h}=1/\lambda_{h}$ is an eigenpair of $\boldsymbol{T}^{\varepsilon}_{h}$,
i.e.,
$\boldsymbol{T}^{\varepsilon}_{h}\boldsymbol{u}_{h}=\frac{1}{\mu_{h}}\boldsymbol{u}_{h}.$

######  Proposition 3.6.

Let $\boldsymbol{f}\in\boldsymbol{P}(\mathcal{V})$ be such that
$\boldsymbol{T}\boldsymbol{f}:=\widetilde{\boldsymbol{u}}$. Then, for any
$\varepsilon\in\\{-1,0,1\\}$ the following best approximation result holds

$\|(\boldsymbol{T}-\boldsymbol{T}^{\varepsilon}_{h})\boldsymbol{f}\|_{\mathrm{DG}}\leq\frac{M_{DG}}{\alpha_{DG}}\inf_{\boldsymbol{v}_{h}\in\mathcal{V}_{h}}\|\boldsymbol{T}\boldsymbol{f}-\boldsymbol{v}_{h}\|_{\mathrm{DG}}^{*}.$

Moreover, for $s$ as in Lemma 2.1, the following estimate holds

$\|(\boldsymbol{T}-\boldsymbol{T}^{\varepsilon}_{h})\boldsymbol{f}\|_{\mathrm{DG}}\lesssim
h^{s}(\|\widetilde{\boldsymbol{u}}\|_{s,\Omega}+\|\mathop{\mathrm{div}}\nolimits\widetilde{\boldsymbol{u}}\|_{1+s,\Omega}),$

where the hidden constant is independent of $h$.

Finally we derive the following result which holds for symmetric and non
symmetric methods.

######  Lemma 3.7.

For any $\boldsymbol{v}_{h}\in\mathcal{V}_{h}$, the following estimate holds

$\|(\boldsymbol{T}-\boldsymbol{T}^{\varepsilon}_{h})\boldsymbol{v}_{h}\|_{\mathrm{DG}}\lesssim
h^{s}\|\boldsymbol{v}_{h}\|_{\mathrm{DG}},$

where $s$ is the exponent provided by Lemma 2.1

###  3.2 Convergence analysis and error estimates

Since the solution operator $\boldsymbol{T}$ is non compact, the analysis of
the spectral convergence and error estimates will be performed under the
approach of the theory of [11, 12]. Let us recall the following definitions,
needed to perform the analysis.

  * •

P1.
$\lVert\boldsymbol{T}-\boldsymbol{T}_{h}^{\varepsilon}\rVert_{\mathcal{L}(\mathcal{V}_{h},\mathcal{V}(h))}\to
0$ as $h\to 0$.

  * •

P2. $\forall\boldsymbol{v}\in\mathcal{V}$, there holds

$\inf_{\boldsymbol{v}_{h}\in\mathcal{V}_{h}}\lVert\boldsymbol{v}-\boldsymbol{v}_{h}\rVert_{\mathrm{DG}}\,\,\to
0\quad\text{as}\quad h\to 0.$

In order to prove convergence between eigenspaces, we introduce the following
definitions: let $\boldsymbol{x}\in\mathcal{V}(h)$ and $\mathbb{E}$ and
$\mathbb{F}$ be closed subspaces of $\mathcal{V}(h)$. We define

$\delta(\boldsymbol{x},\mathbb{E}):=\inf_{\boldsymbol{y}\in\mathbb{E}}\lVert\boldsymbol{x}-\boldsymbol{y}\rVert_{\mathrm{DG}},\quad\delta(\mathbb{E},\mathbb{F}):=\sup_{\boldsymbol{y}\in\mathbb{E}:\,\lVert\boldsymbol{y}\rVert_{\mathrm{DG}}=1}\delta(\boldsymbol{y},\mathbb{F}).$

Hence, the gap between two closed subspaces is defined by

$\widehat{\delta}(\mathbb{E},\mathbb{F}):=\max\\{\delta(\mathbb{E},\mathbb{F}),\delta(\mathbb{F},\mathbb{E})\\}.$

Let $\lambda\in(0,1)$ be an isolated eigenvalue of $\boldsymbol{T}$ and let
$D$ an open disk in the complex plane with boundary $\gamma$ such that
$\lambda$ is the only eigenvalue of $\boldsymbol{T}$ lying in $D$ and
$\gamma\cap\mathop{\mathrm{sp}}\nolimits(\boldsymbol{T})=\emptyset$. We
introduce the spectral projector corresponding to the continuous and discrete
solution operators $\boldsymbol{T}$ and $\boldsymbol{T}_{h}^{\varepsilon}$,
respectively

$\displaystyle\mathcal{E}$ $\displaystyle:=\frac{1}{2\pi
i}\int_{\gamma}\left(z\boldsymbol{I}-\boldsymbol{T}\right)^{-1}\,dz:\mathcal{V}(h)\longrightarrow\mathcal{V}(h),$
$\displaystyle\mathcal{E}_{h}$ $\displaystyle:=\frac{1}{2\pi
i}\int_{\gamma}\left(z\boldsymbol{I}-\boldsymbol{T}^{\varepsilon}_{h}\right)^{-1}\,dz:\mathcal{V}(h)\longrightarrow\mathcal{V}(h)$

where $\mathcal{E}_{h}$ is well-defined and bounded uniformly in $h$.
Moreover, $\mathcal{E}|_{\mathcal{V}}$ is a spectral projection in
$\mathcal{V}$ onto the (finite dimensional) eigenspace
$\mathcal{E}(\mathcal{V})$ corresponding to the eigenvalue $\lambda$ of
$\boldsymbol{T}$. In fact, we have that
$\mathcal{E}(\mathcal{V}(h))=\mathcal{E}(\mathcal{V})$ (see [19] for further
details). Moreover, $\mathcal{E}_{h}|_{\mathcal{V}_{h}}$ is a projector in
$\mathcal{V}_{h}$ onto the eigenspace $\mathcal{E}_{h}(\mathcal{V}_{h})$
corresponding to the eigenvalues of
$\boldsymbol{T}^{\varepsilon}_{h}:\,\mathcal{V}_{h}\to\mathcal{V}_{h}$
contained in $\gamma$. We also have that
$\mathcal{E}_{h}(\mathcal{V}(h))=\mathcal{E}_{h}(\mathcal{V}_{h})$.

Now, we will compare $\mathcal{E}_{h}(\mathcal{V}_{h})$ to
$\mathcal{E}(\mathcal{V})$ in terms of the gap $\widehat{\delta}$. The proof
of the next auxiliary result follows from the definition of $\mathcal{E}$ and
$\mathcal{E}_{h}$.

######  Lemma 3.8.

The following estimate holds

$\displaystyle\lVert\mathcal{E}-\mathcal{E}_{h}\rVert_{\mathcal{L}(\mathcal{V}_{h},\mathcal{V}(h))}\leq
C\lVert\boldsymbol{T}-\boldsymbol{T}^{\varepsilon}_{h}\rVert_{\mathcal{L}(\mathcal{V}_{h},\mathcal{V}(h))},$

where the hidden constant is independent of $h$.

The following result will be used to establish the approximation properties of
the eigenfunctions of the continuous problem by means of those of the discrete
DG discretization (see [20, Lemma 4.5]).

######  Lemma 3.9.

The following estimate holds

$\widehat{\delta}(\mathcal{E}(\mathcal{V}),\mathcal{E}_{h}(\mathcal{V}_{h}))\leq
C\Big{(}\lVert\boldsymbol{T}-\boldsymbol{T}^{\varepsilon}_{h}\rVert_{\mathcal{L}(\mathcal{V}_{h},\mathcal{V}(h))}+\delta(\mathcal{E}(\mathcal{V}),\mathcal{V}_{h})\Big{)},$

where the hidden constant is independent of $h$.

Now, we state the convergence properties of the DG methods.

######  Theorem 3.10.

Let $\lambda\in(0,1)$ be an eigenvalue of $\boldsymbol{T}$ of algebraic
multiplicity $m$ and let $D_{\lambda}$ be a closed disk in the complex plane
centered at $\kappa$ with boundary $\gamma$ such that
$D_{\lambda}\cap\mathop{\mathrm{sp}}\nolimits(\boldsymbol{T})=\\{\kappa\\}$.
Let $\lambda_{1,h},\ldots,\lambda_{m(h),h}$ be the eigenvalues of
$\boldsymbol{T}_{h}^{\varepsilon}$ lying in $D_{\lambda}$ and repeated
according to their algebraic multiplicity. Then, for any DG method defined by
$\varepsilon\in\\{-1,0,1\\}$, we have that $m(h)=m$ for $h$ sufficiently small
and

$\lim_{h\to 0}\max_{1\leq i\leq m}|\lambda-\lambda_{i,h}|=0.$

Moreover, if $\mathcal{E}(\mathcal{V})$ is the eigenspace corresponding to
$\lambda$ and $\mathcal{E}_{h}(\mathcal{V}_{h})$ is the
$\boldsymbol{T}_{h}^{\varepsilon}$-invariant subspace of $\mathcal{V}_{h}$
spanned by the eigenspaces corresponding to $\\{\lambda_{i,h},\hskip
2.84544pti=1,\ldots,m\\}$ then

$\lim_{h\to
0}\widehat{\delta}(\mathcal{E}(\mathcal{V}),\mathcal{E}_{h}(\mathcal{V}_{h}))=0.$

######  Proof 3.11.

See proof of [19, Theorem 5.2].

Let us introduce the following distance

$\delta^{*}(\mathcal{E}(\mathcal{V}),\mathcal{V}_{h}):=\sup_{{\boldsymbol{v}\in\mathcal{E}(\mathcal{V})}\atop{\lVert\boldsymbol{v}\rVert_{\mathrm{DG}}=1}}\inf_{\boldsymbol{v}_{h}\in\mathcal{V}_{h}}\lVert\boldsymbol{v}-\boldsymbol{v}_{h}\rVert^{*}_{\mathrm{DG}}.$

With this distance at hand, we derive the following results, which have been
established in [19, Theorem 6.1] for a fixed eigenvalue $\lambda\in(0,1)$ of
$\boldsymbol{T}$.

######  Theorem 3.12.

For $h$ small enough, the following estimate holds

$\widehat{\delta}\big{(}\mathcal{E}(\mathcal{V}),\mathcal{E}_{h}(\mathcal{V}_{h})\big{)}\lesssim\delta^{*}(\mathcal{E}(\mathcal{V}),\mathcal{V}_{h}),$

where the hidden constant is independent of $h$.

Finally, with the aid of Proposition 2.1, we present the rates of convergence
of the proposed DG methods. The proof of the following result is identical to
the one contained in [20, Theorem 4.4] and holds for both, symmetric and
nonsymmetric DG methods.

######  Theorem 3.13.

Let $r>0$ be such that
$\mathcal{E}(\mathcal{V})\subset\\{\boldsymbol{v}\in\mathrm{H}^{r}(\Omega):\,\,\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\in\mathrm{H}^{1+r}(\Omega)^{n}\\}$
(cf. Proposition 2.1). Then, for $\varepsilon\in\\{-1,0,1\\}$ we have

$\displaystyle\widehat{\delta}(\mathcal{E}_{h}(\mathcal{V}_{h}),\mathcal{E}(\mathcal{V}))\lesssim
h^{\min\\{r,k\\}}\quad\text{ and }\quad\max_{1\leq i\leq
m}|\lambda-\lambda_{i,h}|\lesssim\,h^{\sigma\min\\{r,k\\}},$

where $\sigma:=\frac{(\varepsilon+2)^{\varepsilon}-1}{2}+1$ and the hidden
constants are independent of $h$.

##  4 The pressure formulation

An alternative approach to analyze the acoustic problem (2.1) is to consider
the pressure formulation which is obtained eliminating the displacement in
(2.1), using the second equation, in order to write

(4.10) $\left\\{\begin{aligned}
c^{2}\mathop{\mathrm{div}}\nolimits\left(\frac{1}{\rho(x)}\nabla
p\right)+\frac{\omega^{2}}{\rho(x)}p&=&0&\text{ in }\,\Omega,\\\
\frac{\partial p}{\partial\boldsymbol{n}}&=&0&\text{ on
}\,\partial\Omega,\end{aligned}\right.$

where (4.10) is now a scalar problem. Clearly the displacement is recovered
with a postprocess of the pressure, using the first equation on (2.1). Let us
remark that the computed solutions, namely the eigenvalues and eigenfunctions,
are exactly the same when the pure displacement formulation is considered. A
variational formulation for this problem is: Find $\omega\in\mathbb{R}$ and
$0\neq p\in\mathrm{H}^{1}(\Omega)$ such that

$\displaystyle\int_{\Omega}\frac{1}{\rho(x)}\nabla p\cdot\nabla
v=\omega^{2}\int_{\Omega}\frac{1}{\rho(x)}pv\quad\forall
v\in\mathrm{H}^{1}(\Omega).$

With a shift argument we arrive to the following problem: Find
$\omega\in\mathbb{R}$ and $0\neq p\in\mathrm{H}^{1}(\Omega)$ such that

(4.11) $a^{p}(p,v)=(\omega^{2}+1)b^{p}(p,v)\quad\forall
v\in\mathrm{H}^{1}(\Omega),$

where the bilinear forms
$a:\mathrm{H}^{1}(\Omega)\times\mathrm{H}^{1}(\Omega)\rightarrow\mathbb{R}$
and
$b:\mathrm{H}^{1}(\Omega)\times\mathrm{H}^{1}(\Omega)\rightarrow\mathbb{R}$
are defined for all $q,v\in\mathrm{H}^{1}(\Omega)$ by

$a^{p}(q,v):=c^{2}\int_{\Omega}\frac{1}{\rho(x)}\nabla q\cdot\nabla
v+\int_{\Omega}\frac{1}{\rho(x)}qv\quad\text{and}\quad
b^{p}(q,v):=\int_{\Omega}\frac{1}{\rho(x)}qv.$

Note that $a^{p}(\cdot,\cdot)$ is coercive in $\mathrm{H}^{1}(\Omega)$. Also,
the solution operator
$T^{p}:\mathrm{L}^{2}(\Omega)\rightarrow\mathrm{H}^{1}(\Omega)\hookrightarrow\mathrm{L}^{2}(\Omega)$
is defined by $T^{p}f=\widetilde{p}$ where
$\widetilde{p}\in\mathrm{H}^{1}(\Omega)$ is the solution of the corresponding
associated source problem. This solution operator results to be compact due
the compact inclusion of $\mathrm{H}^{1}(\Omega)$ onto
$\mathrm{L}^{2}(\Omega)$. On the other hand, from Lemma 2.2, we have that
$\boldsymbol{u}\in\mathrm{H}^{s}(\Omega)$ with $s\in(1/2,1]$. Hence, from the
first equation of (2.1) we have that $\nabla p\in\mathrm{H}^{s}(\Omega)$ and
hence $p\in\mathrm{H}^{1+s}(\Omega)$, providing the regularity requirement for
the eigenfunctions when problem (4.11) is considered. For this case, and for
any $k\geq 1$, the finite element spaces are
$\mathcal{V}_{h}^{p}:=\mathcal{P}_{k}(\mathcal{T}_{h})$ and
$\mathcal{V}_{h}^{p,c}:=\mathcal{V}_{h}^{p}\cap\mathrm{H}^{1}(\Omega)$, which
now are scalar. Inspired by the analysis of [1], let us define
$\mathcal{V}^{p}(h):=\mathrm{H}^{1}(\Omega)+\mathcal{V}^{p}_{h}$ which we
endow with the following seminorm
$\|v\|_{\mathcal{V}(h)}^{2}=\|\nabla_{h}v\|_{0,\Omega}^{2}+\|h_{\mathcal{F}}^{-1/2}\llbracket
v\rrbracket\|_{0,\mathcal{F}_{h}}^{2}.$ Then, we define the norm
$\|v_{h}\|_{DG}^{2}=\|v_{h}\|_{0,\Omega}^{2}+\|v_{h}\|_{\mathcal{V}(h)}^{2}$.

The DG discretization for problem (4.11) is as follows: Find
$\omega_{h}\in\mathbb{C}$ and $0\neq p_{h}\in\mathcal{V}_{h}^{p}$ such that

(4.12)
$a_{h}^{p}(p_{h},v_{h})=(\omega_{h}^{2}+1)b^{p}(p_{h},v_{h})\quad\forall
v_{h}\in\mathcal{V}_{h}^{p},$

where

$a_{h}^{p}(p_{h},v_{h}):=c^{2}\int_{\Omega}\frac{1}{\rho(x)}\nabla_{h}p_{h}\cdot\nabla_{h}v_{h}+\int_{\Omega}\frac{1}{\rho(x)}p_{h}v_{h}+\int_{\mathcal{F}^{*}_{h}}\frac{\texttt{a}_{S,p}}{h_{\mathcal{F}}}\llbracket
p_{h}\rrbracket\cdot\llbracket v_{h}\rrbracket\\\
-\int_{\mathcal{F}^{*}_{h}}\frac{1}{\rho(x)}\\{c^{2}\nabla_{h}p_{h}\\}\cdot\llbracket
v_{h}\rrbracket-\varepsilon\int_{\mathcal{F}^{*}_{h}}\frac{1}{\rho(x)}\\{c^{2}\nabla_{h}v_{h}\\}\cdot\llbracket
p_{h}\rrbracket.$

Let us prove that $a_{h}^{p}(\cdot,\cdot)$ is coercive in
$\mathcal{V}_{h}^{p}$.

######  Lemma 4.1 (ellipticity of $a_{h}^{p}(\cdot,\cdot)$).

For any $\varepsilon\in\\{-1,0,1\\}$, there exists a positive parameter
$\texttt{a}_{p}^{*}$ such that for all
$\texttt{a}_{S,p}\geq\texttt{a}_{p}^{*}$ there holds

$a_{h}^{p}(v_{h},v_{h})\geq\alpha_{p}\|v_{h}\|_{DG}^{2}\quad\forall
v_{h}\in\mathcal{V}_{h}^{p},$

where $\alpha_{p}>0$ is independent of $h$.

######  Proof 4.2.

Let $v_{h}\in\mathcal{V}_{h}^{p}$. Then, from the definition of
$a_{h}^{p}(\cdot,\cdot)$ we have

$a_{h}^{p}(v_{h},v_{h})=c^{2}\int_{\Omega}\frac{1}{\rho(x)}\nabla_{h}v_{h}\cdot\nabla_{h}v_{h}+\int_{\Omega}\frac{1}{\rho(x)}v_{h}v_{h}+\int_{\mathcal{F}^{*}_{h}}\frac{\texttt{a}_{S,p}}{h_{\mathcal{F}}}\llbracket
v_{h}\rrbracket\cdot\llbracket v_{h}\rrbracket\\\
-\int_{\mathcal{F}^{*}_{h}}\frac{1}{\rho(x)}\\{c^{2}\nabla_{h}v_{h}\\}\cdot\llbracket
v_{h}\rrbracket-\varepsilon\int_{\mathcal{F}^{*}_{h}}\frac{1}{\rho(x)}\\{c^{2}\nabla_{h}v_{h}\\}\cdot\llbracket
v_{h}\rrbracket\\\
\geq\frac{c^{2}}{\overline{\rho}}\|\nabla_{h}v_{h}\|_{0,\Omega}^{2}+\frac{1}{\overline{\rho}}\|v_{h}\|_{0,\Omega}^{2}+\texttt{a}_{S,p}c^{2}\|h_{\mathcal{F}}^{-1/2}\llbracket
v_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}^{2}+\frac{c^{2}(1-\varepsilon)}{\overline{\rho}}\int_{\mathcal{F}^{*}_{h}}\\{\nabla_{h}v_{h}\\}\cdot\llbracket
v_{h}\rrbracket\\\
\geq\frac{c^{2}}{\overline{\rho}}\|\nabla_{h}v_{h}\|_{0,\Omega}^{2}+\frac{1}{\overline{\rho}}\|v_{h}\|_{0,\Omega}^{2}+\texttt{a}_{S,p}c^{2}\|h_{\mathcal{F}}^{-1/2}\llbracket
v_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}^{2}\\\
+\frac{c^{2}(1-\varepsilon)}{2\overline{\rho}}(-\|h_{e}^{1/2}\\{\nabla
v_{h}\\}\|_{0,\mathcal{F}_{h}^{*}}^{2}-\|h_{\mathcal{F}}^{-1/2}\llbracket
v_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}^{2})\\\
\geq\underbrace{\min\left\\{\frac{1}{\overline{\rho}}\left(c^{2}-C\left(\frac{1+\varepsilon}{2}\right)\right),\frac{1}{\overline{\rho}}\right\\}}_{C_{1}}(\|\nabla_{h}v_{h}\|_{0,\Omega}^{2}+\|v_{h}\|_{0,\Omega}^{2})\\\
+c^{2}\underbrace{\left(\texttt{a}_{S,p}-\left(\frac{1+\varepsilon}{2\overline{\rho}}\right)\right)}_{C_{2}}\|h_{\mathcal{F}}^{-1/2}\llbracket
v_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}^{2}\geq\alpha_{p}\|v_{h}\|_{DG}^{2},$

where the constant $C$ is the one provided by (3.5). Observe that, since
$\nu>0$, the ellipticity holds for $\nu>C(1+\varepsilon)/2$ and
$\texttt{a}_{S,p}>(1+\varepsilon)/2\overline{\rho}$, with
$\varepsilon\in\\{-1,0,1\\}$. Hence, defining
$\alpha_{p}:=\min\\{C_{1},C_{2}\\}$ and choosing $\texttt{a}_{S}$ such that
$\texttt{a}_{S,p}>\texttt{a}_{p}^{*}:=(1+\varepsilon)/2\overline{\rho}$ we
conclude the proof.

Now, the discrete solution operator
$T_{h}^{p,\varepsilon}:\mathrm{H}^{1}(\Omega)\rightarrow\mathcal{V}_{h}^{p}$
is such that $T_{h}^{p,\varepsilon}f:=\widetilde{p}_{h}$ where
$\widetilde{p}_{h}$ is the solution of the source problem
$a_{h}(\widetilde{p}_{h},v_{h})=b^{p}(f,v_{h})$, for all
$v_{h}\in\mathcal{V}_{h}^{p}$. In this formulation,
$\varepsilon\in\\{-1,0,1\\}$ plays the same role as in the displacement
formulation, providing symmetric or nonsymmetric DG methods.

Let us remark that, despite to the fact that the solution operator
$T^{p,\varepsilon}$ is continuous and compact, the classic theory of [2] is
not enough to our IPDG methods, since this theory requires that the numerical
method must be conforming, implying that the noncompact theory of operators
must be used for the pressure formulation. Then, properties P1 and P2 in this
context reads as follows

  * •

P1. $\lVert
T^{p}-T_{h}^{p,\varepsilon}\rVert_{\mathcal{L}(\mathcal{V}^{p}_{h},\mathcal{V}^{p}(h))}\to
0$ as $h\to 0$.

  * •

P2. $\forall\tau\in\mathcal{V}^{p}_{h}$, there holds

$\inf_{\tau\in\mathcal{V}^{p}_{h}}\lVert\tau-\tau_{h}\rVert_{\mathcal{V}^{p}(h)}\,\,\to
0\quad\text{as}\quad h\to 0.$

The following convergence result holds for the continuous and discrete
solution operators.

######  Lemma 4.3.

For all $f\in\mathcal{V}^{p}$, the following estimate holds

$\|(T^{p}-T^{p,\varepsilon}_{h})f\|_{\mathcal{V}^{p}(h)}\lesssim
h^{s}\|f\|_{0,\Omega},$

with $s\in(1/2,1]$ and the hidden constant is independent of $h$.

######  Proof 4.4.

From the definition of the continuous and discrete solutions operators, we
have $\widehat{p}:=T^{p}f$ and $\widehat{p}_{h}:=T^{p,\varepsilon}_{h}f$. Now,
since $\widehat{p}_{h}$ is the solution of the corresponding source problem,
then $\|\widehat{p}-\widehat{p}_{h}\|_{\mathcal{V}^{p}(h)}$ is precisely the
error of the DG method applied on the source problem. Then, the result follow
with $s\in(1/2,1]$.

If we consider discrete sources, we obtain the analogous of the previous
lemma, in the sense that For all $f_{h}\in\mathcal{V}_{h}$, the following
estimate holds

$\|(T^{p}-T^{p,\varepsilon}_{h})f_{h}\|_{\mathcal{V}^{p}(h)}\lesssim
h^{s}\|f_{h}\|_{\mathcal{V}^{p}(h)},$

with $s\in(1/2,1]$ and the hidden constant is independent of $h$. Now we are
in position to establish P1.

######  Lemma 4.5.

There following estimate holds

$\lVert
T^{p}-T^{p,\varepsilon}_{h}\rVert_{\mathcal{L}(\mathcal{V}_{h},\mathcal{V}^{p}(h))}\lesssim
h^{s},$

where the hidden constant is independent of $h$.

######  Proof 4.6.

Given $f_{h}\in\mathcal{V}^{p}_{h}$, we have

$\lVert
T-T^{p,\varepsilon}_{h}\rVert_{\mathcal{L}(\mathcal{V}^{p}_{h},\mathcal{V}^{p}(h)}:=\sup_{0\neq
f_{h}\in\mathcal{V}^{p}_{h}}\frac{\lVert(T^{p}-T^{p,\varepsilon}_{h})f_{h}\rVert_{\mathcal{V}^{p}(h)}}{{\lVert
f_{h}\rVert_{\mathcal{V}^{p}(h)}}}\lesssim h^{s}.$

This concludes the proof.

If $\mathcal{E}(\mathcal{V}^{p})$ represents the invariant space for the
eigenfunctions of the pressure formulation, such that if
$p\in\mathcal{E}(\mathcal{V}^{p})$ then $p\in\mathrm{H}^{1+r}(\Omega)$ for
$r>0$ and $\mathcal{E}(\mathcal{V}^{p}_{h})$ is its discrete counterpart, we
are able to provide an error estimate for the eigenfunctions. To do this task,
and taking into account that the pressure formulation is nothing else that the
Laplacian operator, the arguments of [1] can be easily followed for our
purposes.

######  Lemma 4.7.

For $h$ small enough it holds

$\widehat{\delta}_{h}(\mathcal{E}(\mathcal{V}^{p}),\mathcal{E}_{h}(\mathcal{V}^{p}_{h}))\lesssim
h^{\min\\{r,k\\}},$

where the hidden constant is independent of $h$.

Finally we present error estimates for the eigenvalues for the symmetric and
nonsymmetric methods, where optimal and suboptimal order of convergence are
attained, respectively.

######  Theorem 4.8.

There exists a strictly positive constant $h_{0}$ such that, for $h<h_{0}$
there holds

  1. 1.

If the symmetric IPDG method is considered $(\varepsilon=1)$, then there holds

$|\lambda-\lambda_{h}|\lesssim h^{2\min\\{r,k\\}},$

  2. 2.

If any of the nonsymmetric IPDG methods are considered
$(\varepsilon\in\\{-1,0\\})$, then there holds

$|\lambda-\lambda_{h}|\lesssim h^{\min\\{r,k\\}},$

where in each estimate the hidden constant is independent of $h$.

##  5 Numerical experiments

Now we report a series of numerical experiments in order to corroborate the
robustness of the proposed method with $\varepsilon\in\\{-1,0,1\\}$. All the
results have been obtained with a FEniCS script [16]. The computed order of
convergence for the eigenvalues were obtained by means of a standard least
square fitting. When no analytical solution is available, we compare the
computed eigenvalues for each mesh with extrapolated values. Let us denote by
$N$ the refinement level of the mesh and $\lambda$ the eigenvalues. For the
experiments we consider three cases. In the first case, we study the influence
of the stabilization parameter a when the spectrum is computed, in order to
observe the arising of spurious modes. This is an important analysis since,
according to Lemma 3.3, the method is stable for a certain threshold of a. Let
us remark that this threshold changes when the configuration of the method
changes, namely, physical constants, geometry of the domain, boundary
conditions, etc. The second situation studies the convergence of the methods,
where for $\varepsilon=1$ the optimal order is expected, whereas for
$\varepsilon\in\\{-1,0\\}$ suboptimal order of convergence must be attained
according to, for instance, [8, 20]. Finally, we consider domains with
singularities and non-polygonal domains to observe the behavior of the method
in these geometries.

Throughout this section, the stabilization parameter will be chosen
proportional to the square of the polynomial degree $k$ as
$\texttt{a}_{S}=\texttt{a}k^{2}$, with $\texttt{a}>0$. We begin our report of
numerical results considering first the discrete displacement formulation
(3.6). Then, we will study the pure pressure numerical scheme (4.12).

From the results reported in [19, 20] for the elasticity and Stokes equations,
respectively it follows that for $\texttt{a}<1$, the spurious eigenvalues
appear more strongly compared with those computed when $\texttt{a}>1$.
Therefore, considering that the coercivity constant depends on the density, in
the spurious analysis experiments, we will present different combinations of
stabilization parameters that depend on the lower and upper bounds of the
density.

###  5.1 Test 1: Influence of the stabilization

In this experiment we test the effect of the stabilization parameter on each
discontinuous scheme, namely SIP ($\varepsilon=1$), IIP ($\varepsilon=0$) and
NIP ($\varepsilon=-1$). The importance of this study lies in determining which
is the safest parameter to avoid pollution on the spectrum. To detect spurious
modes, we compare eigenvalues and eigenfunctions.

The domain is the rectangle $\Omega:=(0,a)\times(0,b)$. We choose
$\rho_{1}(x,y)=\frac{1}{x^{2}+y^{2}+1}$ and $\rho_{2}(x,y)=e^{xy+1}$, which
correspond to density functions such that $|\rho_{1}(x,y)|\leq 1$ and
$|\rho_{2}(x,y)|>1$, respectively. For the experiment, we choose $a=1$ and
$b=1.1$, whereas $c=1$.

####  5.1.1 SIP method

In this test, we take $\varepsilon=1$ in (3.7). The first goal is to study the
influence of the stabilization parameter $\texttt{a}_{S}$ on the computation
of the spectrum. More precisely, the intention is to analyze the appearance of
spurious eigenvalues for certain values and localize a threshold in which the
method is safe to compute the eigenvalues with physical meaning. In the
forthcoming results, we report the computed eigenvalues fixing the refinement
parameter $N$ and taking different values of a, and polynomial degrees. The
numbers inside boxes represent spurious eigenvalues.

Table 1: Test 1. Computed eigenvalues $\lambda_{h,i}$ for polynomial degrees $k=1,2,3$, mesh level $N=8$ and different values of a and $\rho(x)$, for the SIP method ($\varepsilon=1$). | $\rho_{1}(x,y)$ | $\rho_{2}(x,y)$
---|---|---
k a | $2(\overline{\rho}+\underline{\rho})$ | $4\overline{\rho}$ | $4(\overline{\rho}+\underline{\rho})$ | $8\overline{\rho}$ | $2(\overline{\rho}+\underline{\rho})$ | $4\overline{\rho}$ | $4(\overline{\rho}+\underline{\rho})$ | $8\overline{\rho}$
1 | 7.83005 | 7.85010 | 7.86139 | 7.87395 | 7.92006 | 7.93951 | 7.94676 | 7.95428
9.58484 | 9.28119 | 9.64013 | 9.65956 | 10.29239 | 10.32611 | 10.33702 | 10.35267
17.38317 | 10.02391 | 17.53053 | 17.60191 | 18.16105 | 18.28682 | 18.326917 | 18.36967
26.83626 | 17.65229 | 32.78925 | 33.04550 | 32.91697 | 33.22771 | 33.41198 | 33.58185
32.20629 | 32.68393 | 39.74527 | 40.03675 | 39.80118 | 40.23012 | 40.43040 | 40.62000
39.05536 | 39.66805 | 42.43583 | 42.88156 | 41.68459 | 43.10124 | 43.38814 | 43.63388
41.58004 | 42.41955 | 47.72058 | 48.27411 | 42.68574 | 48.74642 | 49.15382 | 49.52661
46.71721 | 47.77799 | 72.84119 | 74.20859 | 48.66419 | 73.83369 | 74.61207 | 75.31965
49.80943 | 72.54681 | 75.13436 | 76.53154 | 72.59414 | 75.71586 | 76.64642 | 77.55341
70.17544 | 73.49931 | 84.38202 | 86.64709 | 74.59481 | 85.64420 | 87.01155 | 88.42254
2 | 7.83258 | 7.83261 | 7.83262 | 7.83263 | 7.91306 | 7.913092 | 7.913102 | 7.91312
9.59294 | 9.59299 | 9.59300 | 9.59302 | 10.27004 | 10.27012 | 10.27016 | 10.27018
17.42091 | 17.42146 | 17.42162 | 17.42179 | 18.16579 | 18.16631 | 18.16649 | 18.16666
32.34693 | 32.34999 | 32.35098 | 32.35200 | 32.75789 | 32.76038 | 32.76132 | 32.76219
39.15838 | 39.16262 | 39.16390 | 39.16518 | 39.60587 | 39.60917 | 39.61035 | 39.61145
41.87547 | 41.88456 | 41.88704 | 41.88956 | 42.53659 | 42.54257 | 42.54459 | 42.54646
47.03714 | 47.04933 | 47.05232 | 47.05533 | 47.98522 | 47.99580 | 47.99878 | 48.00148
71.48122 | 71.53121 | 71.54302 | 71.55498 | 72.23229 | 72.26248 | 72.27162 | 72.28008
73.16055 | 73.19808 | 73.20927 | 73.22068 | 73.57347 | 73.59884 | 73.60805 | 73.61669
82.66712 | 82.75860 | 82.77748 | 82.79638 | 83.59049 | 83.65245 | 83.66898 | 83.68400
3 | 7.83253 | 7.83253 | 7.83253 | 7.83253 | 7.91295 | 7.91295 | 7.91295 | 7.91295
9.59288 | 9.59288 | 9.59288 | 9.59288 | 10.26991 | 10.26991 | 10.26991 | 10.26992
17.42060 | 17.42060 | 17.42060 | 17.42060 | 18.16508 | 18.16508 | 18.16508 | 18.16509
32.34437 | 32.34437 | 32.34437 | 32.34438 | 32.75262 | 32.75263 | 32.75263 | 32.75264
39.15584 | 39.15584 | 9.155855 | 39.15585 | 39.59998 | 39.59998 | 39.59998 | 39.59998
41.87258 | 41.87260 | 41.87261 | 41.87262 | 42.52780 | 42.52781 | 42.52782 | 42.52784
47.03552 | 47.03554 | 47.03555 | 47.03556 | 47.97646 | 47.97648 | 47.97649 | 47.97650
71.47553 | 71.47572 | 71.47580 | 71.47587 | 72.19433 | 72.19446 | 72.19452 | 72.19457
73.13334 | 73.13345 | 73.13350 | 73.13354 | 73.51695 | 73.51703 | 73.51707 | 73.51710
82.66821 | 82.66852 | 82.66863 | 82.66875 | 83.53586 | 83.53609 | 83.53618 | 83.53628
Table 2: Test 1. Computed eigenvalues $\lambda_{h,i}$ for polynomial degrees $k=1,2,3$, mesh level $N=8$ and different values of a and $\rho(x)$, for the SIP method ($\varepsilon=1$). | $\rho_{1}(x,y)$ | $\rho_{2}(x,y)$
---|---|---
k a | 4 | $4(\overline{\rho}+1)$ | 8 | $8(\overline{\rho}+1)$ | 4 | $4(\overline{\rho}+1)$ | 8 | $8(\overline{\rho}+1)$
1 | 7.85010 | 7.87395 | 7.87395 | 7.88547 | 4.94719 | 7.94189 | 7.85203 | 7.95588
9.28119 | 9.65956 | 9.65956 | 9.67768 | 7.73864 | 10.32434 | 10.15531 | 10.35597
10.02391 | 17.60191 | 17.60191 | 17.66392 | 9.91112 | 18.29696 | 17.80452 | 18.37871
17.65229 | 33.04550 | 33.04550 | 33.27389 | 15.58674 | 33.31167 | 29.18053 | 33.61829
32.68393 | 40.03675 | 40.03675 | 40.29402 | 17.28392 | 40.31934 | 32.13420 | 40.66076
39.66805 | 42.88156 | 42.88156 | 43.24050 | 27.92958 | 43.21682 | 37.94840 | 43.68507
42.41955 | 48.27411 | 48.27411 | 48.73874 | 28.51991 | 48.86380 | 40.57725 | 49.60571
47.77799 | 74.20859 | 74.20859 | 75.17274 | 33.30603 | 59.44931 | 44.65770 | 75.46235
72.54681 | 76.53154 | 76.53154 | 77.73133 | 34.96146 | 74.52412 | 56.35742 | 77.74420
73.49931 | 86.64709 | 86.64709 | 88.35366 | 36.65087 | 76.40702 | 65.81340 | 88.70701
2 | 7.83261 | 7.83263 | 7.83263 | 7.83264 | 7.91487 | 7.91309 | 7.91191 | 7.91312
9.59299 | 9.59302 | 9.59302 | 9.59304 | 10.27009 | 10.27014 | 8.87486 | 10.27019
17.42146 | 17.42179 | 17.42179 | 17.42194 | 18.15392 | 18.16639 | 10.26904 | 18.16669
32.34999 | 32.35200 | 32.35200 | 32.35292 | 26.39878 | 32.76080 | 13.21076 | 32.76238
39.16262 | 39.16518 | 39.16518 | 39.16633 | 32.75376 | 39.60970 | 17.28802 | 39.61168
41.88456 | 41.88956 | 41.88956 | 41.89181 | 39.61667 | 42.54348 | 18.17024 | 42.54685
47.04933 | 47.05533 | 47.05533 | 47.05800 | 42.52010 | 47.99715 | 32.75270 | 48.00205
71.53121 | 71.55498 | 71.55498 | 71.56563 | 47.97542 | 72.26659 | 39.60443 | 72.28185
73.19808 | 73.22068 | 73.22068 | 73.23088 | 53.28165 | 73.60296 | 42.53533 | 73.61851
82.75860 | 82.79638 | 82.79638 | 82.81314 | 72.17575 | 83.65992 | 47.97276 | 83.68714
3 | 7.83253 | 7.83253 | 7.83253 | 7.83253 | 7.91295 | 7.91295 | 7.91295 | 7.91295
9.59288 | 9.59288 | 9.59288 | 9.59288 | 10.26991 | 10.26991 | 10.26991 | 10.26991
17.42060 | 17.42060 | 17.42060 | 17.42060 | 18.16508 | 18.16508 | 18.16508 | 18.16508
32.34437 | 32.34438 | 32.34438 | 32.34438 | 31.51368 | 32.75263 | 32.75260 | 32.75263
39.15584 | 39.15585 | 39.15585 | 39.15585 | 32.75255 | 39.59998 | 39.59995 | 39.59998
41.87260 | 41.87262 | 41.87262 | 41.87263 | 39.59988 | 42.52782 | 42.52768 | 42.52783
47.03554 | 47.03556 | 47.03556 | 47.03557 | 42.52740 | 47.97649 | 47.97635 | 47.97651
71.47572 | 71.47587 | 71.47587 | 71.47594 | 47.97424 | 72.19448 | 72.19350 | 72.19458
73.13345 | 73.13354 | 73.13354 | 73.13359 | 50.98591 | 73.51705 | 73.51650 | 73.51711
82.66852 | 82.66875 | 82.66875 | 82.66885 | 72.19381 | 83.53613 | 83.53511 | 83.53630
Table 3: Test 1. Convergence analysis for $\texttt{a}=10(\overline{\rho}+1)$, $\rho(x,y)=\frac{1}{x^{2}+y^{2}+1}$, and $k=1,2,3$ for the SIP method ($\varepsilon=1$).  k | $N=10$ | $N=20$ | $N=30$ | $N=40$ | Order | $\lambda_{i}$
---|---|---|---|---|---|---
1 | 7.86806 | 7.84149 | 7.83652 | 7.83478 | 1.98 | 7.83248
9.64947 | 9.60705 | 9.59919 | 9.59643 | 2.00 | 9.59290
17.58538 | 17.46235 | 17.43923 | 17.43109 | 1.98 | 17.42052
32.96979 | 32.50146 | 32.41428 | 32.38371 | 1.99 | 32.34399
2 | 7.83258 | 7.83253 | 7.83253 | 7.83253 | 3.99 | 7.83253
9.59295 | 9.59289 | 9.59288 | 9.59288 | 3.99 | 9.59288
17.42116 | 17.42063 | 17.42060 | 17.42060 | 3.98 | 17.42060
32.34797 | 32.34457 | 32.34438 | 32.34435 | 3.97 | 32.34434
3 | 7.83253 | 7.83253 | 7.83253 | 7.83253 | 5.85 | 7.83253
9.59288 | 9.59288 | 9.59288 | 9.59288 | 5.88 | 9.59288
17.42060 | 17.42060 | 17.42060 | 17.42060 | 5.98 | 17.42060
32.34435 | 32.34434 | 32.34434 | 32.34434 | 5.99 | 32.34434
Table 4: Test 1. Convergence analysis for $\texttt{a}=10(\overline{\rho}+1)$, $\rho(x,y)=e^{xy+1}$, and $k=1,2,3$ for the SIP method ($\varepsilon=1$).  k | $N=10$ | $N=20$ | $N=30$ | $N=40$ | Order | $\lambda_{i}$
---|---|---|---|---|---|---
1 | 7.94234 | 7.92041 | 7.91628 | 7.91483 | 1.97 | 7.91291
10.32857 | 10.28466 | 10.27648 | 10.27361 | 1.99 | 10.26989
18.31292 | 18.20280 | 18.18194 | 18.17459 | 1.96 | 18.16476
33.34798 | 32.90287 | 32.81956 | 32.79031 | 1.98 | 32.75189
2 | 7.91302 | 7.91296 | 7.912956 | 7.912956 | 3.99 | 7.91295
10.27003 | 10.26992 | 10.26991 | 10.26991 | 3.99 | 10.26991
18.16577 | 18.16512 | 18.16508 | 18.16508 | 3.98 | 18.16507
32.75676 | 32.75285 | 32.75264 | 32.75260 | 3.97 | 32.75258
3 | 7.91295 | 7.91295 | 7.91295 | 7.91295 | 5.92 | 7.91295
10.26991 | 10.26991 | 10.26991 | 10.26991 | 5.95 | 10.26991
18.16508 | 18.16507 | 18.16507 | 18.16507 | 5.98 | 18.16507
32.75260 | 32.75258 | 32.75258 | 32.75258 | 5.98 | 32.75258

In Table 1 we note that spurious eigenvalues appear in the lowest order
approximation for $\rho_{1}$ since $2(\overline{\rho}+\underline{\rho})<4$,
which is expected from the results in [20]. However, for $\rho_{2}$ we have a
different scenario since $2(\overline{\rho}+\underline{\rho})\approx 21.7689$,
but a spurious eigenvalue appears in the spectrum. Also, for this
configuration of the problem, when the polynomial degree is increased, the
appearance of spurious eigenvalues is less frequent. This fact allows us to
infer that when the polynomial approximation increases, the spurious
eigenvalues begin to vanish, at least when the lowest frequencies are
considered.

From the above, one may be tempted to infer that the safe parameter for the
SIP scheme is $\texttt{a}\overline{\rho}$, for $\texttt{a}\geq 4$. However, in
Table 2 we observe an alternative spurious analysis, where we consider, in
particular, stabilization parameters of the references [19, 20], together with
a particular combination of the form $\texttt{a}(\overline{\rho}+1)$. Note the
remarkable difference in the cleanliness of the spectrum between the schemes
with $\rho_{1}$ and $\rho_{2}$. In particular, the scheme with $\rho_{2}$ and
$4(\overline{\rho}+1)$, whose stabilized is larger compared to
$4\overline{\rho}$, presented in the previous table, gives a spurious
eigenvalue for $k=1$. Moreover, we note that the safe parameters in the
references are not sufficient if $|\rho(x,y)|>1$. Therefore, we infer that the
safe parameter choice for the SIP with variable density lies between the
combinations $\texttt{a}(\overline{\rho}+\underline{\rho})/2$ and
$\texttt{a}(\overline{\rho}+1)$, for $\texttt{a}\geq 8$, at least for the
cases under study.

On the other hand, Tables 3–4 reveal that with the SIP method the optimal
order of convergence is attained. More precisely, the order is exactly
$\mathcal{O}(h^{2k})$. This order is obtained since the domain in which we
state the problem is convex, leading to sufficiently smooth eigenfunctions for
our spectral problem.

####  5.1.2 NIP and IIP method

Now our aim is to study computationally the presence of spurious eigenvalues
for the nonsymmetric methods, more precisely the NIP ($\varepsilon=-1$) and
SIP ($\varepsilon=0$). The combination of stabilization parameters are the
same as those of SIP. In Table 5 we present the computed eigenvalues for
different stabilizations and polynomial degrees when the NIP method is
considered.

Table 5: Test 1. Computed eigenvalues $\lambda_{h,i}$ for polynomial degrees $k=1,2,3$, mesh level $N=8$ and different values of a and $\rho(x)$, for the NIP ($\varepsilon=-1$) method. | $\rho_{1}(x,y)$ | $\rho_{2}(x,y)$
---|---|---
k a | $\overline{\rho}+\underline{\rho}$ | $2\overline{\rho}$ | $2(\overline{\rho}+\underline{\rho})$ | $4\overline{\rho}$ | $\overline{\rho}+\underline{\rho}$ | $2\overline{\rho}$ | $2(\overline{\rho}+\underline{\rho})$ | $4\overline{\rho}$
1 | 7.74084 | 7.79239 | 7.81612 | 7.84287 | 7.89912 | 7.92612 | 7.93896 | 7.95100
9.42728 | 9.50518 | 9.54304 | 9.58812 | 10.16304 | 10.22912 | 10.26408 | 10.30103
16.84491 | 17.11841 | 17.24596 | 17.39284 | 17.90370 | 18.07698 | 18.16537 | 18.25622
30.37398 | 31.33856 | 31.79499 | 32.32444 | 31.60803 | 32.30451 | 32.66855 | 33.05071
36.67343 | 37.91183 | 38.46944 | 39.10970 | 38.24068 | 39.06090 | 39.49054 | 39.94433
38.21148 | 39.99929 | 40.82734 | 41.74630 | 40.64388 | 41.81255 | 42.38055 | 42.94331
43.61388 | 45.13923 | 45.90425 | 46.84143 | 45.58253 | 46.85679 | 47.56141 | 48.33365
59.21923 | 65.54764 | 68.22675 | 70.98606 | 65.44416 | 69.58263 | 71.41674 | 73.13639
63.44700 | 68.01181 | 70.23721 | 72.89916 | 67.67393 | 71.02756 | 72.85119 | 74.81097
70.50419 | 75.99282 | 78.72432 | 82.06519 | 75.21822 | 79.64399 | 82.07951 | 84.72267
2 | 7.84646 | 7.84285 | 7.84091 | 7.83847 | 7.92746 | 7.92356 | 7.92134 | 7.91889
9.61147 | 9.60681 | 9.60426 | 9.60101 | 10.29693 | 10.29000 | 10.28596 | 10.28140
17.49628 | 17.47764 | 17.46744 | 17.45441 | 18.24717 | 18.22656 | 18.21449 | 18.20084
32.64374 | 32.56841 | 32.52784 | 32.47677 | 33.02521 | 32.95571 | 32.91572 | 32.87114
39.53666 | 39.44163 | 39.39020 | 39.32520 | 39.94652 | 39.85876 | 39.80801 | 39.75122
42.36415 | 42.24515 | 42.18022 | 42.09761 | 42.95950 | 42.85318 | 42.79122 | 42.72144
47.62310 | 47.48150 | 47.40405 | 47.30527 | 48.54451 | 48.40647 | 48.32566 | 48.23419
72.84790 | 72.52868 | 72.35369 | 72.13023 | 73.37488 | 73.09611 | 72.93325 | 72.74945
74.68225 | 74.30735 | 74.10584 | 73.85231 | 74.84721 | 74.52393 | 74.33815 | 74.13114
84.65230 | 84.19038 | 83.93873 | 83.61867 | 85.32188 | 84.90646 | 84.66393 | 84.39006
3 | 7.83256 | 7.83255 | 7.83255 | 7.83254 | 7.91298 | 7.91297 | 7.91296 | 7.91296
9.59293 | 9.59291 | 9.59291 | 9.59290 | 10.26997 | 10.26995 | 10.26994 | 10.26993
17.42105 | 17.42093 | 17.42087 | 17.42079 | 18.16550 | 18.16538 | 18.16531 | 18.16524
32.34726 | 32.34647 | 32.34607 | 32.34556 | 32.75509 | 32.75440 | 32.75402 | 32.75360
39.15956 | 39.15856 | 39.15803 | 39.15738 | 39.60313 | 39.60225 | 39.60176 | 39.60122
41.87933 | 41.87756 | 41.87662 | 41.87545 | 42.53310 | 42.53165 | 42.53083 | 42.52994
47.04380 | 47.04165 | 47.04051 | 47.03908 | 47.98422 | 47.98215 | 47.98097 | 47.97967
71.50754 | 71.49928 | 71.49490 | 71.48939 | 72.21964 | 72.21292 | 72.20910 | 72.20487
73.16741 | 73.15827 | 73.15349 | 73.14759 | 73.54462 | 73.53705 | 73.53280 | 73.52813
82.72109 | 82.70733 | 82.70006 | 82.69098 | 83.58124 | 83.56921 | 83.56237 | 83.55478
Table 6: Test 1. Computed eigenvalues $\lambda_{h,i}$ for polynomial degrees $k=1,2,3$, mesh level $N=8$ and different values of a and $\rho(x)$, for the NIP ($\varepsilon=-1$) method. | $\rho_{1}(x,y)$ | $\rho_{2}(x,y)$
---|---|---
k a | 2 | $2(\overline{\rho}+1)$ | 4 | $4(\overline{\rho}+1)$ | 2 | $2(\overline{\rho}+1)$ | 4 | $4(\overline{\rho}+1)$
1 | 7.79239 | 7.84287 | 7.84287 | 7.86916 | 7.50961 | 7.93180 | 7.751918 | 7.95346
9.50518 | 9.58812 | 9.58812 | 9.63647 | 9.41399 | 10.24419 | 9.85755 | 10.30940
17.11841 | 17.3928 | 17.3928 | 17.54436 | 15.58637 | 18.11537 | 17.02452 | 18.27644
31.33856 | 32.32444 | 32.32444 | 32.87046 | 20.20981 | 32.46180 | 28.30139 | 33.13710
37.91183 | 39.10970 | 39.10970 | 39.77174 | 22.88114 | 39.24633 | 33.68424 | 40.04751
39.99929 | 41.74630 | 41.74630 | 42.63479 | 24.96308 | 42.06174 | 34.82433 | 43.06573
45.13923 | 46.84143 | 46.84143 | 47.87919 | 26.94440 | 47.15775 | 38.22918 | 48.51337
65.54764 | 70.98606 | 70.98606 | 73.47962 | 28.69336 | 70.40197 | 40.38647 | 73.50271
68.01180 | 72.89916 | 72.89916 | 75.69995 | 29.71622 | 71.80870 | 45.80849 | 75.25708
75.99282 | 82.06519 | 82.06519 | 85.68356 | 30.54102 | 80.68556 | 52.33948 | 85.33071
2 | 7.84285 | 7.83847 | 7.83847 | 7.83579 | 7.11555 | 7.92261 | 7.94064 | 7.94064
9.60680 | 9.60101 | 9.60101 | 9.59739 | 7.93392 | 10.28829 | 10.31911 | 10.28033
17.47764 | 17.45441 | 17.45441 | 17.43976 | 10.27103 | 18.22145 | 18.31268 | 18.19761
32.56841 | 32.47677 | 32.47677 | 32.42035 | 18.17500 | 32.93874 | 33.25550 | 32.86073
39.44163 | 39.32520 | 39.32520 | 39.25305 | 32.71864 | 39.83724 | 40.23387 | 39.73792
42.24515 | 42.09761 | 42.09761 | 42.00507 | 39.60958 | 42.82695 | 43.30082 | 42.70503
47.48150 | 47.30527 | 47.30527 | 47.19430 | 42.28333 | 48.37230 | 48.98499 | 48.21260
72.52868 | 72.13023 | 72.13023 | 71.87892 | 47.90767 | 73.02722 | 74.26323 | 72.70616
74.30735 | 73.85231 | 73.85231 | 73.57203 | 50.66857 | 74.44508 | 75.92709 | 74.08280
84.19038 | 83.61866 | 83.61866 | 83.26009 | 71.03522 | 84.80387 | 86.65687 | 84.32550
3 | 7.83255 | 7.83254 | 7.83254 | 7.83253 | 7.91308 | 7.91297 | 7.91302 | 7.91296
9.59291 | 9.59290 | 9.59290 | 9.59289 | 10.27009 | 10.26995 | 10.27003 | 10.26993
17.42093 | 17.42079 | 17.42079 | 17.42070 | 18.16644 | 18.16535 | 18.16597 | 18.16523
32.34647 | 32.34556 | 32.34556 | 32.34501 | 32.76030 | 32.75424 | 32.75768 | 32.75350
39.15855 | 39.15738 | 39.15738 | 39.15668 | 39.60958 | 39.60204 | 39.60637 | 39.60110
41.87755 | 41.87545 | 41.87545 | 41.87417 | 42.54384 | 42.53130 | 42.53849 | 42.52973
47.04165 | 47.03908 | 47.03908 | 47.03749 | 47.99840 | 47.98165 | 47.99156 | 47.97936
71.49928 | 71.48939 | 71.48939 | 71.48330 | 72.26535 | 72.21129 | 72.24347 | 72.20388
73.15827 | 73.14759 | 73.14759 | 73.14115 | 73.59971 | 73.53524 | 73.57238 | 73.52705
82.70733 | 82.69098 | 82.69098 | 82.68094 | 83.66459 | 83.56630 | 83.62420 | 83.55301

We observe that when the NIP method is performing the approximation of the
spectrum with the selected stabilization parameters, there are no spurious
eigenvalues arising for small values of a. For $\rho_{1}$, since
$\overline{\rho}\geq 1$, we fall in the results presented in [20], whereas the
scheme with density $\rho_{2}$ remains clean with $\overline{\rho}>2$.

In contrast, we observe the results in Table 6, where a completely clean
spectrum is observed for $\rho_{1}$, while the spectrum with $\rho_{2}$ and
$\texttt{a}=2,4,$ result in the computation of spurious eigenvalues. When we
compare this behavior with the SIP method, we notice that the NIP method
introduces more spurious modes compared with the SIP method, for the same
stabilization parameter. However, when the polynomial degree is increased, the
NIP methods behaves better than the SIP method where for $\texttt{a}\geq 4$

Let us remark that when the NIP method is considered, Table 6 shows that for
$k=1$ the number of spurious eigenvalues is clearly large compared for $k=2$.
Moreover, for $k=3$ there is no visible spurious eigenvalues for any
stabilization parameter. This has been also observed for $k=4,5,6$. From these
results, we infer that a safe parameter in which we can operate with the NIP
method is for stabilizations such that $\texttt{a}\geq 4$.

Now in Tables 7–8 we present convergence rates for the NIP method when
different polynomial degrees are considered. From these tables, we observe
that a suboptimal order of convergence is attained. More precisely, when $k$
is odd, the order of convergence is $\mathcal{O}(h^{k})$, whereas when $k$ is
even, the order of convergence is $\mathcal{O}(h^{k+1})$, which clearly is not
optimal compared with the SIP method.

We remark that for the IIP method ($\varepsilon=0$), the results related to
the spurious analysis and order of convergence are similar to the NIP method.
For the sake of brevity, we skip the results.

Table 7: Test 1. Convergence analysis for $\texttt{a}=10(\overline{\rho}+1)$, $\rho(x,y)=\frac{1}{x^{2}+y^{2}+1}$, and $k=1,2,3$ for the NIP method ($\varepsilon=-1$).  k | $N=10$ | $N=20$ | $N=30$ | $N=40$ | Order | $\lambda_{i}$
---|---|---|---|---|---|---
1 | 7.86660 | 7.84112 | 7.83636 | 7.83468 | 1.99 | 7.83253
9.64235 | 9.60531 | 9.59842 | 9.59600 | 1.99 | 9.59286
17.56622 | 17.45746 | 17.43704 | 17.42986 | 1.98 | 17.42050
32.91484 | 32.48800 | 32.40836 | 32.38040 | 1.99 | 32.34428
2 | 7.83341 | 7.83273 | 7.83262 | 7.83258 | 2.16 | 7.83253
9.59410 | 9.59316 | 9.59300 | 9.59295 | 2.16 | 9.59289
17.42594 | 17.42177 | 17.42110 | 17.42087 | 2.23 | 17.42064
32.36612 | 32.34895 | 32.34630 | 32.34542 | 2.29 | 32.34454
3 | 7.83253 | 7.83253 | 7.83253 | 7.83253 | 4.10 | 7.83253
9.59288 | 9.59288 | 9.59288 | 9.59288 | 4.08 | 9.59288
17.42061 | 17.42060 | 17.42060 | 17.42060 | 4.14 | 17.42060
32.34446 | 32.34434 | 32.34434 | 32.34434 | 4.15 | 32.34434
Table 8: Test 1. Convergence analysis for $\texttt{a}=10(\overline{\rho}+1)$, $\rho(x,y)=e^{xy+1}$, and $k=1,2,3$ for the NIP method ($\varepsilon=-1$).  k | $N=10$ | $N=20$ | $N=30$ | $N=40$ | Order | $\lambda_{i}$
---|---|---|---|---|---|---
1 | 7.94590 | 7.92131 | 7.91668 | 7.91505 | 1.98 | 7.91294
10.32299 | 10.28334 | 10.27591 | 10.27329 | 1.98 | 10.26988
18.30356 | 18.20050 | 18.18092 | 18.17402 | 1.96 | 18.16484
33.28860 | 32.88857 | 32.81330 | 32.78681 | 1.97 | 32.75168
2 | 7.91440 | 7.91328 | 7.91309 | 7.91303 | 2.16 | 7.91296
10.27278 | 10.27057 | 10.27020 | 10.27007 | 2.14 | 10.26993
18.17420 | 18.16711 | 18.16595 | 18.16556 | 2.19 | 18.16513
32.78360 | 32.75927 | 32.75544 | 32.75417 | 2.25 | 32.75282
3 | 7.91295 | 7.91295 | 7.91295 | 7.91295 | 4.25 | 7.91295
10.26992 | 10.26991 | 10.26991 | 10.26991 | 4.20 | 10.26991
18.16510 | 18.16508 | 18.16507 | 18.16507 | 4.22 | 18.16507
32.75274 | 32.75259 | 32.75258 | 32.75258 | 4.17 | 32.75258

###  5.2 Test 2. A 2D reactor

This experiment aims to show the performance of the scheme when there is a
non-polygonal domain and the density $\rho(x)$ is non constant. This implies
that a priori we have zones within the domain where the fluid density is
higher or lower. The considered domain is given by
$\Omega:=\widetilde{\Omega}\backslash\bigcup_{i=1}^{4}\Omega_{i}$, where

$\displaystyle\widetilde{\Omega}$
$\displaystyle:=\\{(x,y)\in\mathbb{R}^{2}\,:\,x^{2}+y^{2}\leq 1\\},$
$\displaystyle\Omega_{i}$
$\displaystyle:=\left\\{(x,y)\in\mathbb{R}^{2}\,:\,\left(x+\frac{1}{3}a\right)^{2}+\left(y+\frac{1}{3}b\right)^{2}<\frac{1}{8}\right\\},$

with $a=(-1)^{i}$ and $b=(-1)^{\frac{i(i-1)}{2}+1}$, for $i=1,2,3,4$. This
domain resembles a cross-section of a reactor. The domain is meshed such that
$h\approx 1/2^{N}$, with $N=3,4,5,6$. For a clean spectrum, we choose the
stabilization parameter $\texttt{a}=20$. It notes that the expected order of
convergence of the eigenvalues is $\mathcal{O}(h^{2})$. This is due to the
variational crime that is committed when a curved domain is approximated by
polygons. This fact has been also observed in [20].

In Figure 1 we report the experimental convergence obtained by least squares
fitting, solving the problem using the SIP and NIP scheme. Here, dof denotes
the number of degrees of freedom, whereas the error on the $i$-th eigenvalue
is denoted by $\texttt{err}(\lambda_{i})$, with

$\texttt{err}(\lambda_{i}):=|\lambda_{h_{i}}-\lambda_{i}|.$

For the IIP method, the values obtained were similar to those obtained with
NIP. Note that in all cases, the order of convergence remains as predicted
theoretically, for all selected fluid densities. This is followed by Figure 2,
where we show the eigenfunctions associated to the first eigenvalue for each
density function. Note that the fluid particles tend to move faster in the
lower density zones, as expected. It is worth noting that, when using the NIP
method, the spectrum may contain complex values, so we have only reported the
real part of the spectrum. These complex contributions can appear at any place
and any refinement. For example, for $\rho(x,y)=1/[\cos(\pi x)\cos(\pi y)+2]$
and $N=3$, the first two computed eigenvalues are of the form
$\lambda_{h}=3.47649\pm 0.00009i$.

$\rho(x,y)=e^{\sin\pi x\cos\pi y}$

$\rho(x,y)=e^{\sin\left(\frac{\pi x}{2}\right)\cos\left(\frac{\pi
y}{2}\right)}$

$\rho(x,y)=\frac{1}{\cos(\pi x)\cos(\pi y)+2}$

$\rho(x,y)=e^{xy+1}$

Figure 1: Test 2. Error curves for the first four lowest computed eigenvalues
when using different density functions.

$\boldsymbol{u}_{1h}$

$p_{1h}$

$\rho(x,y)=e^{\sin\pi x\cos\pi y}$

$\boldsymbol{u}_{1h}$

$p_{1h}$

$\rho(x,y)=e^{\sin\left(\frac{\pi x}{2}\right)\cos\left(\frac{\pi
y}{2}\right)}$

$\boldsymbol{u}_{1h}$

$p_{1h}$

$\rho(x,y)=\frac{1}{\cos(\pi x)\cos(\pi y)+2}$

$\boldsymbol{u}_{1h}$

$p_{1h}$

$\rho(x,y)=e^{xy+1}$

Figure 2: Test 2. Lowest computed eigenmodes and postprocessed pressure using
different density functions.

###  5.3 Test 3: A 3D torus

This test aims to assess the DG method on a three dimensional non-polygonal
domain. We consider a torus azimuthally symmetric about the $z-$axis, whose
domain is characterized by

$\Omega:=\left\\{(x,y,z)\in\mathbb{R}^{3}\;:\;\left(\sqrt{x^{2}+y^{2}}-R\right)^{2}+z^{2}=r^{2}\right\\},$

where $R=1/2$ and $r=1/4$. For simplicity, the physical parameters are taken
equal to $1$. In Table 9 we observe that for the lowest order, the variants of
our discontinuous scheme converge with quadratic order, as the theory
predicts. This is due to the fact that the eigenfunctions in this domain are
regular. However, a variational crime is also being committed when
approximating this domain with polyhedrons, so this order will not to improve
for $k>1$. Finally, in Figures 3 and 4 we depict the velocity field along with
the post-processed pressure for the third and fifth eigenmodes.

Table 9: Test 3. Lowest computed eigenvalues for $k$$\,=\,$$1$, and $\texttt{a}=10$ for the SIP, NIP and IIP method.  method | $N=15$ | $N=20$ | $N=25$ | $N=30$ | Order | $\lambda_{extr}$
---|---|---|---|---|---|---
| 4.20546 | 4.20096 | 4.19911 | 4.19778 | 1.97 | 4.19504
| 4.20588 | 4.20127 | 4.19923 | 4.19794 | 1.94 | 4.19497
SIP | 15.90982 | 15.87822 | 15.86392 | 15.85563 | 2.00 | 15.83625
15.91094 | 15.87873 | 15.86493 | 15.85597 | 2.00 | 15.83658
| 33.31118 | 33.20592 | 33.15811 | 33.12984 | 1.96 | 33.06265
| 4.20528 | 4.20086 | 4.19904 | 4.19774 | 1.98 | 4.19506
| 4.20571 | 4.20117 | 4.19917 | 4.19789 | 1.94 | 4.19497
NIP | 15.90732 | 15.87681 | 15.86295 | 15.85497 | 1.99 | 15.83604
15.90853 | 15.87728 | 15.86398 | 15.85527 | 2.01 | 15.83669
| 33.30065 | 33.19987 | 33.15391 | 33.12693 | 1.96 | 33.06252
| 4.20535 | 4.20090 | 4.19907 | 4.19775 | 1.97 | 4.19503
| 4.20577 | 4.20121 | 4.19920 | 4.19791 | 1.94 | 4.19497
IIP | 15.90892 | 15.87736 | 15.86333 | 15.85523 | 2.00 | 15.83624
15.90946 | 15.87784 | 15.86435 | 15.85554 | 2.00 | 15.83653
| 33.30471 | 33.20221 | 33.15554 | 33.12806 | 1.96 | 33.06259

Figure 3: Test 3. Third lowest computed eigenmode $\boldsymbol{u}_{h}$
streamlines (left) and postprocessed pressure (right).

Figure 4: Test 3. Fifth lowest computed eigenmode $\boldsymbol{u}_{h}$
streamlines (left) and postprocessed pressure (right).

##  6 Test 4. The pressure formulation

Since $\varepsilon\in\\{-1,0,1\\}$ as in the displacement formulation, for the
pressure formulation we analyze the spurious behavior when the DG method is
considered. The domain, physical parameters and eigenvalues are the same as
those of Test 1. In Tables 10–11 we report the computed eigenvalues obtained
with the SIP method, fixing the mesh refinement in $N=8$ and considering
different polynomial degrees and values of a. We observe that the results are
similar compared with those on Tables 1–2. However, for the smallest
stabilization parameters, the pressure formulation gives more spurious
eigenvalues when the polynomial degree is $k=1$.

Table 10: Test 4. Computed eigenvalues $\lambda_{h,i}$ for polynomial degrees $k=1,2,3$, mesh level $N=8$ and different values of a and $\rho(x)$, for the SIP method ($\varepsilon=1$). | $\rho_{1}(x,y)$ | $\rho_{2}(x,y)$
---|---|---
k a | $2(\overline{\rho}+\underline{\rho})$ | $4\overline{\rho}$ | $4(\overline{\rho}+\underline{\rho})$ | $8\overline{\rho}$ | $2(\overline{\rho}+\underline{\rho})$ | $4\overline{\rho}$ | $4(\overline{\rho}+\underline{\rho})$ | $8\overline{\rho}$
1 | 7.85831 | 7.87585 | 6.28236 | 7.89707 | 8.00791 | 8.00878 | 8.00921 | 8.00965
9.61673 | 9.64631 | 7.99171 | 9.66568 | 10.39198 | 10.39333 | 10.39401 | 10.39470
11.12628 | 14.56096 | 9.65834 | 17.78279 | 18.81280 | 18.82297 | 18.82815 | 18.83341
17.58573 | 17.70166 | 17.72448 | 33.38125 | 34.38199 | 34.39806 | 34.40617 | 34.41436
32.74643 | 33.05492 | 33.12913 | 40.25775 | 41.54168 | 41.56437 | 41.57583 | 41.58741
35.66546 | 37.81773 | 40.04185 | 43.71383 | 45.60117 | 45.64630 | 45.66926 | 45.69258
39.50152 | 40.23103 | 42.42326 | 49.10155 | 52.15709 | 52.22394 | 52.25799 | 52.29260
42.56402 | 43.36561 | 46.34274 | 76.38572 | 81.30052 | 81.44241 | 81.51355 | 81.58492
47.59880 | 48.28572 | 51.09714 | 78.06806 | 81.90845 | 82.01209 | 82.06557 | 82.12053
57.83388 | 56.48285 | 72.73759 | 89.46875 | 96.68894 | 96.89340 | 96.99754 | 97.10336
2 | 7.83258 | 7.83261 | 7.83264 | 7.83269 | 7.91319 | 7.91319 | 7.91320 | 7.91320
9.59304 | 9.59295 | 9.59300 | 9.59306 | 10.27022 | 10.27022 | 10.27023 | 10.27023
17.42250 | 17.42202 | 17.42253 | 17.42303 | 18.16847 | 18.16851 | 18.16853 | 18.16854
32.34518 | 32.34930 | 32.35180 | 32.35447 | 32.76848 | 32.76863 | 32.76870 | 32.76877
39.14741 | 39.15855 | 39.16265 | 39.16658 | 39.61886 | 39.61907 | 39.61917 | 39.61928
41.88102 | 41.89014 | 41.89664 | 41.90275 | 42.56796 | 42.56834 | 42.56853 | 42.56872
47.02938 | 47.05119 | 47.06061 | 47.06864 | 48.03134 | 48.03189 | 48.03216 | 48.03244
65.11345 | 71.55640 | 71.59582 | 71.62566 | 72.39744 | 72.39937 | 72.40035 | 72.40134
71.73920 | 73.19101 | 73.21856 | 73.24738 | 73.69022 | 73.69170 | 73.69245 | 73.69320
73.17229 | 82.78775 | 82.83178 | 82.87428 | 83.84312 | 83.84562 | 83.84688 | 83.84816
3 | 7.83253 | 7.83253 | 7.83253 | 7.83253 | 7.91295 | 7.91295 | 7.91295 | 7.91295
9.59288 | 9.59288 | 9.59288 | 9.59288 | 10.26991 | 10.26991 | 10.26991 | 10.26991
17.42060 | 17.42060 | 17.42060 | 17.42060 | 18.16508 | 18.16508 | 18.16508 | 18.16508
32.34439 | 32.34439 | 32.34440 | 32.34440 | 32.75265 | 32.75265 | 32.75265 | 32.75265
39.15585 | 39.15586 | 39.15587 | 39.15587 | 39.60000 | 39.60000 | 39.60000 | 39.60000
41.87267 | 41.87270 | 41.87272 | 41.87274 | 42.52793 | 42.52793 | 42.52793 | 42.52793
47.03560 | 47.03564 | 47.03566 | 47.03568 | 47.97662 | 47.97662 | 47.97662 | 47.97662
71.47621 | 71.47651 | 71.47664 | 71.47679 | 72.19550 | 72.19551 | 72.19551 | 72.19551
73.13362 | 73.13384 | 73.13394 | 73.13406 | 73.51765 | 73.51766 | 73.51766 | 73.51766
82.66900 | 82.66951 | 82.66975 | 82.67002 | 83.53781 | 83.53783 | 83.53784 | 83.53784
Table 11: Test 4. Computed eigenvalues $\lambda_{h,i}$ for polynomial degrees $k=1,2,3$, mesh level $N=8$ and different values of a and $\rho(x)$, for the SIP method ($\varepsilon=1$). | $\rho_{1}(x,y)$ | $\rho_{2}(x,y)$
---|---|---
k a | 4 | $4(\overline{\rho}+1)$ | 8 | $8(\overline{\rho}+1)$ | 4 | $4(\overline{\rho}+1)$ | 8 | $8(\overline{\rho}+1)$
1 | 7.87585 | 7.89707 | 7.89707 | 7.91503 | 7.99747 | 8.00897 | 8.00365 | 8.00975
9.64631 | 9.66568 | 9.66568 | 9.68920 | 10.37575 | 10.39362 | 10.38534 | 10.39485
14.56096 | 17.78279 | 17.78279 | 17.89578 | 18.69974 | 18.82522 | 18.76446 | 18.83457
17.70166 | 33.38125 | 33.38125 | 33.66847 | 34.18685 | 34.40160 | 34.30248 | 34.41615
33.05492 | 40.25775 | 40.25775 | 40.66849 | 41.26919 | 41.56936 | 41.43004 | 41.58995
37.81773 | 43.71383 | 43.71383 | 44.32293 | 45.09283 | 45.65629 | 45.38552 | 45.69770
40.23103 | 49.10155 | 49.10155 | 49.89565 | 51.41133 | 52.23875 | 51.83909 | 52.30021
43.36561 | 76.38572 | 76.38572 | 78.37028 | 79.63021 | 81.47347 | 80.59982 | 81.60048
48.28572 | 78.06806 | 78.06806 | 79.75233 | 80.75383 | 82.03530 | 81.42428 | 82.13270
56.48285 | 89.46875 | 89.46875 | 92.09170 | 94.40090 | 96.93872 | 95.71525 | 97.12663
2 | 7.83261 | 7.83269 | 7.83269 | 7.83273 | 7.91316 | 7.91320 | 7.91318 | 7.91320
9.59295 | 9.59306 | 9.59306 | 9.59312 | 10.27018 | 10.27022 | 10.27020 | 10.27023
17.42202 | 17.42303 | 17.42303 | 17.42357 | 18.16806 | 18.16852 | 18.16830 | 18.16855
32.34930 | 32.35447 | 32.35447 | 32.35720 | 32.76667 | 32.76866 | 32.76775 | 32.76879
39.15855 | 39.16658 | 39.16658 | 39.17049 | 39.61629 | 39.61911 | 39.61782 | 39.61930
41.89014 | 41.90275 | 41.90275 | 41.90920 | 42.56350 | 42.56842 | 42.56612 | 42.56876
47.05119 | 47.06864 | 47.06864 | 47.07696 | 48.02489 | 48.03201 | 48.02868 | 48.03250
71.55640 | 71.62566 | 71.62566 | 71.65676 | 72.37484 | 72.39979 | 72.38805 | 72.40155
73.19101 | 73.24738 | 73.24738 | 73.27627 | 73.67168 | 73.69203 | 73.68279 | 73.69336
82.78775 | 82.87428 | 82.87428 | 82.91733 | 83.81269 | 83.84617 | 83.83073 | 83.84844
3 | 7.83253 | 7.83253 | 7.83253 | 7.83253 | 7.91295 | 7.91295 | 7.91295 | 7.91295
9.59288 | 9.59288 | 9.59288 | 9.59288 | 10.26991 | 10.26991 | 10.26991 | 10.26991
17.42060 | 17.42060 | 17.42060 | 17.42060 | 18.16508 | 18.16508 | 18.16508 | 18.16508
32.34439 | 32.34440 | 32.34440 | 32.34441 | 32.75264 | 32.75265 | 32.75264 | 32.75265
39.15586 | 39.15587 | 39.15587 | 39.15588 | 39.60000 | 39.60000 | 39.60000 | 39.60000
41.87270 | 41.87274 | 41.87274 | 41.87275 | 42.52792 | 42.52793 | 42.52793 | 42.52793
47.03564 | 47.03568 | 47.03568 | 47.03570 | 47.97660 | 47.97662 | 47.97661 | 47.97662
71.47651 | 71.47679 | 71.47679 | 71.47695 | 72.19541 | 72.19551 | 72.19546 | 72.19551
73.13384 | 73.13406 | 73.13406 | 73.13419 | 73.51757 | 73.51766 | 73.51762 | 73.51766
82.66951 | 82.67002 | 82.67002 | 82.67031 | 83.53763 | 83.53783 | 83.53774 | 83.53784

##  7 Test 5. A 3D benchmark

In order to obtain more information of the discrete formulations (3.6) and
(4.12), now we compare some aspects of the computational cost of these
schemes. To achieve this goal, we resort to the debug report delivered by
FENICS through the command set_log_level(LogLevel.DEBUG), where we can recover
several debug information, among them, the assembly time, the non-zero entries
of the assembled matrices and the sparsity percentage. We solve the
eigenproblems in the unit cube $\Omega:=(0,1)^{2}$ using uniform refinements,
where the number of cells scales as $6(N+1)^{3}$. It is clear that for the
cases studied, the generalized eigenvalue system is of the form
$\boldsymbol{K_{u}}\boldsymbol{u}=\lambda\boldsymbol{M_{u}}\boldsymbol{u}$ for
(3.6) and $\boldsymbol{K_{p}}p=\lambda\boldsymbol{M_{p}}p$ for (4.12). In this
experiment, we denote by $K_{0,\boldsymbol{u}}$ and $M_{0,\boldsymbol{u}}$ the
number of non-zero entries of $\boldsymbol{K_{u}}$ and $\boldsymbol{M_{u}}$ in
each refinement. Similarly, we denote by $K_{0,p}$ and $M_{0,p}$ the number of
non-zero entries of $\boldsymbol{K_{p}}$ and $\boldsymbol{M_{p}}$. In Figure 5
can observe the behavior of the assembly and solution times, in seconds, of
the eigenvalue system in both cases, taken as the average of 5 continuous
runs. It is clear that solving (4.12) is considerably cheaper, mainly because
of the number of degrees of freedom (dof) in the problem. This is followed by
the results depicted in Figure 6, where we observe the number of non-zero
inputs with respect to the number of degrees of freedom, and also the ratio
between the non-zero inputs and dof2. The $y-$axis scale is measured with the
logarithm of the computed values. Note that the ratio between non-zero entries
and $\texttt{dof}^{2}$ of the matrices associated with (4.12) are lower than
those of (3.6), which is expected since (4.12) is associated with fewer
degrees of freedoms.

Figure 5: Test 5. Comparison of time employed by the velocity and pressure
schemes when assembling and solving the eigenproblems in a 3D unit cube with
uniform refinements. Left: velocity formulation. Right: pure pressure
formulation.

Figure 6: Test 5. Comparison of the sparsity in the the velocity and pressure
schemes on each mesh level for the 3D unit cube with uniform refinements.
Left: velocity formulation. Right: pure pressure formulation.

##  8 Conclusions

We have presented DG methods to approximate the eigenvalues and eigenfunctions
of the acoustic eigenvalue problem. The proposed symmetric and nonsymmetric
methods lead to an accurate computation and approximation of the corresponding
eigenfunctions. The acoustic problem has as unknowns the displacement and
pressure of the fluid. These DG methods were considered for two formulations:
a formulation where the displacement is the only unknown and other where the
pressure is the unknown. For both approaches, the remaining unknown is always
recovered with a postprocess of the solution operator under study. A spurious
analysis is presented in two dimensions, where we compare the effects of the
stabilization parameter in the symmetric and nonsymmetric methods when the
spectrum es computed, performing both the displacement and pressure
formulations. The numerical tests also reveal that, for a correct choice of
the aforementioned stabilization, the convergence order is the expected in two
and three dimensions. On the other hand, a benchmark test shows that the
pressure formulation is less expensive in terms of computational costs
compared with the displacement formulation. This conclusion is expectable
since the pure pressure problem is only scalar. However, in both formulations,
the eigenvalues are exactly the same, together with the computed convergence
orders for both, the symmetric or nonsymmetric approaches.

## References

  * [1] P. F. Antonietti, A. Buffa, and I. Perugia, Discontinuous Galerkin approximation of the Laplace eigenproblem, Comput. Methods Appl. Mech. Engrg., 195 (2006), pp. 3483–3503, https://doi.org/10.1016/j.cma.2005.06.023.
  * [2] I. Babuška and J. Osborn, Eigenvalue problems, in Handbook of numerical analysis, Vol. II, Handb. Numer. Anal., II, North-Holland, Amsterdam, 1991, pp. 641–787.
  * [3] Y. Bai and Q. Bai, Subsea pipeline integrity and risk management, Gulf Professional Publishing, 2014.
  * [4] L. Beirão da Veiga, D. Mora, G. Rivera, and R. Rodríguez, A virtual element method for the acoustic vibration problem, Numer. Math., 136 (2017), pp. 725–763, https://doi.org/10.1007/s00211-016-0855-5.
  * [5] A. Bermúdez, R. Durán, M. A. Muschietti, R. Rodríguez, and J. Solomin, Finite element vibration analysis of fluid-solid systems without spurious modes, SIAM J. Numer. Anal., 32 (1995), pp. 1280–1295, https://doi.org/10.1137/0732059.
  * [6] A. Bermúdez, P. Gamallo, and R. Rodríguez, Finite element methods in local active control of sound, SIAM J. Control Optim., 43 (2004), pp. 437–465, https://doi.org/10.1137/S0363012903431785.
  * [7] A. Bermúdez, R. Rodríguez, and D. Santamarina, Two discretization schemes for a time-domain dissipative acoustics problem, Math. Models Methods Appl. Sci., 16 (2006), pp. 1559–1598, https://doi.org/10.1142/S0218202506001637.
  * [8] A. Buffa, P. Houston, and I. Perugia, Discontinuous Galerkin computation of the Maxwell eigenvalues on simplicial meshes, J. Comput. Appl. Math., 204 (2007), pp. 317–333, https://doi.org/10.1016/j.cam.2006.01.042.
  * [9] A. Buffa and I. Perugia, Discontinuous Galerkin approximation of the Maxwell eigenproblem, SIAM J. Numer. Anal., 44 (2006), pp. 2198–2226, https://doi.org/10.1137/050636887.
  * [10] A. Buffa, I. Perugia, and T. Warburton, The mortar-discontinuous Galerkin method for the 2D Maxwell eigenproblem, J. Sci. Comput., 40 (2009), pp. 86–114, https://doi.org/10.1007/s10915-008-9238-0.
  * [11] J. Descloux, N. Nassif, and J. Rappaz, On spectral approximation. I. The problem of convergence, RAIRO Anal. Numér., 12 (1978), pp. 97–112, iii, https://doi.org/10.1051/m2an/1978120200971.
  * [12] J. Descloux, N. Nassif, and J. Rappaz, On spectral approximation. II. Error estimates for the Galerkin method, RAIRO Anal. Numér., 12 (1978), pp. 113–119, iii, https://doi.org/10.1051/m2an/1978120201131.
  * [13] D. A. Di Pietro and A. Ern, Mathematical aspects of discontinuous Galerkin methods, vol. 69 of Mathématiques & Applications (Berlin) [Mathematics & Applications], Springer, Heidelberg, 2012, https://doi.org/10.1007/978-3-642-22980-0.
  * [14] E. Hernández and R. Rodriguez, Finite element approximation of spectral acoustic problems on curved domains, Numer. Math., 97 (2004), pp. 131–158, https://doi.org/10.1007/s00211-003-0501-x.
  * [15] T. Kirstein, Multidisciplinary know-how for smart-textiles developers, Elsevier, 2013.
  * [16] H. P. Langtangen and A. Logg, Solving PDEs in Python, vol. 3 of Simula SpringerBriefs on Computing, Springer, Cham, 2016, https://doi.org/10.1007/978-3-319-52462-7.  The FEniCS tutorial I.
  * [17] W. Larbi, J.-F. Deü, and R. Ohayon, A new finite element formulation for internal acoustic problems with dissipative walls, Internat. J. Numer. Methods Engrg., 68 (2006), pp. 381–399, https://doi.org/10.1002/nme.1727.
  * [18] F. Lepe, S. Meddahi, D. Mora, and R. Rodríguez, Acoustic vibration problem for dissipative fluids, Math. Comp., 88 (2019), pp. 45–71, https://doi.org/10.1090/mcom/3336.
  * [19] F. Lepe, S. Meddahi, D. Mora, and R. Rodríguez, Mixed discontinuous Galerkin approximation of the elasticity eigenproblem, Numer. Math., 142 (2019), pp. 749–786, https://doi.org/10.1007/s00211-019-01035-9.
  * [20] F. Lepe and D. Mora, Symmetric and nonsymmetric discontinuous Galerkin methods for a pseudostress formulation of the Stokes spectral problem, SIAM J. Sci. Comput., 42 (2020), pp. A698–A722, https://doi.org/10.1137/19M1259535.
  * [21] F. Lepe and G. Rivera, A virtual element approximation for the pseudostress formulation of the Stokes eigenvalue problem, Comput. Methods Appl. Mech. Engrg., 379 (2021), pp. Paper No. 113753, 21, https://doi.org/10.1016/j.cma.2021.113753.
  * [22] A. Márquez, S. Meddahi, and T. Tran, Analyses of mixed continuous and discontinuous Galerkin methods for the time harmonic elasticity problem with reduced symmetry, SIAM J. Sci. Comput., 37 (2015), pp. A1909–A1933, https://doi.org/10.1137/14099022X.
  * [23] D. Mora, G. Rivera, and R. Rodríguez, A virtual element method for the Steklov eigenvalue problem, Math. Models Methods Appl. Sci., 25 (2015), pp. 1421–1445, https://doi.org/10.1142/S0218202515500372.
  * [24] R. Osiander and A. Darrin, Wireless mems for space applications, in Wireless MEMS Networks and Applications, Elsevier, 2017, pp. 197–214.
  * [25] J. Querales, R. Rodríguez, and P. Venegas, Numerical approximation of the displacement formulation of the axisymmetric acoustic vibration problem, SIAM J. Sci. Comput., 43 (2021), pp. A1583–A1606, https://doi.org/10.1137/20M1346225.
  * [26] J. Querales and P. Venegas, Mixed approximation of the axisymmetric acoustic eigenvalue problem, Comput. Math. Appl., 108 (2022), pp. 1–11, https://doi.org/10.1016/j.camwa.2021.12.013.