The Leray-Gårding method for finite difference schemes
TThe Leray-G˚arding methodfor finite difference schemes
Jean-Fran¸cois
Coulombel ∗ September 20, 2018
Abstract
In [Ler53] and [G˚ar56],
Leray and
G˚arding have developed a multiplier technique for derivinga priori estimates for solutions to scalar hyperbolic equations in either the whole space or the torus.In particular, the arguments in [Ler53, G˚ar56] 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 isthe starting point of the argument by
Rauch in [Rau72] for the derivation of semigroup estimates forhyperbolic initial boundary value problems. In this article, we explain how this multiplier technique canbe adapted to the framework of finite difference approximations of transport equations. The techniqueapplies to numerical schemes with arbitrarily many time levels, and encompasses a somehow magicaltrick that has been known for a long time for the leap-frog scheme. More importantly, the existenceand properties of the local multiplier enable us to derive optimal semigroup estimates for fully discretehyperbolic initial boundary value problems, which answers a problem raised by
Trefethen , Kreiss and Wu [Tre84, KW93]. AMS classification:
Keywords: hyperbolic equations, difference approximations, stability, boundary conditions, semigroup.
Throughout this article, we use the notation U := { ζ ∈ C , | ζ | > } , U := { ζ ∈ C , | ζ | ≥ } , D := { ζ ∈ C , | ζ | < } , S := { ζ ∈ C , | ζ | = 1 } , D := D ∪ S . 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 ∗ denotesthe conjugate transpose of M . We let I denote the identity matrix or the identity operator when it actson 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 ) / . The induced matrix norm on M n ( C ) is denoted (cid:107) · (cid:107) .The letter C denotes a constant that may vary from line to line or within the same line. The dependenceof the constants on the various parameters is made precise throughout the text. ∗ CNRS and Universit´e de Nantes, Laboratoire de Math´ematiques Jean Leray (UMR CNRS 6629), 2 rue de la Houssini`ere,BP 92208, 44322 Nantes Cedex 3, France. Email: [email protected] . Research of the author wassupported by ANR project BoND, ANR-13-BS01-0009-01. a r X i v : . [ m a t h . NA ] M a y n what follows, we let d ≥ (cid:96) of square integrable sequences. Sequences maybe valued in C k for some integer k . Some sequences will be indexed by Z d − 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 > i = 1 , . . . , d be d space steps. We shall make use of the (cid:96) ( Z d − )-norm that we define as follows: for all v ∈ (cid:96) ( Z d − ), (cid:107) v (cid:107) (cid:96) ( Z d − ) := (cid:32) d (cid:89) k =2 ∆ x k (cid:33) d (cid:88) i =2 (cid:88) j i ∈ Z | v j ,...,j d | . The corresponding scalar product is denoted (cid:104)· , ·(cid:105) (cid:96) ( Z d − ) . Then for all integers m ≤ m , we set ||| u ||| m ,m := ∆ x m (cid:88) j = m (cid:107) u j , · (cid:107) (cid:96) ( Z d − ) , to denote the (cid:96) -norm on the set [ m , m ] × Z d − ( m may equal −∞ and m may equal + ∞ ). Thecorresponding scalar product is denoted (cid:104)· , ·(cid:105) m ,m . Other notation throughout the text is meant to beself-explanatory. The ultimate goal of this article is to derive semigroup estimates for finite difference approximationsof hyperbolic initial boundary value problems. Up to now, the only available general stability theoryfor such numerical schemes is due to
Gustafsson , Kreiss and
Sundstr¨om [GKS72]. It relies ona Laplace transform with respect to the time variable, and the corresponding stability estimates arethereby restricted to zero initial data. A long standing problem in this line of research is, starting fromthe GKS stability estimates, which are resolvent type estimates, to incorporate nonzero initial data andto derive semigroup estimates, see, e.g., the discussion in [Tre84, section 4]. This problem is delicate forthe following reason: the validity of the GKS stability estimate is known to be equivalent to a slightlystronger version of the resolvent estimatesup z ∈ U ( | z | − (cid:107) ( z I − T ) − (cid:107) L ( (cid:96) ( N )) < + ∞ , (1)where T is some bounded operator on (cid:96) ( N ) that incorporates both the discretization of the hyperbolicequation and the numerical boundary conditions. Deriving an optimal semigroup estimate amounts toshowing 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 equivalenceis known to fail in general in infinite dimension. Worse, even the strong resolvent conditionsup n ≥ sup z ∈ U ( | z | − n (cid:107) ( z I − T ) − n (cid:107) L ( (cid:96) ( N )) < + ∞ , does not imply in general that T is power bounded, see, e.g., the review [SW97] or [TE05] for details andhistorical comments.Optimal semigroup estimates have nevertheless been derived for some discretized hyperbolic initialboundary value problems. More specifically, the first general derivation of semigroup estimates is due2o 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 arbitraryspace dimension, but the arguments in [CG11] are still restricted to numerical schemes with two timelevels. The present article gives, as far as we are aware of, the first systematic derivation of semigroupestimates for fully discrete hyperbolic initial boundary value problems in the case of numerical schemeswith arbitrarily many time levels. It generalizes the arguments of [Wu95, CG11] and provides new insightfor 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 withby
Thomas [Tho72] by using a multiplier technique. It is precisely this technique which we aim atdeveloping 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 nj for the leap-frog scheme, see, e.g., [RM67],follows from a general theory that is the analogue of the Leray - G˚arding method for partial differentialequations, which we briefly recall now.The method by
Leray and
G˚arding [Ler53, G˚ar56] provides with suitable multipliers for scalarhyperbolic operators of arbitrary order. Namely, given an integer m ≥
0, we consider a partial differentialoperator of the form L := ∂ m +1 t + m +1 (cid:88) k =1 P k ( ∂ x ) ∂ m +1 − kt , where t ∈ R stands for the time variable, x ∈ R d stands for the space variable , and each operator P k ( ∂ x )is a linear combination of spatial partial derivatives of order k : P k ( ∂ x ) = (cid:88) | α | = k p k,α ∂ αx , ∂ αx := ∂ α x · · · ∂ α d x d , | α | := α + · · · + α d . In the above formula, the p k,α ’s are real numbers . Well-posedness of the Cauchy problem L u = 0 , ( u, ∂ t u, . . . , ∂ mt u ) | t =0 = ( u , u , . . . , u m ) , (2)in Sobolev spaces is known to be linked with hyperbolicity of L . Namely, if L is strictly hyperbolic,meaning that for all ξ ∈ R d \ { } , the (homogeneous) polynomial P ( τ, ξ ) := τ m +1 + m +1 (cid:88) k =1 P k ( i ξ ) τ m +1 − k , P k ( i ξ ) := i k (cid:88) | α | = k p k,α ξ α , (3)has m + 1 simple purely imaginary roots with respect to τ , then the Cauchy problem (2) is well-posed in H m ( R d ) × · · · × L ( R d ). In particular, there exists a constant C >
0, that is independent of the solution u and the initial data u , u , . . . , u m , such that there holds:sup t ∈ R m (cid:88) k =0 (cid:107) ∂ kt u ( t ) (cid:107) H m − k ( R d ) ≤ C m (cid:88) k =0 (cid:107) u k (cid:107) H m − k ( R d ) . (4)The method by Leray and
G˚arding gives a quick and elegant way to derive the estimate (4) assumingthat the solution u to (2) is sufficiently smooth. By standard duality arguments, the validity of the a The periodic case x ∈ T d can be dealt with in a similar way and is actually the one considered in [G˚ar56]. We restrict here for simplicity to linear operators with constant coefficients. u is sufficiently smooth and decayingat infinity so that all integration by parts arising in the computations are legitimate. The main idea isto find a suitable quantity M u , which we call a multiplier and that will be linear with respect to u , suchthat 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˚ar56,Section 3], one possible choice of a multiplier is given by L (cid:48) u where L (cid:48) stands for the partial differentialoperator of order m whose symbol is ∂ τ P , with P given in (3). Why L (cid:48) u is a good multiplier is justified in[Ler53, G˚ar56]. A well-known particular case is the choice of 2 ∂ t u as a multiplier for the wave equation.Here P ( τ, ξ ) = τ + | ξ | and therefore ∂ τ P = 2 τ , hence the choice 2 ∂ t u . The latter quantity is indeed asuitable multiplier for the wave operator because of the formula :2 ∂ t u ( ∂ t u − ∆ x u ) = ∂ t (cid:16) ( ∂ t u ) + d (cid:88) j =1 ( ∂ x j u ) (cid:17) − x (cid:0) ∂ t u ∇ x u (cid:1) . The important fact here is that the energy :( ∂ t u ) + d (cid:88) j =1 ( ∂ x j u ) , is a positive definite quadratic form of the first order partial derivatives of u . Let us observe that themultiplier L (cid:48) 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˚ar56] showthat this property is not specific to the wave operator. The fact that both the multiplier and the energyare 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 thearguments of [Rau72] which use time reversibility, but rather construct dissipative boundary conditionswhich will yield the optimal semigroup estimate we are aiming at. We first set a few notations. We let ∆ x , . . . , ∆ x d , ∆ t > Courant-Friedrichs-Lewy parameters, λ i := ∆ t/ ∆ x i , i = 1 , . . . , d , are fixedpositive constants. We keep ∆ t ∈ (0 ,
1] as a small parameter and let the space steps ∆ x , . . . , ∆ x d varyaccordingly. The (cid:96) -norms with respect to the space variables have been previously defined and thusdepend on ∆ t and the CFL parameters through the mesh volume (∆ x · · · ∆ x d on Z d − , and ∆ x · · · ∆ x d on Z d ). We always identify a sequence w indexed by either N (for time), Z d − or Z d (for space), withthe corresponding step function. In particular, we shall feel free to take Fourier or Laplace transforms ofsuch sequences.For all j ∈ Z d , we set j = ( j , j (cid:48) ) with j (cid:48) := ( j , . . . , j d ) ∈ Z d − . We let p, q, r ∈ N d denote some fixedmulti-integers, and define p , q , r , p (cid:48) , q (cid:48) , r (cid:48) according to the above notation. We also let s ∈ N denote We refer to [G˚ar56, page 74] for the generalization of such ”integration by parts” formula to partial derivatives of higherorder. s +1 (cid:88) σ =0 Q σ u n + σj = ∆ t F n + s +1 j , j ∈ Z d , j ≥ , n ≥ ,u n + s +1 j + s +1 (cid:88) σ =0 B j ,σ u n + σ ,j (cid:48) = g n + s +1 j , j ∈ Z d , j = 1 − r , . . . , , n ≥ ,u nj = f nj , j ∈ Z d , j ≥ − r , n = 0 , . . . , s , (5)where the operators Q σ and B j ,σ are given by: Q σ := p (cid:88) (cid:96) = − r p (cid:48) (cid:88) (cid:96) (cid:48) = − r (cid:48) a (cid:96),σ S (cid:96) , B j ,σ := q (cid:88) (cid:96) =0 q (cid:48) (cid:88) (cid:96) (cid:48) = − q (cid:48) b (cid:96),j ,σ S (cid:96) . (6)In (6), the a (cid:96),σ , b (cid:96),j ,σ are real numbers and are independent of the small parameter ∆ t (they may dependon the CFL parameters though), while S denotes the shift operator on the space grid: ( S (cid:96) v ) j := v j + (cid:96) for j, (cid:96) ∈ Z d . We have also used the short notation p (cid:48) (cid:88) (cid:96) (cid:48) = − r (cid:48) := d (cid:88) i =2 p i (cid:88) (cid:96) i = − r i , q (cid:48) (cid:88) (cid:96) (cid:48) = − q (cid:48) := d (cid:88) i =2 q i (cid:88) (cid:96) i = − q i . The numerical scheme (5) is understood as follows: one starts with (cid:96) initial data ( f j ), ..., ( f sj ) definedfor j ≥ − r . Assuming that the solution has been defined up to some time index n + s , n ≥
0, then thefirst and second equations in (5) should uniquely determine u n + s +1 j for j ≥ − r , j (cid:48) ∈ Z d − . The meshesassociated with j ≥ interior domain while those associated with j = 1 − r , . . . , discrete boundary . We wish to deal here simultaneously with explicit and implicit schemesand therefore make the following solvability assumption. Assumption 1 (Solvability of (5)) . The operator Q s +1 is an isomorphism on (cid:96) ( Z d ) . Moreover, for all ( F j ) ∈ (cid:96) ( N ∗ × Z d − ) and for all g − r , · , . . . , g , · ∈ (cid:96) ( Z d − ) , there exists a unique solution ( u j ) j ≥ − r ∈ (cid:96) to the equations (cid:40) Q s +1 u j = F j , j ∈ Z d , j ≥ ,u j + B j ,s +1 u ,j (cid:48) = g j , j ∈ Z d , j = 1 − r , . . . , . In particular, Assumption 1 is trivially satisfied in the case of explicit schemes for which Q s +1 is theidentity ( a (cid:96),s +1 = δ (cid:96) , · · · δ (cid:96) d , in (6), with δ the Kronecker symbol).The first and second equations in (5) therefore uniquely determine u n + s +1 j for j ≥ − r , and onethen proceeds to the following time index n + s + 2. Existence and uniqueness of a solution ( u nj ) to (5)follows from Assumption 1, so the last requirement for well-posedness is continuous dependence of thesolution on the three possible source terms ( F nj ), ( g nj ), ( f nj ). This is a stability problem for which severaldefinitions 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 isspecifically relevant when the boundary conditions are non-homogeneous (( g nj ) (cid:54)≡ 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 γ > and all ∆ t ∈ (0 , , the solution ( u nj ) to (5) with f j ) = · · · = ( f sj ) = 0 satisfies the estimate: γγ ∆ t + 1 (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t ||| u n ||| − r , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t p (cid:88) j =1 − r (cid:107) u nj , · (cid:107) (cid:96) ( Z d − ) ≤ C γ ∆ t + 1 γ (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t ||| F n ||| , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j =1 − r (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) . (7)The main contributions in [GKS72, Mic83] are to show that strong stability can be characterized by acertain 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 assumefrom the start that (5) is strongly stable. We can thus control, with zero initial data, (cid:96) type norms of thesolution 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 j ) , . . . , ( f sj ) in (cid:96) . Our main assumption is the following. Assumption 2 (Stability for the discrete Cauchy problem) . For all ξ ∈ R d , the dispersion relation s +1 (cid:88) σ =0 (cid:99) Q σ (e i ξ , . . . , e i ξ d ) z σ = 0 , (cid:99) Q σ ( κ ) := p (cid:88) (cid:96) = − r κ (cid:96) a (cid:96),σ , (8) 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 κ (cid:96) := κ (cid:96) · · · κ (cid:96) d d , for κ ∈ ( C \ { } ) d and (cid:96) ∈ Z d . From Assumption 1, we know that Q s +1 is an isomorphism on (cid:96) , which implies by Fourier analysisthat (cid:91) Q s +1 (e i ξ , . . . , e i ξ d ) does not vanish for any ξ ∈ R d . In particular, the dispersion relation (8) is apolynomial equation of degree s + 1 in z for any ξ ∈ R d . We now make the following assumption, whichalready appeared in [GKS72, Mic83] and several other works on the same topic. Assumption 3 (Noncharacteristic discrete boundary) . For (cid:96) = − r , . . . , p , z ∈ C and η ∈ R d − , let usdefine a (cid:96) ( z, η ) := s +1 (cid:88) σ =0 z σ p (cid:48) (cid:88) (cid:96) (cid:48) = − r (cid:48) a (cid:96),σ e i (cid:96) (cid:48) · η . (9) Then a − r and a p do not vanish on U × R d − , and they have nonzero degree with respect to z for all η ∈ R d − . Our main result is comparable with [Wu95, Theorem 3.3] and [CG11, Theorems 2.4 and 3.5] and showsthat strong stability (or ”GKS stability”) is a sufficient condition for incorporating (cid:96) initial conditionsin (5) and proving optimal semigroup estimates. The main price to pay in Assumption 2 is that theroots of the dispersion relation (8), which are nothing but the eigenvalues of the so-called amplificationmatrix for the Cauchy problem, need to be simple . This property is satisfied for instance by the leap-frogand modified leap-frog schemes in several space dimensions, under an appropriate CFL condition, seeParagraph 1.3. Our main result reads as follows. 6 heorem 1. Let Assumptions 1, 2 and 3 be satisfied, and assume that the scheme (5) is strongly stablein the sense of Definition 1. Then there exists a constant C such that for all γ > and all ∆ t ∈ (0 , ,the solution to (5) satisfies the estimate: sup n ≥ e − γ n ∆ t ||| u n ||| − r , + ∞ + γγ ∆ t + 1 (cid:88) n ≥ ∆ t e − γ n ∆ t ||| u n ||| − r , + ∞ + (cid:88) n ≥ ∆ t e − γ n ∆ t p (cid:88) j =1 − r (cid:107) u nj , · (cid:107) (cid:96) ( Z d − ) ≤ C s (cid:88) σ =0 ||| f σ ||| − r , + ∞ + γ ∆ t + 1 γ (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t ||| F n ||| , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j =1 − r (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) . (10) In particular, the scheme (5) is ”semigroup stable” in the sense that there exists a constant C such thatfor all ∆ t ∈ (0 , , the solution ( u nj ) to (5) with ( F nj ) = ( g nj ) = 0 satisfies the estimate sup n ≥ ||| u n ||| − r , + ∞ ≤ C s (cid:88) σ =0 ||| f σ ||| − r , + ∞ . (11) The scheme (5) is also (cid:96) -stable with respect to boundary data, see [Tre84, Definition 4.5], in the sensethat there exists a constant C such that for all ∆ t ∈ (0 , , the solution ( u nj ) to (5) with ( F nj ) = ( f nj ) = 0 satisfies the estimate sup n ≥ ||| u n ||| − r , + ∞ ≤ C (cid:88) n ≥ s +1 ∆ t (cid:88) j =1 − r (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) . Theorem 1 gives the optimal semigroup estimate (11), and is therefore an improvement with respect toour earlier work [Cou14] where in one space dimension, and under an appropriate non-glancing condition ,we were able to derive the estimate (here r = r , p = p since d = 1): γγ ∆ t + 1 (cid:88) n ≥ ∆ t e − γ n ∆ t ||| u n ||| − r, + ∞ + (cid:88) n ≥ ∆ t e − γ n ∆ t p (cid:88) j =1 − r | u nj | ≤ C s (cid:88) σ =0 ||| f σ ||| − r, + ∞ + γ ∆ t + 1 γ (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t ||| F n ||| , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j =1 − r | g nj | . The latter estimate does not incorporate on the left hand side the quantity:sup n ≥ e − γ n ∆ t ||| u n ||| − r, + ∞ , and was unfortunately still not sufficient for deriving the semigroup estimate (11). Our main contributionin this article is to exhibit a suitable multiplier for the multistep recurrence relation in (5). With thismultiplier, we can readily show that, for zero initial data, the (discrete) derivative of an energy can becontrolled, as in [Rau72], by the trace estimate of ( u nj ) and this is where strong stability comes into play. The non-glancing condition is unfortunately not met by the leap-frog scheme. . Hence we must find another argument for dealing with nonzero initialdata. Hopefully, the properties of our multiplier enable us to construct an auxiliary problem, where wemodify the boundary conditions of (5), and for which we can prove optimal semigroup and trace estimatesby ”hand-made” calculations. In other words, we exhibit an alternative set of boundary conditions thatyields strict dissipativity . Using these auxiliary numerical boundary conditions, the proof of Theorem1 follows from a standard superposition argument, see, e.g., [BGS07, Section 4.5] for partial differentialequations or [Wu95, CG11] for numerical schemes. Remark 1.
Assumption 3 excludes the case of explicit two level schemes for which s = 0 and Q = I ,for in that case a − r and/or a p 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 and a p are not constant isinvolved, and why the proof is actually simpler in the case s = 0 and Q = I . Our goal is to approximate the outgoing transport equation ( d = 1 here): ∂ t u + a ∂ x u = 0 , u | t =0 = u , (12)with t, x > 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 oneconsiders an upwind type scheme with a space stencil ”on the right” (meaning r = 0 in (5)). We nowdetail two possible multistep schemes for discretizing (12). Both are obtained by the so-called methodof lines, which amounts to first discretizing the space derivative ∂ x u and then choosing an integrationtechnique for discretizing the time evolution, see [GKO95]. The leap-frog scheme.
It is obtained by approximating the space derivative ∂ x u by the centereddifference ( u j +1 − u j − ) / (2 ∆ x ), and by then applying the so-called Nystr¨om method of order 2, see[HNW93, Chapter III.1]. The resulting approximation reads u n +2 j + λ a ( u n +1 j +1 − u n +1 j − ) − u nj = 0 , which corresponds to s = p = r = 1. Recall that λ > t/ ∆ x . Even though (12)does not require any boundary condition at x = 0, the leap-frog scheme stencil includes one point to theleft, and we therefore need to prescribe some numerical boundary condition at j = 0. One possibility is to prescribe the homogeneous or inhomogeneous Dirichlet boundary condition. With general sourceterms, the corresponding scheme reads u n +2 j + λ a ( u n +1 j +1 − u n +1 j − ) − u nj = ∆ t F n +2 j , j ≥ , n ≥ ,u n +20 = g n +20 , n ≥ , ( u j , u j ) = ( f j , f j ) , j ≥ . (13) With the notable exception of the leap-frog scheme that is indeed time reversible ! This is of course not the only possibility and we refer to [GKO95, Oli74, Slo83, Tre84] for some other possible choiceswhich might be more meaningful from a consistency and accuracy point of view. Our main concern here is a discussion onstability for (5) and the Dirichlet boundary conditions are a good illustration for this aspect. λ | a | <
1. In that case, the two roots to the dispersion relation z + 2 i λ a sin ξ z − , 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 ( z ) = − a − ( 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 propertyis given in the numerical simulation of a bump function, propagating at speed a = − j = 0. The reflection of the bump generatesa highly oscillatory wave packet that propagates with velocity +1 towards the right. The envelope of thiswave packet coincides with the profile of the initial condition, which indicates that the (cid:96) -norm is roughlypreserved 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 extensionare discussed, e.g., in [Oli74, Slo83, Tho72, Tre84]. The main discussion in [Oli74, Slo83, Tre84] is to verifystrong stability for a wide choice of numerical boundary conditions, and if strong stability holds, thenTheorem 1 automatically gives semigroup boundedness, which was not achieved in these earlier works.Figure 1: Reflection of a bump by the leap-frog scheme with homogeneous Dirichlet condition at foursuccessive times. 9 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 − ) / (2 ∆ x ), andthen apply the backwards differentiation formula of order 2, see [HNW93, Chapter III.1]. The resultingscheme reads: 32 u n +2 j + λ a u n +2 j +1 − u n +2 j − ) − u n +1 j + 12 u nj = 0 . This corresponds to s = 1 and Q : ( u j ) j ∈ Z (cid:55)−→ (cid:18) u j + λ a u j +1 − u j − ) (cid:19) j ∈ Z . The operator Q is an isomorphism on (cid:96) ( Z ) since Q is an isomorphism for any small λ a (as a perturbationof 3 / I ), Q depends continuously on λ a , and there holds (uniformly with respect to λ a ):32 ||| u ||| −∞ , + ∞ ≤ ||| Q u ||| −∞ , + ∞ . The operator Q is therefore an isomorphism on (cid:96) ( Z ) for any λ a > (cid:18)
32 + i λ a sin ξ (cid:19) z − z + 12 = 0 . It is clear that the latter equation has two simple roots in z for any ξ ∈ R . Moreover, if sin ξ = 0, the rootsare 1 and 1 / D . In the case sin ξ (cid:54) = 0, none of the roots belongs to S and examiningthe case λ a sin ξ = 1, we find that for sin ξ (cid:54) = 0, both roots belong to D (which is consistent with theshape of the stability region for the backwards differentiation formula of order 2, see [HW96, ChapterV.1]). Assumption 2 is therefore satisfied. Assumption 3 is satisfied as long as a is nonzero since thereholds p = r = 1 and a ( z ) = a − ( z ) = λ a z / λ a small enough) and strongstability holds. Here we wish to approximate the two-dimensional transport equation ( d = 2): ∂ t u + a ∂ x u + a ∂ x u = 0 , u | t =0 = u , in the space domain { x > , x ∈ R } . When a is negative, the latter problem does not necessitate anyboundary condition at x = 0. Following [AG76], we use one of the following two-dimensional versions ofthe leap-frog scheme, either u n +2 j,k + λ a ( u n +1 j +1 ,k − u n +1 j − ,k ) + λ a ( u n +1 j,k +1 − u n +1 j,k − ) − u nj,k = 0 , (14)or u n +2 j,k + λ a (cid:32) u n +1 j +1 ,k +1 + u n +1 j +1 ,k − − u n +1 j − ,k +1 + u n +1 j − ,k − (cid:33) + λ a (cid:32) u n +1 j +1 ,k +1 + u n +1 j − ,k +1 − u n +1 j +1 ,k − + u n +1 j − ,k − (cid:33) − u nj,k = 0 . (15)10ssumption 1 is trivially satisfied because (14) and (15) are explicit schemes. The scheme (14) satisfiesAssumption 2 if and only if λ | a | + λ | a | <
1, while the scheme (15) satisfies Assumption 2 if and onlyif max( λ | a | , λ | a | ) <
1. Let us now study when Assumption 3 is valid. For the scheme (14), we have r = p = 1, and a ( z, η ) = λ a z , a − ( z, η ) = − a ( z, η ) , so Assumption 3 is valid as long as a (cid:54) = 0. For the scheme (15), we have again r = p = 1, and a ( z, η ) = z ( λ a cos η + i λ a sin η ) , a − ( z, η ) = z ( − λ a cos η + i λ a sin η ) , so Assumption 3 is valid as long as both a and a are nonzero. We refer to [AG79] for the verificationof strong stability depending on the choice of some numerical boundary conditions for (14) or (15). Onceagain, if strong stability holds, then Theorem 1 yields semigroup boundedness and (cid:96) -stability with respectto boundary data. This section is devoted to proving stability estimates for discretized Cauchy problems, which is the firststep before considering the discretized initial boundary value problem (5). More precisely, we considerthe simpler case of the whole space j ∈ Z d , and the recurrence relation: s +1 (cid:88) σ =0 Q σ u n + σj = 0 , j ∈ Z d , n ≥ ,u nj = f nj , j ∈ Z d , n = 0 , . . . , s , (16)where the operators Q σ are given by (6). We recall that in (6), the a (cid:96),σ are real numbers and areindependent of the small parameter ∆ t (they may depend on the CFL parameters λ , . . . , λ d ), while S denotes the shift operator on the space grid: ( S (cid:96) v ) j := v j + (cid:96) for j, (cid:96) ∈ Z d . Stability of (16) is defined asfollows. Definition 2 (Stability for the discrete Cauchy problem) . The numerical scheme defined by (16) is ( (cid:96) -)stable if Q s +1 is an isomorphism from (cid:96) ( Z d ) onto itself, and if furthermore there exists a constant C > such that for all ∆ t ∈ (0 , , for all initial conditions ( f j ) j ∈ Z d , . . . , ( f sj ) j ∈ Z d in (cid:96) ( Z d ) , there holds sup n ∈ N ||| u n ||| −∞ , + ∞ ≤ C s (cid:88) σ =0 ||| f σ ||| −∞ , + ∞ . (17)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 (cid:81) i ∆ x i ), and can be characterizedin terms of the uniform power boundedness of the so-called amplification matrix A ( κ ) := − (cid:99) Q s ( κ ) / (cid:91) Q s +1 ( κ ) . . . . . . − (cid:99) Q ( κ ) / (cid:91) Q s +1 ( κ )1 0 . . .
00 . . . . . . ...0 0 1 0 ∈ M s +1 ( C ) , (18)where the (cid:99) Q σ ( κ )’s are defined in (8) and where it is understood that A is defined on the largest openset of C d on which (cid:91) Q s +1 does not vanish. Let us also recall that if Q s +1 is an isomorphism from (cid:96) ( Z d )11nto itself, then (cid:91) Q s +1 does not vanish on ( S ) d , and therefore does not vanish on an open neighborhoodof ( S ) 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 (cid:96) ( Z d ) onto itself. Then the scheme (16) is stable in the sense of Definition 2 ifand only if there exists a constant C > such that the amplification matrix A in (18) satisfies ∀ n ∈ N , ∀ ξ ∈ R d , (cid:12)(cid:12)(cid:12) A (e i ξ , . . . , e i ξ d ) n (cid:12)(cid:12)(cid:12) ≤ C . (19) In particular, the spectral radius of A (e i ξ , . . . , e i ξ d ) should not be larger than (the so-called von Neu-mann condition). The eigenvalues of A (e i ξ , . . . , e i ξ d ) are the roots to the dispersion relation (8). When these roots aresimple 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 discreteCauchy 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) bya suitable local multiplier. The analysis relies first on the simpler case where one only considers the timeevolution and no additional space variable. In this Paragraph, we consider sequences ( v n ) n ∈ N with values in C . The index n should be thought of asthe discrete time variable, and we therefore introduce the new notation T for the shift operator on thetime 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˚ar56, Lemme 1.1]. Lemma 1 (The energy-dissipation balance law) . Let P ∈ C [ X ] be a polynomial of degree s + 1 whoseroots are simple and located in D . Then there exists a positive definite Hermitian form q e on C s +1 , and anonnegative 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 ∀ n ∈ N , (cid:16) T ( P (cid:48) ( T ) v n ) P ( T ) v n (cid:17) = ( s + 1) | P ( T ) v n | + ( T − I ) ( q e ( v n , . . . , v n + s )) + q d ( v n , . . . , v n + s ) . In particular, for all sequence ( v n ) n ∈ N that satisfies the recurrence relation ∀ n ∈ N , P ( T ) v n = 0 , the sequence ( q e ( v n , . . . , v n + s )) n ∈ N is nonincreasing. The fact that there exists a Hermitian norm on C s +1 that is nonincreasing along solutions to therecurrence 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 directlythe ”energy boundedness” (or decay). The fact that the coefficients of this multiplier are integer multiplesof 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˚ar56, Lemme 1.1] and introduce the interpolation polynomials: ∀ k = 1 , . . . , s + 1 , P k ( X ) := a (cid:89) j (cid:54) = k ( X − x j ) , x , . . . , x s +1 denote the roots of P , and a (cid:54) = 0 its dominant coefficient. Since the roots of P arepairwise 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 P (cid:48) = s +1 (cid:88) k =1 P k . We then consider a sequence ( v n ) n ∈ N with values in C and compute2 Re (cid:16) T ( P (cid:48) ( T ) v n ) P ( T ) v n (cid:17) − ( s + 1) | P ( T ) v n | = s +1 (cid:88) k =1 T ( P k ( T )) v n ( T − x k ) P k ( T ) v n + T ( P k ( T ) v n ) ( T − x k ) P k ( T ) v n − s +1 (cid:88) k =1 ( T − x k )( P k ( T ) v n ) ( T − x k )( P k ( T ) v n )= s +1 (cid:88) k =1 ( T − | x k | ) | P k ( T ) v n | . The conclusion follows by defining: ∀ ( w , . . . , w s ) ∈ C s +1 , q e ( w , . . . , w s ) := s +1 (cid:88) k =1 | P k ( T ) w | , (20) q d ( w , . . . , w s ) := s +1 (cid:88) k =1 (1 − | x k | ) | P k ( T ) w | . (21)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 becausethe roots of P are located in D . Both forms depend in a C ∞ way on the coefficients of P because theroots of P are simple.Lemma 1 shows that the polynomial P (cid:48) yields the good multiplier T P (cid:48) ( T ) v n for the recurrencerelation P ( T ) v n = 0. Of course, P (cid:48) is not the only possible choice, though it will be our favorite one inwhat follows. As in [G˚ar56, Lemme 1.1], any polynomial of the form Q := s +1 (cid:88) k =1 α k P k , α , . . . , α s +1 > , provides with an energy balance of the form2 Re (cid:16) T ( Q ( T ) v n ) P ( T ) v n (cid:17) = ( α + · · · + α s +1 ) | P ( T ) v n | + ( T − I ) ( q e ( v n , . . . , v n + s )) + q d ( v n , . . . , v n + s ) , with suitable Hermitian forms q e , q d that have the same properties as stated in Lemma 1. The sign condition here on the coefficients α k is the analogue of the separation condition for the roots in [Ler53, G˚ar56]. .2 The energy-dissipation balance for finite difference schemes In this Paragraph, we consider the numerical scheme (16). We introduce the following notation: L := s +1 (cid:88) σ =0 T σ Q σ , M := s +1 (cid:88) σ =0 σ T σ Q σ . (22)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. Thenthere exist a continuous coercive quadratic form E and a continuous nonnegative quadratic form D on (cid:96) ( Z d ; R ) s +1 such that for all sequences ( v n ) n ∈ N with values in (cid:96) ( Z d ; R ) and for all n ∈ N , there holds (cid:104) M v n , L v n (cid:105) −∞ , + ∞ = ( s + 1) ||| L v n ||| −∞ , + ∞ + ( T − I ) E ( v n , . . . , v n + s ) + D ( v n , . . . , v n + s ) . In particular, for all initial data f , . . . , f s ∈ (cid:96) ( Z d ; R ) , the solution to (16) satisfies sup n ∈ N E ( v n , . . . , v n + s ) ≤ E ( f , . . . , f s ) , and (16) is ( (cid:96) -)stable.Proof. We use the same notation v n for the sequence ( v nj ) j ∈ Z d and the corresponding step function on R d whose value on the cell [ j ∆ x , ( j + 1) ∆ x ) × · · · × [ j d ∆ x d , ( j d + 1) ∆ x d ) equals v nj . Then PlancherelTheorem gives2 (cid:104) M v n , L v n (cid:105) −∞ , + ∞ − ( s + 1) ||| L v n ||| −∞ , + ∞ = (cid:90) R d (cid:16) T ( P (cid:48) ζ ( T ) (cid:99) v n ( ξ )) P ζ ( T ) (cid:99) v n ( ξ ) (cid:17) − ( s + 1) | P ζ ( T ) (cid:99) v n ( ξ ) | d ξ (2 π ) d , where (cid:99) v n denotes the Fourier transform of v n , and where we have let P ζ ( z ) := s +1 (cid:88) σ =0 (cid:99) Q σ (cid:0) e i ζ , . . . , e i ζ d (cid:1) z σ , ζ j := ξ j ∆ x j , and P (cid:48) ζ ( z ) denotes the derivative of P ζ with respect to z .From Assumption 2, we know that for all ζ ∈ R d , P ζ has degree s + 1 and has s + 1 simple roots in D . We can apply Lemma 1 and get2 (cid:104) M v n , L v n (cid:105) −∞ , + ∞ − ( s + 1) ||| L v n ||| −∞ , + ∞ = (cid:90) R d ( T − I ) q e,ζ ( (cid:99) v n ( ξ ) , . . . , (cid:100) v n + s ( ξ )) + q d,ζ ( (cid:99) v n ( ξ ) , . . . , (cid:100) v n + s ( ξ )) d ξ (2 π ) d , where q e,ζ , q d,ζ depend in a C ∞ way on ζ ∈ R d and are 2 π -periodic in each ζ j . Furthermore, q e,ζ is positivedefinite and q d,ζ is nonnegative. The conclusion of Proposition 2 follows by a standard compactnessargument for showing coercivity of E . 14 .3 Examples The first basic example corresponds to the case s = 0 for which the multiplier provided by Proposition 2is Q v n +1 j . In that case, the energy E reads ||| Q v ||| −∞ , + ∞ (recall that Q is an isomorphism) and theenergy-dissipation balance law is nothing but the trivial identity2 (cid:104) Q v n +1 , Q v n +1 + Q v n (cid:105) −∞ , + ∞ = ||| Q v n +1 + Q v n ||| −∞ , + ∞ + ||| Q v n +1 ||| −∞ , + ∞ − ||| Q v n ||| −∞ , + ∞ . The second line of this algebraic identity can be rewritten as ||| Q v n +1 ||| −∞ , + ∞ − ||| Q v n ||| −∞ , + ∞ + ||| Q v n ||| −∞ , + ∞ − ||| Q v n ||| −∞ , + ∞ , and (cid:96) -stability for the Cauchy problem amounts to assuming that the operator norm of Q − Q is notlarger than 1. Hence the dissipation term ||| Q v n ||| −∞ , + ∞ − ||| Q v n ||| −∞ , + ∞ is nonnegative.Let us now consider the leap-frog scheme in one space dimension, for which we have s = 1 and L = T + λ a T ( S − S − ) − I .
The corresponding dispersion relation (8) reduces to z + 2 i λ a sin ξ z − . For λ | a | <
1, the latter equation has two simple roots x ( ξ ) , x ( ξ ) of modulus 1. Following the previousanalysis, see (20)-(21), the form q e,ζ is given by q e,ζ ( w , w ) = | w − x ( ζ ) w | + | w − x ( ζ ) w | = 2 | w | + 2 | w | + 4 λ a Re (cid:0) i sin ζ w w (cid:1) , and q d,ζ is zero. The associated forms in Proposition 2 are D ≡ d = 1): E ( v , v ) = 2 (cid:88) j ∈ Z ∆ x ( v j ) + 2 (cid:88) j ∈ Z ∆ x ( v j ) + 2 λ a (cid:88) j ∈ Z ∆ x ( v j +1 − v j − ) v j . The latter energy functional E is coercive under the condition λ | a | <
1, which is the necessary andsufficient condition of stability for the leap-frog scheme, and E is conserved for solutions to the leap-frogscheme. The conservation of E is usually proved by starting from the recurrence relation ∀ j ∈ Z , ∀ n ∈ N , u n +2 j + λ a ( u n +1 j +1 − u n +1 j − ) − u nj = 0 , using the multiplier u n +2 j + u nj , and summing with respect to j . This is equivalent, for solutions to theleap-frog scheme, to what we propose here, since our multiplier reads M u nj = 2 u n +2 j + λ a ( u n +1 j +1 − u n +1 j − ) = u n +2 j + u nj + L u nj (cid:124)(cid:123)(cid:122)(cid:125) =0 . However, it will appear more clearly in Section 3 why our choice for
M u nj has a major advantage whenconsidering initial boundary value problems.Let us observe here that the energy functional E is associated with a local energy density E ,j ( v , v ) := 2 ( v j ) + 2 ( v j ) + 2 λ a ( v j +1 − v j − ) v j . 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 , D do not necessarily admit local densities. Thisis one main difference with [Ler53, G˚ar56]. 15 Semigroup estimates for fully discrete initial boundary value prob-lems
We now turn to the proof of Theorem 1 for which we shall use the results of Section 2 as a toolbox. Bylinearity of (5), it is sufficient to prove Theorem 1 separately in the case ( f j ) = · · · = ( f sj ) = 0, and in thecase ( F nj ) = 0, ( g nj ) = 0. The latter case is the most difficult and requires the introduction of an auxiliaryset of “dissipative” boundary conditions. Solutions to (5) are always assumed to be real valued, whichmeans that the data are real valued. For complex valued initial data and/or forcing terms, one just usesthe linearity of (5). We first assume ( f j ) = · · · = ( f sj ) = 0. By strong stability, we already know that (7) holds with a constant C that is independent of γ > t ∈ (0 , C , that is independent of γ > t ∈ (0 ,
1] such that the solution to (5) with( f j ) = · · · = ( f sj ) = 0 satisfiessup n ≥ e − γ n ∆ t ||| u n ||| − r , + ∞ ≤ C γ ∆ t + 1 γ (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t ||| F n ||| , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j =1 − r (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) . (23)We thus consider a parameter γ > t ∈ (0 , f j ) = · · · = ( f sj ) = 0). For all n ∈ N , we extend the sequence ( u nj ) by zerofor j ≤ − r : v nj := (cid:40) u nj if j ≥ − r , T − I ) E ( v n , . . . , v n + s ) + D ( v n , . . . , v n + s ) = 2 (cid:104) M v n , L v n (cid:105) −∞ , + ∞ − ( s + 1) ||| L v n ||| −∞ , + ∞ . Due to the form of the operator L , see (22), and the fact that v nj vanishes for j ≤ − r , there holds: L v nj = (cid:40) ∆ t F n + s +1 j if j ≥ , j ≤ − r − p , and we thus get( T − I ) E ( v n , . . . , v n + s ) + D ( v n , . . . , v n + s )= (cid:32) d (cid:89) k =1 ∆ x k (cid:33) (cid:88) j ≥ (cid:88) j (cid:48) ∈ Z d − t ( M v nj ) F n + s +1 j − ( s + 1) ∆ t ( F n + s +1 j ) + (cid:32) d (cid:89) k =1 ∆ x k (cid:33) (cid:88) j =1 − r − p (cid:88) j (cid:48) ∈ Z d − M v nj ) L v nj − ( s + 1) ( L v nj ) .
16e multiply the latter equality by exp( − γ ( n + s + 1) ∆ t ), sum with respect to n from 0 to some N anduse the fact that D is nonnegative. Recalling that the initial data in (5) vanish, we gete − γ ( N + s +1) ∆ t E (cid:0) v N +1 , . . . , v N + s +1 (cid:1) + (cid:0) − e − γ ∆ t (cid:1) N (cid:88) n =1 e − γ ( n + s ) ∆ t E ( v n , . . . , v n + s ) (cid:124) (cid:123)(cid:122) (cid:125) ≥ ≤ S ,N + S ,N , (24)with S ,N := N (cid:88) n =0 e − γ ( n + s +1) ∆ t (cid:16) t (cid:104) M v n , F n + s +1 (cid:105) , + ∞ − ( s + 1) ∆ t ||| F n + s +1 ||| , + ∞ (cid:17) , (25)and S ,N := (cid:32) d (cid:89) k =1 ∆ x k (cid:33) N (cid:88) n =0 e − γ ( n + s +1) ∆ t (cid:88) j =1 − r − p (cid:88) j (cid:48) ∈ Z d − M v nj ) L v nj − ( s + 1) ( L v nj ) . (26)Let us now estimate the two source terms S ,N , S ,N in (24). We begin with the term S ,N defined in(26). Let us recall that the ratio ∆ t/ ∆ x is fixed. Furthermore, the form of the operators L and M in(22) gives the estimate (recall that v nj vanishes for j ≤ − r ): S ,N ≤ C ∆ t (cid:32) d (cid:89) k =2 ∆ x k (cid:33) N (cid:88) n =0 e − γ ( n + s +1) ∆ t p (cid:88) j =1 − r (cid:88) j (cid:48) ∈ Z d − ( u nj ) + · · · + ( u n + s +1 j ) , for a constant C that does not depend on N , γ nor on ∆ t . We thus have, uniformly with respect to N ∈ N , γ > t ∈ (0 , S ,N ≤ C N + s +1 (cid:88) n = s +1 ∆ t e − γ n ∆ t p (cid:88) j =1 − r (cid:107) u nj , · (cid:107) (cid:96) ( Z d − ) ≤ C (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t p (cid:88) j =1 − r (cid:107) u nj , · (cid:107) (cid:96) ( Z d − ) ≤ C γ ∆ t + 1 γ (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t ||| F n ||| , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j =1 − r (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) , where we have used the trace estimate (7) that follows from the strong stability assumption.Let us now focus on the term S ,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 stability17ssumption): S ,N ≤ N (cid:88) n =0 ∆ t e − γ ( n + s +1) ∆ t ||| M v n ||| , + ∞ ||| F n + s +1 ||| , + ∞ ≤ C N (cid:88) n =0 ∆ t e − γ ( n + s +1) ∆ t (cid:16) ||| v n +1 ||| − r , + ∞ + · · · + ||| v n + s +1 ||| − r , + ∞ (cid:17) ||| F n + s +1 ||| , + ∞ ≤ C γγ ∆ t + 1 N + s +1 (cid:88) n = s +1 ∆ t e − γ n ∆ t ||| u n ||| − r , + ∞ + C γ ∆ t + 1 γ N + s +1 (cid:88) n = s +1 ∆ t e − γ n ∆ t ||| F n ||| , + ∞ ≤ C γ ∆ t + 1 γ (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t ||| F n ||| , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j =1 − r (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) . Ignoring the nonnegative term on the left hand-side of (24) and using the coercivity of E , we have provedthat there exists a constant C >
N, γ, ∆ t such that:e − γ ( N + s +1) ∆ t ||| v N + s +1 ||| −∞ , + ∞ ≤ C γ ∆ t + 1 γ (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t ||| F n ||| , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j =1 − r (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) , which yields (23) and therefore the validity of Theorem 1 in the case of zero initial data. In this paragraph, we consider an auxiliary problem for which we shall be able to prove simultaneouslyan optimal semigroup estimate and a trace estimate for the solution. More precisely, we shall prove thefollowing result.
Theorem 2.
Let Assumptions 1, 2 and 3 be satisfied. Then for all P ∈ N , there exists a constant C P > such that, for all initial data ( f j ) , . . . , ( f sj ) ∈ (cid:96) ( Z d ) and for all source term ( g nj ) j ≤ ,n ≥ s +1 thatsatisfies ∀ Γ > , (cid:88) n ≥ s +1 e − n (cid:88) j ≤ (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) < + ∞ , there exists a unique sequence ( u nj ) j ∈ Z d ,n ∈ N solution to L u nj = 0 , j ∈ Z d , j ≥ , n ≥ ,M u nj = g n + s +1 j , j ∈ Z d , j ≤ , n ≥ ,u nj = f nj , j ∈ Z d , n = 0 , . . . , s . (27)18 oreover for all γ > and ∆ t ∈ (0 , , this solution satisfies sup n ≥ e − γ n ∆ t ||| u n ||| −∞ , + ∞ + γγ ∆ t + 1 (cid:88) n ≥ ∆ t e − γ n ∆ t ||| u n ||| −∞ , + ∞ + (cid:88) n ≥ ∆ t e − γ n ∆ t P (cid:88) j =1 − r (cid:107) u nj , · (cid:107) (cid:96) ( Z d − ) ≤ C P s (cid:88) σ =0 ||| f σ ||| −∞ , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j ≤ (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) . (28)Theorem 2 justifies why we advocate the choice M u nj = 2 u n +2 j + λ a ( u n +1 j +1 − u n +1 j − ) rather than themore standard u n +2 j + u nj as a multiplier for the leap-frog scheme. Despite repeated efforts, we have notbeen able to prove the estimate (28) when using the numerical boundary condition u n +2 j + u nj on j ≤ j ≥ Proof.
Let us first quickly observe that the solution to (27) is well-defined since, as long as we havedetermined the solution up to a time index n + s , n ≥
0, then u n + s +1 is sought as a solution to anequation of the form Q s +1 u n + s +1 = F , where F belongs to (cid:96) ( Z d ) (this is due to the form of L and M , see (22)). Hence u n is uniquely definedand belongs to (cid:96) ( Z d ) for all n ∈ N .The proof of Theorem 2 starts again with the application of Proposition 2. Using the nonnegativityof the dissipation form D , we get ( T − I ) E ( u n , . . . , u n + s ) + ( s + 1) ||| L u n ||| −∞ , + ∞ ≤ (cid:104) M u n , L u n (cid:105) −∞ , + ∞ = 2 (cid:104) g n + s +1 , L u n (cid:105) −∞ , . By the Young inequality, we get( T − I ) E ( u n , . . . , u n + s ) + s + 12 ||| L u n ||| −∞ , + ∞ ≤ s + 1 ||| g n + s +1 ||| −∞ , . We multiply the latter inequality by exp( − γ ( n + s + 1) ∆ t ), sum from n = 0 to some arbitrary N andalready derive the estimate (here we use again the fact that ∆ t/ ∆ x is a fixed positive constant):sup n ≥ e − γ ( n + s ) ∆ t E ( u n , . . . , u n + s ) + (cid:0) − e − γ ∆ t (cid:1) (cid:88) n ≥ e − γ ( n + s ) ∆ t E ( u n , . . . , u n + s )+ (cid:88) n ≥ ∆ t e − γ ( n + s +1) ∆ t (cid:88) j ∈ Z (cid:107) L u nj , · (cid:107) (cid:96) ( Z d − ) ≤ C e − γ s ∆ t E ( f , . . . , f s ) + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j ≤ (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) . Using the coercivity of E and the inequality1 − e − γ ∆ t ≥ γ ∆ tγ ∆ t + 1 , Since
L u nj = 0 for j ≥
1, one could also write |||
L u n ||| −∞ , rather than ||| L u n ||| −∞ , + ∞ on the left hand-side of theinequality.
19e have therefore derived the estimatesup n ≥ e − γ n ∆ t ||| u n ||| −∞ , + ∞ + γγ ∆ t + 1 (cid:88) n ≥ ∆ t e − γ n ∆ t ||| u n ||| −∞ , + ∞ + (cid:88) n ≥ ∆ t e − γ ( n + s +1) ∆ t (cid:88) j ∈ Z (cid:107) L u nj , · (cid:107) (cid:96) ( Z d − ) ≤ C s (cid:88) σ =0 ||| f σ ||| −∞ , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j ≤ (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) , (29)where the constant C is independent of γ , ∆ t and on the solution ( u nj ). In order to prove (28), the mainremaining task is to derive the trace estimate for ( u nj ). 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 > J suchthat ( L u nj ) ≥
12 ( Q s +1 u n + s +1 j ) − C s (cid:88) σ =0 (cid:88) | (cid:96) |≤ J ( u n + σj + (cid:96) ) . Since Q s +1 is an isomorphism, there exists a constant c > (cid:88) j ∈ Z d ( L u nj ) ≥ c (cid:88) j ∈ Z d ( u n + s +1 j ) − c s (cid:88) σ =0 (cid:88) j ∈ Z d ( u n + σj ) . Multiplying by exp( − γ ( n + s + 1) ∆ t ) and summing with respect to n ∈ N , we get (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j ∈ Z (cid:107) u nj , · (cid:107) (cid:96) ( Z d − ) ≤ C (cid:88) n ≥ ∆ t e − γ ( n + s +1) ∆ t (cid:88) j ∈ Z (cid:107) L u nj , · (cid:107) (cid:96) ( Z d − ) +e − γ ∆ t (cid:88) n ≥ ∆ t e − γ n ∆ t (cid:88) j ∈ Z (cid:107) u nj , · (cid:107) (cid:96) ( Z d − ) . Choosing γ ∆ t large enough, that is γ ∆ t ≥ ln R for some numerical constant R > L , we have derived the estimate (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j ∈ Z (cid:107) u nj , · (cid:107) (cid:96) ( Z d − ) ≤ C (cid:88) n ≥ ∆ t e − γ ( n + s +1) ∆ t (cid:88) j ∈ Z (cid:107) L u nj , · (cid:107) (cid:96) ( Z d − ) +e − γ ∆ t s (cid:88) σ =0 ||| f σ ||| −∞ , + ∞ (cid:41) . It remains to use (29) and we get an even better estimate than (28) which we were originally aiming at: (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j ∈ Z (cid:107) u nj , · (cid:107) (cid:96) ( Z d − ) ≤ C s (cid:88) σ =0 ||| f σ ||| −∞ , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j ≤ (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) . γ ∆ t can be small). • From now on, we have fixed a constant R > γ ∆ t ≥ ln R and we thusassume γ ∆ t ∈ (0 , ln R ]. 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 nj ) for finitely many values of j .We first observe from (29) that for all γ > t ∈ (0 , C γ, ∆ t such that ∀ n ∈ N , e − γ n ∆ t (cid:88) j ∈ Z d ( u nj ) ≤ C γ, ∆ t . In particular, for any j ∈ Z , the Laplace-Fourier transforms (cid:99) u j of the step functions u j : ( t, y ) ∈ R + × R d − (cid:55)→ u nj if ( t, y ) ∈ [ n ∆ t, ( n + 1) ∆ t ) × d (cid:89) k =2 [ j k ∆ x k , ( j k + 1) ∆ x k ) , is well-defined on { τ ∈ C , Re τ > } × R d − . The dual variables are denoted τ = γ + i θ , γ >
0, and η = ( η , . . . , η d ) ∈ R d − . It will also be convenient to introduce the notation η ∆ := ( η ∆ x , . . . , η d ∆ x d ).Given Γ >
0, the sequence ( (cid:99) u j (Γ + i θ, η )) j ∈ Z belongs to (cid:96) ( Z ) for almost every ( θ, η ) ∈ R × R d − .We first show the following estimate. Lemma 2.
With R > fixed as above, there exists a constant C > such that for all γ > and ∆ t ∈ (0 , satisfying γ ∆ t ∈ (0 , ln R ] , there holds (cid:88) j ∈ Z (cid:90) R × R d − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p (cid:88) (cid:96) = − r a (cid:96) (cid:0) e ( γ + i θ ) ∆ t , η ∆ (cid:1) (cid:92) u j + (cid:96) ( γ + i θ, η ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d θ d η + (cid:88) j ≤ (cid:90) R × R d − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p (cid:88) (cid:96) = − r e ( γ + i θ ) ∆ t ∂ z a (cid:96) (cid:0) e ( γ + i θ ) ∆ t , η ∆ (cid:1) (cid:92) u j + (cid:96) ( γ + i θ, η ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d θ d η ≤ C s (cid:88) σ =0 ||| f σ ||| −∞ , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j ≤ (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) . (30) Proof of Lemma 2.
Given τ = γ + i θ and η , we compute (here j ∈ Z is fixed): p (cid:88) (cid:96) = − r a (cid:96) (cid:0) e τ ∆ t , η ∆ (cid:1) (cid:92) u j + (cid:96) ( τ, η ) = (cid:92) L u j , · ( τ, η ) + 1 − e − τ ∆ t τ s +1 (cid:88) σ =1 σ − (cid:88) σ (cid:48) =0 e ( σ − σ (cid:48) ) τ ∆ t F σ,σ (cid:48) j ( η ) , (31) p (cid:88) (cid:96) = − r e τ ∆ t ∂ z a (cid:96) (cid:0) e τ ∆ t , η ∆ (cid:1) (cid:92) u j + (cid:96) ( τ, η ) = (cid:92) M u j , · ( τ, η ) + 1 − e − τ ∆ t τ s +1 (cid:88) σ =1 σ − (cid:88) σ (cid:48) =0 σ e ( σ − σ (cid:48) ) τ ∆ t F σ,σ (cid:48) j ( η ) . (32)where, in (31) and (32), we have set F σ,σ (cid:48) j ( η ) = p (cid:88) (cid:96) = − r p (cid:48) (cid:88) (cid:96) (cid:48) = − r (cid:48) a (cid:96),σ e i (cid:96) (cid:48) · η ∆ (cid:92) f σ (cid:48) j + (cid:96) , · ( η ) , y = ( x , . . . , x d ) ∈ R d − , of the stepfunction associated with the sequence ( Q σ f σ (cid:48) j ) (no Laplace transform here).We need to estimate integrals with respect to ( θ, η ) of the right hand side of (31) and (32). The firstterm on the right of (31) and (32) are easy. For instance, we have (applying Plancherel Theorem): (cid:88) j ∈ Z (cid:90) R × R d − (cid:12)(cid:12) (cid:92) L u j , · ( τ, η ) (cid:12)(cid:12) d θ d η = (2 π ) d (cid:88) j ∈ Z (cid:88) n ≥ (cid:90) ( n +1) ∆ tn ∆ t e − γ s (cid:107) L u nj , · (cid:107) (cid:96) ( Z d − ) d s = (2 π ) d − e − γ ∆ t γ ∆ t (cid:88) n ≥ ∆ t e − γ n ∆ t (cid:88) j ∈ Z (cid:107) L u nj , · (cid:107) (cid:96) ( Z d − ) . We now recall that γ ∆ t is restricted to the interval (0 , ln R ], and we use (29) to derive (cid:88) j ∈ Z (cid:90) R × R d − (cid:12)(cid:12) (cid:92) L u j , · ( τ, η ) (cid:12)(cid:12) d θ d η ≤ C s (cid:88) σ =0 ||| f σ ||| −∞ , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j ≤ (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) . Similarly, we have (cid:88) j ≤ (cid:90) R × R d − (cid:12)(cid:12) (cid:92) M u j , · ( τ, η ) (cid:12)(cid:12) d θ d η = (2 π ) d − e − γ ∆ t γ ∆ t (cid:88) n ≥ ∆ t e − γ n ∆ t (cid:88) j ≤ (cid:107) M u nj , · (cid:107) (cid:96) ( Z d − ) , 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 σ (cid:48) ): (cid:88) j ∈ Z (cid:90) R × R d − (cid:12)(cid:12)(cid:12)(cid:12) − e − τ ∆ t τ (cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12) F σ,σ (cid:48) j ( η ) (cid:12)(cid:12) d θ d η = (cid:32)(cid:90) R (cid:12)(cid:12)(cid:12)(cid:12) − e − τ ∆ t τ (cid:12)(cid:12)(cid:12)(cid:12) d θ (cid:33) (cid:88) j ∈ Z (cid:90) R d − (cid:12)(cid:12) F σ,σ (cid:48) j ( η ) (cid:12)(cid:12) d η , where we have applied Fubini Theorem. Applying first Plancherel Theorem with respect to the d − (cid:88) j ∈ Z (cid:90) R d − (cid:12)(cid:12) F σ,σ (cid:48) j ( η ) (cid:12)(cid:12) d η ≤ C (cid:88) j ∈ Z (cid:88) j (cid:48) ∈ Z d − (cid:32) d (cid:89) k =2 ∆ x k (cid:33) ( f σ (cid:48) j ) ≤ C ∆ t s (cid:88) σ =0 ||| f σ ||| −∞ , + ∞ . The conclusion then follows by computing (cid:90) R (cid:12)(cid:12)(cid:12)(cid:12) − e − τ ∆ t τ (cid:12)(cid:12)(cid:12)(cid:12) d θ = 2 π ∆ t − e − γ ∆ t γ ∆ t , and by recalling that γ ∆ t belongs to (0 , ln R ]. We can eventually bound the integrals on the left handside 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 > be fixed as aboveand let P ∈ N . Then there exists a constant C P > such that for all z ∈ U with | z | ≤ R , for all η ∈ R d − and for all sequence ( w j ) j ∈ Z ∈ (cid:96) ( Z ; C ) , there holds P (cid:88) j = − r − p | w j | ≤ C P (cid:88) j ∈ Z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p (cid:88) (cid:96) = − r a (cid:96) ( z, η ) w j + (cid:96) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:88) j ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p (cid:88) (cid:96) = − r z ∂ z a (cid:96) ( z, η ) w j + (cid:96) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (33) Recall that the functions a (cid:96) , (cid:96) = − r , . . . , p , are defined in (9) . z = exp( τ ∆ t ), τ = γ + i θ with γ ∆ t ∈ (0 , ln R ], η ∆ ∈ R d − and the sequence ( (cid:99) u j ( τ, η )) j ∈ Z . We then integrate (33) with respect to ( θ, η ) and use Lemma2 to derive P (cid:88) j = − r − p (cid:90) R × R d − | (cid:99) u j ( γ + i θ, η ) | d θ d η ≤ C s (cid:88) σ =0 ||| f σ ||| −∞ , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j ≤ (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) . It remains to apply Plancherel Theorem and we get P (cid:88) j = − r − p (cid:88) n ∈ N − e − γ ∆ t γ ∆ t ∆ t e − γ n ∆ t (cid:107) u nj , · (cid:107) (cid:96) ( Z d − ) ≤ C s (cid:88) σ =0 ||| f σ ||| −∞ , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j ≤ (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) . Recalling that γ ∆ t is restricted to the interval (0 , ln R ], we have thus derived the trace estimate (cid:88) n ∈ N ∆ t e − γ n ∆ t P (cid:88) j = − r − p (cid:107) u nj , · (cid:107) (cid:96) ( Z d − ) ≤ C s (cid:88) σ =0 ||| f σ ||| −∞ , + ∞ + (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j ≤ (cid:107) g nj , · (cid:107) (cid:96) ( Z d − ) . Combined with the semigroup and interior estimate (29), this gives the estimate (28) of Theorem 2 for γ ∆ t ∈ (0 , ln R ]. Proof of Lemma 3.
Let us recall that the functions a (cid:96) are 2 π -periodic with respect to each coordinate of η . We can therefore restrict to η ∈ [0 , π ] d − rather than considering η ∈ R d − . We argue by contradictionand assume that the conclusion to Lemma 3 does not hold. This means the following, up to normalizingand extracting subsequences; there exist three sequences (indexed by k ∈ N ): • a sequence ( w k ) k ∈ N with values in (cid:96) ( Z ; C ) such that ( w k − r − p , . . . , w kP ) belongs to the unit sphereof C P + r + p +1 for all k , and ( w k − r − p , . . . , w kP ) converges towards ( w − r − p , . . . , w P ) as k tends toinfinity, • a sequence ( z k ) k ∈ N with values in U ∩ { ζ ∈ C , | ζ | ≤ R } , which converges towards z ∈ U , • a sequence ( η k ) k ∈ N with values in [0 , π ] d − , which converges towards η ∈ [0 , π ] d − ,and these sequences satisfy:lim k → + ∞ (cid:88) j ∈ Z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p (cid:88) (cid:96) = − r a (cid:96) ( z k , η k ) w kj + (cid:96) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:88) j ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p (cid:88) (cid:96) = − r z k ∂ z a (cid:96) ( z k , η k ) w kj + (cid:96) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 0 . (34)We are going to show that (34) implies that ( w − r − p , . . . , w P ) must be zero, which will yield a contra-diction since this vector must have norm 1. 23 Let us first show that each component ( w kj ) k ∈ N , j ∈ Z , has a limit as k tends to infinity. This isalready clear for j = − r − p , . . . , P . For j > P , we argue by induction. From (34), we havelim k → + ∞ p (cid:88) (cid:96) = − r a (cid:96) ( z k , η k ) w kP − p +1+ (cid:96) = 0 , and by Assumption 3, we know that a p ( z, η ) is nonzero. Hence ( w kP +1 ) k ∈ N converges towards − a p ( z, η ) p − (cid:88) (cid:96) = − r a (cid:96) ( z, η ) w P − p +1+ (cid:96) , which we define as w P +1 . We can argue by induction in the same way for all indices j > P + 1, butalso for indices j < − r − p because the function a − r also does not vanish on U × R d − .Using (34), we have thus shown that for each j ∈ Z , ( w kj ) k ∈ N tends towards some limit w j as k tends to infinity, and the sequence w , which does not necessarily belong to (cid:96) ( Z ; C ), satisfies the inductionrelations: ∀ j ∈ Z , p (cid:88) (cid:96) = − r a (cid:96) ( z, η ) w j + (cid:96) = 0 , (35) ∀ j ≤ , p (cid:88) (cid:96) = − r z ∂ z a (cid:96) ( z, η ) w j + (cid:96) = 0 . (36) • The induction relation (35) is the one that arises in [GKS72, Mic83] and all the works that dealwith strong stability. The main novelty here is to use simultaneously (35) for controlling the unstablecomponents of ( w − r − p , . . . , w − ) and (36) for controlling the stable components of ( w − r − p , . . . , w − ).The fact that w satisfies simultaneously (35) and (36) for j ≤ F kj := p (cid:88) (cid:96) = − r a (cid:96) ( z k , η k ) w kj + (cid:96) , G kj := p (cid:88) (cid:96) = − r z k ∂ z a (cid:96) ( z k , η k ) w kj + (cid:96) , which, according to (34), satisfylim k → (cid:88) j ∈ Z | F kj | = 0 , lim k → (cid:88) j ≤ | G kj | = 0 . (37)We also introduce the vectors (here T denotes transposition) W kj := (cid:16) w kj + p , . . . , w kj +1 − r (cid:17) T , W j := (cid:16) w j + p , . . . , w j +1 − r (cid:17) T , L ( z, η ) := − a p − ( z, η ) /a p ( z, η ) . . . . . . − a − r ( z, η ) /a p ( z, η )1 0 . . .
00 . . . . . . ...0 0 1 0 ∈ M p + r ( C ) , (38) M ( z, η ) := − ∂ z a p − ( z, η ) /∂ z a p ( z, η ) . . . . . . − ∂ z a − r ( z, η ) /∂ z a p ( z, η )1 0 . . .
00 . . . . . . ...0 0 1 0 ∈ M p + r ( C ) . (39)The matrix L is well-defined on U × R d − according to Assumption 3. The matrix M is also well-definedon U × R d − because for any η ∈ R d − , Assumption 3 asserts that a p ( · , η ) is a nonconstant polynomialwhose roots lie in D . From the Gauss-Lucas Theorem, the roots of ∂ z a p ( · , η ) lie in the convex hull ofthose of a p ( · , η ). Therefore ∂ z a p ( · , η ) does not vanish on U . In the same way, ∂ z a − r ( · , η ) does notvanish on U .With our above notation, the vectors W kj , W j , satisfy the one step induction relations: ∀ j ∈ Z , W kj +1 = L ( z k , η k ) W kj + (cid:16) F kj +1 /a p ( z k , η k ) , , . . . , (cid:17) T , W j +1 = L ( z, η ) W j , (40) ∀ j ≤ − , W kj +1 = M ( z k , η k ) W kj + (cid:16) G kj +1 / ( z k ∂ z a p ( z k , η k )) , , . . . , (cid:17) T , W j +1 = M ( z, η ) W j . (41) • From Assumption 3 and the above application of the Gauss-Lucas Theorem, we already know thatboth matrices L ( z, η ) and M ( z, η ) are invertible for ( z, η ) ∈ U × R d − . Furthermore, Assumption 2 showsthat L ( z, η ) has no eigenvalue on S for ( z, η ) ∈ U × R d − . This property dates back at least to [Kre68].However, central eigenvalues on S may occur for L when z belongs to S . The crucial point for provingLemma 3 is that Assumption 2 precludes central eigenvalues of M for all z ∈ U . Namely, for all z ∈ U and all η ∈ R d − , M ( z, η ) has no eigenvalue on S . This property holds because otherwise, for some( z, η ) ∈ U × R d − , there would exist a solution κ ∈ S to the dispersion relation p (cid:88) (cid:96) = − r z ∂ z a (cid:96) ( z, η ) κ (cid:96) = 0 . For convenience, the coordinates of η are denoted ( η , . . . , η d ). Using the definition (9) of a (cid:96) , and defining κ := ( κ , e i η , . . . , e i η d ), we have found a root z ∈ U to the relation s +1 (cid:88) σ =1 σ (cid:99) Q σ ( κ ) z σ − = 0 , (42)but this is not possible because the s + 1 roots (in z ) to the dispersion relation (8) are simple and belongto D . The Gauss-Lucas Theorem thus shows that the roots to the relation (42) belong to D (and thereforenot to U ).At this stage, we know that the eigenvalues of M ( z, η ), ( z, η ) ∈ U × R d − , split into two groups: thosein U , which we call the unstable ones, and those in D , which we call the stable ones. For ( z, η ) ∈ U × R d − ,25e then introduce the spectral projector Π s M ( z, η ), resp. Π u M ( z, η ), of M ( z, η ) on the generalized eigenspaceassociated with eigenvalues in D , resp. U . We can then integrate the first induction relation in (41) andget Π s M ( z k , η k ) W k = 1 z k ∂ z a p ( z k , η k ) (cid:88) j ≤ M ( z k , η k ) | j | Π s M ( z k , η k ) (cid:16) G kj , , . . . , (cid:17) T . The projector Π s M depends continuously on ( z, η ) ∈ U × R d − . Furthermore, since the spectrum of M does not meet S even for z ∈ S , there exists a constant C > δ ∈ (0 ,
1) that are independent of k ∈ N and such that ∀ j ≤ , (cid:107) M ( z k , η k ) | j | Π s M ( z k , η k ) (cid:107) ≤ C δ | j | . We thus get a uniform estimate with respect to k : | Π s M ( z k , η k ) W k | ≤ C (cid:88) j ≤ | G kj | . Passing to the limit and using (37), we get Π s M ( z, η ) W = 0, or in other words W = Π u M ( z, η ) W . • The sequence ( W j ) j ≤ satisfies both induction relations (40) and (41). Due to the form of thecompanion matrices L and M , see (38)-(39), we can conclude that the vector W belongs to the generalizedeigenspace (of either L or M ) associated with the common eigenvalues of M ( z, η ) and L ( z, η ). We havealready seen that M ( z, η ) has no eigenvalue on S and W = Π u M ( z, η ) W , so we can conclude that W belongs to the generalized eigenspace of L associated with those common eigenvalues of M ( z, η ) and L ( z, η ) in U .The matrix L ( z, η ) has N u eigenvalues in U , N s in D and N c on S . (Since z may belong to S , N c is not necessarily zero.) With obvious notations, we let Π u,s,c L ( z, η ) denote the correspondingspectral projectors of L for ( z, η ) sufficiently close to ( z, η ). In particular, the eigenvalues correspondingto Π u L ( z, η ) lie in U uniformly away from S for ( z, η ) sufficiently close to ( z, η ). We can then integratethe first induction relation in (40) and derive (for k sufficiently large):Π u L ( z k , η k ) W k = − a p ( z k , η k ) (cid:88) j ≥ L ( z k , η k ) − j − Π u L ( z k , η k ) (cid:16) F kj , , . . . , (cid:17) T . Using the uniform exponential decay of L ( z k , η k ) − j − Π u L ( z k , η k ) and (37), we finally end up withΠ u L ( z, η ) W = 0 . Since W belongs to the generalized eigenspace of L associated with those common eigenvalues of M ( z, η )and L ( z, η ) in U , we can conclude that W equals zero. Applying the induction relation (40), the wholesequence ( W j ) j ∈ Z is zero which yields the expected contradiction.The crucial property that we use in the proof of Lemma 3 is the fact that up to z ∈ S , the eigenvaluesof M ( z, η ) lie either in D or U . For the leap-frog scheme, this property would not be true if we hadimposed the auxiliary numerical boundary condition u n +2 j + u nj rather than 2 u n +2 j + λ a ( u n +1 j +1 − u n +1 j − ).Let us also observe that we have used the fact that a p and a − r are nonconstant in order to studythe induction relation (36). There might be some schemes for which a p and/or a − r are constant but forwhich one can still apply similar arguments as in the previous proof, even though (36) is an inductionrelation with fewer steps than (35). In this respect, Assumption 3 might be relaxed in specific applications.26 emark 2. The auxiliary problem (27) is in general not of the same form as (5) because in (27) onehas to impose infinitely many numerical boundary conditions. This is due to the fact that the stencilof M incorporates points “on the left” with respect to the first space variable. A remarkable exceptionoccurs for explicit schemes with s = 0 , for in that case the multiplier M v nj reads v n +1 j and (27) is exactlythe auxiliary problem considered (and labeled (2.7)) in [CG11] where one imposes Dirichlet boundaryconditions on finitely many boundary meshes (just use g nj = 0 for j ≤ − r ). The reader can then checkthat the energy-dissipation balance law of Proposition 2 in that case ( s = 0 , Q = I ) coincides exactlywith the algebra involved in the derivation of the estimate (2.12) in [CG11]. The reader can also checkthat for s = 0 , and Q = I , Lemma 3 becomes a rather trivial exercise...There still remains the problem of constructing a set of dissipative numerical boundary conditions ofthe same form as (5) with s ≥ , that is with finitely many numerical boundary conditions, and for whichone can prove by hand both a semigroup and a trace estimate as in Theorem 2. As explained in the introduction of Section 3, the linearity of (5) reduces the proof of Theorem 1 to thecase ( F nj ) = 0, ( g nj ) = 0, since we have already dealt with the case of zero initial data. We thus focus on(5) with ( F nj ) = 0 and ( g nj ) = 0, and write the corresponding solution ( u nj ) as u nj = v nj + w nj , where thesequence ( v nj ) solves: L v nj = 0 , j ∈ Z d , j ≥ , n ≥ ,M v nj = 0 , j ∈ Z d , j ≤ , n ≥ ,v nj = f nj , j ∈ Z d , n = 0 , . . . , s , (43)and ( w nj ) solves: L w nj = 0 , j ∈ Z d , j ≥ , n ≥ ,w n + s +1 j + s +1 (cid:88) σ =0 B j ,σ w n + σ ,j (cid:48) = ˜ g n + s +1 j , j ∈ Z d , j = 1 − r , . . . , , n ≥ ,w nj = 0 , j ∈ Z d , n = 0 , . . . , s . (44)For v nj + w nj to coincide with the solution ( u nj ) to (5), it is sufficient to extend the initial data f j , . . . , f sj by zero for j ≤ − r , which provides with the initial data in (43) on all Z d , and to define the boundarysource term in (44) by: ˜ g n + s +1 j := − v n + s +1 j − s +1 (cid:88) σ =0 B j ,σ v n + σ ,j (cid:48) . (45)We can estimate the solution ( v nj ) to (43) by applying Theorem 2. In particular, the trace estimate: (cid:88) n ≥ ∆ t e − γ n ∆ t P (cid:88) j =1 − r (cid:107) v nj , · (cid:107) (cid:96) ( Z d − ) ≤ C s (cid:88) σ =0 ||| f σ ||| − r , + ∞ , P = max( p , q + 1) gives (recall the definition (45) of ˜ g n + s +1 j ): (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j =1 − r (cid:107) ˜ g nj , · (cid:107) (cid:96) ( Z d − ) ≤ C (cid:88) n ≥ ∆ t e − γ n ∆ t max( p ,q +1) (cid:88) j =1 − r (cid:107) v nj , · (cid:107) (cid:96) ( Z d − ) ≤ C s (cid:88) σ =0 ||| f σ ||| − r , + ∞ . We can apply Theorem 1 to the solution ( w nj ) to (44) because the initial data in (44) vanish. We get:sup n ≥ e − γ n ∆ t ||| w n ||| − r , + ∞ + γγ ∆ t + 1 (cid:88) n ≥ ∆ t e − γ n ∆ t ||| w n ||| − r , + ∞ + (cid:88) n ≥ ∆ t e − γ n ∆ t p (cid:88) j =1 − r (cid:107) w nj , · (cid:107) (cid:96) ( Z d − ) ≤ C (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t (cid:88) j =1 − r (cid:107) ˜ g nj , · (cid:107) (cid:96) ( Z d − ) ≤ C s (cid:88) σ =0 ||| f σ ||| − r , + ∞ . Combining with the similar estimate provided by Theorem 2 for ( v nj ), we end up with the expectedestimate:sup n ≥ e − γ n ∆ t ||| u n ||| − r , + ∞ + γγ ∆ t + 1 (cid:88) n ≥ ∆ t e − γ n ∆ t ||| u n ||| − r , + ∞ + (cid:88) n ≥ ∆ t e − γ n ∆ t p (cid:88) j =1 − r (cid:107) u nj , · (cid:107) (cid:96) ( Z d − ) ≤ C s (cid:88) σ =0 ||| f σ ||| − r , + ∞ , which completes the proof of Theorem 1. Let us first observe that in [Wad90],
Wade has constructed symmetrizers for deriving stability esti-mates for multistep schemes, even in the case of variable coefficients. His conditions for constructing asymmetrizer are less restrictive than Assumption 2. However, the symmetrizer in [Wad90] is genuinelynonlocal and it is therefore not clear that it may be useful for boundary value problems. The main nov-elty 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 . Forthe approximation of parabolic equations, this term may give some extra dissipation, but a crucial pointto keep in mind is that the coefficients of the numerical scheme are assumed to be constant (which mayin turn yield rather severe CFL conditions for implicit approximations of parabolic equations). Hence itdoes not seem very clear that our approach will yield stability estimates with “optimal” CFL conditionswhen 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 are simple. Here we have assumed that all the roots, including those in D are simple. If28e could manage to deal with multiple roots in D , then Theorem 1 would be applicable to numericalapproximations of the transport equation (12) that are based on Adams-Bashforth methods of order 3 orhigher (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 strongstability, in the sense of Definition 1, implies semigroup stability. However, an even stronger assumptionwas made in [KW93], namely that the sole fulfillment of the interior estimate γγ ∆ t + 1 (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t ||| u n ||| − r , + ∞ ≤ C γ ∆ t + 1 γ (cid:88) n ≥ s +1 ∆ t e − γ n ∆ t ||| F n ||| , + ∞ , when both the initial and boundary data for (5) vanish, does imply semigroup stability. The analogousconjecture for partial differential equations seems to be still open so far, but we do hope that our multipliertechnique may yield some insight for dealing with the strong form of the conjecture in [KW93]. Acknowledgments
This article was completed while the author was visiting the Institute of Math-ematical Science at Nanjing University. The author warmly thanks Professor Yin Huicheng and theInstitute for their hospitality during this visit.
References [AG76] S. Abarbanel and D. Gottlieb. A note on the leap-frog scheme in two and three space dimensions.
J. Computational Phys. , 21(3):351–355, 1976.[AG79] S. Abarbanel and D. Gottlieb. Stability of two-dimensional initial boundary value problemsusing leap-frog type schemes.
Math. Comp. , 33(148):1145–1155, 1979.[BGS07] S. Benzoni-Gavage and D. Serre.
Multidimensional hyperbolic partial differential equations .Oxford University Press, 2007. First-order systems and applications.[CG11] J.-F. Coulombel and A. Gloria. Semigroup stability of finite difference schemes for multidimen-sional hyperbolic initial boundary value problems.
Math. Comp. , 80(273):165–203, 2011.[Cou09] J.-F. Coulombel. Stability of finite difference schemes for hyperbolic initial boundary valueproblems.
SIAM J. Numer. Anal. , 47(4):2844–2871, 2009.[Cou13] J.-F. Coulombel. Stability of finite difference schemes for hyperbolic initial boundary valueproblems. In
HCDTE Lecture Notes. Part I. Nonlinear Hyperbolic PDEs, Dispersive and Trans-port Equations , pages 97–225. American Institute of Mathematical Sciences, 2013.[Cou14] J.-F. Coulombel. Fully discrete hyperbolic initial boundary value problems with nonzero initialdata.
Preprint , http://arxiv.org/abs/1412.0851 , 2014.[G˚ar56] L. G˚arding. Solution directe du probl`eme de Cauchy pour les ´equations hyperboliques. In Lath´eorie des ´equations aux d´eriv´ees partielles , Colloques Internationaux du C. N. R. S., pages71–90. C. N. R. S., Paris, 1956.[GKO95] B. Gustafsson, H.-O. Kreiss, and J. Oliger.
Time dependent problems and difference methods .John Wiley & Sons, 1995. 29GKS72] B. Gustafsson, H.-O. Kreiss, and A. Sundstr¨om. Stability theory of difference approximationsfor mixed initial boundary value problems. II.
Math. Comp. , 26(119):649–686, 1972.[GT81] M. Goldberg and E. Tadmor. Scheme-independent stability criteria for difference approxima-tions of hyperbolic initial-boundary value problems. II.
Math. Comp. , 36(154):603–626, 1981.[HNW93] E. Hairer, S. P. Nørsett, and G. Wanner.
Solving ordinary differential equations. I . Springer-Verlag, second edition, 1993. Nonstiff problems.[HW96] E. Hairer and G. Wanner.
Solving ordinary differential equations. II . Springer-Verlag, secondedition, 1996. Stiff and differential-algebraic problems.[Kre68] H.-O. Kreiss. Stability theory for difference approximations of mixed initial boundary valueproblems. I.
Math. Comp. , 22:703–714, 1968.[KW93] H.-O. Kreiss and L. Wu. On the stability definition of difference approximations for the initial-boundary value problem.
Appl. Numer. Math. , 12(1-3):213–227, 1993.[Ler53] J. Leray.
Hyperbolic differential equations . The Institute for Advanced Study, Princeton, N. J.,1953.[Mic83] D. Michelson. Stability theory of difference approximations for multidimensional initial-boundary value problems.
Math. Comp. , 40(161):1–45, 1983.[Oli74] J. Oliger. Fourth order difference methods for the initial boundary-value problem for hyperbolicequations.
Math. Comp. , 28:15–25, 1974.[Rau72] J. Rauch. L is a continuable initial condition for Kreiss’ mixed problems. Comm. Pure Appl.Math. , 25:265–285, 1972.[RM67] R. D. Richtmyer and K. W. Morton.
Difference methods for initial value problems . Grad-uate Texts in Mathematics. Interscience Publishers John Wiley & Sons, 1967. Theory andapplications.[Slo83] D. M. Sloan. Boundary conditions for a fourth order hyperbolic difference scheme.
Math.Comp. , 41:1–11, 1983.[SW97] J. C. Strikwerda and B. A. Wade. A survey of the Kreiss matrix theorem for power boundedfamilies of matrices and its extensions. In
Linear operators (Warsaw, 1994) , volume 38 of
Banach Center Publ. , pages 339–360. Polish Acad. Sci., 1997.[TE05] L. N. Trefethen and M. Embree.
Spectra and pseudospectra . Princeton University Press, 2005.The behavior of nonnormal matrices and operators.[Tho72] J. M. Thomas. Discr´etisation des conditions aux limites dans les sch´emas saute-mouton.
Rev.Fran¸caise Automat. Informat. Recherche Op´erationnelle , 6(Ser. R-2):31–44, 1972.[Tre84] L. N. Trefethen. Instability of difference models for hyperbolic initial boundary value problems.
Comm. Pure Appl. Math. , 37:329–367, 1984.[Wad90] B. A. Wade. Symmetrizable finite difference operators.
Math. Comp. , 54(190):525–543, 1990.30Wu95] L. Wu. The semigroup stability of the difference approximations for initial-boundary valueproblems.