STABILITY AND SPACE/TIME CONVERGENCE OF ST ¨ORMER-VERLET TIME INTEGRATION OF THE MIXED FORMULATION OF LINEAR WAVE EQUATIONS

. 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¨ormer-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.


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 ≤  ≤  , ) The source terms and the initial conditions are supposed regular enough so that the following hypothesis holds, in three Hilbert spaces  ,  = { ∈  |ℬ ∈ } ⊂  and : Hypothesis 1.2 (Stability of the continuous system).The source terms and the initial conditions are regular enough such that there exists a constant  > 0 such that ‖‖  3 (0, ; ) + ‖‖  0 (0, ; ) + ‖‖  2 (0, ;) ≤ . (1.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 an interleaved scheme usually referred to as "centered explicit scheme" and which is formally equivalent to the Störmer-Verlet scheme, which is of primary importance for solving linear wave equations models because of its good mathematical properties, efficiency and ease of implementation [14].
We therefore suppose that the spatial discretisation is done with usual methods such as Finite Differences [18], Finite Elements [3], Finite Volumes [10], 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 ( ℎ ,  ℎ ) is sought in finite dimensional spaces  ℎ ⊂  ⊂  and  ℎ ⊂ : Find  ℎ ∈  ℎ and  ℎ ∈  ℎ such that ) where  ℎ :  ℎ →  ℎ is a discrete approximation of the operator ℬ, and  * ℎ :  ℎ →  ℎ is its adjoint, and  ℎ,0 ,  ℎ,0 ,  ℎ and  ℎ are discrete representations of  0 ,  0 ,  and  in  ℎ and  ℎ .We also denote  ℎ the identity operator of  ℎ .
Example 1.3.A possible choice of spatial discretisation for the acoustic wave equation given in example 1.1 is to partition the domain Ω into  disjoint intervals {  } 1≤≤ of length at most ℎ and to approximate  and  as  ℎ and  ℎ which are piecewise polynomial functions of degree  on each interval   , with  ℎ being continuous at the junction between intervals.This ensures that the approximation is conformal ( ℎ ⊂  ).In this example, the operator  ℎ is the restriction of ℬ to  ℎ .
(1.5) Time is discretised using a centered Störmer-Verlet algorithm (see [8], Chap.5, Eq. (5.12)) on staggered time grids with a time step ∆, which is formally equivalent to using the leap frog scheme on the second order system obtained after elimination of one equation, as developed further later.The unknown  ℎ is sought for on the grid { + 1 2 = ( + 1 2 )∆} 1≤≤ while the unknown  ℎ is sought for on the grid {  = ∆} 0≤≤ .Let us introduce the discrete operators  and  defined as They satisfy the useful following properties (1.7b) The Störmer Verlet scheme reads: Find For  ∈ [0,  ], where The motivations of the present paper are to provide proof of space/time convergence of the Störmer Verlet scheme in the mixed formulation of wave equations.In simple cases such as the one presented here, eliminating one unknown by applying the   (resp.) operator to (1.3c) (resp.(1.8c)) and using (1.3d) (resp.(1.8d)) to replace  ℎ yields a second order linear wave equation, if the source terms are regular enough: and is equivalent to applying the leap frog scheme [7,14]: So, for a field  ℎ such that  ℎ (0) =  ℎ,0 , uℎ =  ℎ , and its discrete counterpart {   } 0≤≤ such that  0 ℎ =  ℎ,0 , , the Störmer Verlet scheme on the mixed formulation amounts to a leap frog scheme on the second order wave equation üℎ Hence we expect similar estimates, at least on this remaining field, as those found in [5,6] for the leap frog scheme applied to the second order wave equation: If the CFL stability condition is respected, e.g.
where  stands for the spectral radius: then, there exists a constant  < 0 independent of ℎ and  such that Existing proofs of space/time convergence of the Störmer Verlet scheme on the mixed formulation of wave equations are found in [13], [11], [15], [1], [9].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  ≤ 1 2 or  < 1.Indeed, those proofs turn out to not be uniform as ∆ approaches its greatest admissible value, in other words the stability and convergence bounding constants depend on  and blow up as  → 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 [1] or dissipation [2].The paper is organized as follows.In Section 2, 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 (Props.3.12 and 3.15 of Sect.3) and convergence (Thm.4.1 of Sect.4) are proposed.

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 , and even blow up as  → 1.Let us manipulate the discrete equations in the following way.Let (, , ) ∈ (R + ) 3 .We take the scalar product of (1.8c)centered at time   with   ℎ in  ℎ .We then apply the operator  to (1.8d) at times  + 1 2 and  − 1 2 , and take the scalar product with   ℎ in  ℎ .We take the scalar product of (1.8d) centered at in  ℎ .Finally we take the scalar product of (1.8d) centered at  + 1 2 with ( + ) in  ℎ .We add up and get This is a discrete energy variation identity with Note that this identity is centered around  − 1 2 if and only if  ≡ 0, which leads to interpret  as an off-centering parameter.

Choice 𝛽 = 1 and 𝛼 = 𝛾 = 0
The choice  = 1 and  =  = 0 is done in [14] and proof of Lemma 5.3 of [9].It leads to a cross term in the energy by using (1.7a).

Special case
If  ℎ ≡ 0, we have that , therefore with ̃︀  ℎ =  ℎ − Δ 2 4  * ℎ  ℎ a modified identity operator.It is positive under the usual (CFL-condition).Then if  ≤ 1, If moreover  < 1 (strict CFL condition), we have A discrete summation from  = 1 to  =  yields, after telescopic elimination, which grants the stability of the energy, under the strict CFL condition.It then provides the stability of ‖ by using (2.6).The stability constants blow up as  → 1.
Remark 2.1.By introducing auxiliary functions on top of more usual elliptic projections, [9] write a system of equations on the error terms which are of the form of (1.8) with  ℎ = 0 (see Eq. (5.5)-(5.6)).They suppose  ≤ 1 2 in the second part of their hypothesis (A4).

General case
If  ℎ ̸ = 0, a Young's inequality with any  > 0 leads to Hence Recall the definition of  from (CFL-condition).The -CFL condition is equivalent to supposing that  < 1 (strict CFL condition), since it suffices to pose  = 1− 2 1+ 2 > 0 to satisfy (2.11).Therefore (2.12) We then have The first RHS term can not be bounded by the square root of the energy since ‖  .We cannot conclude.

Choice
The choice  = 1 2 ,  = 0 and  = 0 is done in [1].This leads to which is positive by definition, and satisfies The two first terms of the right-hand side of the previous equation have no defined sign but appear to behave as perturbations of the positive energy (2.14).This prompts to denoting = 0, it is immediately positive under the usual (CFL-condition) since . Otherwise, the analysis of the positivity (or more precisely the boundedness-by-below) of )  with the same techniques as above (Cauchy-Schwarz and Young's inequalities): with ̃︀  ℎ, =  ℎ − Δ 2 4 (1 + ) * ℎ  ℎ a modified identity operator which is positive under the -CFL condition (2.11), i.e. under the strict CFL condition  < 1.We get that (2.22) The modified energy satisfies (2.23) A discrete summation from 0 to  yields the stability of the energy.Under the strict CFL condition, it yields the stability of ‖  by using (2.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  → 1.

Choice 𝛼 = 𝛽 = 0 and 𝛾 = 1
The choice  = 0,  = 0 and  = 1 is done in the convergence proof of [13] and is remarkable since it breaks the symmetry of the manipulation.This leads to with We denote which has no definite sign.However, a Young's inequality with  > 0 yields (2.28) The largest positivity condition is obtained with  = 1 and reads ∆ ||| ℎ ||| ≤ 2. Under the usual (CFL-condition)  ≤ 1, we have Note that the modified energy is non-negative even in the equality case.It satisfies which is bounded by the square root of a the modified energy with constants that blow up as  → 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 ‖ ‖  and ‖  ℎ ‖  directly.However, the stability constants still blow up as  → 1.

Stability at the CFL
We provide here two strategies to obtain stability even as  → 1.This will not be possible directly on the unknowns of the system, but on post averaged or post processed quantities.

Case 𝛿𝑔 ℎ ≡ 0
We follow the ideas developed for the second order equation in [14].The leading principle is to use the second equation of (2.6) on two consecutive time steps to compensate the crossing term appearing in a new energy obtained on the system (1.8) averaged between two time steps.More precisely, Let us pose because 2( If satisfies Adding the telescopic sum from  = 1 to  =  we get which yields uniform stability estimates on the averaged fields: and ‖ However, trying to bound directly  ℎ only gets that which blows up as  → 1.

General case
In this section, we apply the methodology developed in [5,6] and adapt it to the mixed formulation and the Störmer-Verlet scheme.It is based on the manipulation described above, with  = 1 and  =  = 0 as in [14] and proof of Lemma 5.3 [9].A naive manipulation would be to take the scalar product of (1.8c) with   ℎ and of (1.8d) with  + 1 2 ℎ .One would get that Transferring the usual concepts for energy analysis directly here would imply to use (1.7a) and (1.8d) to replace the term ∆ 2 ( )  /4.However it induces a difficulty due to the presence of the source term  ℎ , as mentionned in Section 2.1.2.To mitigate this issue, let us proceed to a change of variables and suppose that  ℎ is regular enough uniformly with ℎ, which occurs for the classical discretization described in Example 1.3 as long as the source term is spatially distributed in the domain with a regularity related to the polynomial order : Hypothesis 3.1.In the sequel we will suppose that there exists a constant   > 0 such that for all ℎ > 0, The couple ({ ,ℎ also writes where m is the bilinear form defined on  ℎ such that for all (, ) Proof.Using (1.7a) we get where  is defined in (CFL-condition).
Corollary 3.5.Suppose (3.20) is satisfied.Then, m( To show the stability uniformly as ∆ approaches its greatest allowed value, let us exploit the spectral properties of the operator  * ℎ  ℎ as in [5]. We call ( ℎ, ,  ℎ, ) its eigenpairs which are chosen orthonormal in  .
Proof.In finite dimensional spaces, any symmetric real operator is diagonalizable in an orthonormal basis.
Following [5], we introduce the polynomial which is non-negative on the interval [0, 4].Then ̃︀ ℐ ℎ can be expressed as a polynomial of the operator ∆ 2  * ℎ  ℎ : Note that this polynomial appears naturally in the definition of the bilinear form m which refined the kinetic energy part of the discrete energy, hence the subscript "k" in the polynomial notation.For more intricate integration schemes, another polynomial   is introduced to treat the potential part of the discrete energy, see [5].Also note that the polynomial is used with ∆ 2  * ℎ  ℎ as argument to define ̃︀ ℐ ℎ , which relates the upper bound of the interval [0, 4] to the (CFL-condition).If ∆ reaches its greatest admissible value, the operator ̃︀ ℐ ℎ has a kernel, which prevents from using the kinetic part of the energy as a norm on the solution.This is why a uniform bound is sought for through a partitioning of the interval [0, 4], that will allow a uniform control on the different spectral components of the discrete solution.

Proposition 3.11 (Stability of the energy). The energy satisfies
Proof.Apply the Cauchy-Schwarz inequality on (3.12).Since  ℎ satisfies Hypothesis 3.1, we have Add this telescopic sum from  = 1 to  to get the first expected result: ,ℎ by using the expression (3.13): ℎ , ℎ , where Hence ‖ Since (CFL-condition) is satisfied, we have that We write that Since (CFL-condition) is satisfied, we also have that Proposition 3.12 (Stability of the averaged unknowns). ‖ Proof.Use jointly Prop.3.10 and 3.11 to get the two first equations.Moreover, the initial system was posed on the unknown   ℎ =   ℎ +   ℎ .Using the previous equations, one can only hope a direct control on Remark 3.13.This result generalizes the result recalled in Section 3.1 to the case of two source terms.It confirms that the uniform stability is only obtained on the averaged unknowns.
Remark 3.14.One can not hope for more control on the unknowns and   ℎ , unfortunately, at least without supposing more restrictive bounds on ∆.Only the Π  -projection of  + 1 2 ℎ is controlled directly.However, a global control is achieved on a post-processed field.Proposition 3.15 (Stability of the underlying field).Let {  ℎ } 0≤≤ +1 be the series of elements of  ℎ defined as Then,

.56)
Proof.From the orthogonality between the projection spaces, we get that Here we treat the two terms differently.For the first one we use (3.25a) with  ℎ = From the definition of  ℎ in (3.55) and using (3.11d), we get that  ℎ . Hence Both telescopic sums are canceled up from 0 to , and then added up.

Convergence at the CFL
Subtracting directly equations (1.3) to1 (1.8) yields truncation errors on each appearance of the operator  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  ℎ ∈  ℎ be defined as The second line can be integrated in time to get Along with the other equations, we get üℎ +  * ℎ  ℎ =  ℎ , (4.4a) On the fully discrete level, as in the previous section, let us define  ℎ ∈  ℎ as 0 ℎ = 0. (4.6b) Inserting this field in equation (1.8) yields The second line can be summed up in time to get ) Let ū ℎ =  ℎ (  ) and v ℎ =  ℎ (  ), where  ℎ is solution to the semi-discrete equation (1.3),  ℎ is the solution to the semi-discrete equation (4.1), therefore solutions to the semi-discrete system (4.4).where  depends on the final time  as a quadratic polynomial.
Proof.Let us subtract the two first lines of systems (4.4) and (4.8).
Let us define the error terms    = ū ℎ −   ℎ and    = v ℎ −   ℎ , which satisfy ) where the right hand side of the second line vanishes because of the initial condition  0 ℎ is defined as (1.8a), and a Taylor expansion shows that there exists   1 ∈ [ − 1 2 ,  + The result follows from the stability analysis performed in Section 3.2.Indeed, the two first lines of system (4.19) are of the form (1.8-3.55), with the substitutions where ℰ is defined as ).(4.23) Therefore, using (4.18), where there exists ). (4.27) Similarly, where there exists Next, let us show that ℰ is small.Using the naive formulation of the energy (3.12), we write that ,  ℎ are defined as (1.8a) and (1.8b), we directly have that  0  = 0.It remains to estimate  A Taylor expansion of pℎ yields that there exists  0 5 ∈ [0,  Along with Hypothesis 1.5, this proves the uniform space/time convergence of the Störmer-Verlet integration scheme of the average unknowns and a post processed field.Performing the same eliminations for the variable  ℎ yields a similar second order equation and leap-frog scheme but with differentiated sources in time.This distinction is of importance when using the mixed formulation model coupled with other models, in the context of multi-physics phenomena for instance.

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

2 𝑢 . The value of 𝑒 1 2 𝑝
can be found by subtracting both sides of (1.8b) to p