Compactness estimates for difference schemes for conservation laws with discontinuous flux

We establish quantitative compactness estimates for finite difference schemes used to solve nonlinear conservation laws. These equations involve a flux function $f(k(x,t),u)$, where the coefficient $k(x,t$ is $BV$-regular and may exhibit discontinuities along curves in the $(x,t)$ plane. Our approach, which is technically elementary, relies on a discrete interaction estimate and the existence of one entropy function. While the details are specifically outlined for the Lax-Friedrichs scheme, the same framework can be applied to other difference schemes. Notably, our compactness estimates are new even in the homogeneous case ($k\equiv 1$).


Introduction
The main part of this paper investigates a finite difference algorithm as it applies to the Cauchy problem for scalar conservation laws with the form (1.1) where (x, t) ∈ R × R + ; u(x, t) is the scalar unknown function; and u 0 (x), k(x, t), f (k, u) are given functions to be detailed later.Here it suffices to say that for the compactness estimates we need k(x, t) ∈ BV (R × R + ), u → f (k(x, t), u) genuinely nonlinear, and u 0 (x) bounded (see Section 2 for the complete list of assumptions).
The special feature of (1.1) is the nonlinear flux function f (k(x, t), u) that depends explicitly on the spatial and temporal variables through a discontinuous coefficient k(x, t).Conservation laws with discontinuous flux functions are encountered in various applications.For example, in oil reservoirs, rock permeability may vary significantly in different locations, resulting in a discontinuous flux function (see, e.g., [11]).Similarly, in traffic flow, abrupt changes in vehicle density may occur at bottlenecks, which leads to a discontinuous flux (see, e.g., [8]).
The study of conservation laws with discontinuous flux has been a heavily investigated area for the past three decades.This is partly due to its links to various applications, but also because it possesses several non-trivial mathematical properties.These properties include the existence of several L 1 stable semigroups, which are based on different entropy conditions, as well as the lack of uniform bounds on the total variation (i.e., BV estimates).As a result, it constitutes a non-trivial generalization of scalar conservation laws [19].While a comprehensive review of the extensive literature is beyond the scope of this paper, we provide a few select references for interested readers and direct them to the reference lists in these papers [1,2,3,5,7,18,16,17].
Proving the existence of solutions for conservation laws is associated with establishing strong a priori estimates for approximate solutions.Classically, this involves bounding the total variation of the approximate solutions, independent of the approximation parameter.However, bounding the total variation when the flux is discontinuous is in general impossible.To address this issue, several alternative convergence approaches have been applied over the years, including singular mapping techniques as well as compensated compactness and other advanced weak convergence methods [9,13,16,17,18,21] (this is just a few examples).Specifically, the weak convergence methods (see, e.g., [9,21]) are profound and involve a significant amount of functional analysis, making them technically challenging to comprehend.
The aim of this paper is to derive quantitative L 1 translation estimates that can be utilized for various applications, such as demonstrating the convergence of finite difference approximations.Let us consider a time step ∆t and a grid size parameter ∆x, which collectively form a parameter pair ∆ = (∆x, ∆t).We use the notation u ∆ ∆>0 to represent the approximate solutions.The translation estimates, which maintain uniformity across all ∆, are defined as follows: (1.2) This estimate holds true for any temporal (τ > 0) and spatial (h > 0) translations.Here, µ t , µ x are parameters in the interval (0, 1) that quantify the "degree of compactness", where µ t = µ x = 1 corresponds to uniformly bounded total variation.In (1.2), χ is a weight function that can be employed to control the growth of solutions as they approach infinity in x.An example of a weight function is χ(x) = (1 + |x| 2 ) −N , where N > 1/2, see (2.15) and also [15].If supp u ∆ is contained within a compact set [0, T ] × [−R, R] ⊂ R 2 , independent of ∆, we may set χ ≡ 1, which is what we do for the rest of the introduction!
We refer to translation estimates like (1.2) as quantitative compactness estimates.They can be used to derive convergence results via the well-known Kolmogorov-Riesz-Fréchet characterization of precompact subsets of L 1 in terms of the uniform continuity of translations in L 1 , see, e.g., [6,Theorem 4.26].Beyond implying convergence, quantitative compactness estimates can be used to derive continuous dependence and error estimates [20,4].
Our approach is technically elementary and relies on discrete interaction estimates and the existence of one uniformly convex entropy.We draw inspiration from previous work of Golse and Perthame [12], who developed a quantitative compensated compactness framework for establishing Besov space regularity of solutions to homogenous conservation laws.We adapt their approach and apply it to establish quantitative compactness estimates for sequences of approximate solutions.A related approach has recently been employed in the study of vanishing viscosity approximations of stochastic conservation laws, as described in [15].While the details are specifically outlined for the Lax-Friedrichs scheme, the same framework can be applied to other difference schemes as well.The obtained estimates are new also in the homogenous case k ≡ 1 (scalar conservation law).
To clarify further, let us consider the Lax-Friedrichs scheme for the Cauchy problem (1.1), which can be represented by the following equation [17]: Here, U n j approximates the exact solution u at the grid point (j∆x, n∆t), where n and j are integers, λ = ∆t/∆x, and f n j = f k n j , U n j .The temporal and spatial discretization parameters are linked through a CFL condition (2.8), so that ∆t ∼ ∆x.
Assuming reasonable conditions on f , k ∈ BV , and u 0 ∈ L ∞ , we can establish the following two a priori estimates (see next section): These estimates are the only bounds that remain robust with respect to the regularity of the coefficient k(x, t).If k ≡ 1, the approximations are of bounded total variation [14].The second bound is referred to as a dissipation (or entropy stability) estimate, which is known to hold for many numerical schemes for conservation laws [10,22].These estimates are sometimes referred to as weak H 1 estimates since they suggest that ∂ x u ∆ 2 dx dt 1/ |∆|.The estimate for the Lax-Friedrichs scheme (1.3) with variable and discontinuous k was derived in [17].In this paper, we present a slightly generalized variant of the estimate.
Given the genuinely nonlinear flux f (k, u), in the sense of (2.7), consider an entropy/entropy flux pair S(k, u), Q(k, u) .Set S n j = S k n j , U n j and Q n j = Q k n j , U n j .The Lax-Friedrichs scheme satisfies the following entropy balance (to be derived later): where, by (1.4) and k ∈ BV , we have ∆x n j Ψ n j 1.Note that the entropy production Ψ n j is a signed measure (even if S = S(u) is convex), which is distinct from the homogenous conservation law case where it is negative.
Let ν be an arbitrary nonnegative integer and introduce the "spatial difference" quantities The appearance of 2ν, as opposed to simply ν, in the above equations is due to that we use a staggered mesh, where the various grid functions are only defined at grid points where j + n is an even integer.Then it is straightforward to deduce from (1.3) and (1.5) that the following 2 × 2 system of finite difference equations hold: where C n A,j ≡ 0 and C n D,j = Ψ n j+2ν − Ψ n j .We are also going to need the following more regular quantities obtained by applying "inverse-difference" operators to A n j , D n j : Indeed, one can prove that sup n,j recalling that ν is a nonnegative integer.
Next, we introduce the interaction functional which is a measure of future potential interaction at time level n of the finite difference solutions.This is a discrete version of a functional referred to as the Varadhan functional by Tartar [23, p. 182] and Golse & Perthame [12].
In this paper, we establish the following discrete interaction identity for the 2 × 2 system (1.6): (1.7) 1 2 ∆t∆x To obtain an evolution difference equation satisfied by the interaction potential I n , we compute the temporal difference . This identity is then multiplied by ∆t and summed over n, resulting in the equation given by (1.7).The derivation of (1.7) employs a discrete chain rule, as well as the difference equations for A n j and D n j .A more complicated form of (1.7) holds when there is a weight function χ present.The interaction identity (1.7) can be interpreted as a discrete form of the interaction identity (7) introduced in [12], as well as the stochastic interaction identity (3.12) in [15] (with σ = 0).It should be noted that the term E appearing on the right-hand side of (4.12) is solely a consequence of the discretization process and does not have a corresponding counterpart in the interaction identities of [12,15].For the detailed form of E, see Lemma 4.1.
The terms on the right-hand side of (1.7) can all be bounded by a ∆-independent constant C times 2ν∆x.The next step is to convert (1.7) into a quantitative compactness estimate, which requires exploiting the genuine nonlinearity of f (k, u) in u.The precise assumption can be found in (2.7), and a corresponding assumption for homogeneous conservation laws can be found in [12,Section 5].A relevant special case occurs if u → f (k, u) is uniformly convex and S = f .Then and therefore A similar estimate can be deduced for the temporal differences.Along with the compact support assumption (or if we use a weight function) and Hölder's inequality, this translates into the L 1 translation estimate given by (1.2) with h = τ = 2ν∆x and µ t = µ x = 1/4, which is the main result of this paper.In summary, the above outline provides an overview of the quantitative compactness approach.To delve deeper, we will present detailed proofs incorporating a weight function in the following sections.We demonstrate the effectiveness of the quantitative compactness approach on the Lax-Friedrichs scheme.However, with slight modifications to the discrete interaction identity given by (1.7), the same approach can be applied to other classical schemes, provided they are uniformly bounded in L ∞ and satisfy a weak H 1 (dissipation) estimate, that is to say, they satisfy modified forms of the a priori estimates described in (1.4).This allows for elementary convergence proofs, especially in cases where obtaining total variation estimates is challenging.
Our results provide an existence theorem for the Cauchy problem (1.1).The uniqueness question for conservation laws with discontinuous flux is another problem entirely, even in the most basic setting where the flux has a single spatial discontinuity and no temporal discontinuity (the so-called two-flux problem).For example, it turns out that the two-flux problem generally has infinitely many definitions of entropy solution, each one of which generates its own distinct L 1 contraction semigroup [1].In [17] we used results from [16] to prove a uniqueness result applicable to the Lax-Friedrichs scheme of this paper in the special case where k is piecewise Lipschitz continuous, meaning that all of the discontinuities of k occur along Lipschitz continuous curves in the (x, t)plane.This uniqueness result was proven under the additional assumption that all of the flux discontinuities satisfy a certain "crossing condition".At least for the simple two-flux version of the problem (1.1), the solution generated by the Lax-Friedrichs scheme of this paper corresponds to the so-called vanishing viscosity solution [2].In the absence of the crossing condition it is not known whether the (subsequential) limit of the Lax-Friedrichs scheme is the vanishing viscosity solution.Finally, given (1.8), it is worth noting that the solution u derived as the limit of the Lax-Friedrichs scheme, exhibits Besov space regularity as given by This regularity aligns with the known regularization effect established in [12,Theorem 5.1] for homogeneous equations with a single convex entropy.
The structure of this paper is as follows: Section 2 presents the assumptions related to the data of the problem and precisely defines the Lax-Friedrich scheme.This section also introduces our main result.In Section 3, we establish preliminary L ∞ and weak H 1 estimates.Section 4 proves the spatial translation estimate, while Section 5 details the temporal estimate.By bringing together the spatial and temporal estimates, we provide the proof of our main result.

Lax-Friedrichs scheme and main result
We begin by listing some assumptions on u 0 , k, f which will be needed.Regarding the initial function we assume For the discontinuous coefficient k : R × R + → R we assume that

Regarding the flux function
We need also an assumption on f that guarantees that the Lax-Friedrichs approximations stay uniformly bounded.For example, we can require which in fact implies that the interval [a, b] becomes an invariant region.We use the notation ∂ k G and ∂ u G to denote the first order partial derivatives of G(•, •) with respect to the first and second variables.Letting U := [α, β] × [a, b], we will use the following abbreviations: We are given functions S and Q (referred to earlier) that are assumed to form an entropy/entropy flux pair (S(k, u), Q(k, u)), meaning that We assume that We make the following genuine nonlinearity assumptions about f and S.
Next we describe the Lax-Friedrichs scheme.Let ∆x > 0 and ∆t > 0 denote the spatial and temporal discretization parameters, which are chosen so that they always obey the CFL condition Here κ is a positive parameter which we can choose to be very small so that the allowable time step is reduced only negligibly.We will work under the standing assumption that the space step ∆x and the time step ∆t are comparable, i.e., there are constants c 1 , c 2 > 0 such that c 1 ≤ ∆t ∆x ≤ c 2 .Therefore, when we declare that a constant C is independent of ∆x (or ∆t), it implies that C is also independent of ∆ = (∆x, ∆t).
The time domain [0, ∞) is discretized via t n = n∆t for n ∈ Z 0 + := {0, 1, . ..} (Z + := {1, 2, . ..}), resulting in time strips [t n , t n+1 ).The spatial domain R is divided into cells [x j−1 , x j+1 ) with centers at the points x j = j∆x for j ∈ Z.Let ρ j (x) be the characteristic function for the interval [x j−1 , x j+1 ) and ρ n j the characteristic function for the rectangle [x j−1 , x j+1 ) × [t n , t n+1 ).The finite difference scheme then generates, for each mesh size ∆ = (∆x, ∆t), with ∆x and ∆t taking values in sequences tending to zero, a piecewise constant approximation (2.9) where the values U n j : (j, n) ∈ Z × Z 0 + , j + n = even remain to be defined.We define U 0 j : j = even by (2.10) with Riemann initial data where the coefficient k(x, t) has been discretized via the piecewise constant approximation (2.12) Recall that k(x, t) is a BV function.For every t, the function k(t, •) can be considered as a precise representative, being defined everywhere and normalized through right-continuity.Consequently, we are justified in setting k n j = k(x j , tn ) in (2.12), where tn is an arbitrary point lying in the interval [t n , t n+1 ) (for example, tn = t n ).We then define Integrating the weak formulation of (2.11) over the control volume [x j−1 , x j+1 ) × (0, ∆t) gives After a direct evaluation of the integrals for ∆t small, we obtain the staggered Lax-Friedrichs scheme Notice that in this paper we restrict our attention to the sublattice which means that U 0 j : j = even , U 1 j : j = odd , U 2 j : j = even etc.are calculated.The following abbreviations can be used to shorten some of the expressions that will arise.For fixed n ∈ {0, 1, . . ., N }, Note that U n j , U n j±2 , . . .and U n+1 j±1 , U n+1 j±1 , . . .are defined for j ∈ Ω n , while U n j±1 , U n j±3 , . . .and U n+1 j , U n+1 j±2 , . . .are defined for j ∈ Ω ′ n .We work with the so-called "even" sublattice described above in the interest of conceptual simplicity.Note that when one applies the Lax-Friedrichs scheme on the standard lattice, what is generated is two uncoupled numerical solutions, one solution on the even sublattice, and one solution on the odd sublattice.Our analysis for the even sublattice would then apply to each of those solutions separately.
The validation of (2.14) is directly derived from the results in Section 4 (Proposition 4.5) and Section 5 (Proposition 5.7).
Remark 2.2.Assuming that the entropy u → S(k, u) in (2.7) is uniformly convex, (2.16) If we choose S = f as the entropy, then µ = 1 4 , see also Remark 3.5.Remark 2.3.In this paper, we employ estimates denoted as "a b", signifying that there exists a constant C such that a ≤ Cb.Notably, C may depend on the specific constants associated with the assumptions of the problem.However, C does not depend on the grid parameters ∆.
We will conclude this section by outlining the key concepts underpinning the proof of the spatial part of Theorem 2.1.Our discussion will be focused on a homogenous equation, augmented with artificial viscosity, to provide a clear outline of the underlying ideas.
For any fixed ε > 0, let u ε ∈ C 2 satisfy the equations where f, η ∈ C 2 are (for example) uniformly convex and To prevent ambiguity and ensure clarity, we denote the entropy/entropy flux pair as (η, q), distinguishing them from the functions (S, Q) used for (1.1).We denote by ∆ h W (x, t) := W (t, x + h) − W (x, t) the spatial difference operator with step size h.Set Then the system (2.17) takes the form where, applying Lemma 3.4 and Remark 3.5 with S = η and Q We need the (spatial) anti-derivatives of a and d: A(t, y) := Denote by I(t) the (spatial) interaction functional: A straightforward calculation will confirm that the following (spatial) interaction identify holds: Hence, by integrating the above identify in t ∈ [0, T ], T > 0, we arrive at where Utilizing the assumed bounds on u ε , we can estimate these three terms following the approach outlined in [15].The final result is that |J i | h for i = 1, 2, 3, which implies that A similar temporal translation estimate can be derived for u ε .This brings us to the end of the outline detailing the proof for the vanishing viscosity approximation u ε .

Preliminary results
The following result is taken from [17, Lemma 4.1].
In [17] we employed the single entropy S(u) = u 2 /2.We wish to allow for an entropy of the form S(k, u).The first part (inequality (3.1)) of the lemma that follows is a generalization of [17,Lemma 4.3], which we prove by modifying the proof in [17].The second part of the lemma (inequality (3.2)) is new, and is required for the analysis that follows.
Lemma 3.2.Let (S, Q) be an arbitrary entropy/entropy flux pair defined by (2.5), satisfying the uniform convexity condition (2.16).With k n j−1 , U n j−1 and k n j+1 , U n j+1 given, compute U n+1 j by (2.13).Define Then where K 1 is a constant that is independent of ∆x.
In addition, for some constants K 2 , K 3 , and K 4 that are independent of ∆x, The estimate (3.2) holds without the convexity condition (2.16).
Proof.Let us introduce the functions w, v, Φ : [a, b] → R defined by We collect the following elementary facts about these functions in one place before continuing with the proof: We seek an estimate of Φ ′ (s).It is readily verified that where This yields As a consequence of the CFL condition (2.8), F ≥ κ, and Also, G ≥ γ and |G| ≤ ∂ 2 uu S .Thus, Integrating this last inequality from 0 to 1 gives Combining this with the fact that we obtain , the proof of (3.1) will be complete as soon as we show that |Φ(0)| is bounded by a constant times which yields the inequality Recalling that 2), the inequality (3.4) provides the desired bound on |Φ(0)| and the proof of (3.1) is complete.
The proof of (3.2) will be complete if we can obtain a suitable estimate of |Φ(1)|.We have Recalling which completes the proof of (3.

2).
In what follows, we consider an arbitrary weight function χ ∈ C 1 (R) ∩ L 1 (R) satisfying the properties in (2.15).We will use the following abbreviations: Also, the following facts will be helpful.For a fixed integer i, if {Z n j } is bounded, then The following result is the version of [17, Lemma 4.3] that is appropriate for the assumptions of this paper.Lemma 3.3.For T > 0, N = ⌊T /∆t⌋, we have the bounds where C 1 (T ) and C 2 (T ) are independent of ∆.
Proof.For the first inequality of (3.5), we use (3.1) with S = u 2 /2: For the first sum on the right side of (3.6) we use the identity Summing over n and j, the terms −χ j S n+1 j The desired bound for the first sum on the right side of (3.6) then follows from (3.7), using (2.15).For the second sum on the right side of (3.6) we use the identity . Summing over n and j, the contribution from the first part of the right side of (3.8) telescopes, resulting in (3.9) ∆x The desired bound for the second sum then follows from (3.9), using (2.15).The third and fourth sums on the right side of (3.6) are bounded in absolute value due to the fact that χ is bounded and k ∈ BV (R × R + ), and the proof of the first part of (3.5) is complete.
For the second part of (3.5), We prove a bound for the first sum on the right side of (3.10).A bound for the second sum is proven in a similar manner.To this end, (3.11) Here we have used the inequality (a + b) 2 ≤ 2(a 2 + b 2 ).Recalling (2.13), we can replace (3.11) by We use along with an application of Young's inequality, ab ≤ (a 2 + b 2 )/2, to obtain where d 1 and d 2 are constants that are independent of ∆x.
The proof of the second part of (3.5) is completed by substituting (3.12) into (3.10),along with a similar inequality for the second sum on the right side of (3.10).The desired inequality then follows from the first part of (3.5), the fact that k ∈ BV (R × R + ), and which follows from (2.15).
The following lemma is a mild adaptation of [12,Lemma 5.2].It will be used in the upcoming two sections.Lemma 3.4.Suppose the flux f satisfies (2.3), (2.4), and the genuine nonlinearity assumption in (2.7).Suppose also that we are given an entropy/entropy flux pair (S, Q) that satisfies (2.5), (2.6), and (2.7).Then for all v, w ∈ [a, b] and k ∈ [α, β], Proof.A straightforward calculation using (2.5) yields By combining (3.14) with the version that results by swapping ξ and ζ, we obtain Recalling (2.7), for (ξ, ζ) ∈ R 2 we have the inequality By combining (3.15), (3.16) and (3.17), we obtain This completes the proof for the case where v ≤ w.We obtain the desired inequality for w ≤ v by swapping the roles of v and w.Finally, combining the two cases yields (3.13).

Estimate of spatial translates
This section focuses on stating and proving Proposition 4.5, crucial for deriving Theorem 2.1 as mentioned earlier.To facilitate this, we will first lay the groundwork with a series of technical lemmas.A condensed version of the proof, particularly addressing the viscosity approximation, is presented at the end of Section 2 to guide the upcoming calculations.
For now we assume that h = 2ν, where ν is a nonnegative integer.In what follows ∆ h and ∆ −h denote the spatial finite difference operators, We will also use the difference operators ∆ u h and ∆ h h defined by The immediate goal is to prove Lemma 4.4 below.We will use the following key fact, which is the content of Lemma 3.4: We write the Lax-Friedrichs scheme and entropy inequality in the form (4.5) where, by (3.2) of Lemma 3.2, With the assumption that k ∈ BV (R × R + ), along with (3.5) of Lemma 3.
where C 3 (T ) is a constant that is independent of ∆x.We multiply both equations of (4.5) by χ j ∆ h , and after some algebra, arrive at the following 2 × 2 discrete system: where

Define
A n ℓ = ∆x Ωn:j≤ℓ where The following identities hold: Ωn:ℓ≥j where we are using the notational convention Ωn:ℓ≥j P j,ℓ = {(j,ℓ)∈Ωn×Ωn:ℓ≥j} with a similar definition for sums where Ω n is replaced by Ω ′ n .The interaction identity (4.12) below can be viewed as a discrete, Lax-Friedrichs version of the interaction identity (7) in [12], and also the stochastic interaction identity (3.12) of [15] (with σ = 0).The term E appearing on the right side of (4.12) is purely an artifact of discretization, and does not have a counterpart in the interaction identities of [12] and [15].
Then the following interaction identity holds: where Proof.Using the identity in combination with (4.11), yields (4.13) Substituting the identity After substituting (4.15) into (4.13) the result is Due to the identities and the identities (4.10), we can express (4.16) in the form (4.17) Next we use the identities ) Substituting (4.18), (4.19), and (4.20) into (4.17),we obtain The proof is completed by summing over n, recalling λ = ∆t/∆x, and then solving for the sum containing A n j E n j − D n j B n j .
Lemma 4.2.Assuming that h = 2ν∆x, with ν a nonnegative integer, the following estimates hold: Proof.We prove that |A n ℓ | h.The proof that D n j h is similar.Recalling (4.9), This yields via (2.15) Proof.For the proof of (4.22), Here we have used the Cauchy-Schwarz inequality, (3.5) of Lemma 3.3, and ∆t = λ∆x.For the proof of (4.23), we repeat the calculation in (4.24),with G n j replacing U n j , resulting in Summing over n and j, invoking (4.22) and the fact that k ∈ BV (R × R + ) then yields (4.23) Lemma 4.4.Let h = 2ν∆x, where ν is a nonnegative integer.We have the following estimates for the spatial translates: Under the assumptions of Remark 3.5, we obtain Proof.We start with the proof of (4.25).Recalling (4.8), (4.28) In the last step we used that fact that k ∈ BV (R × R + ).Similarly, using h.
We are now in a position to prove the main result of this section, specifically, the spatial part of Theorem 2.1.Proposition 4.5.For t ∈ [0, T ] and h > 0, the spatial translates satisfy the following estimate: Proof.In the special case where h = 2ν∆x, with ν a nonnegative integer, The first term on the right side of (4.32) is h µ , according to Lemma 4.4.The second term is O(∆x), and therefore also h µ .For more general h > 0, we write h = 2ν∆x + 2α∆x, where ν is a nonnegative integer and α ∈ (0, 1).Then  The first integral on the right side of (4.33) is h µ , according to the previous part of the proof.The second integral is equal to  In the last step we used the fact that u ∆ is piecewise constant, and that Using the estimate from the first part of the proof, the sum on the right side of (4.34) is α∆x µ , and therefore also h µ .
Applying the Cauchy-Schwarz inequality to the right side of (5.17 .Recalling (3.5) we have (5.15).For the proof of (5.16), we use (5.18) The proof of (5.16) is completed by estimating each of the sums on the right side of (5.18).The first sum is estimated using the Cauchy-Schwarz inequality as in the proof of (5.15).For the second sum we use the fact that k ∈ BV (R × R + ).
Lemma 5.5.Assuming that θ is an even nonnegative integer, we have the following estimate for the temporal translates: In view of the interaction identity (5.13), (5.14), the proof then reduces to proving that each term on the right side of (5.13) is τ .
Estimate of S3 .For the estimate of S3 , we have via (5.1) and (2.15) Cn A,j ∆xχ j .Combining this with (5.9), we obtain To estimate R4 , we use the fact that Ẽn j−1 is bounded, along with the estimate Cn A,j ∆xχ j , from which it is clear that (5.22) R4 ∆t.