# Two implicit parametrized discretization strategies for the mixed formulation of linear wave equations

Juliette Chabassier

► To cite this version:

Juliette Chabassier. Two implicit parametrized discretization strategies for the mixed formulation of linear wave equations. RR-9527, Inria. 2023, pp.18. hal-04262905

HAL Id: hal-04262905

https://inria.hal.science/hal-04262905v1

Submitted on 6 Nov 2023

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Ínria

Two implicit
parametrized
discretization strategies
for the mixed formulation
of linear wave equations

Juliette Chabassier

RESEARCH
REPORT

Nº 9527

October 2023

Project-Team MAKUTU
<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Two implicit parametrized discretization
strategies for the mixed formulation of linear
wave equations

Juliette Chabassier *

Project-Team MAKUTU

Research Report n° 9527 — October 2023 — 17 pages

**Abstract:** This report presents two parametric implicit strategies for discretizing the mixed formulation of linear wave equations, based on the first on the Störmer-Verlet scheme and for the second on the Crank-Nicolson scheme. If the initial conditions are well chosen, these two strategies are equivalent, and equivalence formulas allow the variables of one scheme to be reconstructed from the other. These two schemes are also equivalent to the θ-scheme discretization of the wave equation formulated at order 2. However, the first scheme has a lower complexity than the second one. Numerical results illustrate these equivalences and these efficiency comparisons.

**Key-words:** numerical analysis, implicit discretization, mixed formulation

RESEARCH CENTRE
BORDEAUX - SUD-OUEST

200 avenue de la Vieille Tour
33405 Talence Cedex

<footnote>* Inria Project-Team MAKUTU, E2S-UPPA, CNRS, Pau, France.</footnote>
# Deux stratégies de discrétisation implicites paramétriques de la formulation mixte des équations d'ondes linéaires

**Résumé :** Ce rapport présente deux stratégies implicites paramétriques de discrétisation de la formulation mixte des équations d'onde linéaire, basées pour la première sur le schéma de Störmer-Verlet et pour la seconde sur le schéma de Crank-Nicolson. Si les conditions initiales sont bien choisies, ces deux stratégies sont équivalentes et des formules d'équivalence permettent de reconstruire les variables d'un schéma à partir de l'autre. Ces deux schémas sont également équivalents à la discrétisation par $\theta$-schéma de l'équation des ondes formulée à l'ordre 2. Cependant, la complexité du premier est inférieure à celle du second. Des résultats numériques illustrent ces équivalences et ces comparaisons d'efficacité.

**Mots-clés :** analyse numérique, discrétisation implicite, formulation mixte
# Contents

<table><tr><td>1</td><td>General model, space discretization and usual time integration</td><td>4</td></tr><tr><td>1.1</td><td>Störmer-Verlet (SV)</td><td>4</td></tr><tr><td>1.2</td><td>Crank-Nicolson (CN)</td><td>5</td></tr><tr><td>2</td><td>Implicit parametrisation of the SV and CN schemes</td><td>5</td></tr><tr><td>2.1</td><td>Parametrized Störmer-Verlet (θ-SV)</td><td>5</td></tr><tr><td>2.2</td><td>Parametrized Crank-Nicolson (θ-CN)</td><td>6</td></tr><tr><td>2.3</td><td>Relation between θ-SV and θ-CN</td><td>6</td></tr><tr><td>3</td><td>Equivalence to second order θ-scheme</td><td>8</td></tr><tr><td>3.1</td><td>Field associated to the parametrized Störmer-Verlet (θ-SV)</td><td>8</td></tr><tr><td>3.2</td><td>Field associated to the parametrized Crank-Nicolson (θ-CN)</td><td>9</td></tr><tr><td>3.3</td><td>Equivalence of the fields</td><td>10</td></tr><tr><td>4</td><td>Numerical results</td><td>11</td></tr><tr><td>4.1</td><td>Test case</td><td>11</td></tr><tr><td>4.2</td><td>Verification of the reconstruction formulas</td><td>12</td></tr><tr><td>4.3</td><td>Computational efficiency</td><td>15</td></tr><tr><td>5</td><td>Conclusion and prospects</td><td>15</td></tr></table>
# 1 General model, space discretization and usual time integration

The mixed formulation of linear wave equations is one possible modeling of wave propagation phenomena. It is chosen over its second order formulation counterpart in situations where the unknowns are more relevant to the physical context [Makridakis and Monk, 1995, Lanteri and Scheid, 2013], or where these unknowns are more natural for modeling purposes (as for instance coupling with other parts [Banjai et al., 2015]), or even where it is not possible to formulate the equations as a second order equation (as for instance in presence of intricate dissipative [Bilbao and Harrison, 2016] or nonlinear phenomena [He and Sun, 2020]).

Numerical methods to solve this system are numerous and can rely on several analysis tools. In this work, we want to introduce and compare two parametrized time discretisation strategies, and we therefore suppose that the spatial discretisation is done with usual methods such as Finite Differences [Zuazua, 2005], Finite Elements [Brezzi and Fortin, 1991], Finite Volumes [Eymard et al., 2000], or any other method that provides the following semi discrete system, along with some necessary bounds on the space discretisation error (typically those supposed in [Chabassier, 2023]). More precisely, we assume that, after introducing a small space discretization parameter $h$, and following the steps that lead to the semi discrete system, the semi discrete solution $(p_h, v_h)$ is sought in finite dimensional spaces $U_h \subset U \subset P$ and $D_h \subset D$, where $P$ and $D$ are two Hilbert spaces naturally prompted by the equation as the variational spaces for both variables:

$$ \left\{ \begin{array}{l} p_h(0) = p_{h,0} \end{array} \right. \qquad (1.1a) $$

$$ \left\{ \begin{array}{l} v_h(0) = v_{h,0} \end{array} \right. \qquad (1.1b) $$

$$ \left\{ \begin{array}{l} \dot{p}_h + B_h^* v_h = f_h \end{array} \right. \qquad (1.1c) $$

$$ \left\{ \begin{array}{l} \dot{v}_h - B_h p_h = \dot{g}_h \end{array} \right. \qquad (1.1d) $$

where $B_h: U_h \to D_h$ is a discrete operator and $B_h^*: D_h \to U_h$ is its adjoint, and $p_{h,0}, v_{h,0}, f_h$ and $g_h$ are discrete representations of the initial conditions and sources in $U_h$ and $D_h$ respectively. We also denote $I_h$ the identity operator of $U_h$.

For the time discretization, we introduce the time step $\Delta t > 0$. Two usual discretization strategies are particularly spread and used in the community. The usual Störmer-Verlet (leap-frog) scheme [Hairer et al., 2003] is an explicit and conditionally stable strategy, and Crank-Nicolson scheme [Hairer et al., 2006] is an implicit and unconditionally stable strategy. They are recalled below.

## 1.1 Störmer-Verlet (SV)

Time is discretized on staggered time grids with the time step $\Delta t$. The unknown $p_h$ is sought for on a regular time grid $t^n = n\Delta t$ for $n \in [0, N]$ while the unknown $v_h$ is sought for on the interleaved regular time grid $t^{n+1/2} = (n+1/2)\Delta t$ for $n \in [0, N]$. Let us introduce the discrete operators $\delta$ and $\mu$ defined as

$$ \delta v_h^n = \frac{v_h^{n+1/2} - v_h^{n-1/2}}{\Delta t}, \quad \mu v_h^n = \frac{v_h^{n+1/2} + v_h^{n-1/2}}{2}, \quad \delta p_h^{n+1/2} = \frac{p_h^{n+1} - p_h^n}{\Delta t}, \quad \mu p_h^{n+1/2} = \frac{p_h^{n+1} + p_h^n}{2} \qquad (1.2) $$

The usual Störmer-Verlet scheme is defined as

$$ \left\{ \begin{array}{l} p_h^0 = p_h(0) \end{array} \right. \qquad (1.3a) $$

$$ \left\{ \begin{array}{l} v_h^{1/2} = v_h(0) + \frac{\Delta t}{2} [\delta g_h^{1/2} + B_h p_h(0)] \end{array} \right. \qquad (1.3b) $$

$$ \delta p_h^{n+1/2} + B_h^* v_h^{n+1/2} = f_h^{n+1/2}, \quad 0 \le n \le N \qquad (1.3c) $$

$$ \delta v_h^n - B_h p_h^n = \mu \delta g_h^n, \quad 1 \le n \le N \qquad (1.3d) $$

Note that the initial conditions are chosen in order to guarantee second order accuracy in time even in the presence of initial data and nonzero source terms at initial time.
## 1.2 Crank-Nicolson (CN)

We seek both unknowns $q_h$ and $w_h$ a regular time grid $t^n = n\Delta t$ for $n \in [0, N]$. Let us introduce the discrete operators $\delta$ and $\mu$ defined as

$$ \delta q_h^{n+\frac{1}{2}} = \frac{q_h^{n+1} - q_h^n}{\Delta t}, \quad \mu q_h^{n+\frac{1}{2}} = \frac{q_h^{n+1} + q_h^n}{2}, \quad \delta w_h^{n+\frac{1}{2}} = \frac{w_h^{n+1} - w_h^n}{\Delta t}, \quad \mu w_h^{n+\frac{1}{2}} = \frac{w_h^{n+1} + w_h^n}{2} \quad (1.4) $$

The usual Crank-Nicolson scheme reads

$$ q_h^0 = p_h(0) \qquad (1.5a) $$

$$ w_h^0 = v_h(0) \qquad (1.5b) $$

$$ \delta q_h^{n+\frac{1}{2}} + B_h^* \mu w_h^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}}, \quad 0 \le n \le N \qquad (1.5c) $$

$$ \delta w_h^{n+\frac{1}{2}} - B_h \mu q_h^{n+\frac{1}{2}} = \delta g_h^{n+\frac{1}{2}}, \quad 0 \le n \le N \qquad (1.5d) $$

# 2 Implicit parametrisation of the SV and CN schemes

Let $\theta \ge 0$ be a real number. We now consider two implicit parametrized perturbations of the Störmer-Verlet (leap-frog) and Crank-Nicolson schemes. The treatment of the source terms is tailored so that the comparison of the resulting discrete series is as straightforward as possible.

## 2.1 Parametrized Störmer-Verlet ($\theta$-SV)

The unknown $p_h$ is sought for on a regular time grid $t^n = n\Delta t$ for $n \in [0, N]$ while the unknown $v_h$ is sought for on the interleaved regular time grid $t^{n+\frac{1}{2}} = (n + \frac{1}{2})\Delta t$ for $n \in [0, N]$. The proposed scheme reads

$$ v_h^{\frac{1}{2}} = v_h(0) + \frac{\Delta t}{2} [\delta g_h^{\frac{1}{2}} + B_h p_h(0)] \qquad (2.1a) $$

$$ p_h^0 = p_h(0) \qquad (2.1b) $$

$$ (I_h + \theta \Delta t^2 B_h^* B_h) \delta p_h^{n+\frac{1}{2}} + B_h^* v_h^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}}, \quad 0 \le n \le N \qquad (2.1c) $$

$$ \delta v_h^n - B_h p_h^n = \mu \delta g_h^n, \quad 1 \le n \le N \qquad (2.1d) $$

We see that the usual SV scheme (1.3) is recovered from the choice $\theta = 0$. This is the only choice that makes the scheme explicit. This scheme satisfies an energy identity which prompts the stability condition for some values of $\theta$.

**PROPOSITION 2.1.** Any solution to (2.1) satisfies, for $n \in [0, N]$,

$$ \delta \mathcal{E}_h^{n+\frac{1}{2}} = \left( f_h^{n+\frac{1}{2}} + B_h^* \mu g_h^{n+\frac{1}{2}} \right) \cdot \mu p_h^{n+\frac{1}{2}} \quad \text{where} \quad \mathcal{E}_h^n = \frac{1}{2} \|p_h^n\|_{\tilde{I}_h}^2 + \frac{1}{2} \|\mu v_h^n - \mu^2 g_h^n\|^2 \quad (2.2) $$

where

$$ \tilde{I}_h = I_h + (\theta - \frac{1}{4}) \Delta t^2 B_h^* B_h \qquad (2.3) $$

*Proof.* Take the scalar product of (2.1c) with $\mu p_h^{n+\frac{1}{2}}$ and of $\mu$ (2.1d) with $\mu v_h^n - \mu^2 g_h^n$. $\square$

**PROPOSITION 2.2 (Stability condition).** If $\theta \ge \frac{1}{4}$, the discrete energy $\mathcal{E}_h^n$ is always positive. If $\theta < \frac{1}{4}$, the discrete energy is positive if

$$ \eta \le 1, \quad \text{where} \quad \eta = \frac{\Delta t}{2} \sqrt{1 - 4\theta} \sqrt{\rho(B_h^* B_h)} \qquad (2.4) $$

*Proof.* This directly follows from the definition of $\tilde{I}_h$. $\square$
## 2.2 Parametrized Crank-Nicolson ($\theta$-CN)

We seek both unknowns $q_h$ and $w_h$ a regular time grid $t^n = n\Delta t$ for $n \in [0, N]$. The proposed scheme reads

$$ w_h^0 = v_h(0) \qquad (2.5a) $$

$$ q_h^0 = p_h(0) \qquad (2.5b) $$

$$ \left( I_h + \left(\theta - \frac{1}{4}\right) \Delta t^2 B_h^* B_h \right) \delta q_h^{n+\frac{1}{2}} + B_h^* \mu w_h^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}}, \quad 0 \le n \le N \qquad (2.5c) $$

$$ \delta w_h^{n+\frac{1}{2}} - B_h \mu q_h^{n+\frac{1}{2}} = \delta g_h^{n+\frac{1}{2}}, \quad 0 \le n \le N \qquad (2.5d) $$

We see that the usual CV scheme (1.5) is recovered from the choice $\theta = \frac{1}{4}$. No choice of $\theta$ makes the scheme explicit. This scheme satisfies an energy identity which prompts the stability condition for some values of $\theta$.

**PROPOSITION 2.3.** Any solution to (2.5) satisfies, for $n \in [0, N]$,

$$ \delta \mathcal{E}_h^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}} \cdot \mu p_h^{n+\frac{1}{2}} + \delta g_h^{n+\frac{1}{2}} \cdot \mu w_h^{n+\frac{1}{2}} \quad \text{where} \quad \mathcal{E}_h^n = \frac{1}{2} \|q_h^n\|_{\tilde{I}_h}^2 + \frac{1}{2} \|w_h^n\|^2 \qquad (2.6) $$

where

$$ \tilde{I}_h = I_h + \left(\theta - \frac{1}{4}\right) \Delta t^2 B_h^* B_h \qquad (2.7) $$

*Proof.* Take the scalar product of (2.5c) with $\mu q_h^{n+\frac{1}{2}}$ and of (2.5d) with $\mu w_h^{n+\frac{1}{2}}$. $\square$

**PROPOSITION 2.4 (Stability condition).** If $\theta \ge \frac{1}{4}$, the discrete energy $\mathcal{E}_h^n$ is always positive. If $\theta < \frac{1}{4}$, the discrete energy is positive if

$$ \eta \le 1, \quad \text{where} \quad \eta = \frac{\Delta t}{2} \sqrt{1 - 4\theta} \sqrt{\rho(B_h^* B_h)} \qquad (2.8) $$

*Proof.* This directly follows from the definition of $\tilde{I}_h$. $\square$

## 2.3 Relation between $\theta$-SV and $\theta$-CN

The following proposition allows to reconstruct the fields of the $\theta$-SV scheme (2.1) from the fields obtained by the $\theta$-CN scheme (2.5).

**PROPOSITION 2.5.** Let $(w_h, q_h)$ be the solution to (2.5). Then, $(v_h, p_h)$ solution to (2.1) satisfy

$$ p_h^n = q_h^n, \quad n \ge 0 \qquad (2.9a) $$

$$ v_h^{n+\frac{1}{2}} = \mu w_h^{n+\frac{1}{2}} - \frac{\Delta t^2}{4} B_h \delta q_h^{n+\frac{1}{2}}, \quad n \ge 0 \qquad (2.9b) $$

*Proof.* Define the field $x_h^{n+\frac{1}{2}}$ for all $n \ge 0$ as

$$ x_h^{n+\frac{1}{2}} = \mu w_h^{n+\frac{1}{2}} - \frac{\Delta t^2}{4} B_h \delta q_h^{n+\frac{1}{2}} \qquad (2.10) $$

From the definition of $x_h$ we can derive the value of $\delta x_h^n$:

$$ \delta x_h^n = \delta \mu w_h^n - \frac{\Delta t^2}{4} B_h \delta^2 q_h^n \qquad (2.11) $$

$$ = B_h \left( \mu^2 - \frac{\Delta t^2}{4} \delta^2 \right) q_h^n + \delta \mu g_h^n \quad \text{using (2.5d)} \qquad (2.12) $$

$$ = B_h q_h^n + \delta \mu g_h^n \qquad (2.13) $$
If we use the definition of $x_h$ in (2.5c), we get the system

$$ \left\{ \begin{array}{l} q_h^0 = p_h(0) \\ x_h^{\frac{1}{2}} = \mu w_h^{\frac{1}{2}} - \frac{\Delta t^2}{4} B_h \delta q_h^{\frac{1}{2}} \end{array} \right. \qquad (2.14a) $$

$$ (I_h + \theta \Delta t^2 B_h^* B_h) \delta q_h^{n+\frac{1}{2}} + B_h^* x_h^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}}, \quad 0 \le n \le N \qquad (2.14c) $$

$$ \delta x_h^{n+\frac{1}{2}} - B_h q_h^{n+\frac{1}{2}} = \delta \mu g_h^{n+\frac{1}{2}}, \quad 0 \le n \le N \qquad (2.14d) $$

Finally, from (2.5d) taken for $n = 0$, we replace $\delta w_h^{\frac{1}{2}}$ and $\mu q_h^{\frac{1}{2}}$ with $\mu w_h^{\frac{1}{2}}$ and $\delta q_h^{\frac{1}{2}}$ and remaining terms to get

$$ \frac{2}{\Delta t} \mu w_h^{\frac{1}{2}} - \frac{2}{\Delta t} w_h^0 - B_h \left[ \frac{\Delta t}{2} \delta q_h^{\frac{1}{2}} + q_h^0 \right] = \delta g_h^{\frac{1}{2}} \qquad (2.15) $$

We recognize the terms forming the expression of $x_h^{\frac{1}{2}}$ to get that

$$ x_h^{\frac{1}{2}} = w_h^0 + \frac{\Delta t}{2} \left[ \delta g_h^{\frac{1}{2}} + B_h q_h^0 \right] \qquad (2.16) $$

In the end, the couple $(x_h, q_h)$ is solution to

$$ \left\{ \begin{array}{l} q_h^0 = p_h(0) \\ x_h^{\frac{1}{2}} = w_h^0 + \frac{\Delta t}{2} \left[ \delta g_h^{\frac{1}{2}} + B_h q_h^0 \right] \end{array} \right. \qquad (2.17a) $$

$$ (I_h + \theta \Delta t^2 B_h^* B_h) \delta q_h^{n+\frac{1}{2}} + B_h^* x_h^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}}, \quad 0 \le n \le N \qquad (2.17c) $$

$$ \delta x_h^{n+\frac{1}{2}} - B_h q_h^{n+\frac{1}{2}} = \delta \mu g_h^{n+\frac{1}{2}}, \quad 0 \le n \le N \qquad (2.17d) $$

which is strictly identical to (2.1). □

The following proposition gives the formula to reconstruct the fields of the $\theta$-CN scheme (2.5) from the fields of the $\theta$-SV scheme (2.1).

**PROPOSITION 2.6.** Let $(v_h, p_h)$ be the solution to (2.1). Then, $(w_h, q_h)$ solution to (2.5) satisfy

$$ \left\{ \begin{array}{l} q_h^n = p_h^n, \quad \forall n \ge 0 \\ w_h^n = \mu v_h^n - \frac{\Delta t^2}{4} \delta^2 g_h^n, \quad \forall n \ge 1 \end{array} \right. \qquad (2.18a) $$

$$ w_h^n = \mu v_h^n - \frac{\Delta t^2}{4} \delta^2 g_h^n, \quad \forall n \ge 1 \qquad (2.18b) $$

*Proof.* It is also possible to prove the proposition without using the field $u$. Indeed, let $(v_h, p_h)$ be the solution to the $\theta$-SV scheme (2.1), and define for all $n \ge 1$,

$$ y_h^n = \mu v_h^n - \frac{\Delta t^2}{4} \delta^2 g_h^n \qquad (2.19) $$

Let us compute

$$ \delta y_h^{n+\frac{1}{2}} - B_h \mu p_h^{n+\frac{1}{2}} = \delta \left( \mu v_h - \frac{\Delta t^2}{4} \delta^2 g_h \right)^{n+\frac{1}{2}} - B_h \mu p_h^{n+\frac{1}{2}} \qquad (2.20) $$

$$ = \mu (\delta v_h)^{n+\frac{1}{2}} - \frac{\Delta t^2}{4} \delta^3 g_h^n - B_h \mu p_h^{n+\frac{1}{2}} \qquad (2.21) $$

$$ =(2.1d) \mu (\mu \delta g_h + B_h p_h)^{n+\frac{1}{2}} - \frac{\Delta t^2}{4} \delta^3 g_h^n - B_h \mu p_h^{n+\frac{1}{2}} \qquad (2.22) $$

$$ = \delta \left( \mu^2 - \frac{\Delta t^2}{4} \delta^2 \right) g_h^{n+\frac{1}{2}} \qquad (2.23) $$

$$ = \delta g_h^{n+\frac{1}{2}} \qquad (2.24) $$
On the other hand, let $\tilde{I}_h = (I_h + (\theta - \frac{1}{4})\Delta t^2 B_h^* B_h)$. We have

$$ \tilde{I}_h \delta p_h^{n+\frac{1}{2}} + B_h^* \mu y_h^{n+\frac{1}{2}} = \tilde{I}_h \delta p_h^{n+\frac{1}{2}} + B_h^* \mu \left[ \mu v_h - \frac{\Delta t^2}{4} \delta^2 g_h \right]^{n+\frac{1}{2}} \quad (2.25) $$

$$ = \tilde{I}_h \delta p_h^n + B_h^* \mu^2 v_h^{n+\frac{1}{2}} - B_h^* \frac{\Delta t^2}{4} \mu \delta^2 g_h^{n+\frac{1}{2}} \quad (2.26) $$

$$ = (2.31) \quad \tilde{I}_h \delta p_h^n + B_h^* \left[ v_h^{n+\frac{1}{2}} + \frac{\Delta t^2}{4} \delta^2 v_h^{n+\frac{1}{2}} \right] - B_h^* \frac{\Delta t^2}{4} \mu \delta^2 g_h^{n+\frac{1}{2}} \quad (2.27) $$

$$ = \delta_{(2.1d)} \tilde{I}_h \delta p_h^n + B_h^* v_h^{n+\frac{1}{2}} + \frac{\Delta t^2}{4} B_h^* \left[ \mu \delta^2 g_h^{n+\frac{1}{2}} + B_h \delta p_h^{n+\frac{1}{2}} \right] - B_h^* \frac{\Delta t^2}{4} \mu \delta^2 g_h^{n+\frac{1}{2}} \quad (2.28) $$

$$ = (I_h + \theta \Delta t^2 B_h^* B_h) \delta p_h^n + B_h^* v_h^{n+\frac{1}{2}} \quad (2.29) $$

$$ = (2.1c) \quad f_h^{n+\frac{1}{2}} \quad (2.30) $$

where we used the identity

$$ X^{n+\frac{1}{2}} = \mu^2 X^{n+\frac{1}{2}} - \frac{\Delta t^2}{4} \delta^2 X^{n+\frac{1}{2}} \quad (2.31) $$

Therefore we recover that $(y_h, p_h)$ is solution to a system similar to (2.5). It remains to prove that the initial conditions match. It is the case for $p_h$ and it suffices to pose $y_h^0 = v_h(0)$ to ensure the full equivalence, hence $y_h^n = w_h^n$. $\square$

In other words, the two formulations are equivalent, and we have formulas to recover one from the other.

## 3 Equivalence to second order $\theta$-scheme

In this section, the solutions to (2.1) and to (2.5) are written as post processed quantities based on two different fields $u_1$ and $u_2$ which are in fact solutions to $\theta$-schemes with appropriate initial data. Then, these two fields are shown to be equal, which provides an alternate proof to Prop. 2.6.

### 3.1 Field associated to the parametrized Störmer-Verlet ($\theta$-SV)

**PROPOSITION 3.1.** The unknowns $(p_h^n, v_h^{n+\frac{1}{2}})$ solution to (2.1) can be written as:

$$ p_h^n = \delta u_1^n, \quad v_h^{n+\frac{1}{2}} = B_h u_1^{n+\frac{1}{2}} + \mu g_h^{n+\frac{1}{2}} \quad (3.1) $$

where $u_1^n$ is defined on the interleaved grid $t^{n+\frac{1}{2}}$ such that for all $0 \le n \le N$,

$$ u_1^{-\frac{1}{2}} = u_1^{\frac{1}{2}} - \Delta t p_h^0 \quad (3.2a) $$

$$ u_1^{\frac{1}{2}} = (B_h^* B_h)^{-1} B_h^* \left[ v_h^{\frac{1}{2}} - \mu g_h^{\frac{1}{2}} \right] \quad (3.2b) $$

$$ \delta^2 u_1^{n+\frac{1}{2}} + B_h^* B_h \{u_1\}_\theta^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}} - B_h^* \mu g_h^{n+\frac{1}{2}}, \quad 0 \le n \le N \quad (3.2c) $$

*Proof.* Let us introduce a field $u_1$ defined on the interleaved grid $t^{n+\frac{1}{2}}$ such that for all $0 \le n \le N$,

$$ \delta u_1^n = p_h^n, \quad u_1^{\frac{1}{2}} = (B_h^* B_h)^{-1} B_h^* \left[ v_h^{\frac{1}{2}} - \mu g_h^{\frac{1}{2}} \right] \quad (3.3) $$

Then

$$ (I_h + \theta \Delta t^2 B_h^* B_h) \delta^2 u_1^{n+\frac{1}{2}} + B_h^* v_h^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}}, \quad 1 \le n \le N \quad (3.4a) $$

$$ \delta v_h^n - B_h \delta u_1^n = \mu \delta g_h^n, \quad 1 \le n \le N \quad (3.4b) $$

We get from summation of $B_h^*$ applied to the second equation from 1 to $n$ that

$$ B_h^* v_h^{n+\frac{1}{2}} = B_h^* B_h u_1^{n+\frac{1}{2}} + B_h^* \mu g_h^{n+\frac{1}{2}} + \underbrace{B_h^* v_h^{\frac{1}{2}} - B_h^* B_h u_1^{\frac{1}{2}} - B_h^* \mu g_h^{\frac{1}{2}}}_{=0} \quad (3.5) $$
Hence

$$ (I_h + \theta \Delta t^2 B_h^* B_h) \delta^2 u_1^{n+\frac{1}{2}} + B_h^* B_h u_1^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}} - B_h^* \mu g_h^{n+\frac{1}{2}} \quad (3.6) $$

Which also can be written as the usual $\theta$-scheme on $u_1$

$$ u_1^{-\frac{1}{2}} = u_1^{\frac{1}{2}} - \Delta t p_h^0 \quad (3.7a) $$

$$ u_1^{\frac{1}{2}} = (B_h^* B_h)^{-1} B_h^* \left[ v_h^{\frac{1}{2}} - \mu g_h^{\frac{1}{2}} \right] \quad (3.7b) $$

$$ \delta^2 u_1^{n+\frac{1}{2}} + B_h^* B_h \{u_1\}_\theta^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}} - B_h^* \mu g_h^{n+\frac{1}{2}}, \quad 0 \le n \le N \quad (3.7c) $$

Moreover,

$$ v_h^{n+\frac{1}{2}} = B_h u_1^{n+\frac{1}{2}} + \mu g_h^{n+\frac{1}{2}} + e_v^{\frac{1}{2}}, \quad \text{where } e_v^{\frac{1}{2}} = v_h^{\frac{1}{2}} - B_h u_1^{\frac{1}{2}} - \mu g_h^{\frac{1}{2}} \quad (3.8) $$

But

$$ B_h^* B_h u_1^{\frac{1}{2}} = B_h^* \left[ v_h^{\frac{1}{2}} - \mu g_h^{\frac{1}{2}} \right] = B_h^* e_v^{\frac{1}{2}} + B_h^* B_h u_1^{\frac{1}{2}} \Rightarrow B_h^* e_v^{\frac{1}{2}} = 0 \Rightarrow B_h B_h^* e_v^{\frac{1}{2}} = 0 \Rightarrow e_v^{\frac{1}{2}} = 0 \quad (3.9) $$

□

## 3.2 Field associated to the parametrized Crank-Nicolson ($\theta$-CN)

**PROPOSITION 3.2.** The unknowns $w_h^n$ and $q_h^n$ solution to (2.5) can be written as:

$$ q_h^n = \delta u_2^n, \quad w_h^n = B_h \mu u_2^n + g_h^n \quad (3.10) $$

where the field $u_2$ is defined on the interleaved grid $t^{n+\frac{1}{2}}$ such that for all $1 \le n \le N$,

$$ u_2^{-\frac{1}{2}} = u_2^{\frac{1}{2}} - \Delta t q_h^0 \quad (3.11a) $$

$$ u_2^{\frac{1}{2}} = (B_h^* B_h)^{-1} [B_h^* w_h^0 - B_h^* g_h^0] + \frac{\Delta t}{2} q_h^0 \quad (3.11b) $$

$$ \delta^2 u_2^{n+\frac{1}{2}} + B_h^* B_h \{u_2\}_\theta^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}} - B_h^* \mu g_h^{n+\frac{1}{2}}, \quad 0 \le n \le N \quad (3.11c) $$

*Proof.* Let us introduce a field $u_2$ defined on the interleaved grid $t^{n+\frac{1}{2}}$ such that for all $1 \le n \le N$,

$$ \delta u_2^n = q_h^n, \quad u_2^{\frac{1}{2}} = (B_h^* B_h)^{-1} [B_h^* w_h^0 - B_h^* g_h^0] + \frac{\Delta t}{2} q_h^0 \quad (3.12) $$

$$ \left\{ \left( I_h + \left(\theta - \frac{1}{4}\right) \Delta t^2 B_h^* B_h \right) \delta^2 u_2^{n+\frac{1}{2}} + B_h^* \mu w_h^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}} \right. \quad (3.13a) $$

$$ \left. \delta w_h^{n+\frac{1}{2}} - B_h \mu \delta u_2^{n+\frac{1}{2}} = \delta g_h^{n+\frac{1}{2}} \right. \quad (3.13b) $$

We get from summation of $B_h^*$ applied to the second equation from 0 to n that

$$ B_h^* w_h^n = B_h^* B_h \mu u_2^n + B_h^* g_h^n + \underbrace{B_h^* w_h^0 - B_h^* B_h \mu u_2^0 - B_h^* g_h^0}_{=0} \quad (3.14) $$

Hence,

$$ \left( I_h + \left(\theta - \frac{1}{4}\right) \Delta t^2 B_h^* B_h \right) \delta^2 u_2^{n+\frac{1}{2}} + B_h^* B_h \mu^2 u_2^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}} - B_h^* \mu g_h^{n+\frac{1}{2}} \quad (3.15) $$

Which also can be written as the usual $\theta$-scheme on $u_2$

$$ u_2^{-\frac{1}{2}} = u_2^{\frac{1}{2}} - \Delta t q_h^0 \quad (3.16a) $$

$$ u_2^{\frac{1}{2}} = (B_h^* B_h)^{-1} [B_h^* w_h^0 - B_h^* g_h^0] + \frac{\Delta t}{2} q_h^0 \quad (3.16b) $$

$$ \delta^2 u_2^{n+\frac{1}{2}} + B_h^* B_h \{u_2\}_\theta^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}} - B_h^* \mu g_h^{n+\frac{1}{2}}, \quad 0 \le n \le N \quad (3.16c) $$
Moreover,

$$w_h^n = B_h \mu u_2^n + g_h^n + e_w^0, \quad \text{where } e_w^0 = w_h^0 - B_h \mu u_2^0 - g_h^0 \qquad (3.17)$$

But

$$B_h^* w_h^0 - B_h^* B_h \mu u_2^0 - B_h^* g_h^0 = 0 \Rightarrow B_h [B_h^* w_h^0 - B_h^* B_h \mu u_2^0 - B_h^* g_h^0] = 0 \qquad (3.18)$$

$$\Rightarrow w_h^0 - B_h \mu u_2^0 - g_h^0 = 0 = e_w^0 \quad \text{by inverting } B_h B_h^* \qquad (3.19)$$

which proves the relation between $u_2$ and $(q_h^n, w_h^n)$. $\square$

## 3.3 Equivalence of the fields

The similarity between (3.2) and (3.11) is striking. Only the initial conditions seem to differ, but the following proposition shows that they are in fact the same.

**PROPOSITION 3.3.** The fields $u_1$ defined by (3.2) and $u_2$ defined by (3.11) are equal.

*Proof.* The difference $\Delta_u = u_1 - u_2$ is solution to the following system of equations:

$$\Delta_u^{-\frac{1}{2}} = u_1^{-\frac{1}{2}} - u_2^{-\frac{1}{2}} \qquad (3.20a)$$

$$\Delta_u^{\frac{1}{2}} = u_1^{\frac{1}{2}} - u_2^{\frac{1}{2}} \qquad (3.20b)$$

$$\delta^2 \Delta_u^{n+\frac{1}{2}} + B_h^* B_h \{\Delta_u\}_\theta^{n+\frac{1}{2}} = 0, \quad 0 \le n \le N \qquad (3.20c)$$

where

$$u_1^{\frac{1}{2}} = (B_h^* B_h)^{-1} B_h^* [v_h^{\frac{1}{2}} - \mu g_h^{\frac{1}{2}}], \quad u_1^{-\frac{1}{2}} = u_1^{\frac{1}{2}} - \Delta t p_h^0 \qquad (3.21)$$

and

$$u_2^{\frac{1}{2}} = (B_h^* B_h)^{-1} B_h^* [w_h^0 - g_h^0] + \frac{\Delta t}{2} q_h^0, \quad u_2^{-\frac{1}{2}} = u_2^{\frac{1}{2}} - \Delta t q_h^0 \qquad (3.22)$$

From discrete energy considerations [Chabassier and Imperiale, 2017], and given that the source term is null in (3.20c), we get that

$$\|\Delta_u^{n+\frac{1}{2}}\| \le \|\Delta_u^{-\frac{1}{2}}\| + t^{n+\frac{1}{2}} \sqrt{2\epsilon_\Delta^0} \qquad (3.23)$$

where

$$\mathcal{E}_{\Delta}^{0} = \frac{1}{2} \|\delta \Delta_{u}^{0}\|^{2} + \frac{1}{2} \|\mu \Delta_{u}^{0}\|_{\tilde{I}_{h,\theta}}^{2}, \quad \text{with } \tilde{I}_{h,\theta} = I_{h} + \Delta t^{2}(\theta - \frac{1}{4}) B_{h}^{*} B_{h} \qquad (3.24)$$

Since $p_h^0 = q_h^0$ (see eqs. (2.1b) and (2.5b)), we get that

$$\Delta_u^{\frac{1}{2}} = u_1^{\frac{1}{2}} - u_2^{\frac{1}{2}} = u_1^{-\frac{1}{2}} - u_2^{-\frac{1}{2}} = \Delta_u^{-\frac{1}{2}} \qquad (3.25)$$

and

$$\Delta_u^{\frac{1}{2}} = (B_h^* B_h)^{-1} B_h^* [v_h^{\frac{1}{2}} - \mu g_h^{\frac{1}{2}} - w_h^0 + g_h^0] - \frac{\Delta t}{2} q_h^0 \qquad (3.26)$$

$$= (B_h^* B_h)^{-1} B_h^* \left[ v_h(0) + \frac{\Delta t}{2} \delta g_h^{\frac{1}{2}} + \frac{\Delta t}{2} B_h p_h(0) - \mu g_h^{\frac{1}{2}} - v_h(0) + g_h^0 - \frac{\Delta t}{2} B_h p_h(0) \right] \qquad (3.27)$$

$$= (B_h^* B_h)^{-1} B_h^* \left[ \frac{g_h^1 - g_h^0}{2} - \frac{g_h^1 + g_h^0}{2} + g_h^0 \right] \qquad (3.28)$$

$$= 0 \qquad (3.29)$$

Hence, for all $n \ge 0$,

$$\|\Delta_u^{n-\frac{1}{2}}\| \le 0 \Rightarrow \Delta_u^{n-\frac{1}{2}} = 0 \qquad \square \qquad (3.30)$$

From this proposition, we can now introduce a unique field $u = u_1 = u_2$. This field allows to derive a reconstruction of the fields of the $\theta$-CN scheme (2.5) from the fields obtained by the $\theta$-SV scheme (2.1). This prompts a second proof to Prop. 2.6:
Proof 2 of Prop 2.6, using the fields. The fields $p_h$, $q_h$, $v_h$, $w_h$ can be reconstructed from $u$ as:

$$p_h^n = \delta u^n, \quad q_h^n = \delta u^n, \quad \forall n \ge 0 \qquad (3.31)$$

$$v_h^{n+\frac{1}{2}} = B_h u^{n+\frac{1}{2}} + \mu g_h^{n+\frac{1}{2}}, \quad w_h^n = B_h \mu u^n + g_h^n, \quad \forall n \ge 0 \qquad (3.32)$$

Hence, $p_h^n = q_h^n$ and

$$\mu v_h^n = B_h \mu u^n + \mu^2 g_h^n = w_h^n - g_h^n + \mu^2 g_h^n = w_h^n + \frac{\Delta t^2}{4} \delta^2 g_h^n \qquad (3.33)$$

□

# 4 Numerical results

## 4.1 Test case

This section presents very simple numerical results that illustrate the reconstruction formulas given in this work. An impulse source is put at the entry of a 1D pipe in which waves propagate. More precisely, the pressure $p = p(x, t)$ and the velocity $v = v(x, t)$ are solution to

$$\left\{
\begin{aligned}
p(0) &= 0 && (4.1a) \\
v(0) &= 0 && (4.1b) \\
m_p(x)\dot{p} + \partial_x v &= 0 && (4.1c) \\
m_v(x)\dot{v} + \partial_x p &= 0 && (4.1d) \\
v(0) &= \lambda(t) && (4.1e) \\
v(L) &= 0 && (4.1f)
\end{aligned}
\right.$$

Following the usual mixed finite elements theory [Brezzi and Fortin, 1991], the variables are sought to in variational spaces of unequal regularity.

$$\mathcal{H} = \{(p, v) \in H^1([0, L]) \times L^2([0, L])\} \qquad (4.2)$$

The considered variational formulation is to seek $(p, v) \in \mathcal{H}$ such that for all $(\tilde{p}, \tilde{v}) \in \mathcal{H}$,

$$\int_{0}^{L} m_{p}(x) \dot{p} \tilde{p} - \int_{0}^{L} v \partial_{x} \tilde{p} - \lambda(t) \tilde{p}(0) = 0 \qquad (4.3a)$$

$$\int_{0}^{L} m_{v}(x) \dot{v} \tilde{v} + \int_{0}^{L} \partial_{x} p \tilde{v} = 0 \qquad (4.3b)$$

and the spatial discretization is done as described in [Tournemenne and Chabassier, 2019] by choosing conformal spectral finite elements of the same order on each element of the mesh, allowing jumps between elements for the variable $v$. This specific choice is not classical but is shown to converge in [Thibault, 2023, Chap 2]. This prompts the matricial formulation

$$p_h(0) = 0 \qquad (4.4a)$$

$$v_h(0) = 0 \qquad (4.4b)$$

$$\dot{p}_h + B_h^* v_h = -\lambda_0 E_{h,0} \qquad (4.4c)$$

$$\dot{v}_h - B_h p_h = 0 \qquad (4.4d)$$

This is a special case of (1.1) with a source term which has a contribution from the edge. In the following numerical example, we take $L = 1$, we use a $C^\infty$ source

$$\lambda(t) = \exp \left( - \frac{1}{1 - \left| \frac{2t - t_0}{t_0} \right|^2} \right), \quad 0 < t < t_0, \quad t_0 = 30 \qquad (4.5)$$
and the parameters are chosen as

$$m_p(x) = \pi R(x)^2, \quad m_v(x) = \frac{1}{\pi R(x)^2}, \quad R(x) = R_0 + x \frac{R_1 - R_0}{L}, \quad R_0 = 3 \times 10^{-3}, \quad R_1 = 8 \times 10^{-3} \quad (4.6)$$

We use 4th order finite elements of length 0.5 so that 2 elements compose the mesh. We set a time step $\Delta t$ chosen according to a chosen value of $\eta \le 1$ defined in (2.8) and (2.4) (which are the same definition) when $\theta < \frac{1}{4}$. When $\theta \ge \frac{1}{4}$, $\Delta t$ is chosen as 10 times the maximal allowed value of the explicit case. In other words,

$$\Delta t = \frac{\eta}{\sqrt{1 - 4\theta}} \frac{2}{\sqrt{\rho(B_h^* B_h)}} \quad \text{if } \theta < \frac{1}{4} \quad (4.7a)$$

$$\Delta t = 10 \frac{2}{\sqrt{\rho(B_h^* B_h)}} \quad \text{if } \theta \ge \frac{1}{4} \quad (4.7b)$$

## 4.2 Verification of the reconstruction formulas

We first run the $\theta$-SV scheme. Then, as a post processing step, we compute the quantities defined in (2.18) and compare them numerically to the ones obtained by running directly $\theta$-CN. The same process is done by running first the $\theta$-CN scheme, then computing the quantities defined in (2.9) and comparing them numerically to the ones obtained by running directly $\theta$-SV.

The relative error between both sides of equations (2.18) and (2.9) are computed and displayed in Fig. 1 for $\theta \in \{0, 0.25\}$, $\eta \in \{0.5, 1\}$. We observe that the relations are indeed respected up to machine precision.

<figure><img src="image_4.png" /><figcaption>Figure 1: Evolution of different ways of computing the velocity along time, $\theta \in \{0, 0.25\}$, $\eta \in \{0.5, 1\}$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 2: Evolution of different ways of computing the velocity along time, $\theta = 0, \eta \in \{0.5, 1\}$.</figcaption></figure>

The results obtained for $\eta \in \{0.5, 1\}$ and $\theta = 0, 1/12, 0.24, 0.25$ and $0.5$ are shown respectively in Fig. 2, 3, 4, 5 and 6. In those figures, the values are zoomed in for better visibility. We see that both relations (2.18) and (2.9) are respected.

<figure><img src="image_6.png" /><figcaption>Figure 3: Evolution of different ways of computing the velocity along time, $\theta = 1/12, \eta \in \{0.5, 1\}$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 4: Evolution of different ways of computing the velocity along time, $\theta = 0.24$, $\eta \in \{0.5, 1\}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5: Evolution of different ways of computing the velocity along time, $\theta = 0.25$, $\eta \in \{0.5, 1\}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6: Evolution of different ways of computing the velocity along time, $\theta = 0.5$, $\eta \in \{0.5, 1\}$.</figcaption></figure>
## 4.3 Computational efficiency

<figure><img src="image_10.png" /><figcaption>Figure 7: Comparison of computation times</figcaption></figure>

Finally, we compare the computation times of $\theta$-SV and $\theta$-CN for the same value of $\theta$. Since we have seen that both schemes are in fact equivalent (up to a post-processing step), the less computationally intensive should be preferred.

**Solving (2.1)** requires at each time step one matrix-vector product with the inverse of the matrix

$$I_h + \theta \Delta t^2 B_h^* B_h \qquad (4.8)$$

which is a sparse matrix of size $n_p \times n_p$. Moreover, at each time step, the scheme also requires one matrix-vector product by $B_h^*$ and one matrix-vector product by $B_h$.

**Solving (2.5)** requires at each time step the inversion of the matrix

$$\left( I_h + \left( \theta - \frac{1}{4} \right) \Delta t^2 B_h^* B_h - \frac{\Delta t}{2} B_h \quad \frac{\Delta t}{2} B_h^* \right) \qquad (4.9)$$

which is a sparse matrix of size $(n_p + n_w) \times (n_p + n_w)$. Moreover, at each time step, the scheme also requires one matrix-vector product by $B_h^*$ and one matrix-vector product by $B_h$.

In both cases, the matrix to invert can be factorized once before the time iterations. The $\theta$-SV scheme should always require less computation.

In Fig. 7, the computation time for the two strategies are compared for the same choice of $\theta$ and $\eta$. As $\theta$ increases, the time step of the computation increases, except for $\theta = \frac{1}{4}$ and $\theta = \frac{1}{2}$ where $\Delta t$ is the same, as described in (4.7). Every simulation is run 5 times in order to average the uncertainties. We observe that the computation time of the $\theta$-SV is always smaller (about 1/3 in this specific 1D case) than the one of $\theta$-CN, which is coherent with the complexity description.

**Remark 4.1.** In some sense, the $\theta$-SV scheme can be seen as an efficient algorithm to solve the $\theta$-CN scheme. In particular, the $\frac{1}{4}$-SV scheme is a cheap way of solving the CN scheme.

# 5 Conclusion and prospects

The two parametrized implicit discretization strategies presented in this work have been shown to be equivalent, if the initial conditions are well chosen. This equivalence has been illustrated by numerical results, as well as the
induced computational cost. A post processing step allows to reconstruct the same field from both schemes, that is also solution to a classical $\theta$-scheme. This work focused on numerical and computational aspects of these schemes, but this equivalence prompts the question of their analysis of stability and convergence. Indeed, one might now expect similar results as those obtained for the $\theta$-scheme applied to the second order equation (see [Chabassier and Imperiale, 2021]). However, the derivation of optimal proofs of stability and space/time convergence for the Störmer-Verlet scheme applied to the mixed formulation is not a straightforward question. Existing stability proofs for the usual Störmer-Verlet scheme (1.3), which are reviewed in [Chabassier, 2023], do not apply when the time step reaches its greatest admissible value for stability. In [Chabassier, 2023] the author proposed alternative proofs of stability and convergence that also apply to this case and that use regularity assumptions similar to those of [Chabassier and Imperiale, 2021]. It is an open question to generalize those results to the implicit parametrized schemes proposed here. The equivalence of the two formulations shown in this work could help harmonize the required hypotheses.
## References

[Banjai et al., 2015] Banjai, L., Lubich, C., and Sayas, F.-J. (2015). Stable numerical coupling of exterior and interior problems for the wave equation. *Numerische Mathematik*, 129(4):611–646.

[Bilbao and Harrison, 2016] Bilbao, S. and Harrison, R. (2016). Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section. *The Journal of the Acoustical Society of America*, 140(1):728–740.

[Brezzi and Fortin, 1991] Brezzi, F. and Fortin, M. (1991). *Mixed and Hybrid Finite Element Methods*, volume 2. New York: Springer-Verlag.

[Chabassier, 2023] Chabassier, J. (2023). Stability and space/time convergence of Störmer-Verlet time integration of the mixed formulation of linear wave equations.

[Chabassier and Imperiale, 2017] Chabassier, J. and Imperiale, S. (2017). Space/time convergence analysis of a class of conservative schemes for linear wave equations. *Comptes Rendus Mathematiques*, 355(3):282-289.

[Chabassier and Imperiale, 2021] Chabassier, J. and Imperiale, S. (2021). Construction and convergence analysis of conservative second order local time discretisation for linear wave equations. *ESAIM: Mathematical Modelling and Numerical Analysis*, 55(4):1507–1543.

[Eymard et al., 2000] Eymard, R., Gallouët, T., and Herbin, R. (2000). Finite volume methods. In *Handbook of Numerical Analysis*, volume 7 of *Solution of Equation in $\mathbb{R}$ (Part 3)*, *Techniques of Scientific Computing (Part 3)*, pages 713–1018. Elsevier.

[Hairer et al., 2003] Hairer, E., Lubich, C., and Wanner, G. (2003). Geometric numerical integration illustrated by the Störmer-Verlet method. In Iserles, A., editor, *Acta Numerica 2003*, pages 399–450. Cambridge University Press, 1 edition.

[Hairer et al., 2006] Hairer, E., Lubich, C., and Wanner, G. (2006). *Geometric Numerical Integration*, volume 31 of *Springer Series in Computational Mathematics*. Springer-Verlag, Berlin/Heidelberg.

[He and Sun, 2020] He, M. and Sun, P. (2020). Energy-preserving finite element methods for a class of nonlinear wave equations. *Applied Numerical Mathematics*, 157:446–469.

[Lanteri and Scheid, 2013] Lanteri, S. and Scheid, C. (2013). Convergence of a discontinuous Galerkin scheme for the mixed time-domain Maxwell's equations in dispersive media. *IMA Journal of Numerical Analysis*, 33(2):432–459.

[Makridakis and Monk, 1995] Makridakis, Ch. G. and Monk, P. (1995). Time-discrete finite element schemes for Maxwell's equations. *ESAIM: Mathematical Modelling and Numerical Analysis*, 29(2):171–197.

[Thibault, 2023] Thibault, A. (2023). *Modélisation, Analyse et Simulation de l'acoustique Dissipative Dans Les Tubes Poreux Ou Rugueux. Application Aux Instruments à Vent*. PhD thesis, Université de Pau et des pays de l'Adour.

[Tournemenne and Chabassier, 2019] Tournemenne, R. and Chabassier, J. (2019). A comparison of a one-dimensional finite element method and the transfer matrix method for the computation of wind music instrument impedance. *Acta Acustica united with Acustica*, 105(5):838–849.

[Zuazua, 2005] Zuazua, E. (2005). Propagation, Observation, and Control of Waves Approximated by Finite Difference Methods. *SIAM Review*, 47(2):197–243.
Inria

RESEARCH CENTRE
BORDEAUX - SUD-OUEST

200 avenue de la Vieille Tour
33405 Talence Cedex

Publisher
Inria
Domaine de Voluceau - Rocquencourt
BP 105 - 78153 Le Chesnay Cedex
inria.fr

ISSN 0249-6399