# Stability and space/time convergence of Störmer-Verlet time integration of the mixed formulation of linear wave equations

Juliette Chabassier

► To cite this version:

Juliette Chabassier. Stability and space/time convergence of Störmer-Verlet time integration of the mixed formulation of linear wave equations. 2023. hal-04189789

HAL Id: hal-04189789

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

Preprint submitted on 29 Aug 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
# STABILITY AND SPACE/TIME CONVERGENCE OF STÖRMER-VERLET TIME INTEGRATION OF THE MIXED FORMULATION OF LINEAR WAVE EQUATIONS

JULIETTE CHABASSIER<fnref n="1" />

**Abstract.** This work focuses on the mixed formulation of linear wave equations. It provides a proof of stability and convergence of time discretisation of a semi discrete linear wave equation in mixed form with Störmer-Verlet time integration, that is uniform as the time step reaches its largest allowed value for stability (Courant-Friedrich-Levy condition), contrary to the proofs recalled here from the literature.

**Résumé.** Ce travail se concentre sur la formulation mixte des équations d'ondes linéaires. Il fournit une preuve de la stabilité et de la convergence de la discrétisation temporelle d'une équation d'onde linéaire semi-discrète sous forme mixte par l'intégration temporelle de Störmer-Verlet, qui est uniforme lorsque le pas de temps atteint sa plus grande valeur autorisée pour la stabilité (condition de Courant-Friedrich-Levy), contrairement aux preuves ici rappelées issues de la littérature.

**2020 Mathematics Subject Classification.** 35L05, 35M30, 65M12.

August 2023.

## INTRODUCTION

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, or where these unknowns are more natural for modeling purposes (as for instance coupling with other parts), or even where it is not possible to formulate the equations as a second order equation (as for instance in presence of intricate dissipative or nonlinear phenomena). This encompasses acoustic waves, elastodynamics, electromagnetic waves, etc. On an abstract level, this system reads, for $0 \le t \le T$,

$$ egin{cases} p(0) = p_0 & (0.1a) \\ v(0) = v_0 & (0.1b) \\ \dot{p} + B^*v = f & (0.1c) \\ \dot{v} - Bp = \dot{g} & (0.1d) \end{cases} $$

<footnote>Keywords and phrases: linear wave equation, mixed formulation, stability, energy, space/time convergence</footnote>

<footnote n="1">¹ Inria Bordeaux University, MAKUTU team, 200 avenue de la vieille Tour, 33400 TALENCE, France</footnote>
where $\mathcal{B}: U \subset P \to D$ is an operator and $\mathcal{B}^*$ its adjoint. The source terms and the initial conditions are supposed regular enough so that the following hypothesis holds, in three Hilbert spaces $P, U = \{p \in P | \mathcal{B}p \in D\} \subset P$ and $D$:

**Hypothesis 0.1** (Stability of the continuous system). The source terms and the initial conditions are regular enough such that there exists a constant $C > 0$ such that

$$ \|p\|_{C^3(0,T;P)} + \|p\|_{C^0(0,T;U)} + \|v\|_{C^2(0,T;D)} \leq C \quad (0.2) $$

Numerical methods to solve this system are numerous and can rely on several analysis tools. We want to focus on this work on the time discretisation with a Störmer-Verlet scheme, 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. More precisely, we assume that, after 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$: Find $p_h \in U_h$ and $v_h \in D_h$ such that

$$ \begin{cases} p_h(0) = p_{h,0} & (0.3a) \\ v_h(0) = v_{h,0} & (0.3b) \\ \dot{p}_h + B_h^* v_h = f_h & (0.3c) \\ \dot{v}_h - B_h p_h = \dot{g}_h & (0.3d) \end{cases} $$

where $B_h: U_h \to D_h$ is a discrete approximation of the operator $\mathcal{B}$, 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 $p_0, v_0, f$ and $g$ in $U_h$ and $D_h$. We also denote $I_h$ the identity operator of $U_h$.

For the sequel, we will suppose that the spatial discretisation satisfies the

**Hypothesis 0.2** (Stability of the semi discrete system). The spatial discretisation is such that there exists a constant $C > 0$ independent of $h$ such that

$$ \|p_h\|_{C^3(0,T;P)} + \|p_h\|_{C^0(0,T;U)} + \|v_h\|_{C^2(0,T;D)} \leq C \quad (0.4) $$

**Hypothesis 0.3** (Convergence of the semi discrete system). The spatial discretisation is such that there exists a function $\delta: \mathbb{R}_*^+ \to \mathbb{R}^+$, such that

$$ \|p - p_h\|_{C^0(0,T;P)} + \|v - v_h\|_{C^0(0,T;D)} \leq \delta(h), \quad \text{with } \delta(h) \xrightarrow{h \to 0} 0 \quad (0.5) $$

Time is discretised using a Störmer-Verlet algorithm on staggered time grids with a time step $\Delta t$. The unknown $p_h$ is sought for on the grid $\{t^{n+\frac{1}{2}} = (n + \frac{1}{2}) \Delta t\}_{1 \leq n \leq N}$ while the unknown $v_h$ is sought for on the grid $\{t^n = n \Delta t\}_{0 \leq n \leq N}$. Let us introduce the discrete operators $\delta$ and $\mu$ defined as

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

They satisfy the useful following properties

$$ v_h^{n-\frac{1}{2} \pm \frac{1}{2}} = \mu v_h^{n-\frac{1}{2}} \pm \frac{\Delta t}{2} \delta v_h^{n-\frac{1}{2}} \quad (0.7a) $$

$$ \|\delta X\| \leq \frac{2}{\Delta t} \mu \|X\| \quad (0.7b) $$
The Störmer Verlet scheme reads: Find {$p_h^{n-\frac{1}{2}}$}$_{1\le n\le N+1}$ $\in U_h^{N+1}$ and {$v_h^n$}$_{0\le n\le N+1}$ $\in D_h^{N+2}$ such that

$$v_h^0 = v_h(0) \tag{0.8a}$$

$$p_h^{\frac{1}{2}} = p_h(0) + \frac{\Delta t}{2} [f_h^0 - B_h^* v_h(0)] \tag{0.8b}$$

$$\text{For } n \in [1, N], \delta p_h^n + B_h^* v_h^n = f_h^n \tag{0.8c}$$

$$\text{For } n \in [0, N], \delta v_h^{n+\frac{1}{2}} - B_h p_h^{n+\frac{1}{2}} = \delta g_h^{n+\frac{1}{2}} \tag{0.8d}$$

where $f_h^n = f_h(t^n)$ and $g_h^n = g_h(t^n)$.

The motivations of the present paper are to provide proof of space/time convergence of the Störmer Verlet scheme on the mixed formulation of wave equations. In simple cases such as the one presented here, eliminating one unknown by applying the $\partial_t$ (resp. $\delta$) operator to (0.3c) (resp. (0.8c)) and using (0.3d) (resp. (0.8d)) to replace $v_h$ yields a second order linear wave equation, if the source terms are regular enough:

$$\ddot{p}_h + B_h^* B_h p_h = \dot{f}_h - B_h^* \dot{g}_h$$

and is equivalent to applying the leap frog scheme [Joly, 2003, Cohen and Fauqueux, 2011]:

$$\delta^2 p_h^{n+\frac{1}{2}} + B_h^* B_h p_h^{n+\frac{1}{2}} = \delta f_h^{n+\frac{1}{2}} - B_h^* \delta g_h^{n+\frac{1}{2}}$$

So, for a field $u_h$ such that $u_h(0) = u_{h,0}$, $\dot{u}_h = p_h$ and its discrete counterpart {$u_h^n$}$_{0\le n\le N}$ such that $u_h^0 = u_{h,0}$, $\delta u_h^{n+\frac{1}{2}} = p_h^{n+\frac{1}{2}}$, the Störmer Verlet scheme on the mixed formulation amounts to a leap frog scheme on the second order wave equation

$$\ddot{u}_h + B_h^* B_h u_h = f_h - B_h^* g_h$$

such that

$$\delta^2 u_h^{n+\frac{1}{2}} + B_h^* B_h u_h^{n+\frac{1}{2}} = f_h^{n+\frac{1}{2}} - B_h^* g_h^{n+\frac{1}{2}}$$

Hence we expect similar estimates, at least on this remaining field, as those found in [Chabassier and Imperiale, 2021, Chabassier and Imperiale, 2017] for the leap frog scheme applied to the second order wave equation: If the CFL stability condition is respected, e.g.

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

where $\rho$ stands for the spectral radius:

$$\rho(B_h^* B_h) = \sup_{p_h \in U_h} \frac{(B_h^* B_h p_h, p_h)_P}{(p_h, p_h)_P} = \sup_{p \in U_h} \frac{(B_h p_h, B_h p_h)_D}{(p_h, p_h)_P} = \|\| B_h \|\| \tag{0.9}$$

then, there exists a constant $C < 0$ independant of $h$ and $\eta$ such that

$$\| \bar{u}_h^n - u_h^n \| \le C \Delta t^2 \| u_h \|_{C^4(0,T;P)} \tag{0.10}$$

Existing proofs of space/time convergence of the Störmer Verlet scheme on the mixed formulation of wave equations are found in [Jenkins et al., 2002], [Fezoui et al., 2005], [Lanteri and Scheid, 2013], [Banjai et al., 2015], [Egger and Radu, 2020]. Some of them are quite generic, while others are performed for a very specific spatial discretisation, but in all those references, the authors either suppose that $\eta \le \frac{1}{2}$ or $\eta < 1$. Indeed, those proofs turn out to not be uniform as $\Delta t$ approaches its greatest admissible value, in other words the stability and convergence bounding constants depend on $\eta$ and blow up as $\eta \to 1$. Numerical evidence and equivalence
to second order formulation suggest that this limitation is only technical, which should be possible to overcome theoretically.

Moreover, we aim at providing a uniform proof that does not eliminate one of the two unknowns, because this procedure could prevent to generalize the result to the addition of terms as coupling [Banjai et al., 2015] or dissipation [Bilbao and Harrison, 2016].

The paper is organized as follows. In Section 1, several identified proofs of stability are put in the same framework in order to exhibit their blow up as the time step approaches its greatest allowed value. Then, uniform proofs of stability (Propositions 2.10 and 2.11 of Section 2) and convergence (Theorem 3.1 of Section 3) are proposed.

# 1. STABILITY FAR FROM THE CFL

In this section, several proofs of the literature are put under the same form on a very simple discrete system, in order to show their underlying mechanisms and exhibit the fact that the bounding stability constants depend on $\eta$, and even blow up as $\eta \to 1$.

Let us manipulate the discrete equations in the following way. Let $(\alpha, \beta, \gamma) \in (\mathbb{R}^+)^3$, we take the scalar product of $(0.8c)^n$ with $\mu p_h^n$ in $U_h$, of $\mu(0.8d)^n$ with $\beta v_h^n$ in $D_h$, of $(0.8d)^{n-\frac{1}{2}}$ with $\alpha\mu v_h^{n-\frac{1}{2}}$ in $D_h$, and of $(0.8d)^{n+\frac{1}{2}}$ with $(\alpha + \gamma)\mu v_h^{n+\frac{1}{2}}$ in $D_h$. We add up and get

$$ \begin{aligned} & \frac{\alpha}{2\Delta t} (\|v_h^n\|_D^2 - \|v_h^{n-1}\|_D^2) + \frac{\alpha + \gamma}{2\Delta t} (\|v_h^{n+1}\|_D^2 - \|v_h^n\|_D^2) + \frac{1}{2\Delta t} (\|p_h^{n+\frac{1}{2}}\|_P^2 - \|p_h^{n-\frac{1}{2}}\|_P^2) + \frac{\beta}{2\Delta t} ((v_h^{n+1}, v_h^n)_D - (v_h^n, v_h^{n-1})_D) \\ & = \alpha(B_h p_h^{n-\frac{1}{2}}, \mu v_h^{n-\frac{1}{2}})_D + (\alpha + \gamma)(B_h p_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D + \beta(B_h \mu p_h^n, v_h^n)_D - (v_h^n, B_h \mu p_h^n)_D \\ & \quad + \alpha(\delta g_h^{n-\frac{1}{2}}, \mu v_h^{n-\frac{1}{2}})_D + (\alpha + \gamma)(\delta g_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D + \beta(\mu \delta g_h^n, v_h^n)_D + (f_h^n, \mu p_h^n)_P \end{aligned} \quad (1.1) $$

This is a discrete energy variation identity

$$ \begin{aligned} \delta \mathcal{E}_{\alpha, \beta, \gamma}^n ={}& \alpha \left[ (B_h p_h^{n-\frac{1}{2}}, \mu v_h^{n-\frac{1}{2}})_D + (B_h p_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D \right] + \gamma (B_h p_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D + (\beta - 1) (B_h \mu p_h^n, v_h^n)_D \\ & + \alpha \left[ (\delta g_h^{n-\frac{1}{2}}, \mu v_h^{n-\frac{1}{2}})_D + (\delta g_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D \right] + \gamma (\delta g_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D + \beta (\mu \delta g_h^n, v_h^n)_D + (f_h^n, \mu p_h^n)_P \end{aligned} \quad (1.2) $$

with

$$ \mathcal{E}_{\alpha, \beta, \gamma}^{n-\frac{1}{2}} = \frac{1}{2} \|p_h^{n-\frac{1}{2}}\|_P^2 + \alpha \left[ \frac{1}{2} \|v_h^{n-1}\|_D^2 + \frac{1}{2} \|v_h^n\|_D^2 \right] + \frac{\gamma}{2} \|v_h^n\|_D^2 + \frac{\beta}{2} (v_h^n, v_h^{n-1})_D \quad (1.3) $$

## 1.1. Choice $\beta = 1$ and $\alpha = \gamma = 0$

The choice $\beta = 1$ and $\alpha = \gamma = 0$ is done in [Joly, 2003] and [Egger and Radu, 2020, proof of Lemma 5.3]. It leads to a crossing term in the energy

$$ \mathcal{E}_{0,1,0}^{n-\frac{1}{2}} = \frac{1}{2} \|p_h^{n-\frac{1}{2}}\|_P^2 + \frac{1}{2}(v_h^n, v_h^{n-1})_D \equiv \frac{1}{2} \|p_h^{n-\frac{1}{2}}\|_P^2 + \frac{1}{2} \|\mu v_h^{n-\frac{1}{2}}\|_D^2 - \frac{\Delta t^2}{8} \|\delta v_h^{n-\frac{1}{2}}\|_D^2 \quad (1.4) $$

by using (0.7a).

### 1.1.1. Special case.

If $\delta g_h \equiv 0$, we have that $\delta v_h^{n-\frac{1}{2}} = B_h p_h^{n-\frac{1}{2}}$, therefore

$$ \mathcal{E}_{0,1,0}^{n-\frac{1}{2}} = \frac{1}{2} (\tilde{I}_h p_h^{n-\frac{1}{2}}, p_h^{n-\frac{1}{2}})_P + \frac{1}{2} \|\mu v_h^{n-\frac{1}{2}}\|_D^2 \quad (1.5) $$
with $\tilde{I}_h = I_h - \frac{\Delta t^2}{4} B_h^* B_h$ a modified identity operator. It is positive under the usual (CFL-condition). Then if $\eta \le 1$,

$$ \mathcal{E}_{0,1,0}^{n-\frac{1}{2}} \ge \frac{1-\eta^2}{2} \|p_h^{n-\frac{1}{2}}\|_P^2 \quad \text{and} \quad \mathcal{E}_{0,1,0}^{n-\frac{1}{2}} \ge \frac{1}{2} \|\mu v_h^{n-\frac{1}{2}}\|_D^2 \quad (1.6) $$

If moreover $\eta < 1$ (strict CFL condition), we have

$$ \delta \mathcal{E}_{0,1,0}^n = (f_h^n, \mu p_h^n)_P \le \|f_h^n\|_P \mu \|p_h\|_P^n \le \|f_h^n\|_P \frac{\sqrt{2}}{\sqrt{1-\eta^2}} \mu(\sqrt{\mathcal{E}_{0,1,0}})^n \quad (1.7) $$

A discrete summation from $n=1$ to $n=N$ yields, after telescopic elimination,

$$ \sqrt{\mathcal{E}_{0,1,0}^{N+\frac{1}{2}}} \le \sqrt{\mathcal{E}_{0,1,0}^{\frac{1}{2}}} + \frac{\sqrt{2}}{\sqrt{1-\eta^2}} \Delta t \sum_{n=1}^{N} \|f_h^n\|_P \quad (1.8) $$

which grants the stability of the energy, under the strict CFL condition. It then provides the stability of $\|p_h^{n-\frac{1}{2}}\|$ and $\|\mu v_h^{n-\frac{1}{2}}\|$ by using (1.6). The stability constants blow up as $\eta \to 1$.

*Remark 1.1.* By introducing auxiliary functions on top of more usual elliptic projections, [Egger and Radu, 2020] write a system of equations on the error terms which are of the form of (0.8) with $\delta g_h = 0$ (see Eq. (5.5)-(5.6)). They suppose $\eta \le \frac{1}{2}$ in the second part of their hypothesis (A4).

### 1.1.2. General case.

If $\delta g_h \ne 0$, a Young's inequality with any $\varepsilon > 0$ leads to

$$ \begin{aligned} \|\delta v_h^{n-\frac{1}{2}}\|_D^2 &= \|B_h p_h^{n-\frac{1}{2}} + \delta g_h^{n-\frac{1}{2}}\|_D^2 = \|B_h p_h^{n-\frac{1}{2}}\|_D^2 + 2(B_h p_h^{n-\frac{1}{2}}, \delta g_h^{n-\frac{1}{2}})_D + \|\delta g_h^{n-\frac{1}{2}}\|_D^2 \\ &\le \|B_h p_h^{n-\frac{1}{2}}\|_D^2 + \|\delta g_h^{n-\frac{1}{2}}\|_D^2 + \varepsilon \|B_h p_h^{n-\frac{1}{2}}\|_D^2 + \frac{1}{\varepsilon} \|\delta g_h^{n-\frac{1}{2}}\|_D^2 \\ &\le (1+\varepsilon) \|B_h p_h^{n-\frac{1}{2}}\|_D^2 + (1+\frac{1}{\varepsilon}) \|\delta g_h^{n-\frac{1}{2}}\|_D^2 \end{aligned} \quad (1.9) $$

Hence

$$ \begin{aligned} \mathcal{E}_{0,1,0}^{n-\frac{1}{2}} &\ge \frac{1}{2} \|p_h^{n-\frac{1}{2}}\|_P^2 + \frac{1}{2} \|\mu v_h^{n-\frac{1}{2}}\|_D^2 - \frac{\Delta t^2}{8} \left[ (1+\varepsilon) \|B_h p_h^{n-\frac{1}{2}}\|_D^2 + (1+\frac{1}{\varepsilon}) \|\delta g_h^{n-\frac{1}{2}}\|_D^2 \right] \\ &\ge \frac{1}{2} (\tilde{I}_{h,\varepsilon} p_h^{n-\frac{1}{2}}, p_h^{n-\frac{1}{2}})_P + \frac{1}{2} \|\mu v_h^{n-\frac{1}{2}}\|_D^2 - \frac{\Delta t^2}{8} (1+\frac{1}{\varepsilon}) \|\delta g_h^{n-\frac{1}{2}}\|_D^2 \end{aligned} \quad (1.10) $$

with $\tilde{I}_{h,\varepsilon} = I_h - \frac{\Delta t^2}{4}(1+\varepsilon)B_h^* B_h$ a modified identity operator. It is positive under the $\varepsilon$-CFL condition

$$ \Delta t^2 \rho(B_h^* B_h) \le \frac{4}{1+\varepsilon} < 4 \quad (1.11) $$

Recall the definition of $\eta$ from (CFL-condition). The $\varepsilon$-CFL condition is equivalent to supposing that $\eta < 1$ (strict CFL condition), since it suffices to pose $\varepsilon = \frac{1-\eta^2}{1+\eta^2} > 0$ to satisfy (1.11). Therefore

$$ \mathcal{E}_{0,1,0}^{n-\frac{1}{2}} + \frac{\Delta t^2}{4} \frac{1}{1-\eta^2} \|\delta g_h^{n-\frac{1}{2}}\|_D^2 \ge \frac{1-\eta^2}{1+\eta^2} \|p_h^{n-\frac{1}{2}}\|_P^2 + \frac{1}{2} \|\mu v_h^{n-\frac{1}{2}}\|_D^2 \ge 0 \quad (1.12) $$

We then have

$$ \delta \mathcal{E}_{0,1,0}^n = (\mu \delta g_h^n, v_h^n)_D + (f_h^n, \mu p_h^n)_P \quad (1.13) $$
The first RHS term can not be bounded by the square root of the energy since $||v_h^n||_D \le ||\mu v_h^{n-\frac{1}{2}}||_D$. We cannot conclude.

1.2. **Choice** $\alpha = \frac{1}{2}$, $\beta = \gamma = 0$

The choice $\alpha = \frac{1}{2}$, $\beta = 0$ and $\gamma = 0$ is done in [Banjai et al., 2015]. This leads to

$$ \mathcal{E}_{\frac{1}{2},0,0}^{n-\frac{1}{2}} = \frac{1}{2} \|p_h^{n-\frac{1}{2}}\|_P^2 + \frac{1}{4} \|v_h^{n-1}\|_D^2 + \frac{1}{4} \|v_h^n\|_D^2 \quad (1.14) $$

which is positive by definition, and satisfies

$$ \begin{aligned} \delta \mathcal{E}_{\frac{1}{2},0,0}^{n} &= \frac{1}{2} (B_h p_h^{n-\frac{1}{2}}, \mu v_h^{n-\frac{1}{2}})_D + \frac{1}{2} (B_h p_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D - (B_h \mu p_h^n, v_h^n)_D \\ & \qquad + \frac{1}{2} \left[ (\delta g_h^{n-\frac{1}{2}}, \mu v_h^{n-\frac{1}{2}})_D + (\delta g_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D \right] + (f_h^n, \mu p_h^n)_P \end{aligned} \quad (1.15) $$

$$ = \frac{1}{2} (B_h p_h^{n-\frac{1}{2}}, \frac{v_h^n + v_h^{n-1}}{2})_D + \frac{1}{2} (B_h p_h^{n+\frac{1}{2}}, \frac{v_h^n + v_h^{n+1}}{2})_D - (B_h \frac{p_h^{n+\frac{1}{2}} + p_h^{n-\frac{1}{2}}}{2}, v_h^n)_D + \dots \quad (1.16) $$

$$ = \frac{1}{4} (B_h p_h^{n-\frac{1}{2}}, v_h^{n-1})_D + \frac{1}{4} (B_h p_h^{n+\frac{1}{2}}, v_h^{n+1})_D - \frac{1}{4} (B_h p_h^{n+\frac{1}{2}}, v_h^n)_D - \frac{1}{4} (B_h p_h^{n-\frac{1}{2}}, v_h^n)_D + \dots \quad (1.17) $$

$$ = \frac{\Delta t}{4} (B_h p_h^{n+\frac{1}{2}}, \delta v_h^{n+\frac{1}{2}})_D - \frac{\Delta t}{4} (B_h p_h^{n-\frac{1}{2}}, \delta v_h^{n-\frac{1}{2}})_D + \frac{1}{2} \left[ (\delta g_h^{n-\frac{1}{2}}, \mu v_h^{n-\frac{1}{2}})_D + (\delta g_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D \right] \\ + (f_h^n, \mu p_h^n)_P \quad (1.18) $$

We denote

$$ \tilde{\mathcal{E}}_{\frac{1}{2},0,0}^{n-\frac{1}{2}} = \mathcal{E}_{\frac{1}{2},0,0}^{n-\frac{1}{2}} - \frac{\Delta t^2}{4} (B_h p_h^{n-\frac{1}{2}}, \delta v_h^{n-\frac{1}{2}})_D = \frac{1}{2} \|p_h^{n-\frac{1}{2}}\|_P^2 + \frac{1}{4} \|v_h^{n-1}\|_D^2 + \frac{1}{4} \|v_h^n\|_D^2 - \frac{\Delta t^2}{4} (B_h p_h^{n-\frac{1}{2}}, \delta v_h^{n-\frac{1}{2}})_D \quad (1.19) $$

which has no defined sign. If $\delta g_h^{n-\frac{1}{2}} = 0$, it is immediately positive under the usual (CFL-condition) since $\delta v_h^{n-\frac{1}{2}} = B_h p_h^{n-\frac{1}{2}}$. Otherwise, the analysis of the positivity (or more precisely the boundedness-by-below) of $\tilde{\mathcal{E}}_{\frac{1}{2},0,0}^n$ relies on bounding $(B_h p_h^{n-\frac{1}{2}}, \delta v_h^{n-\frac{1}{2}})_D$ with the same techniques as above (Cauchy-Schwarz and Young's inequalities):

$$ (B_h p_h^{n-\frac{1}{2}}, \delta v_h^{n-\frac{1}{2}})_D \le (1+\varepsilon) \|B_h p_h^{n-\frac{1}{2}}\|_D^2 + \frac{1}{4\varepsilon} \|\delta g_h^{n-\frac{1}{2}}\|_D^2 \quad (1.20) $$

We get

$$ \tilde{\mathcal{E}}_{\frac{1}{2},0,0}^{n-\frac{1}{2}} \ge \frac{1}{2} (\tilde{I}_{h,\varepsilon} p_h^{n-\frac{1}{2}}, p_h^{n-\frac{1}{2}})_P + \frac{1}{4} \|v_h^{n-1}\|_D^2 + \frac{1}{4} \|v_h^n\|_D^2 - \frac{\Delta t^2}{16\varepsilon} \|\delta g_h^{n-\frac{1}{2}}\|_D^2 \quad (1.21) $$

with $\tilde{I}_{h,\varepsilon} = I_h - \frac{\Delta t^2}{4}(1+\varepsilon)B_h^* B_h$ a modified identity operator which is positive under the $\varepsilon$-CFL condition (1.11), i.e. under the strict CFL condition $\eta < 1$. We get that

$$ \tilde{\mathcal{E}}_{\frac{1}{2},0,0}^{n-\frac{1}{2}} + \frac{\Delta t^2}{16} \frac{1+\eta^2}{1-\eta^2} \|\delta g_h^{n-\frac{1}{2}}\|_D^2 \ge \frac{1-\eta^2}{1+\eta^2} \|p_h^{n-\frac{1}{2}}\|_P^2 + \frac{1}{4} \|v_h^{n-1}\|_D^2 + \frac{1}{4} \|v_h^n\|_D^2 \ge \frac{1-\eta^2}{1+\eta^2} \|p_h^{n-\frac{1}{2}}\|_P^2 + \frac{1}{2} \|\mu v_h^{n-\frac{1}{2}}\|_D^2 \ge 0 \quad (1.22) $$

The modified energy satisfies

$$ \delta \tilde{\mathcal{E}}_{\frac{1}{2},0,0}^{n} = \frac{1}{2} \left[ (\delta g_h^{n-\frac{1}{2}}, \mu v_h^{n-\frac{1}{2}})_D + (\delta g_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D \right] + (f_h^n, \mu p_h^n)_P \quad (1.23) $$

A discrete summation from 0 to N yields the stability of the energy. Under the strict CFL condition, it yields the stability of $\|p_h^{n-\frac{1}{2}}\|_P$ and $\|\mu v_h^{n-\frac{1}{2}}\|_D$ by using (1.22). This method improves the previous one since it is valid with source terms on each equation. However, the stability constants still blow up as $\eta \to 1$.
1.3. **Choice** $\alpha = \beta = 0$ and $\gamma = 1$

The choice $\beta = 0$, $\alpha = 0$ and $\gamma = 1$ is done in the convergence proof of [Jenkins et al., 2002] and is remarkable since it breaks the symmetry of the manipulation. This leads to

$$ \mathcal{E}_{0,0,1}^{n-\frac{1}{2}} = \frac{1}{2} \|p_h^{n-\frac{1}{2}}\|_P^2 + \frac{1}{2} \|v_h^n\|_D^2 \quad (1.24) $$

with

$$ \begin{aligned} \delta\mathcal{E}_{0,0,1}^{n} &= (B_h p_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D - (B_h \mu p_h^n, v_h^n)_D + (\delta g_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D + (f_h^n, \mu p_h^n)_P \\ &= (B_h p_h^{n+\frac{1}{2}}, \frac{v_h^{n+1} + v_h^n}{2})_D - (B_h \frac{p_h^{n+\frac{1}{2}} + p_h^{n-\frac{1}{2}}}{2}, v_h^n)_D + (\delta g_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D + (f_h^n, \mu p_h^n)_P \\ &= \frac{1}{2}(B_h p_h^{n+\frac{1}{2}}, v_h^{n+1})_D - \frac{1}{2}(B_h p_h^{n-\frac{1}{2}}, v_h^n)_D + (\delta g_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D + (f_h^n, \mu p_h^n)_P \end{aligned} \quad (1.25) $$

We denote

$$ \tilde{\mathcal{E}}_{0,0,1}^{n-\frac{1}{2}} = \mathcal{E}_{0,0,1}^{n-\frac{1}{2}} - \frac{\Delta t}{2} (B_h p_h^{n-\frac{1}{2}}, v_h^n)_D = \frac{1}{2} \|p_h^{n-\frac{1}{2}}\|_P^2 + \frac{1}{2} \|v_h^n\|_D^2 - \frac{\Delta t}{2} (B_h p_h^{n-\frac{1}{2}}, v_h^n)_D \quad (1.26) $$

which has no definite sign. However, a Young's inequality with $\varepsilon > 0$ yields

$$ (B_h p_h^{n-\frac{1}{2}}, v_h^n)_D \le \|B_h p_h^{n-\frac{1}{2}}\|_D \|v_h^n\|_D \le \|B_h\| \|p_h^{n-\frac{1}{2}}\|_P \|v_h^n\|_D \le \|B_h\| \left[ \frac{\varepsilon}{2} \|p_h^{n-\frac{1}{2}}\|_P^2 + \frac{1}{2\varepsilon} \|v_h^n\|_D^2 \right] \quad (1.27) $$

Hence

$$ \tilde{\mathcal{E}}_{0,0,1}^{n-\frac{1}{2}} \ge \frac{1}{2} \left( 1 - \frac{\Delta t \|B_h\| \varepsilon}{2} \right) \|p_h^{n-\frac{1}{2}}\|_P^2 + \frac{1}{2} \left( 1 - \frac{\Delta t \|B_h\|}{2\varepsilon} \right) \|v_h^n\|_D^2 \quad (1.28) $$

The largest positivity condition is obtained with $\varepsilon = 1$ and reads $\Delta t \|B_h\| \le 2$. Under the usual (CFL-condition) $\eta \le 1$, we have

$$ \tilde{\mathcal{E}}_{0,0,1}^{n-\frac{1}{2}} \ge \frac{1}{2} (1 - \eta) \|p_h^{n-\frac{1}{2}}\|_P^2 + \frac{1}{2} (1 - \eta) \|v_h^n\|_D^2 \ge 0 \quad (1.29) $$

Note that the modified energy is non-negative even in the equality case. It satisfies

$$ \delta \tilde{\mathcal{E}}_{0,0,1}^{n} = (\delta g_h^{n+\frac{1}{2}}, \mu v_h^{n+\frac{1}{2}})_D + (f_h^n, \mu p_h^n)_P \quad (1.30) $$

which is bounded by the square root of a the modified energy with constants that blow up as $\eta \to 1$. This methods improves the previous ones since a discrete summation yields the stability of a sharper modified energy under the strict CFL condition for source terms on each equation. Moreover, it provides bounds on $\|p_h^{n-\frac{1}{2}}\|_P$ and $\|v_h^n\|_D$ directly. However, the stability constants still blow up as $\eta \to 1$.

## 2. STABILITY AT THE CFL

We provide here two strategies to obtain stability even as $\eta \to 1$. This will not be possible directly on the unknowns of the system, but on post averaged or post processed quantities.
## 2.1. Case $\delta g_h \equiv 0$

We follow the ideas developed for the second order equation in [Joly, 2003]. The leading principle is to use the second equation of (1.6) on two consecutive time steps to compensate the crossing term appearing in a new energy obtained on the system (0.8) averaged between two time steps. More precisely,

$$ \mathcal{H}_{0,1,0}^{n} = \frac{1}{2} \|\mu p_h^n\|_P^2 + \frac{1}{2} (\mu v_h^{n+\frac{1}{2}}, \mu v_h^{n-\frac{1}{2}})_D \quad (2.1) $$

satisfies

$$ \delta \mathcal{H}_{0,1,0}^{n-\frac{1}{2}} = (\mu f_h^{n-\frac{1}{2}}, \mu^2 p_h^{n-\frac{1}{2}})_P \quad (2.2) $$

Let us pose

$$ \begin{aligned} \mathcal{F}_{0,1,0}^n &= \mathcal{H}_{0,1,0}^n + \mathcal{E}_{0,1,0}^{n+\frac{1}{2}} + \mathcal{E}_{0,1,0}^{n-\frac{1}{2}} = \frac{1}{2} \|\mu p_h^n\|_P^2 + \frac{1}{2}(\mu v_h^{n+\frac{1}{2}}, \mu v_h^{n-\frac{1}{2}})_D + \mathcal{E}_{0,1,0}^{n+\frac{1}{2}} + \mathcal{E}_{0,1,0}^{n-\frac{1}{2}} \\ &\geq (1.6) \frac{1}{2} \|\mu p_h^n\|_P^2 + \frac{1}{2}(\mu v_h^{n+\frac{1}{2}}, \mu v_h^{n-\frac{1}{2}})_D + \frac{1}{2} \|\mu v_h^{n+\frac{1}{2}}\|_D^2 + \frac{1}{2} \|\mu v_h^{n-\frac{1}{2}}\|_D^2 \quad \text{if } \eta \leq 1 \\ &\geq \frac{1}{2} \|\mu p_h^n\|_P^2 \end{aligned} \quad (2.3) $$

because

$$ 2(\mu v_h^{n+\frac{1}{2}}, \mu v_h^{n-\frac{1}{2}})_D = \|\mu v_h^{n+\frac{1}{2}} + \mu v_h^{n-\frac{1}{2}}\|_D^2 - \|\mu v_h^{n+\frac{1}{2}}\|_D^2 - \|\mu v_h^{n-\frac{1}{2}}\|_D^2 \quad (2.4) $$

If satisfies

$$ \begin{aligned} \delta \mathcal{F}_{0,1,0}^{n-\frac{1}{2}} &= (\mu f_h^{n-\frac{1}{2}}, \mu^2 p_h^{n-\frac{1}{2}})_P + (f_h^n, \mu p_h^n)_P + (f_h^{n-1}, \mu p_h^{n-1})_P \\ &\leq \|\mu f_h^{n-\frac{1}{2}}\|_P \|\mu^2 p_h^{n-\frac{1}{2}}\|_P + \|f_h^n\|_P \|\mu p_h^n\|_P + \|f_h^{n-1}\|_P \|\mu p_h^{n-1}\|_P \\ &\leq \left[ \frac{1}{2} \|\mu f_h^{n-\frac{1}{2}}\|_P + \|f_h^n\|_P \right] \|\mu p_h^n\|_P + \left[ \frac{1}{2} \|\mu f_h^{n-\frac{1}{2}}\|_P + \|f_h^{n-1}\|_P \right] \|\mu p_h^{n-1}\|_P \\ &\leq \left[ \frac{1}{2} \|\mu f_h^{n-\frac{1}{2}}\|_P + \|f_h^n\|_P \right] \sqrt{2 \mathcal{F}_{0,1,0}^n} + \left[ \frac{1}{2} \|\mu f_h^{n-\frac{1}{2}}\|_P + \|f_h^{n-1}\|_P \right] \sqrt{2 \mathcal{F}_{0,1,0}^{n-1}} \\ &\leq \left[ \|\mu f_h^{n-\frac{1}{2}}\|_P + \|f_h^n\|_P + \|f_h^{n-1}\|_P \right] \sqrt{2} \left[ \sqrt{\mathcal{F}_{0,1,0}^n} + \sqrt{\mathcal{F}_{0,1,0}^{n-1}} \right] \end{aligned} \quad (2.5) $$

Adding the telescopic sum from $n=1$ to $n=N$ we get

$$ \sqrt{\mathcal{F}_{0,1,0}^N} \leq \sqrt{\mathcal{F}_{0,1,0}^{\frac{1}{2}}} + \Delta t \sum_{n=1}^{N} \|\mu f_h^{n-\frac{1}{2}}\|_P + \|f_h^n\|_P + \|f_h^{n-1}\|_P \quad (2.6) $$

which yields uniform stability estimates on the averaged fields:

$$ \|\mu p_h^n\|_P \leq \sqrt{2} \sqrt{\mathcal{F}_{0,1,0}^n} \quad \text{and} \quad \|\mu v_h^{n+\frac{1}{2}}\|_D^2 + \|\mu v_h^{n-\frac{1}{2}}\|_D^2 \leq 2 \mathcal{F}_{0,1,0}^n \quad (2.7) $$

However, trying to bound directly $p_h$ only gets that

$$ \|p_h^{n-\frac{1}{2}}\|_P^2 + \|p_h^{n+\frac{1}{2}}\|_P^2 \leq \frac{2}{1-\eta^2} \mathcal{F}_{0,1,0}^n \quad (2.8) $$

which blows up as $\eta \to 1$.
## 2.2. General case

In this section, we apply the methodology developed in [Chabassier and Imperiale, 2017, Chabassier and Imperiale, 2021] and adapt it to the mixed formulation and the Störmer-Verlet scheme. It is based on the manipulation described above, with $\beta = 1$ and $\alpha = \gamma = 0$ as in [Joly, 2003] and [Egger and Radu, 2020, proof of Lemma 5.3]. A naive manipulation would be to take the scalar product of (0.8c) with $\mu p_h^n$ and of $\mu(0.8d)$ with $v_h^{n+\frac{1}{2}}$. One would get that

$$ \delta\mathcal{E}_{pv,h}^n = (f_h^n, \mu p_h^n)_P + (\mu\delta g_h^n, v_h^n)_D \quad \text{with} \quad \mathcal{E}_{pv,h}^{n+\frac{1}{2}} = \frac{1}{2}(p_h^{n+\frac{1}{2}}, p_h^{n+\frac{1}{2}})_P + \frac{1}{2}(v_h^{n+1}, v_h^n)_D \quad (2.9) $$

Transferring the usual concepts for energy analysis directly here would imply to use (0.7a) and (0.8d) to replace the term $\Delta t^2(\delta v_h^{n+\frac{1}{2}}, \delta v_h^{n+\frac{1}{2}})_D/4$. However it induces a difficulty due to the presence of the source term $\delta g_h$, as mentioned in Sec. 1.1.2. To mitigate this issue, let us proceed to a change of variables and suppose that $g_h$ is regular enough uniformly with $h$:

**Hypothesis 2.1.** In the sequel we will suppose that there exists a constant $C_g > 0$ such that for all $h > 0$,

$$ \|B_h^* g_h\|_{C^0(0,T,P)} \leq C_g \quad (2.10) $$

The couple $(\{p_h^{n+\frac{1}{2}}\}_n, \{w_h^n\}_n)$ defined as $w_h^n = v_h^n - g_h^n \in D_h$ is solution to

$$ w_h^0 = v_h^0 - g_h^0 \quad (2.11a) $$

$$ p_h^{\frac{1}{2}} = p_h(0) + \frac{\Delta t}{2} [f_h^0 - B_h^* v_h(0)] \quad (2.11b) $$

$$ \delta p_h^n + B_h^* w_h^n = f_h^n - B_h^* g_h^n \quad (2.11c) $$

$$ \delta w_h^{n+\frac{1}{2}} - B_h p_h^{n+\frac{1}{2}} = 0 \quad (2.11d) $$

**Proposition 2.1 (Discrete power balance).** Any solution to (0.8) satisfies, for $n \in [1, N]$,

$$ \delta \mathcal{E}_{pw,h}^n = (f_h^n, \mu p_h^n)_P - (B_h^* g_h^n, \mu p_h^n)_P \quad \text{with} \quad \mathcal{E}_{pw,h}^{n+\frac{1}{2}} = \frac{1}{2}(p_h^{n+\frac{1}{2}}, p_h^{n+\frac{1}{2}})_P + \frac{1}{2}(w_h^{n+1}, w_h^n)_D \quad (2.12) $$

where $w_h^n = v_h^n - g_h^n$.

*Proof.* Take the scalar product of (2.11c) with $\mu p_h^n \in U_h$ and of $\mu(2.11d)^n$ with $w_h^n \in D_h$. $\square$

**Proposition 2.2 (Reformulation of the energy).** The discrete energy $\mathcal{E}_{pw,h}^{n+\frac{1}{2}}$ also writes

$$ \mathcal{E}_{pw,h}^{n+\frac{1}{2}} = \frac{1}{2} \tilde{m}(p_h^{n+\frac{1}{2}}, p_h^{n+\frac{1}{2}}) + \frac{1}{2} (\mu w_h^{n+\frac{1}{2}}, \mu w_h^{n+\frac{1}{2}})_D \quad (2.13) $$

where $\tilde{m}$ is the bilinear form defined on $U_h$ such that for all $(p, q) \in U_h \times U_h$,

$$ \tilde{m}(p, q) = (p, q)_P - \frac{\Delta t^2}{4} (B_h p, B_h q)_D = (\tilde{I}_h p, q)_P \quad \text{where} \quad \tilde{I}_h = \mathcal{I}_h - \frac{\Delta t^2}{4} B_h^* B_h \quad (2.14) $$
Proof. Using (0.7a) we get

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

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

$$ = (2.11d) \quad (\mu w_h^{n+\frac{1}{2}}, \mu w_h^{n+\frac{1}{2}})_D - \frac{\Delta t^2}{4} (B_h p_h^{n+\frac{1}{2}}, B_h p_h^{n+\frac{1}{2}})_D \quad (2.17) $$

Hence

$$ \mathcal{E}_{pw,h}^{n+\frac{1}{2}} = \frac{1}{2}(p_h^{n+\frac{1}{2}}, p_h^{n+\frac{1}{2}})_P + \frac{1}{2}(\mu w_h^{n+\frac{1}{2}}, \mu w_h^{n+\frac{1}{2}})_D - \frac{\Delta t^2}{4}(B_h p_h^{n+\frac{1}{2}}, B_h p_h^{n+\frac{1}{2}})_D \quad (2.18) $$

$$ = \frac{1}{2} \tilde{m}(p_h^{n+\frac{1}{2}}, p_h^{n+\frac{1}{2}}) + \frac{1}{2}(\mu w_h^{n+\frac{1}{2}}, \mu w_h^{n+\frac{1}{2}})_D \quad (2.19) $$

□

**Proposition 2.3** (CFL condition). The energy (2.13) is positive if

$$ \eta \le 1 \quad (2.20) $$

where $\eta$ is defined in (CFL-condition).

**Corollary 2.1.** Suppose (2.20) is satisfied. Then,

$$ \tilde{m}(p_h^{n+\frac{1}{2}}, p_h^{n+\frac{1}{2}}) \le 2\mathcal{E}_{pw,h}^{n+\frac{1}{2}}, \quad (\mu w_h^{n+\frac{1}{2}}, \mu w_h^{n+\frac{1}{2}})_D \le 2\mathcal{E}_{pw,h}^{n+\frac{1}{2}} \quad (2.21) $$

To show the stability uniformly as $\Delta t$ approaches its greatest allowed value, let us exploit the spectral properties of the operator $B_h^*B_h$ as in [Chabassier and Imperiale, 2017].

**Proposition 2.4** (Spectral decomposition of $B_h^*B_h$). The operator $B_h^*B_h: U_h \to U_h$ is diagonalizable in $\mathbb{R}$. We call $(\lambda_{h,i}, e_{h,i})$ its eigenpairs which are chosen orthonormal in $P$.

Proof. In finite dimensional spaces, any symmetric real operator is diagonalizable in an orthonormal basis. □

Following [Chabassier and Imperiale, 2017], we introduce the polynomial

$$ P_k(x) = 1 - \frac{1}{4}x \quad (2.22) $$

which is non-negative on the interval $[0, 4]$. Then $\tilde{\mathcal{I}}_h$ can be expressed as a polynomial of the operator $\Delta t^2 B_h^*B_h$:

$$ \tilde{\mathcal{I}}_h = \mathcal{I}_h - \frac{\Delta t^2}{4} B_h^* B_h = P_k(\Delta t^2 B_h^* B_h) \quad (2.23) $$

**Proposition 2.5** (Partitionning). The interval $[0, 4]$ can be partitionned as $J_k \cup J_p$ with $J_k \cap J_p = \emptyset$ such that there exist $C_k > 0$ and $C_p > 0$ such that for all $x \in J_k$, $P_k(x) \ge C_k$ and for all $x \in J_p$, $x \ge C_p$.

Proof. See [Chabassier and Imperiale, 2017, appendix] for the scheme called “TS” with $\theta = 0$. We get $C_k = \frac{1}{2}$ and $C_p = 2$. □
**Proposition 2.6.** Let us define the two projectors $\Pi_k$ and $\Pi_p$ such that for all $q_h \in U_h$

$$ \Pi_k q_h = \sum_{\substack{\Delta t^2 \lambda_{h,i} \in J_k \\ \lambda_{h,i} \in Sp(\mathcal{K}_h)}} (q_h, e_{h,i})_P e_{h,i}, \quad \Pi_p q_h = \sum_{\substack{\Delta t^2 \lambda_{h,i} \in J_p \\ \lambda_{h,i} \in Sp(\mathcal{K}_h)}} (q_h, e_{h,i})_P e_{h,i} \quad (2.24) $$

Then, for all $q_h \in U_h$

$$ \|\Pi_k q_h\|_P^2 \le C_k^{-1} (\tilde{\mathcal{I}}_h q_h, q_h)_P \quad (2.25a) $$

$$ \|\Pi_p q_h\|_P^2 \le \Delta t^2 C_p^{-1} (B_h q_h, B_h q_h)_D \quad (2.25b) $$

*Proof.* See [Chabassier and Imperiale, 2017]. □

**Proposition 2.7.** Any $(p_h, w_h) \in U_h \times D_h$ solution to (2.11) satisfies

$$ \|\Pi_k \mu p_h^n\|_P \le \frac{C_k^{-1/2}}{2} \left( \sqrt{2\mathcal{E}_{pw,h}^{n+\frac{1}{2}}} + \sqrt{2\mathcal{E}_{pw,h}^{n-\frac{1}{2}}} \right) \quad (2.26a) $$

$$ \|\Pi_p \mu p_h^n\|_P \le C_p^{-1/2} \left( \sqrt{2\mathcal{E}_{pw,h}^{n+\frac{1}{2}}} + \sqrt{2\mathcal{E}_{pw,h}^{n-\frac{1}{2}}} \right) \quad (2.26b) $$

*Proof.* Take $q_h = \mu p_h^n$ in (2.25a). Using the triangular inequality and (2.21), we directly have

$$ (\tilde{\mathcal{I}}_h \mu p_h^n, \mu p_h^n)_D^{1/2} \le \frac{1}{2} (\tilde{\mathcal{I}}_h p_h^{n+\frac{1}{2}}, p_h^{n+\frac{1}{2}})_D^{1/2} + \frac{1}{2} (\tilde{\mathcal{I}}_h p_h^{n-\frac{1}{2}}, p_h^{n-\frac{1}{2}})_D^{1/2} \le \frac{1}{2} \sqrt{2\mathcal{E}_{pw,h}^{n+\frac{1}{2}}} + \frac{1}{2} \sqrt{2\mathcal{E}_{pw,h}^{n-\frac{1}{2}}} \quad (2.27) $$

Take $q_h = \mu p_h^n$ in (2.25b). Since $\delta\mu w_h^n = B_h\mu p_h^n$ from $\mu(2.11d)$,

$$ \Delta t (B_h \mu p_h^n, B_h \mu p_h^n)_D^{1/2} = \Delta t (\delta\mu w_h^n, \delta\mu w_h^n)_D^{1/2} \le \left( (\mu w_h^{n+\frac{1}{2}}, \mu w_h^{n+\frac{1}{2}})_D^{1/2} + (\mu w_h^{n-\frac{1}{2}}, \mu w_h^{n-\frac{1}{2}})_D^{1/2} \right) \quad (2.28) $$

From (2.21) we then have

$$ \Delta t (B_h \mu p_h^n, B_h \mu p_h^n)_D^{1/2} \le \sqrt{2\mathcal{E}_{pw,h}^{n+\frac{1}{2}}} + \sqrt{2\mathcal{E}_{pw,h}^{n-\frac{1}{2}}} \quad (2.29) $$

□

**Proposition 2.8** (Majoration of the averaged unknowns).

$$ \|\mu w_h^{n+\frac{1}{2}}\|_D \le \sqrt{2\mathcal{E}_{pw,h}^{n+\frac{1}{2}}}, \quad \|\mu p_h^n\|_P \le \gamma \left( \sqrt{2\mathcal{E}_{pw,h}^{n+\frac{1}{2}}} + \sqrt{2\mathcal{E}_{pw,h}^{n-\frac{1}{2}}} \right) \quad (2.30) $$

where

$$ \gamma = \frac{C_k^{-1/2}}{2} + C_p^{-1/2} \quad (2.31) $$

*Proof.* For $\mu w_h^{n+\frac{1}{2}}$ the result is immediate using (2.21). For $\mu p_h^n$ we decompose it as $\Pi_k \mu p_h^n + \Pi_p \mu p_h^n$, and use the orthogonality between the projection spaces. We then use (2.26) to get the expected result. □

**Proposition 2.9** (Stability of the energy). The energy satisfies

$$ \sqrt{\mathcal{E}_{pw,h}^{n+\frac{1}{2}}} \le \sqrt{\mathcal{E}_{pw,h}^{\frac{1}{2}}} + \sqrt{2\gamma\Delta t} \sum_{k=1}^{n} [\|f_h^k\|_P + \|B_h^* g_h^k\|_P] \quad (2.32) $$
where

$$ \sqrt{\mathcal{E}_{pw,h}^{\frac{1}{2}}} \le \frac{1 + \sqrt{2}}{2} \|p_h(0)\|_P + \frac{1 + \sqrt{2}}{2} \frac{\Delta t}{2} \|f_h^0\|_P + \frac{1 + 2\sqrt{2}}{2} \|v_h(0)\|_D + \frac{1}{\sqrt{2}} \|g_h^0\|_D \quad (2.33) $$

*Proof.* Apply the Cauchy-Schwarz inequality on (2.12). Since $g_h$ satisfies Hyp. 2.1, we have

$$ \frac{\mathcal{E}_{pw,h}^{n+\frac{1}{2}} - \mathcal{E}_{pw,h}^{n-\frac{1}{2}}}{\Delta t} \le (\|f_h^n\|_P + \|B_h^* g_h^n\|_P) \|\mu p_h^n\|_P \quad (2.34) $$

Use (2.30) to get

$$ \frac{\left(\sqrt{\mathcal{E}_{pw,h}^{n+\frac{1}{2}}} + \sqrt{\mathcal{E}_{pw,h}^{n-\frac{1}{2}}}\right) \left(\sqrt{\mathcal{E}_{pw,h}^{n+\frac{1}{2}}} - \sqrt{\mathcal{E}_{pw,h}^{n-\frac{1}{2}}}\right)}{\Delta t} \le (\|f_h^n\|_P + \|B_h^* g_h^n\|_P) \gamma \sqrt{2} \left(\sqrt{\mathcal{E}_{pw,h}^{n+\frac{1}{2}}} + \sqrt{\mathcal{E}_{pw,h}^{n-\frac{1}{2}}}\right) \quad (2.35) $$

Cancel out to get

$$ \sqrt{\mathcal{E}_{pw,h}^{n+\frac{1}{2}}} \le \sqrt{\mathcal{E}_{pw,h}^{n-\frac{1}{2}}} + \Delta t (\|f_h^n\|_P + \|B_h^* g_h^n\|_P) \|\gamma \sqrt{2}\| \quad (2.36) $$

Add this telescopic sum from $k=1$ to $n$ to get the first expected result:

$$ \sqrt{\mathcal{E}_{pw,h}^{n+\frac{1}{2}}} \le \sqrt{\mathcal{E}_{pw,h}^{\frac{1}{2}}} + \sqrt{2\gamma\Delta t} \sum_{k=1}^{n} [\|f_h^k\|_P + \|B_h^* g_h^k\|_P] \quad (2.37) $$

Finally, let us express $\sqrt{\mathcal{E}_{pw,h}^{\frac{1}{2}}}$ by using the expression (2.13):

$$ \sqrt{\mathcal{E}_{pw,h}^{\frac{1}{2}}} = \left[ \frac{1}{2} \tilde{m}(p_h^{\frac{1}{2}}, p_h^{\frac{1}{2}}) + \frac{1}{2} (\mu w_h^{\frac{1}{2}}, \mu w_h^{\frac{1}{2}})_D \right]^{1/2} \le \frac{1}{\sqrt{2}} \tilde{m}(p_h^{\frac{1}{2}}, p_h^{\frac{1}{2}})^{1/2} + \frac{1}{\sqrt{2}} \|\mu w_h^{\frac{1}{2}}\|_D \quad (2.38) $$

First,

$$ \tilde{m}(p_h^{\frac{1}{2}}, p_h^{\frac{1}{2}}) = \frac{1}{2} \|p_h^{\frac{1}{2}}\|_P^2 - \frac{\Delta t^2}{4} \|B_h p_h^{\frac{1}{2}}\|_D^2 \le \frac{1}{2} \|p_h^{\frac{1}{2}}\|_P^2 \quad (2.39) $$

Moreover,

$$ \|\mu w_h^{\frac{1}{2}}\|_D \le \frac{1}{2} \|w_h^0\|_D + \frac{1}{2} \|w_h^1\|_D \quad (2.40) $$

where

$$ w_h^0 = v_h(0) - g_h^0 \quad \Rightarrow \quad \|w_h^0\|_D = \|v_h(0)\|_D + \|g_h^0\|_D \quad \text{and} \quad w_h^1 = w_h^0 + \Delta t B_h p_h^{\frac{1}{2}} \quad (2.41) $$

Hence

$$ \|\mu w_h^{\frac{1}{2}}\|_D \le \frac{1}{2} \|w_h^0\|_D + \frac{1}{2} \left[ \|w_h^0\|_D + \Delta t \|B_h p_h^{\frac{1}{2}}\|_D \right] \quad (2.42) $$

Since (CFL-condition) is satisfied, we have that

$$ \eta = \frac{\Delta t}{2} \|B_h\| = \frac{\Delta t}{2} \sup_{p_h \in U_h} \frac{\|B_h p_h\|_D}{\|p_h\|_P} \le 1 \quad \Rightarrow \quad \frac{\Delta t}{2} \|B_h p_h\|_D \le \eta \|p_h\|_P \le \|p_h\|_P, \forall p_h \in U_h \quad (2.43) $$

So

$$ \|\mu w_h^{\frac{1}{2}}\|_D \le \|w_h^0\|_D + \|p_h^{\frac{1}{2}}\|_P \quad (2.44) $$

Summing up,

$$ \sqrt{\mathcal{E}_{pw,h}^{\frac{1}{2}}} \le \frac{1}{2} \|p_h^{\frac{1}{2}}\|_P + \frac{1}{\sqrt{2}} \|w_h^0\|_D + \frac{1}{\sqrt{2}} \|p_h^{\frac{1}{2}}\|_P \le \frac{1 + \sqrt{2}}{2} \|p_h^{\frac{1}{2}}\|_P + \frac{1}{\sqrt{2}} \|v_h(0)\|_D + \frac{1}{\sqrt{2}} \|g_h^0\|_D \quad (2.45) $$
We write that

$$p_h^{\frac{1}{2}} = p_h(0) + \frac{\Delta t}{2} [f_h^0 - B_h^* v_h(0)] \Rightarrow \|p_h^{\frac{1}{2}}\|_P \le \|p_h(0)\|_P + \frac{\Delta t}{2} \|f_h^0\|_P + \frac{\Delta t}{2} \|B_h^* v_h(0)\|_P \quad (2.46)$$

Since (CFL-condition) is satisfied, we also have that

$$\eta = \frac{\Delta t}{2} \|B_h\| = \frac{\Delta t}{2} \|B_h^*\| = \frac{\Delta t}{2} \sup_{v_h \in D_h} \frac{\|B_h^* v_h\|_P}{\|v_h\|_D} \le 1 \Rightarrow \frac{\Delta t}{2} \|B_h^* v_h\|_P \le \eta \|v_h\|_D \le \|v_h\|_D \quad \forall v_h \in D_h \quad (2.47)$$

which yields

$$\|p_h^{\frac{1}{2}}\|_P \le \|p_h(0)\|_P + \frac{\Delta t}{2} \|f_h^0\|_P + \|v_h(0)\|_D \quad (2.48)$$

Hence

$$\sqrt{\varepsilon_{pw,h}^{\frac{1}{2}}} \le \frac{1 + \sqrt{2}}{2} \left[ \|p_h(0)\|_P + \frac{\Delta t}{2} \|f_h^0\|_P + \|v_h(0)\|_D \right] + \frac{1}{\sqrt{2}} \|v_h(0)\|_D + \frac{1}{\sqrt{2}} \|g_h^0\|_D \quad (2.49)$$

$$\le \frac{1 + \sqrt{2}}{2} \|p_h(0)\|_P + \frac{1 + \sqrt{2}}{2} \frac{\Delta t}{2} \|f_h^0\|_P + \frac{1 + 2\sqrt{2}}{2} \|v_h(0)\|_D + \frac{1}{\sqrt{2}} \|g_h^0\|_D \quad (2.50)$$

□

**Proposition 2.10** (Stability of the averaged unknowns).

$$\|\mu w_h^{n+\frac{1}{2}}\|_D \le \sqrt{2} \sqrt{\varepsilon_{pw,h}^{\frac{1}{2}}} + \sqrt{2\gamma} \Delta t \sum_{k=1}^{n} \left[ \|f_h^k\|_P + \|B_h^* g_h^k\|_P \right], \quad (2.51)$$

$$\|\mu p_h^n\| \le 2\sqrt{2\gamma} \left( \sqrt{\varepsilon_{pw,h}^{\frac{1}{2}}} + \sqrt{2\gamma} \Delta t \sum_{k=1}^{n} \left[ \|f_h^k\|_P + \|B_h^* g_h^k\|_P \right] \right) \quad (2.52)$$

$$\|\mu v_h^{n+\frac{1}{2}}\|_D \le \sqrt{2} \sqrt{\varepsilon_{pw,h}^{\frac{1}{2}}} + \sqrt{2\gamma} \Delta t \sum_{k=1}^{n} \left[ \|f_h^k\|_P + \|B_h^* g_h^k\|_P \right] + \frac{\|g_h^{n+1}\|_D + \|g_h^n\|_D}{2} \quad (2.53)$$

*Proof.* Use jointly Prop. 2.8 and 2.9 to get the two first equations. Moreover, the initial system was posed on the unknown $v_h^n = w_h^n + g_h^n$. Using the previous equations, one can only hope a direct control on $\mu v_h^{n+\frac{1}{2}}$ as

$$\|\mu v_h^{n+\frac{1}{2}}\|_D \le \|\mu w_h^{n+\frac{1}{2}}\|_D + \|\mu g_h^{n+\frac{1}{2}}\|_D \quad (2.54)$$

□

*Remark 2.1.* This result generalizes the result recalled in Sec. 2.1 to the case of two source terms. It confirms that the uniform stability is only obtained on the averaged unknowns.

*Remark 2.2.* One can not hope for more control on the unknowns $p_h^{n+\frac{1}{2}}$ and $w_h^n$, unfortunately, at least without supposing more restrictive bounds on $\Delta t$. Only the $\Pi_k$-projection of $p_h^{n+\frac{1}{2}}$ is controlled directly. However, a global control is achieved on a post-processed field.

**Proposition 2.11** (Stability of the underlying field). Let $\{u_h^n\}_{0 \le n \le N+1}$ be the series of elements of $U_h$ defined as

$$u_h^{n+1} = u_h^n + \Delta t p_h^{n+\frac{1}{2}}, \quad u_h^0 \text{ such that } B_h u_h^0 = w_h^0 \quad (2.55)$$
Then,

$$ \|u_h^{n+1}\| \le \sqrt{2}\|u_h^0\| + t^{n+1}2\gamma\sqrt{2\epsilon_{pw,h}^{\frac{1}{2}}} + 4\gamma^2\Delta t^2 \sum_{j=1}^{n+1} \sum_{k=1}^{j} [\|f^k\|_P + \|B_h^*g_h^k\|_P] \quad (2.56) $$

*Proof.* From the orthogonality between the projection spaces, we get that

$$ \|u_h^{n+1}\|_P = \|\Pi_k u_h^{n+1}\|_P + \|\Pi_p u_h^{n+1}\|_P \quad (2.57) $$

Here we treat the two terms differently. For the first one we use (2.25a) with $q_h = p_h^{n+\frac{1}{2}}$ to get

$$ \|\Pi_k u_h^{n+1}\|_P \le \|\Pi_k u_h^n\|_P + \Delta t \|\Pi_k \delta u_h^{n+\frac{1}{2}}\|_P \le \|\Pi_k u_h^n\|_P + \Delta t \|\Pi_k p_h^{n+\frac{1}{2}}\|_P \le \|\Pi_k u_h^n\|_P + \Delta t C_k^{-1/2} \sqrt{2\epsilon_{pw,h}^{n+\frac{1}{2}}} \quad (2.58) $$

For the second one, we use (2.25b) with $q_h = \mu u_h^{n+\frac{1}{2}}$ to get

$$ \|\Pi_p u_h^{n+1}\|_P \le \|\Pi_p u_h^n\|_P + 2\|\Pi_p \mu u_h^{n+\frac{1}{2}}\|_P \le \|\Pi_p u_h^n\|_P + 2\Delta t C_p^{-1/2} (B_h \mu u_h^{n+\frac{1}{2}}, B_h \mu u_h^{n+\frac{1}{2}})_D \quad (2.59) $$

From the definition of $u_h$ in (2.55) and using (2.11d), we get that $B_h \mu u_h^{n+\frac{1}{2}} = \mu w_h^{n+\frac{1}{2}}$. Hence

$$ \|\Pi_p u_h^{n+1}\|_P \le \|\Pi_p u_h^n\|_P + 2\Delta t C_p^{-1/2} (\mu w_h^{n+\frac{1}{2}}, \mu w_h^{n+\frac{1}{2}})_D \le \|\Pi_p u_h^n\|_P + 2\Delta t C_p^{-1/2} \sqrt{2\epsilon_{pw,h}^{n+\frac{1}{2}}} \quad (2.60) $$

Both telescopic sums are canceled up from 0 to $n$, and then added up.

$$ \|u_h^{n+1}\|_P \le \sqrt{2}\|u_h^0\|_P + 2\gamma\Delta t \sum_{k=0}^{n} \sqrt{2\epsilon_{pw,h}^{k+\frac{1}{2}}} \quad (2.61) $$

Use (2.37) to conclude.

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

*Remark 2.3.* In practice, $\gamma = \sqrt{2}$, see [Chabassier and Imperiale, 2017, appendix] for the scheme called “TS” with $\theta = 0$.

### 3. CONVERGENCE AT THE CFL

Subtracting directly Eq. (0.3) to (0.8) yields truncation errors on each appearance of the operator $\delta$ and needs additional manipulation to use the results of the previous section, necessitating non optimal assumptions on the semi discrete field. Instead, we propose to introduce the underlying field at both the fully and semi discrete level, and to manipulate the resulting equations before analysis. Let $u_h \in U_h$ be defined as

$$ \begin{cases} \dot{u}_h = p_h \\ u_h(0) = u_{h,0} \end{cases} \quad \text{such that} \quad B_h u_{h,0} + g_h^0 = v_h(0) \quad (3.1b) $$

Inserting this field in Eq. (0.3) yields

$$ \begin{cases} \ddot{u}_h + B_h^* v_h = f_h \\ \dot{v}_h - B_h \dot{u}_h = \dot{g}_h \end{cases} \quad (3.2a) $$

The second line can be integrated in time to get

$$ v_h - v_{h,0} - B_h u_h + B_h u_{h,0} = g_h - g_h^0 \quad (3.3) $$
Since the initial conditions respect (3.1), we get

$$ \ddot{u}_h + B_h^* v_h = f_h \qquad (3.4a) $$

$$ v_h - B_h u_h = g_h \qquad (3.4b) $$

$$ \dot{u}_h = p_h \qquad (3.4c) $$

Note that Hyp. 0.2 amounts to supposing that

$$ \|u_h\|_{C^4(0,T;P)} \leq C \qquad (3.5) $$

On the fully discrete level, as in the previous section, let us define $u_h \in U_h$ as

$$ \delta u_h^{n+\frac{1}{2}} = p_h^{n+\frac{1}{2}} \qquad (3.6a) $$

$$ u_h^0 = u_h(0) \qquad (3.6b) $$

Inserting this field in Eq. (0.8) yields

$$ \delta^2 u_h^n + B_h^* v_h^n = f_h^n \qquad (3.7a) $$

$$ \delta v_h^{n+\frac{1}{2}} - B_h \delta u_h^{n+\frac{1}{2}} = \delta g_h^{n+\frac{1}{2}} \qquad (3.7b) $$

The second line can be summed up in time to get

$$ \delta^2 u_h^n + B_h^* v_h^n = f_h^n \qquad (3.8a) $$

$$ v_h^n - B_h^* u_h^n = g_h^n \qquad (3.8b) $$

$$ \delta u_h^{n+\frac{1}{2}} = p_h^{n+\frac{1}{2}} \qquad (3.8c) $$

Let $\bar{u}_h^n = u_h(t^n)$ and $\bar{v}_h^n = v_h(t^n)$, where $v_h$ is solution to the semi-discrete equation (0.3), $u_h$ is the solution to the semi-discrete equation (3.1), therefore solutions to the semi-discrete system (3.4).

**Theorem 3.1.** Suppose that (CFL-condition) is satisfied. Then there exists a constant $C > 0$ independent of $h$ and $\eta$, such that

$$ \|\bar{p}_h^n - \mu p_h^n\|_P \leq C \Delta t^2 \|p_h\|_{C^3(0,T,P)} \qquad (3.9) $$

$$ \|\bar{v}_h^{n+\frac{1}{2}} - \mu v_h^{n+\frac{1}{2}}\|_D \leq C \Delta t^2 (\|v_h\|_{C^2(0,T;D)} + \|p_h\|_{C^3(0,T,P)}) \qquad (3.10) $$

$$ \|\bar{u}_h^{n+1} - u_h^{n+1}\|_P \leq C \Delta t^2 \|p_h\|_{C^3(0,T;P)} \qquad (3.11) $$

*Proof.* Let us subtract the two first lines of systems (3.4) and (3.8).

$$ \ddot{\bar{u}}_h^n - \delta^2 u_h^n + B_h^* \bar{v}_h^n - B_h^* v_h^n = f_h^n - f_h^n \qquad (3.12a) $$

$$ \bar{v}_h^n - v_h^n - B_h \bar{u}_h^n + B_h u_h^n = g_h^n - g_h^n \qquad (3.12b) $$

Let us define the error terms $e_u^n = \bar{u}_h^n - u_h^n$ and $e_v^n = \bar{v}_h^n - v_h^n$, which satisfy

$$ \delta^2 e_u^n + B_h^* e_v^n = \varepsilon_u^n \qquad (3.13a) $$

$$ e_v^n - B_h e_u^n = 0 \qquad (3.13b) $$
where a Taylor expansion shows that there exists $t_1^n \in [t^{n-\frac{1}{2}}, t^{n+\frac{1}{2}}]$ such that

$$ \epsilon_u^n = \frac{\Delta t^2}{12} \bar{u}_h^{(4)}(t_1^n) \qquad (3.14) $$

We now introduce the sequence $\eta_p^{n+\frac{1}{2}}$ such that

$$ \delta e_u^{n+\frac{1}{2}} = \eta_p^{n+\frac{1}{2}} \qquad (3.15) $$

Notice that $\eta_p^{n+\frac{1}{2}}$ is not equal to $e_p^{n+\frac{1}{2}} = \bar{p}_h^{n+\frac{1}{2}} - p_h^{n+\frac{1}{2}}$ since this quantity satisfies

$$ \delta e_u^{n+\frac{1}{2}} = e_p^{n+\frac{1}{2}} + \dot{\epsilon}_u^{n+\frac{1}{2}} \qquad (3.16) $$

where there exists $t_2^{n+\frac{1}{2}} \in [t^n, t^{n+1}]$ such that

$$ \dot{\epsilon}_u^{n+\frac{1}{2}} = \frac{\Delta t^2}{12} \bar{u}_h^{(3)}(t_2^{n+\frac{1}{2}}) \qquad (3.17) $$

Hence, in $U_h$,

$$ \eta_p^{n+\frac{1}{2}} = e_p^{n+\frac{1}{2}} + \dot{\epsilon}_u^{n+\frac{1}{2}} \qquad (3.18) $$

Finally, the error terms $e_u, e_v$ and $\eta_p$ are solution to

$$ \delta \eta_p^n + B_h^* e_v^n = \epsilon_u^n \qquad (3.19a) $$

$$ \delta e_v^{n+\frac{1}{2}} - B_h \eta_p^{n+\frac{1}{2}} = 0 \qquad (3.19b) $$

$$ \delta e_u^{n+\frac{1}{2}} = \eta_p^{n+\frac{1}{2}} \qquad (3.19c) $$

The result follows from the stability analysis performed in Sec. 2.2. Indeed, the two first lines of system (3.19) are of the form (0.8-2.55), with the substitutions

$$ \begin{aligned} p_h^{n+\frac{1}{2}} &\leftarrow \eta_p^{n+\frac{1}{2}} \\ v_h^n &\equiv w_h^n \leftarrow e_v^n \\ u_h^n &\leftarrow e_u^n \\ f^n &\leftarrow \epsilon_u^n \\ g_h^n &\leftarrow 0 \end{aligned} $$

Hence,

$$ \|\mu e_v^{n+\frac{1}{2}}\|_D \le \sqrt{2} \sqrt{\mathcal{E}_e^{\frac{1}{2}}} + \sqrt{2\gamma} \Delta t \sum_{k=1}^{n} [\|\epsilon_u^k\|_P], \qquad (3.20) $$

$$ \|\mu \eta_p^n\| \le 2\sqrt{2\gamma} \left( \sqrt{\mathcal{E}_e^{\frac{1}{2}}} + \sqrt{2\gamma} \Delta t \sum_{k=1}^{n} [\|\epsilon_u^k\|_P] \right) \qquad (3.21) $$

$$ \|e_u^{n+1}\| \le \sqrt{2}\|e_u^0\| + t^{n+1}2\gamma\sqrt{2\mathcal{E}_e^{\frac{1}{2}}} + 4\gamma^2\Delta t^2 \sum_{j=1}^{n+1} \sum_{k=1}^{j} [\|\epsilon_u^k\|_P] \qquad (3.22) $$
where $\mathcal{E}_e^{\frac{1}{2}}$ is defined as

$$ \mathcal{E}_e^{\frac{1}{2}} = \frac{1}{2}\tilde{m}(\eta_p^{\frac{1}{2}}, \eta_p^{\frac{1}{2}}) + \frac{1}{2}(\mu e_v^{\frac{1}{2}}, \mu e_v^{\frac{1}{2}}) \qquad (3.23) $$

Therefore, using (3.18),

$$ \|\mu e_p^n\|_P \le 2\sqrt{2\gamma} \left( \sqrt{\mathcal{E}_e^{\frac{1}{2}}} + \sqrt{2\gamma}\Delta t \sum_{k=1}^{n} \left[ \|\varepsilon_u^k\|_P \right] \right) + \frac{\|\dot{\varepsilon}_u^{n+\frac{1}{2}}\|_P + \|\dot{\varepsilon}_u^{n-\frac{1}{2}}\|_P}{2} \qquad (3.24) $$

So,

$$ \|\bar{v}_h^{n+\frac{1}{2}} - \mu v_h^{n+\frac{1}{2}}\|_D \le \|\bar{v}_h^{n+\frac{1}{2}} - \mu \bar{v}_h^{n+\frac{1}{2}}\|_D + \|\mu \bar{v}_h^{n+\frac{1}{2}} - \mu v_h^{n+\frac{1}{2}}\|_D \qquad (3.25) $$

$$ \le \|\varepsilon_v^{n+\frac{1}{2}}\|_D + \|\mu \bar{v}_h^{n+\frac{1}{2}} - \mu v_h^{n+\frac{1}{2}}\|_D \qquad (3.26) $$

where there exists $t_3^{n+\frac{1}{2}} \in [t^n, t^{n+1}]$ such that

$$ \varepsilon_v^{n+\frac{1}{2}} = \mu \bar{v}_h^{n+\frac{1}{2}} - \bar{v}_h^{n+\frac{1}{2}} = \frac{\Delta t^2}{3} \bar{v}_h^{(2)}(t_3^{n+\frac{1}{2}}) \qquad (3.27) $$

Similarly,

$$ \|\bar{p}_h^n - \mu p_h^n\|_P \le \|\bar{p}_h^n - \mu \bar{p}_h^n\|_P + \|\mu \bar{p}_h^n - \mu p_h^n\|_P \qquad (3.28) $$

$$ \le \|\varepsilon_p^n\|_P + \|\mu \bar{p}_h^n - \mu p_h^n\|_P \qquad (3.29) $$

where there exists $t_4^n \in [t^{n-\frac{1}{2}}, t^{n+\frac{1}{2}}]$ such that

$$ \varepsilon_p^n = \mu \bar{p}_h^n - \bar{p}_h^n = \frac{\Delta t^2}{3} \bar{p}_h^{(2)}(t_4^n) \qquad (3.30) $$

Hence,

$$ \|\bar{v}_h^{n+\frac{1}{2}} - \mu v_h^{n+\frac{1}{2}}\|_D \le \|\varepsilon_v^{n+\frac{1}{2}}\|_D + \sqrt{2\sqrt{\mathcal{E}_e^{\frac{1}{2}}}} + \sqrt{2\gamma}\Delta t \sum_{k=1}^{n} \left[ \|\varepsilon_u^k\|_P \right], \qquad (3.31) $$

$$ \|\bar{p}_h^n - \mu p_h^n\|_P \le \|\varepsilon_p^n\|_P + 2\sqrt{2\gamma} \left( \sqrt{\mathcal{E}_e^{\frac{1}{2}}} + \sqrt{2\gamma}\Delta t \sum_{k=1}^{n} \left[ \|\varepsilon_u^k\|_P \right] \right) + \frac{\|\dot{\varepsilon}_u^{n+\frac{1}{2}}\|_P + \|\dot{\varepsilon}_u^{n-\frac{1}{2}}\|_P}{2} \qquad (3.32) $$

$$ \|\bar{u}_h^{n+1} - u_h^{n+1}\|_P \le \sqrt{2}\|e_u^0\|_P + t^{n+1}2\gamma\sqrt{2\mathcal{E}_e^{\frac{1}{2}}} + 4\gamma^2\Delta t^2 \sum_{j=1}^{n+1} \sum_{k=1}^{j} \left[ \|\varepsilon_u^k\|_P \right] \qquad (3.33) $$

Next, let us show that $\mathcal{E}_e^{\frac{1}{2}}$ is small. Using the naive formulation of the energy (2.12), we write that

$$ \mathcal{E}_e^{\frac{1}{2}} = \frac{1}{2}(\eta_p^{\frac{1}{2}}, \eta_p^{\frac{1}{2}})_P + \frac{1}{2}(e_v^{1}, e_v^{0})_D \le \frac{1}{2}\|\eta_p^{\frac{1}{2}}\|_P^2 + \frac{1}{2}\|e_v^{1}\|_D \|e_v^{0}\|_D \qquad (3.34) $$

Since the initial conditions $v_h^0$ and $p_h^{\frac{1}{2}}$ are defined as (0.8a) and (0.8b), we directly have that $e_v^0 = 0$. It remains to estimate $\eta_p^{\frac{1}{2}} = e_p^{\frac{1}{2}} + \dot{\varepsilon}_u^{\frac{1}{2}}$. The value of $e_p^{\frac{1}{2}}$ can be found by subtracting both sides of (0.8b) to $\bar{p}_h^{\frac{1}{2}} - \frac{\Delta t}{2} B_h^* v_h^0$
in order to make also appear $e_v^0$.

$$ \underbrace{\bar{p}_h^{\frac{1}{2}} - p_h^{\frac{1}{2}}}_{e_p^{\frac{1}{2}}} + \frac{\Delta t}{2} B_h^* \underbrace{(\bar{v}_h^0 - v_h^0)}_{e_v^0=0} = \frac{\Delta t}{2} \underbrace{[B_h^* \bar{v}_h^0 - f_h^0]}_{= (0.3c)} + \bar{p}_h^{\frac{1}{2}} - p_h(0) \quad (3.35) $$

A Taylor expansion of $\bar{p}_h$ yields that there exists $t_5^0 \in [0, t^{\frac{1}{2}}]$ such that

$$ \bar{p}_h(t^{\frac{1}{2}}) = \bar{p}_h(0) + \frac{\Delta t}{2} \bar{p}_h^{(1)}(0) + \frac{1}{2} \left( \frac{\Delta t}{2} \right)^2 \bar{p}_h^{(2)}(t_5^0) \quad (3.36) $$

Hence we get that

$$ e_p^{\frac{1}{2}} = \frac{\Delta t^2}{8} \bar{p}_h^{(2)}(t_5^0) \quad (3.37) $$

Take $\tilde{p}_h = e_p^{\frac{1}{2}}$ and apply Cauchy Schwarz inequality to get that

$$ \|e_p^{\frac{1}{2}}\|_P \le \frac{\Delta t^2}{8} \|p_h\|_{C^2(0,T;P)} \quad (3.38) $$

So

$$ \sqrt{2\mathcal{E}_e^{\frac{1}{2}}} \le \|e_p^{\frac{1}{2}}\|_P + \|\dot{\epsilon}_u^{\frac{1}{2}}\|_P \quad (3.39) $$

$$ \le \frac{5\Delta t^2}{24} \|p_h\|_{C^2(0,T;P)} \quad (3.40) $$

Finally, we have that

$$ \epsilon_u^n = \frac{\Delta t^2}{12} \bar{p}_h^{(3)}(t_1^n) \quad \text{and} \quad \dot{\epsilon}_u^{n+\frac{1}{2}} = \frac{\Delta t^2}{12} \bar{p}_h^{(2)}(t_2^{n+\frac{1}{2}}) \quad \text{and} \quad e_u^0 = 0 \quad (3.41) $$

where the last equation comes from (3.6b). Altogether, we get

$$ \|\bar{p}_h^n - \mu p_h^n\|_P \le \frac{\Delta t^2}{12} (5 + 5\gamma + 4\gamma^2 T) \|p_h\|_{C^3(0,T,P)} \quad (3.42) $$

$$ \|\bar{v}_h^{n+\frac{1}{2}} - \mu v_h^{n+\frac{1}{2}}\|_D \le \frac{\Delta t^2}{3} \|v_h\|_{C^2(0,T;D)} + (5 + 2\sqrt{2\gamma T}) \frac{\Delta t^2}{24} \|p_h\|_{C^3(0,T,P)} \quad (3.43) $$

$$ \|\bar{u}_h^{n+1} - u_h^{n+1}\|_P \le (5T\gamma + 4\gamma^2 T^2) \frac{\Delta t^2}{12} \|p_h\|_{C^3(0,T;P)} \quad (3.44) $$

□

Along with Hypothesis 0.3, this proves the uniform space/time convergence of the Störmer-Verlet integration scheme of the average unknowns and a post processed field.

## CONCLUSION AND PROSPECTS

This work proposes a proof of stability and convergence of the Störmer-Verlet integration scheme, towards the semi-discrete solution of wave equations in mixed formulation, uniform as the time step tends towards its largest admissible value. The error constants obtained depend neither on the spatial discretization parameters, nor on the distance of the time step from its largest admissible value, thus providing space/time convergence, if certain stability and convergence assumptions are satisfied by the spatial discretization. This makes it possible
to generalize results from the literature in which, for simplicity's sake, the stability condition is supposed to be
strictly satisfied. The assumptions of the theorem are consistent with existing uniform results for the leap-frop
scheme applied to the second-order formulation of wave equations. It appears that the natural variables that
converge uniformly in space/time are not the unknowns directly calculated by the scheme, but their consecutive
average between two time steps. However, it is possible to reconstruct a field that converges at each discretization
instant to a semi-discrete field. With an adequate change of variable or manipulation of the source terms, the
unknowns may be controlled directly. Note that this trick only works for linear wave equations with no external
coupling, while the present result could be extended to the case of dissipative equations, couplings between
two domains, addition of nonlinear terms, and so on. It could also be generalized in order to analyze other
integration schemes, such as implicit or dissipative schemes.

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, Stefan and Harrison, Reginald (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 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.

[Cohen and Fauqueux, 2011] Cohen, Gary. and Fauqueux, Sandrine. (2011). Mixed finite elements with mass-lumping for the transient wave equation. *Journal of Computational Acoustics*.

[Egger and Radu, 2020] Egger, H. and Radu, B. (2020). A mass-lumped mixed finite element method for acoustic wave propagation. *Numerische Mathematik*, 145(2):239–269.

[Eymard et al, 2000] Eymard, R., Gallouët, T., and Herbin, R. (2000). Finite volume methods. *Handbook of numerical analysis*, 7, 713-1018.

[Fezoui et al., 2005] Fezoui, L., Lanteri, S., Lohrengel, S., and Piperno, S. (2005). Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes. *ESAIM: Mathematical Modelling and Numerical Analysis*, 39(6):1149–1176.

[Geveci, 1988] Geveci, T. (1988). On the application of mixed finite element methods to the wave equations. *ESAIM: Mathematical Modelling and Numerical Analysis*, 22(2):243-250.

[Jenkins et al., 2002] Jenkins, E. W., Rivia're, B., and Wheeler, M. F. (2002). A Priori Error Estimates for Mixed Finite Element Approximations of the Acoustic Wave Equation. *SIAM Journal on Numerical Analysis*, 40(5):1698-1715.

[Joly, 2003] Joly, P. (2003). Variational methods for time-dependent wave propagation problems. *Topics in computational wave propagation: direct and inverse problems*, pages 201–264.

[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.

[Liu and Li, 2022] Liu, Z. and Li, X. (2022). Step-by-step solving schemes based on scalar auxiliary variable and invariant energy quadratization approaches for gradient flows. *Numerical Algorithms*, pages 1–22.

[Makridakis, 1992] Makridakis, C. G. (1992). On mixed finite element methods for linear elastodynamics. *Numerische Mathematik*, 61(1):235–260.

[Zuazua, 2005] Zuazua, E. (2005). Propagation, observation, and control of waves approximated by finite difference methods. *SIAM review*, 47(2), 197-243.