The Leray-G{\aa}rding method for finite difference schemes

Leray and G{\aa}rding have developed a multiplier technique for deriving a priori estimates for solutions to scalar hyperbolic equations in either the whole space or the torus. In particular, the arguments in Leray and G{\aa}rding's work provide with at least one local multiplier and one local energy functional that is controlled along the evolution. The existence of such a local multiplier is the starting point of the argument by Rauch for the derivation of semigroup estimates for hyperbolic initial boundary value problems. In this article, we explain how this multiplier technique can be adapted to the framework of finite difference approximations of transport equations. The technique applies to numerical schemes with arbitrarily many time levels, and encompasses a somehow magical trick that has been known for a long time for the leapfrog scheme. More importantly, the existence and properties of the local multiplier enable us to derive optimal semigroup estimates for fully discrete hyperbolic initial boundary value problems, which answers a problem raised by Trefethen, Kreiss and Wu.

Throughout this article, we use the notation We let M n (K) denote the set of n × N matrices with entries in K = R or C. If M ∈ M n (C), M * denotes the conjugate transpose of M .We let I denote the identity matrix or the identity operator when it acts on an infinite dimensional space.We use the same notation x * y for the Hermitian product of two vectors x, y ∈ C n and for the Euclidean product of two vectors x, y ∈ R n .The norm of a vector x ∈ C n is |x| := (x * x) 1/2 .The induced matrix norm on M n (C) is denoted • .
The letter C denotes a constant that may vary from line to line or within the same line.The dependence of the constants on the various parameters is made precise throughout the text.
In what follows, we let d ≥ 1 denote a fixed integer, which will stand for the dimension of the space domain we are considering.We shall also use the space 2 of square integrable sequences.Sequences may be valued in C k for some integer k.Some sequences will be indexed by Z d−1 while some will be indexed by Z d or a subset of Z d .We thus introduce some specific notation for the norms.Let ∆x i > 0 for i = 1, . . ., d be d space steps.We shall make use of the 2 (Z d−1 )-norm that we define as follows: for all v ∈ 2 (Z d−1 ), The corresponding scalar product is denoted •, • 2 (Z d−1 ) .Then for all integers m 1 ≤ m 2 , we set

Some motivations and a brief reminder
The ultimate goal of this article is to derive semigroup estimates for finite difference approximations of hyperbolic initial boundary value problems.Up to now, the only available general stability theory for such numerical schemes is due to Gustafsson, Kreiss and Sundström [GKS72].It relies on a Laplace transform with respect to the time variable, and the corresponding stability estimates are thereby restricted to zero initial data.A long standing problem in this line of research is, starting from the GKS stability estimates, which are resolvent type estimates, to incorporate nonzero initial data and to derive semigroup estimates, see, e.g., the discussion in [Tre84,section 4].This problem is delicate for the following reason: the validity of the GKS stability estimate is known to be equivalent to a slightly stronger version of the resolvent estimate sup z∈U where T is some bounded operator on 2 (N) that incorporates both the discretization of the hyperbolic equation and the numerical boundary conditions.Deriving an optimal semigroup estimate amounts to showing that T is power bounded.In finite dimension, the equivalence between power boundedness of T and the resolvent condition (1) is known as the Kreiss matrix Theorem, but the analogous equivalence is known to fail in general in infinite dimension.Worse, even the strong resolvent condition does not imply in general that T is power bounded, see, e.g., the review [SW97] or [TE05] for details and historical comments.Optimal semigroup estimates have nevertheless been derived for some discretized hyperbolic initial boundary value problems.More specifically, the first general derivation of semigroup estimates is due to Wu [Wu95], whose analysis deals with numerical schemes with two time levels and scalar equations.The results in [Wu95] were extended by Gloria and the author in [CG11] to systems in arbitrary space dimension, but the arguments in [CG11] are still restricted to numerical schemes with two time levels.The present article gives, as far as we are aware of, the first systematic derivation of semigroup estimates for fully discrete hyperbolic initial boundary value problems in the case of numerical schemes with arbitrarily many time levels.It generalizes the arguments of [Wu95,CG11] and provides new insight for the construction of "dissipative" numerical boundary conditions for discretized evolution equations.Let us observe that the leap-frog scheme, with some specific boundary conditions, has been dealt with by Thomas [Tho72] by using a multiplier technique.It is precisely this technique which we aim at developing in a systematic fashion for numerical schemes with arbitrarily many time levels.In particular, we shall explain why the somehow magical multiplier u n+2 j +u n j for the leap-frog scheme, see, e.g., [RM67], follows from a general theory that is the analogue of the Leray-Gårding method for partial differential equations, which we briefly recall now.
The method by Leray and Gårding [Ler53,Går56] provides with suitable multipliers for scalar hyperbolic operators of arbitrary order.Namely, given an integer m ≥ 0, we consider a partial differential operator of the form where t ∈ R stands for the time variable, x ∈ R d stands for the space variable1 , and each operator P k (∂ x ) is a linear combination of spatial partial derivatives of order k: In the above formula, the p k,α 's are real numbers2 .Well-posedness of the Cauchy problem in Sobolev spaces is known to be linked with hyperbolicity of L. Namely, if L is strictly hyperbolic, meaning that for all ξ ∈ R d \ {0}, the (homogeneous) polynomial has m + 1 simple purely imaginary roots with respect to τ , then the Cauchy problem (2) is well-posed in In particular, there exists a constant C > 0, that is independent of the solution u and the initial data u 0 , u 1 , . . ., u m , such that there holds: The method by Leray and Gårding gives a quick and elegant way to derive the estimate (4) assuming that the solution u to (2) is sufficiently smooth.By standard duality arguments, the validity of the a priori estimate (4) yields well-posedness -meaning existence, uniqueness and continuous dependence on the data-for (2).Hence the main point is to prove (4) assuming that u is sufficiently smooth and decaying at infinity so that all integration by parts arising in the computations are legitimate.The main idea is to find a suitable quantity M u, which we call a multiplier and that will be linear with respect to u, such that when integrating the quantity 0 = (M u) (L u) on the slab [0, T ] × R d , one gets the estimate (4) for free (negative times are obtained by changing t → −t).Following [Ler53, Chapter VI] and [Går56, Section 3], one possible choice of a multiplier is given by L u where L stands for the partial differential operator of order m whose symbol is ∂ τ P , with P given in (3).Why L u is a good multiplier is justified in [Ler53,Går56].A well-known particular case is the choice of 2 ∂ t u as a multiplier for the wave equation.Here P (τ, ξ) = τ 2 + |ξ| 2 and therefore ∂ τ P = 2 τ , hence the choice 2 ∂ t u.The latter quantity is indeed a suitable multiplier for the wave operator because of the formula3 : The important fact here is that the energy: is a positive definite quadratic form of the first order partial derivatives of u.Let us observe that the multiplier L u is local, meaning that its pointwise value at (t, x) only depends on u in a neighborhood of (t, x).This is important in view of using this multiplier in the study of initial boundary value problems.
Another important remark is that the above energy is also local, and the arguments in [Ler53,Går56] show that this property is not specific to the wave operator.The fact that both the multiplier and the energy are local is crucial in the arguments of [Rau72, Lemma 1].In our framework of discretized equations, the multiplier will be local but the energy will not necessarily be so.We shall not exactly follow the arguments of [Rau72] which use time reversibility, but rather construct dissipative boundary conditions which will yield the optimal semigroup estimate we are aiming at.

The main result
We first set a few notations.We let ∆x 1 , . . ., ∆x d , ∆t > 0 denote space and time steps where the ratios, the so-called Courant-Friedrichs-Lewy parameters, λ i := ∆t/∆x i , i = 1, . . ., d, are fixed positive constants.We keep ∆t ∈ (0, 1] as a small parameter and let the space steps ∆x 1 , . . ., ∆x d vary accordingly.The 2 -norms with respect to the space variables have been previously defined and thus depend on ∆t and the CFL parameters through the mesh volume (∆x We always identify a sequence w indexed by either N (for time), Z d−1 or Z d (for space), with the corresponding step function.In particular, we shall feel free to take Fourier or Laplace transforms of such sequences.For all j ∈ Z d , we set j = (j 1 , j ) with j := (j 2 , . . ., j d ) ∈ Z d−1 .We let p, q, r ∈ N d denote some fixed multi-integers, and define p 1 , q 1 , r 1 , p , q , r according to the above notation.We also let s ∈ N denote some fixed integer.We consider a recurrence relation of the form: where the operators Q σ and B j 1 ,σ are given by: In ( 6), the a ,σ , b ,j 1 ,σ are real numbers and are independent of the small parameter ∆t (they may depend on the CFL parameters though), while S denotes the shift operator on the space grid: (S v) j := v j+ for j, ∈ Z d .We have also used the short notation The numerical scheme (5) is understood as follows: one starts with 2 initial data (f 0 j ), ..., (f s j ) defined for j 1 ≥ 1 − r 1 .Assuming that the solution has been defined up to some time index n + s, n ≥ 0, then the first and second equations in (5) should uniquely determine u n+s+1 j for j 1 ≥ 1 − r 1 , j ∈ Z d−1 .The meshes associated with j 1 ≥ 1 correspond to the interior domain while those associated with j 1 = 1 − r 1 , . . ., 0 represent the discrete boundary.We wish to deal here simultaneously with explicit and implicit schemes and therefore make the following solvability assumption.
The first and second equations in (5) therefore uniquely determine u n+s+1 j for j 1 ≥ 1 − r 1 , and one then proceeds to the following time index n + s + 2. Existence and uniqueness of a solution (u n j ) to (5) follows from Assumption 1, so the last requirement for well-posedness is continuous dependence of the solution on the three possible source terms (F n j ), (g n j ), (f n j ).This is a stability problem for which several definitions can be chosen according to the functional framework.The following one dates back to [GKS72] in one space dimension and was also considered by Michelson [Mic83] in several space dimensions.It is specifically relevant when the boundary conditions are non-homogeneous ((g n j ) ≡ 0): Definition 1 (Strong stability).The finite difference approximation (5) is said to be "strongly stable" if there exists a constant C such that for all γ > 0 and all ∆t ∈ (0, 1], the solution (u n j ) to (5) with The main contributions in [GKS72,Mic83] are to show that strong stability can be characterized by a certain algebraic condition, which is usually referred to as the Uniform Kreiss-Lopatinskii Condition, see [Cou13] for an overview of such results.We do not pursue such arguments here but rather assume from the start that (5) is strongly stable.We can thus control, with zero initial data, 2 type norms of the solution to (5).Our goal is to understand which kind of stability estimate holds for the solution to (5) when one now considers nonzero initial data (f 0 j ), . . ., (f s j ) in 2 .Our main assumption is the following.
Assumption 2 (Stability for the discrete Cauchy problem).For all ξ ∈ R d , the dispersion relation has s + 1 simple roots in D. (The von Neumann condition is said to hold when the roots are located in D.) In (8), we have used the classical notation From Assumption 1, we know that Q s+1 is an isomorphism on 2 , which implies by Fourier analysis that Q s+1 (e i ξ 1 , . . ., e i ξ d ) does not vanish for any ξ ∈ R d .In particular, the dispersion relation (8) is a polynomial equation of degree s + 1 in z for any ξ ∈ R d .We now make the following assumption, which already appeared in [GKS72,Mic83] and several other works on the same topic.
Then a −r 1 and a p 1 do not vanish on U × R d−1 , and they have nonzero degree with respect to z for all η ∈ R d−1 .
Our main result is comparable with [Wu95, Theorem 3.3] and [CG11, Theorems 2.4 and 3.5] and shows that strong stability (or "GKS stability") is a sufficient condition for incorporating 2 initial conditions in (5) and proving optimal semigroup estimates.The main price to pay in Assumption 2 is that the roots of the dispersion relation (8), which are nothing but the eigenvalues of the so-called amplification matrix for the Cauchy problem, need to be simple.This property is satisfied for instance by the leap-frog and modified leap-frog schemes in several space dimensions, under an appropriate CFL condition, see Paragraph 1.3.Our main result reads as follows.
Theorem 1.Let Assumptions 1, 2 and 3 be satisfied, and assume that the scheme (5) is strongly stable in the sense of Definition 1. Then there exists a constant C such that for all γ > 0 and all ∆t ∈ (0, 1], the solution to (5) satisfies the estimate: In particular, the scheme (5) is "semigroup stable" in the sense that there exists a constant C such that for all ∆t ∈ (0, 1], the solution (u n j ) to (5) with (F n j ) = (g n j ) = 0 satisfies the estimate The scheme (5) is also 2 -stable with respect to boundary data, see [Tre84,Definition 4.5], in the sense that there exists a constant C such that for all ∆t ∈ (0, 1], the solution Theorem 1 gives the optimal semigroup estimate (11), and is therefore an improvement with respect to our earlier work [Cou14] where in one space dimension, and under an appropriate non-glancing condition4 , we were able to derive the estimate (here r 1 = r, p 1 = p since d = 1): The latter estimate does not incorporate on the left hand side the quantity: and was unfortunately still not sufficient for deriving the semigroup estimate (11).Our main contribution in this article is to exhibit a suitable multiplier for the multistep recurrence relation in (5).With this multiplier, we can readily show that, for zero initial data, the (discrete) derivative of an energy can be controlled, as in [Rau72], by the trace estimate of (u n j ) and this is where strong stability comes into play.
This first argument gives Theorem 1 for zero initial data (and even for nonzero initial data if the nonglancing condition of [Cou14] is satisfied).By linearity we can then reduce to the case of zero forcing terms in the interior and on the boundary.The next arguments in [Rau72] use time reversibility, which basically always fails for numerical schemes5 .Hence we must find another argument for dealing with nonzero initial data.Hopefully, the properties of our multiplier enable us to construct an auxiliary problem, where we modify the boundary conditions of (5), and for which we can prove optimal semigroup and trace estimates by "hand-made" calculations.In other words, we exhibit an alternative set of boundary conditions that yields strict dissipativity.Using these auxiliary numerical boundary conditions, the proof of Theorem 1 follows from a standard superposition argument, see, e.g., [BGS07, Section 4.5] for partial differential equations or [Wu95,CG11] for numerical schemes.
Remark 1. Assumption 3 excludes the case of explicit two level schemes for which s = 0 and Q 1 = I, for in that case a −r 1 and/or a p 1 do not depend on z.However, this case has already been dealt with in [Wu95, CG11], and we shall see in Section 3 where the assumption that a −r 1 and a p 1 are not constant is involved, and why the proof is actually simpler in the case s = 0 and Q 1 = I.

One space dimension
Our goal is to approximate the outgoing transport equation (d = 1 here): with t, x > 0 and a < 0. The latter transport equation does not require any boundary condition at x = 0.However, discretizing (12) usually requires prescribing numerical boundary conditions, unless one considers an upwind type scheme with a space stencil "on the right" (meaning r 1 = 0 in (5)).We now detail two possible multistep schemes for discretizing (12).Both are obtained by the so-called method of lines, which amounts to first discretizing the space derivative ∂ x u and then choosing an integration technique for discretizing the time evolution, see [GKO95].
The leap-frog scheme.It is obtained by approximating the space derivative ∂ x u by the centered difference (u j+1 − u j−1 )/(2 ∆x), and by then applying the so-called Nyström method of order 2, see [HNW93, Chapter III.1].The resulting approximation reads Recall that λ > 0 denotes the fixed ratio ∆t/∆x.Even though (12) does not require any boundary condition at x = 0, the leap-frog scheme stencil includes one point to the left, and we therefore need to prescribe some numerical boundary condition at j = 0.One possibility6 is to prescribe the homogeneous or inhomogeneous Dirichlet boundary condition.With general source terms, the corresponding scheme reads Assumption 1 is trivially satisfied because (13) is explicit.The leap-frog scheme satisfies Assumption 2 provided that λ |a| < 1.In that case, the two roots to the dispersion relation are simple and have modulus 1 for all ξ ∈ R. Assumption 3 is satisfied as long as the velocity a is nonzero, for in that case a 1 (z) = −a −1 (z) = λ a z.The scheme ( 13) is known to be strongly stable, see [GT81].
In particular, Theorem 1 shows that (13) is semigroup stable.An illustration of this stability property is given in the numerical simulation of a bump function, propagating at speed a = −1 towards the left.Homogeneous Dirichlet boundary conditions are enforced at j = 0.The reflection of the bump generates a highly oscillatory wave packet that propagates with velocity +1 towards the right.The envelope of this wave packet coincides with the profile of the initial condition, which indicates that the 2 -norm is roughly preserved by the evolution.This numerical observation is in agreement with semigroup boundedness.Other choices of numerical boundary conditions for the leap-frog scheme or its fourth order extension are discussed, e.g., in [Oli74,Slo83,Tho72,Tre84]. The main discussion in [Oli74,Slo83,Tre84] is to verify strong stability for a wide choice of numerical boundary conditions, and if strong stability holds, then Theorem 1 automatically gives semigroup boundedness, which was not achieved in these earlier works.A scheme based on the backwards differentiation rule.We still start from the transport equation (12), approximate the space derivative ∂ x u by the centered finite difference (u j+1 − u j−1 )/(2 ∆x), and then apply the backwards differentiation formula of order 2, see [HNW93, Chapter III.1].The resulting scheme reads: u n j = 0 .This corresponds to s = 1 and The operator Q 2 is an isomorphism on 2 (Z) since Q 2 is an isomorphism for any small λ a (as a perturbation of 3/2 I), Q 2 depends continuously on λ a, and there holds (uniformly with respect to λ a): The operator Q 2 is therefore an isomorphism on 2 (Z) for any λ a > 0 (see, e.g., [Cou09, Lemma 4.3]).
Let us now study the dispersion relation (8), which reads here It is clear that the latter equation has two simple roots in z for any ξ ∈ R.Moreover, if sin ξ = 0, the roots are 1 and 1/3 which belong to D. In the case sin ξ = 0, none of the roots belongs to S 1 and examining the case λ a sin ξ = 1, we find that for sin ξ = 0, both roots belong to D (which is consistent with the shape of the stability region for the backwards differentiation formula of order 2, see [HW96, Chapter V.1]).Assumption 2 is therefore satisfied.Assumption 3 is satisfied as long as a is nonzero since there holds p = r = 1 and a 1 (z) = a −1 (z) = λ a z 2 /2.Theorem 1 therefore yields semigroup boundedness as long as one uses numerical boundary conditions for which the numerical scheme is well-defined (this is at least the case for λ a small enough) and strong stability holds.

Two space dimensions
Here we wish to approximate the two-dimensional transport equation (d = 2): in the space domain {x 1 > 0 , x 2 ∈ R}.When a 1 is negative, the latter problem does not necessitate any boundary condition at x 1 = 0. Following [AG76], we use one of the following two-dimensional versions of the leap-frog scheme, either or so Assumption 3 is valid as long as a 1 = 0.For the scheme (15), we have again r 1 = p 1 = 1, and so Assumption 3 is valid as long as both a 1 and a 2 are nonzero.We refer to [AG79] for the verification of strong stability depending on the choice of some numerical boundary conditions for ( 14) or (15).Once again, if strong stability holds, then Theorem 1 yields semigroup boundedness and 2 -stability with respect to boundary data.

The Leray-Gårding method for fully discrete Cauchy problems
This section is devoted to proving stability estimates for discretized Cauchy problems, which is the first step before considering the discretized initial boundary value problem (5).More precisely, we consider the simpler case of the whole space j ∈ Z d , and the recurrence relation: where the operators Q σ are given by (6).We recall that in (6), the a ,σ are real numbers and are independent of the small parameter ∆t (they may depend on the CFL parameters λ 1 , . . ., λ d ), while S denotes the shift operator on the space grid: (S v) j := v j+ for j, ∈ Z d .Stability of ( 16) is defined as follows.
Definition 2 (Stability for the discrete Cauchy problem).The numerical scheme defined by ( 16) is ( 2 -) stable if Q s+1 is an isomorphism from 2 (Z d ) onto itself, and if furthermore there exists a constant C 0 > 0 such that for all ∆t ∈ (0, 1], for all initial conditions (f 0 j ) j∈Z d , . . ., (f s j ) j∈Z d in 2 (Z d ), there holds Let us quickly recall that stability in the sense of Definition 2 is in fact independent of ∆t ∈ (0, 1] (because (16) does not involve ∆t and (17) can be simplified on either side by i ∆x i ), and can be characterized in terms of the uniform power boundedness of the so-called amplification matrix where the Q σ (κ)'s are defined in (8) and where it is understood that A is defined on the largest open set of C d on which Q s+1 does not vanish.Let us also recall that if Q s+1 is an isomorphism from 2 (Z d ) onto itself, then Q s+1 does not vanish on (S 1 ) d , and therefore does not vanish on an open neighborhood of (S 1 ) d .With the above definition (18) for A , the following well-known result holds: Proposition 1 (Characterization of stability for the fully discrete Cauchy problem).Assume that Q s+1 is an isomorphism from 2 (Z d ) onto itself.Then the scheme (16) is stable in the sense of Definition 2 if and only if there exists a constant C 1 > 0 such that the amplification matrix A in (18) satisfies In particular, the spectral radius of A (e i ξ 1 , . . ., e i ξ d ) should not be larger than 1 (the so-called von Neumann condition).
The eigenvalues of A (e i ξ 1 , . . ., e i ξ d ) are the roots to the dispersion relation (8).When these roots are simple for all ξ ∈ R d , the von Neumann condition is both necessary and sufficient for stability of (16), see, e. g., [Cou13, Proposition 3].Assumption 2 is therefore a way to assume that ( 16) is stable for the discrete Cauchy problem.Our goal is to derive the semigroup estimate (17) not by applying Fourier transform to (16) and using uniform power boundedness of A , but rather by multiplying the first equation in ( 16) by a suitable local multiplier.The analysis relies first on the simpler case where one only considers the time evolution and no additional space variable.

Stable recurrence relations
In this Paragraph, we consider sequences (v n ) n∈N with values in C. The index n should be thought of as the discrete time variable, and we therefore introduce the new notation T for the shift operator on the time grid: (T m v) n := v n+m for all m, n ∈ N. We start with the following elementary but crucial Lemma, which is the analogue of [Går56, Lemme 1.1].
Lemma 1 (The energy-dissipation balance law).Let P ∈ C[X] be a polynomial of degree s + 1 whose roots are simple and located in D. Then there exists a positive definite Hermitian form q e on C s+1 , and a nonnegative Hermitian form q d on C s+1 , that both depend in a C ∞ way on P , such that for any sequence (v n ) n∈N with values in C, there holds In particular, for all sequence (v n ) n∈N that satisfies the recurrence relation The fact that there exists a Hermitian norm on C s+1 that is nonincreasing along solutions to the recurrence relation is not new.In fact, it is easily seen to be a consequence of the Kreiss matrix Theorem, see [SW97].However, the important point here is that we can construct a multiplier that yields directly the "energy boundedness" (or decay).The fact that the coefficients of this multiplier are integer multiples of the coefficients of P will be crucial in the analysis of Section 3, see also Proposition 2 below.
Proof.We borrow some ideas from [Går56, Lemme 1.1] and introduce the interpolation polynomials: where x 1 , . . ., x s+1 denote the roots of P , and a = 0 its dominant coefficient.Since the roots of P are pairwise distinct, the P k 's form a basis of C s [X] and they depend in a C ∞ way on the coefficients of P .We have We then consider a sequence (v n ) n∈N with values in C and compute

The conclusion follows by defining:
∀ (w 0 , . . ., w s ) ∈ C s+1 , q e (w 0 , . . ., w s ) := q d (w 0 , . . ., w s ) := The form q e is positive definite because the P k 's form a basis of C s [X].The form q d is nonnegative because the roots of P are located in D. Both forms depend in a C ∞ way on the coefficients of P because the roots of P are simple.
Lemma 1 shows that the polynomial P yields the good multiplier T P (T) v n for the recurrence relation P (T) v n = 0. Of course, P is not the only possible choice, though it will be our favorite one in what follows.As in [Går56, Lemme 1.1], any polynomial of the form7 provides with an energy balance of the form with suitable Hermitian forms q e , q d that have the same properties as stated in Lemma 1.

The energy-dissipation balance for finite difference schemes
In this Paragraph, we consider the numerical scheme (16).We introduce the following notation: Thanks to Fourier analysis, Lemma 1 easily gives the following result: Proposition 2 (The energy-dissipation balance law).Let Assumptions 1 and 2 be satisfied.Then there exist a continuous coercive quadratic form E 0 and a continuous nonnegative quadratic form D 0 on 2 (Z d ; R) s+1 such that for all sequences (v n ) n∈N with values in 2 (Z d ; R) and for all n ∈ N, there holds In particular, for all initial data f 0 , . . ., f s ∈ 2 (Z d ; R), the solution to (16) satisfies and ( 16) is ( 2 -)stable.
Proof.We use the same notation v n for the sequence (v n j ) j∈Z d and the corresponding step function on R d whose value on the cell [j 1 ∆x 1 , (j where v n denotes the Fourier transform of v n , and where we have let where q e,ζ , q d,ζ depend in a C ∞ way on ζ ∈ R d and are 2 π-periodic in each ζ j .Furthermore, q e,ζ is positive definite and q d,ζ is nonnegative.The conclusion of Proposition 2 follows by a standard compactness argument for showing coercivity of E 0 .

Examples
The first basic example corresponds to the case s = 0 for which the multiplier provided by Proposition 2 is Q 1 v n+1 j .In that case, the energy E 0 reads |||Q 1 v||| 2 −∞,+∞ (recall that Q 1 is an isomorphism) and the energy-dissipation balance law is nothing but the trivial identity The second line of this algebraic identity can be rewritten as and 2 -stability for the Cauchy problem amounts to assuming that the operator norm of Let us now consider the leap-frog scheme in one space dimension, for which we have s = 1 and The corresponding dispersion relation (8) reduces to For λ |a| < 1, the latter equation has two simple roots x 1 (ξ), x 2 (ξ) of modulus 1.Following the previous analysis, see ( 20)-( 21), the form q e,ζ is given by and q d,ζ is zero.The associated forms in Proposition 2 are D 0 ≡ 0 and (recall here d = 1): The latter energy functional E 0 is coercive under the condition λ |a| < 1, which is the necessary and sufficient condition of stability for the leap-frog scheme, and E 0 is conserved for solutions to the leap-frog scheme.The conservation of E 0 is usually proved by starting from the recurrence relation using the multiplier u n+2 j + u n j , and summing with respect to j.This is equivalent, for solutions to the leap-frog scheme, to what we propose here, since our multiplier reads However, it will appear more clearly in Section 3 why our choice for M u n j has a major advantage when considering initial boundary value problems.
Let us observe here that the energy functional E 0 is associated with a local energy density . This is very specific to the leap-frog scheme.In general, the coefficients of the Hermitian forms q e,ζ , q d,ζ are not trigonometric polynomials of ζ and therefore E 0 , D 0 do not necessarily admit local densities.This is one main difference with [Ler53,Går56].

Semigroup estimates for fully discrete initial boundary value problems
We now turn to the proof of Theorem 1 for which we shall use the results of Section 2 as a toolbox.By linearity of (5), it is sufficient to prove Theorem 1 separately in the case (f 0 j ) = • • • = (f s j ) = 0, and in the case (F n j ) = 0, (g n j ) = 0.The latter case is the most difficult and requires the introduction of an auxiliary set of "dissipative" boundary conditions.Solutions to (5) are always assumed to be real valued, which means that the data are real valued.For complex valued initial data and/or forcing terms, one just uses the linearity of (5).

The case with zero initial data
We first assume (f 0 j ) = • • • = (f s j ) = 0.By strong stability, we already know that (7) holds with a constant C that is independent of γ > 0 and ∆t ∈ (0, 1].Therefore, proving Theorem 1 amounts to showing the existence of a constant C, that is independent of γ > 0 and ∆t ∈ (0, 1] such that the solution to (5) with We thus consider a parameter γ > 0 and a time step ∆t ∈ (0, 1], and focus on the numerical scheme (5) with zero initial data (that is, (f 0 j ) = • • • = (f s j ) = 0).For all n ∈ N, we extend the sequence (u n j ) by zero for j 1 ≤ −r 1 : We use Proposition 2 and compute: Due to the form of the operator L, see ( 22), and the fact that v n j vanishes for j 1 ≤ −r 1 , there holds: and we thus get We multiply the latter equality by exp(−2 γ (n + s + 1) ∆t), sum with respect to n from 0 to some N and use the fact that D 0 is nonnegative.Recalling that the initial data in (5) vanish, we get with and Let us now estimate the two source terms S 1,N , S 2,N in (24).We begin with the term S 2,N defined in (26).Let us recall that the ratio ∆t/∆x 1 is fixed.Furthermore, the form of the operators L and M in (22) gives the estimate (recall that v n j vanishes for j 1 ≤ −r 1 ): for a constant C that does not depend on N , γ nor on ∆t.We thus have, uniformly with respect to N ∈ N, γ > 0 and ∆t ∈ (0, 1]: where we have used the trace estimate (7) that follows from the strong stability assumption.
Let us now focus on the term S 1,N in (24), see the defining equation ( 25).We use the Cauchy-Schwarz inequality and derive (using now the interior estimate in (7) that follows from the strong stability assumption): Ignoring the nonnegative term on the left hand-side of (24) and using the coercivity of E 0 , we have proved that there exists a constant C > 0 that is uniform with respect to N, γ, ∆t such that: which yields (23) and therefore the validity of Theorem 1 in the case of zero initial data.

Construction of dissipative boundary conditions
In this paragraph, we consider an auxiliary problem for which we shall be able to prove simultaneously an optimal semigroup estimate and a trace estimate for the solution.More precisely, we shall prove the following result.
Theorem 2. Let Assumptions 1, 2 and 3 be satisfied.Then for all P 1 ∈ N, there exists a constant C P 1 > 0 such that, for all initial data (f 0 j ), . . ., (f s j ) ∈ 2 (Z d ) and for all source term (g n j ) j 1 ≤0,n≥s+1 that satisfies there exists a unique sequence (u n j ) j∈Z d ,n∈N solution to Moreover for all γ > 0 and ∆t ∈ (0, 1], this solution satisfies Theorem 2 justifies why we advocate the choice M u n j = 2 u n+2 j + λ a (u n+1 j+1 − u n+1 j−1 ) rather than the more standard u n+2 j + u n j as a multiplier for the leap-frog scheme.Despite repeated efforts, we have not been able to prove the estimate (28) when using the numerical boundary condition u n+2 j + u n j on j 1 ≤ 0, in conjunction with the leap-frog scheme on j 1 ≥ 1.
Proof.Let us first quickly observe that the solution to ( 27) is well-defined since, as long as we have determined the solution up to a time index n + s, n ≥ 0, then u n+s+1 is sought as a solution to an equation of the form where F belongs to 2 (Z d ) (this is due to the form of L and M , see ( 22)).Hence u n is uniquely defined and belongs to 2 (Z d ) for all n ∈ N.
The proof of Theorem 2 starts again with the application of Proposition 2. Using the nonnegativity of the dissipation form D 0 , we get8 By the Young inequality, we get We multiply the latter inequality by exp(−2 γ (n + s + 1) ∆t), sum from n = 0 to some arbitrary N and already derive the estimate (here we use again the fact that ∆t/∆x 1 is a fixed positive constant): Using the coercivity of E 0 and the inequality we have therefore derived the estimate where the constant C is independent of γ, ∆t and on the solution (u n j ).In order to prove (28), the main remaining task is to derive the trace estimate for (u n j ).This is done by first dealing with the case where γ ∆t is large.
• From the definition of the operator L, see ( 22), there exists a constant C > 0 and an integer J such that Since Q s+1 is an isomorphism, there exists a constant c > 0 such that Multiplying by exp(−2 γ (n + s + 1) ∆t) and summing with respect to n ∈ N, we get n≥s+1 ∆t e −2 γ n ∆t Choosing γ ∆t large enough, that is γ ∆t ≥ ln R 0 for some numerical constant R 0 > 1 that depends only on the (fixed) coefficients of the operator L, we have derived the estimate .
It remains to use (29) and we get an even better estimate than (28) which we were originally aiming at: This gives a control of infinitely many traces and not only finitely many (this restriction to finitely many traces will appear in the regime where γ ∆t can be small).
• From now on, we have fixed a constant R 0 > 1 such that (28) holds for γ ∆t ≥ ln R 0 and we thus assume γ ∆t ∈ (0, ln R 0 ].We also know that the estimate (29) holds, independently of the value of γ ∆t, and we now wish to estimate the traces of the solution (u n j ) for finitely many values of j 1 .We first observe from (29) that for all γ > 0 and ∆t ∈ (0, 1], there exists a constant C γ,∆t such that In particular, for any j 1 ∈ Z, the Laplace-Fourier transforms u j 1 of the step functions The dual variables are denoted τ = γ + i θ, γ > 0, and η = (η 2 , . . ., η d ) ∈ R d−1 .It will also be convenient to introduce the notation η ∆ := (η 2 ∆x 2 , . . ., η d ∆x d ).
We now recall that γ ∆t is restricted to the interval (0, ln R 0 ], and we use (29) to derive Similarly, we have which we can again uniformly estimate by the right hand side of (30).
Going back to the right hand side terms in (31) and (32), we find that there only remains for proving (30) to estimate the integral (here there are finitely many values of σ and σ ): where we have applied Fubini Theorem.Applying first Plancherel Theorem with respect to the d − 1 last space variables, we get The conclusion then follows by computing and by recalling that γ ∆t belongs to (0, ln R 0 ].We can eventually bound the integrals on the left hand side of (30) by estimating separately the integrals of each term on the right hand side of (31) and (32).
The conclusion now relies on the following crucial result.
Lemma 3 (The trace estimate).Let Assumptions 1, 2 and 3 be satisfied.Let R 0 > 1 be fixed as above and let P 1 ∈ N. Then there exists a constant C P 1 > 0 such that for all z ∈ U with |z| ≤ R 0 , for all η ∈ R d−1 and for all sequence (w j 1 ) j 1 ∈Z ∈ 2 (Z; C), there holds Recall that the functions a 1 , 1 = −r 1 , . . ., p 1 , are defined in (9).
The proof of Lemma 3 is rather long.Before giving it in full details, we indicate how Lemma 3 yields the result of Theorem 2. We apply Lemma 3 to z = exp(τ ∆t), τ = γ + i θ with γ ∆t ∈ (0, ln R 0 ], η ∆ ∈ R d−1 and the sequence ( u j 1 (τ, η)) j 1 ∈Z .We then integrate (33) with respect to (θ, η) and use Lemma 2 to derive It remains to apply Plancherel Theorem and we get Recalling that γ ∆t is restricted to the interval (0, ln R 0 ], we have thus derived the trace estimate n∈N ∆t e −2 γ n ∆t Combined with the semigroup and interior estimate (29), this gives the estimate (28) of Theorem 2 for γ ∆t ∈ (0, ln R 0 ]. Proof of Lemma 3. Let us recall that the functions a 1 are 2 π-periodic with respect to each coordinate of η.We can therefore restrict to η ∈ [0, 2 π] d−1 rather than considering η ∈ R d−1 .We argue by contradiction and assume that the conclusion to Lemma 3 does not hold.This means the following, up to normalizing and extracting subsequences; there exist three sequences (indexed by k ∈ N): • a sequence (w k ) k∈N with values in 2 (Z; C) such that (w k −r 1 −p 1 , . . ., w k P 1 ) belongs to the unit sphere of C P 1 +r 1 +p 1 +1 for all k, and (w k −r 1 −p 1 , . . ., w k P 1 ) converges towards (w −r 1 −p 1 , . . ., w P 1 ) as k tends to infinity, and these sequences satisfy: We are going to show that (34) implies that (w −r 1 −p 1 , . . ., w P 1 ) must be zero, which will yield a contradiction since this vector must have norm 1.
• Let us first show that each component (w k j 1 ) k∈N , j 1 ∈ Z, has a limit as k tends to infinity.This is already clear for j 1 = −r 1 − p 1 , . . ., P 1 .For j 1 > P 1 , we argue by induction.From (34), we have and by Assumption 3, we know that a p 1 (z, η) is nonzero.Hence (w k P 1 +1 ) k∈N converges towards which we define as w P 1 +1 .We can argue by induction in the same way for all indices j 1 > P 1 + 1, but also for indices j 1 < −r 1 − p 1 because the function a −r 1 also does not vanish on U × R d−1 .
Using (34), we have thus shown that for each j 1 ∈ Z, (w k j 1 ) k∈N tends towards some limit w j 1 as k tends to infinity, and the sequence w, which does not necessarily belong to 2 (Z; C), satisfies the induction relations: ∀ j 1 ≤ 0 , • The induction relation (35) is the one that arises in [GKS72,Mic83] and all the works that deal with strong stability.The main novelty here is to use simultaneously (35) for controlling the unstable components of (w −r 1 −p 1 , . . ., w −1 ) and (36) for controlling the stable components of (w −r 1 −p 1 , . . ., w −1 ).The fact that w satisfies simultaneously (35) and (36) for j 1 ≤ 0 automatically annihilates the central components.This sketch of proof is made precise below.
We define the source terms: which, according to (34), satisfy lim We also introduce the vectors (here T denotes transposition) and the matrices: The matrix L is well-defined on U × R d−1 according to Assumption 3. The matrix M is also well-defined on U × R d−1 because for any η ∈ R d−1 , Assumption 3 asserts that a p 1 (•, η) is a nonconstant polynomial whose roots lie in D. From the Gauss-Lucas Theorem, the roots of ∂ z a p 1 (•, η) lie in the convex hull of those of a p 1 (•, η).Therefore ∂ z a p 1 (•, η) does not vanish on U .In the same way, ∂ z a −r 1 (•, η) does not vanish on U .With our above notation, the vectors W k j 1 , W j 1 , satisfy the one step induction relations: • From Assumption 3 and the above application of the Gauss-Lucas Theorem, we already know that both matrices L(z, η) and M(z, η) are invertible for (z, η) ∈ U × R d−1 .Furthermore, Assumption 2 shows that L(z, η) has no eigenvalue on S 1 for (z, η) ∈ U × R d−1 .This property dates back at least to [Kre68].However, central eigenvalues on S 1 may occur for L when z belongs to S 1 .The crucial point for proving Lemma 3 is that Assumption 2 precludes central eigenvalues of M for all z ∈ U .Namely, for all z ∈ U and all η ∈ R d−1 , M(z, η) has no eigenvalue on S 1 .This property holds because otherwise, for some (z, η) ∈ U × R d−1 , there would exist a solution κ 1 ∈ S 1 to the dispersion relation For convenience, the coordinates of η are denoted (η 2 , . . ., η d ).Using the definition (9) of a 1 , and defining κ := (κ 1 , e i η 2 , . . ., e i η d ), we have found a root z ∈ U to the relation but this is not possible because the s + 1 roots (in z) to the dispersion relation (8) are simple and belong to D. The Gauss-Lucas Theorem thus shows that the roots to the relation (42) belong to D (and therefore not to U ).At this stage, we know that the eigenvalues of M(z, η), (z, η) ∈ U × R d−1 , split into two groups: those in U , which we call the unstable ones, and those in D, which we call the stable ones.For (z, η) ∈ U ×R d−1 , for P 1 = max(p 1 , q 1 + 1) gives (recall the definition (45) of gn+s+1 We can apply Theorem 1 to the solution (w n j ) to (44) because the initial data in (44) vanish.We get:

Conclusion and perspectives
Let us first observe that in [Wad90], Wade has constructed symmetrizers for deriving stability estimates for multistep schemes, even in the case of variable coefficients.His conditions for constructing a symmetrizer are less restrictive than Assumption 2. However, the symmetrizer in [Wad90] is genuinely nonlocal and it is therefore not clear that it may be useful for boundary value problems.The main novelty here is to construct a local multiplier whose properties allow for the design of an auxiliary dissipative boundary value problem.This is our key to Theorem 1.
In this article we have always discarded the dissipation term provided by the nonnegative form D 0 .For the approximation of parabolic equations, this term may give some extra dissipation, but a crucial point to keep in mind is that the coefficients of the numerical scheme are assumed to be constant (which may in turn yield rather severe CFL conditions for implicit approximations of parabolic equations).Hence it does not seem very clear that our approach will yield stability estimates with "optimal" CFL conditions when approximating parabolic equations.This extension is left to further study in the future.
The main possible improvement of Theorem 1 would consist of assuming that only the roots to (8) that lie on S 1 are simple.Here we have assumed that all the roots, including those in D are simple.If we could manage to deal with multiple roots in D, then Theorem 1 would be applicable to numerical approximations of the transport equation ( 12) that are based on Adams-Bashforth methods of order 3 or higher (such methods have 0 as a root of multiplicity 2 or more at the zero frequency).
The results in this paper achieve the proof of a "weak form" of the conjecture in [KW93] that strong stability, in the sense of Definition 1, implies semigroup stability.However, an even stronger assumption was made in [KW93], namely that the sole fulfillment of the interior estimate when both the initial and boundary data for (5) vanish, does imply semigroup stability.The analogous conjecture for partial differential equations seems to be still open so far, but we do hope that our multiplier technique may yield some insight for dealing with the strong form of the conjecture in [KW93].

Figure 1 :
Figure 1: Reflection of a bump by the leap-frog scheme with homogeneous Dirichlet condition at four successive times.