Diffusive-Dispersive Approximation

# Convergence of Fully discrete schemes for diffusive dispersive conservation laws with discontinuous coefficient

Rajib Dutta
Institut für Mathematik,
Julius-Maximilians-Universität Würzburg,
Campus Hubland Nord, Emil-Fischer-Strasse 30,
97074, Würzburg, Germany.
Ujjwal Koley
Centre For Applicable Mathematics (CAM)
Tata Institute of Fundamental Research
P.O. Box 6503, GKVK post office
Bangalore-560065, India.
and  Deep Ray
Centre For Applicable Mathematics (CAM)
Tata Institute of Fundamental Research
P.O. Box 6503, GKVK post office
Bangalore-560065, India.
August 19, 2019
###### Abstract.

We are concerned with fully-discrete schemes for the numerical approximation of diffusive-dispersive hyperbolic conservation laws with a discontinuous flux function in one-space dimension. More precisely, we show the convergence of approximate solutions, generated by the scheme corresponding to vanishing diffusive-dispersive scalar conservation laws with a discontinuous coefficient, to the corresponding scalar conservation law with discontinuous coefficient. Finally, the convergence is illustrated by several examples. In particular, it is delineated that the limiting solutions generated by the scheme need not coincide, depending on the relation between diffusion and the dispersion coefficients, with the classical Kružkov-Oleĭnik entropy solutions, but contain nonclassical undercompressive shock waves.

###### Key words and phrases:
conservation laws; discontinuous flux; diffusive-dispersive approximation; finite difference scheme; convergence; entropy condition; nonclassical shock

]rajib.ami@gmail.com ]ujjwal@math.tifrbng.res.in ]deep@math.tifrbng.res.in

## 1. Introduction

In this paper, we consider a finite difference method for the vanishing diffusive-dispersive approximations of scalar conservation laws with a discontinuous flux

 (1.1) {uεt+f(k(x),uε)x=R[ε,μ(ε);uε],  x∈R×(0,T),uε(x,0)=u0(x),  x∈R,

when tends to zero with as . Here is a regularization term, depends upon two parameters and referred to as the diffusion and the dispersion coefficients, motivated by the equations of two-phase flow in porous media, is fixed, is the unknown scalar map, the initial data, is a spatially varying (discontinuous) coefficient, and the flux function is a sufficiently smooth scalar function (see Section  2 for the complete list of assumptions).

Motivated by the dynamic capillary pressure [7], we consider in this paper the simplified model

 (1.2) R[ε,μ(ε);uε]=εβuεxx+μ(ε)γuεxxt

with a third-order mixed derivatives term including one time derivative. Here are fixed parameters. The equation (1.1) along with (1.2) serves as a concrete model of two phase flows in a heterogeneous porous medium.

Moreover, drawing preliminary motivation from phase transition dynamics, we also consider the following specific form of the regularization term

 (1.3) R[ε,μ(ε);uε]=εβuεxx+μ(ε)γuεxxx,

with fixed.

Furthermore, for the simplicity in the exposition, we assume that the flux function has the following particular form

 f(k(x),u)=k(x)f(u).

Note that the flux has a possibly discontinuous spatial dependence through the coefficient , which is allowed to have jump discontinuities.

The scalar conservation laws with a discontinuous flux function

 (1.4) ut+f(k(x),u)x=0

is a special example of this type of problems, corresponds to the case . A simple physical model corresponding to (1.4) is the Witham model of car traffic flow on a highway (consult the monograph by Leveque [21]), where the spatially varying coefficient corresponds to changing road conditions. Several other models such as two phase flow in a heterogeneous porous medium that arise in petroleum industry, the modeling of the clarifier thickener unit used in waste water treatment plants are also corresponding to (1.4).

Independently of the smoothness of the initial data and , solutions to (1.4) are not necessarily smooth due to the presence of nonlinear flux term in the equation (1.4). Thus, weak solutions must be sought.

###### Definition 1.1 (Weak solution).

A weak solution of the initial value problem (1.4) is a bounded measurable function satisfying

 (1.5) ∫R∫T0(φtu+φxk(x)f(u))dxdt+∫Rφ(x,0)u0(x)dx=0,

for all .

It is well known that (weak) solutions may be discontinuous and they are not uniquely determined by their initial data. Consequently, an entropy condition must be imposed to single out the physically correct solution. If is “smooth”, a weak solution satisfies the entropy condition if for all convex functions

 η(u)t+(k(x)Q(u))x+k′(x)(η′(u)f(u)−Q(u))≤0,inD(R×[0,T]),

where is defined by .

By standard limiting argument, this implies the Kružkov-type entropy condition

 |u−c|t+sign(u−c)(k(x)(f(u)−f(c)))x+sign(u−c)f(c)k′(x)≤0,inD(R×[0,T]),

holds for all .

However, the notion of entropy solution described above breaks down when is discontinuous. In view of [14], we use the following notion of entropy solution for (one-dimensional) conservation laws with discontinuous flux equations with coefficients that are only spatially dependent. We assume that the spatially varying coefficient is piecewise with finitely many jumps (in and ), located at .

###### Definition 1.2 (Entropy solution).

A weak solution of the initial value problem (1.4) is called an entropy solution, if the following Kružkov-type entropy inequality holds for all and all test functions .

 (1.6) ∫R ∫T0(|u−c|ψt+sign(u−c)k(x)(f(u)−f(c))ψx)dxdt+∫R|u0−c|ψ(x,0)dx +∫R∖{ξm}Mm=1∫T0sign(u−c)k′(x)f(c)ψdxdt+M∑m=1∫T0∣∣f(c)(k+m−k−m)∣∣ψ(ξm,t)dt≥0.

The last couple of decades have witnessed remarkable advances in the studies of conservation laws with discontinuous flux function. However, we will not be able to discuss the whole literature here, but only refer to the parts that are pertinent to the current paper. In case of “smooth” , the notion of entropy solution was introduced independently by Kružkov [18] and Vol’pert [31] (the latter author considered the smaller BV class). These authors also proved general existence, uniqueness, and stability results for the entropy solution, see also Oleĭnik [24] for similar results in the convex case .

### 1.1. Diffusive Dispersive Approximation

It is well known that the conservation law (1.4) is derived by neglecting underlying small scale effects such as diffusion, dispersion, capillarity etc., and may admit physically relevant discontinuous solutions containing shock waves (non-classical shock) that may depend on underlying small-scale mechanisms. It has been successively recognized that a standard entropy inequality (due to Kružkov, Oleĭnik, and others) does not suffice to single out such a physically relevant solution, and it is important to incorporate these small-scale effects in the entropy condition. In other words, additional admissibility criteria (a kinetic relation) are required in order to characterize these small-scale dependent non classical shock waves uniquely. In [9], the authors developed a framework for the existence and uniqueness of the non-classical shock waves that arise as limits of diffusive-dispersive approximations.

Noting that the solutions to (1.4) can be different due to their explicit dependence on the underlying small scale effects, we focus on a concrete model of two phase flow in porous medium (for a brief derivation of this model consult [4]). The relevant small scale effect is a dynamic capillary pressure term, that was introduced in [7]. Compared to the standard capillary pressure models [1], the addition of the new term resulted in a model that contain higher-order mixed spatio-temporal derivatives (cf. (1.2)).

The diffusive-dispersive model has a long tradition, starting with the analysis of linear diffusion-dispersion model (1.1). A pioneering study of the effect of vanishing diffusion and dispersion terms in scalar conservation laws, with -independent flux function, can be found in Schonbek [26]. The technique of compensated compactness was used to prove convergence toward weak solutions. Kondo and LeFloch [17] studied zero diffusion-dispersion limits for -independent fluxes under an optimal balance between the sizes of the diffusion and dispersion parameters. LeFloch and Natalini [19] used the concept of measure-valued solution and established convergence results assuming that the diffusion dominates the dispersion. Subsequently, the approach of kinetic decomposition and velocity averaging [25] was introduced by Hwang and Tzavaras [11] to analyze singular limits including nonclassical shock waves. Furthermore, we also mention related works by Wu [32] and Jacobs, McKinney, and Shearer [12] which provides the first existence result of undercompressive shocks for the modified Korteweg-de Vries-Burgers equation.

It is well known that the relative scaling between and determines the limiting behavior of solutions, and we can distinguish between three cases:

• Diffusion-dominant regime : The qualitative behavior of solution is same as the solution of the conservation laws.

• Dispersion-dominant regime : In this case, high oscillations develop (as ) especially in regions of steep gradients of the solutions and only weak convergence is observed.

• Balanced diffusion-dispersion regime: This typically corresponds to the scenario where . Only mild oscillations are observed near shocks, and the limit solution is a weak solution to conservation laws. Most importantly, in this case, the solution exhibit non-classical behavior, as they contain undercompressive shocks. However, for the regime , the limit solution coincides with the entropy solution determined by Kružkov theory [18].

### 1.2. Numerical Schemes

It is well known that standard finite difference, finite volume and finite element methods have been very successful in computing solutions to hyperbolic conservation laws with discontinuous coefficients, including those containing shock waves. However, we mention that most of these well-established methods are proven to be not good enough to capture nonclassical shock wave solutions numerically. This well known phenomena has been explained by many authers (Hou and LeFloch [10], Hayes and LeFloch [8], and others) in terms of the equivalent equation associated with discrete schemes through a formal Taylor expansion. The key idea behind capturing nonclassical shocks is to design finite difference schemes whose equivalent equation matches, both, the diffusive and the dispersive terms (cf. (1.3)) in the underlying model. However, these schemes fail to approximate nonclassical solutions with large amplitude, especially strong shocks due to lack of control on higher order error terms present in equivalent equation. A recent work by Ernest et al. [6] has overcome such problems by dominating higher order error terms in amplitude by the leading order terms of the equivalent equation.

In another paper by Chalons and Lefloch [2], the authors introduced a fully -discrete scheme for the numerical approximation of diffusive-dispersive hyperbolic conservation laws (cf. (1.1)-(1.3)) in one-space dimension. An important feature of their scheme is that it satisfies a cell entropy inequality and, as a consequence, the space integral of the entropy is a decreasing function of time. Moreover, they showed that the limiting solutions generated by the scheme contains nonclassical undercompressive shock waves.

On the other hand, there is a sparsity of efficient numerical schemes for (1.1)-(1.2) available in the literature. In fact, to the best of our knowledge, this is the first systematic attempt to construct a provably convergent numerical scheme for (1.1)-(1.2). Having said this, there are some numerical experiments available in the final section of the recent paper by Coclite et al. [4] without rigorous analysis of the scheme.

### 1.3. Scope and Outline of the Paper

In view of the above discussion, it is fair to claim that there are no robust and provably stable numerical schemes currently available to simulate the vanishing capillarity approximations of scalar conservation laws equation (1.1)-(1.2). In this context, we consider a fully-discrete (in both space and time) finite difference scheme for (1.1)-(1.2) which is provably convergent and able to capture non classical shocks quite well. Since diffusion-dispersion model for the conservation laws with discontinuous flux has not been studied in detail, we analyze a fully-discrete scheme for (1.1)-(1.3) as well. While there are several numerical methods which perform well in practice, perhaps better than the one presented here, (see [20] for a recent comparison of diferent numerical methods) we emphasize that we prove the convergence of the schemes proposed in this paper. Here, we mention that a detailed analysis of the scheme introduced by Ernest et al. [6] is beyond the scope of this paper, and will be the topic of an upcoming paper.

To sum up, the schemes in the present paper have the following properties:

• Approximate solutions for (1.1)-(1.2), generated by the scheme (3.2), converge to the unique entropy solution of (1.4) as long as . A scheme (cf. (7.2)) has been formulated for (1.1)-(1.3) and the same techniques can be applied, mutatis mutandis,to prove convergence of approximate solutions to the unique entropy solution of (1.4).

• Approximate solutions for (1.1)-(1.2) have been shown to converge to weak solutions of (1.4), when . Moreover, we show numerically that the limiting solutions generated by the schemes (3.2) and (7.2) contain nonclassical undercompressive shock waves.

The rest of the paper is organized as follows: In section 2, we present the mathematical framework used in this paper. In particular, we have used a compensated compactness result in the spirit of Tartar [27] but the proof is based on div-curl lemma and does not rely on the Young measure. Section 3 introduces the fully-discrete finite difference scheme for (1.1)-(1.2) . In section 4, we derive a priori estimates for the approximate solutions and a detailed convergence analysis towards weak solutions of (1.4) has been discussed in section 5. Convergence towards the unique entropy solution has been considered in section 6, while a brief discussion on the results for diffusive-dispersive approximation (1.1)-(1.3) has been addressed in section 7. Finally, numerical results are presented in section 8 to illustrate the performance of the designed schemes.

## 2. Mathematical Framework

In this section, we list all the assumptions on the data for the problem (1.1), and present relevant mathematical tools to be used in the subsequent analysis. Throughout this paper we use the letters etc. to denote various generic constants independent of approximation parameters, which may change line to line, but the notation is kept unchanged so long as it does not impact the central idea.

The basic assumptions on the data of the problem (1.1) are as follows:

1. [label=A.0]

2. For the initial function , we assume that

 u0∈L2(R)∩L∞(R),a≤u0(x)≤b,for a.ex∈R;
3. For the discontinuous coefficient , we assume that

 k∈L∞(R)∩BVloc(R),α≤k(x)≤β,for a.ex∈R;
4. Regarding the flux function , we assume that

 u↦f(k,u) ∈C2([a,b]),for allk∈[α,β], k↦f(k,u) ∈C1([α,β]),for allu∈[a,b];
5. Furthermore, we assume that is genuinely nonlinear a.e. in , i.e., , for a.e. .

###### Remark 2.1.

It is worth mentioning that the Assumption  4 is typically required in the compensated compactness framework. This condition also imposes a condition on the coefficient . In fact, it implies that is genuinely nonlinear (i.e., ) and , for a.e. .

Next, we recapitulate the results required from the compensated compactness method due to Murat and Tartar [23, 27]. For a nice overview of applications of the compensated compactness method to hyperbolic conservation laws, we refer to Chen [3]. Let denote the space of bounded Radon measures on and

 C0(R)={ψ∈C(R)∣∣lim|x|→∞ψ(x)=0}.

If , then

 ⟨μ,ψ⟩=∫Rψdμ,for allψ∈C0(R).

Recall that if and only if for all . We define the norm

 ∥μ∥M(R):=sup{|⟨μ,ψ⟩|:ψ∈C0(R),∥ψ∥L∞(R)≤1}.

The space is a Banach space and it is isometrically isomorphic to the dual space of . Furthermore, we define the space of probablity measures

 Prob(R):={μ∈M(R):μis nonnegative and∥μ∥M(R)=1}.

Before we state the compensated compactness theorem, we shall recall the celebrated div-curl lemma [23].

###### Lemma 2.1 (div-curl lemma).

Let be a bounded open subset of . Suppose

 u1Δx⇀¯¯¯u1,u2Δx⇀¯¯¯u2,v1Δx⇀¯¯¯v1,andv2Δx⇀¯¯¯v2,

in as . Furthermore, assume that the two sequences and lie in a (common) compact subset of , where and . Then along a subsequence

 (u1Δx,u2Δx)⋅(v1Δx,v2Δx)↦(¯¯¯u1,¯¯¯u2)⋅(¯¯¯v1,¯¯¯v2),inD′(Ω),asΔx↓0.

Suitably modified for our purpose, we shall use the following compensated compactness result. For a proof, we refer to the paper by Kenneth and Towers [13, Lemma 3.2].

###### Theorem 2.1.

Assume that 2, 3 and 4 hold. Let be a bounded open set, and assume that is a sequence of uniformly bounded functions such that , for all . Set

 (2.1) (η1(s),q1(k,s)) =(s−c,f(k,s)−f(k,c)), (η2(k,s),q2(k,s)) =(f(k,s)−f(k,c),∫sc(fs(k,θ))2dθ),

where is an arbitrary constant. If the two sequences

 {η1(uΔx)t+q1(k(x),uΔx)x}Δx>0,and{η2(k(x),uΔx)t+q2(k(x),uΔx)x}Δx>0

belong to a compact subset of , then there exists a subsequence of that converges a.e. to a function .

We remark that, a feature of the compensated compactness result above is that it avoids the use of the Young measure by following an approach developed by Chen and Lu [3, 22] for the standard scalar conservation law. This is preferable as the fundamental theorem of Young measures applies most easily to functions that are continuous in all variables.

The following compactness interpolation result (known as Murat’s lemma [23]) is useful in obtaining the compactness needed in Theorem  2.1.

###### Lemma 2.2.

Let be a bounded open subset of . Suppose that the sequence of distributions is bounded in . Suppose also that

 LΔx=L1,Δx+L2,Δx,

where is in a compact subset of and is in a bounded subset of . Then is in a compact subset of .

## 3. A Fully-Discrete Finite Difference Scheme

We begin by introducing some notation needed to define the fully-discrete finite difference scheme. Throughout this paper, we reserve to denote small positive numbers that represent the spatial and temporal discretizations parameter of the numerical scheme. Given , we set for , to denote the spatial mesh points. Similarly, we set for , where for some fixed time horizon . Moreover, for any function admitting point values, we write . Furthermore, let us introduce the spatial and spatial-temporal grid cells

 Ij=[xj−1/2,xj+1/2),Inj=[xj−1/2,xj+1/2)×[tn,tn+1).

where . Let denote the discrete forward and backward differences in space, i.e.,

 D±uj=±uj±1−ujΔx,

The discrete Leibnitz rule is given by

 D±(ujvj)=ujD±vj+vj±1D±uj

while the summation-by-parts formula is given by

 ∑j∈ZujD±vj=−∑j∈ZvjD∓uj.

Furthermore, for any function , using the Taylor expansion on the sequence we obtain

 D±f(uj)=f′(uj)D±uj±Δx2f′′(ξj±12)(D±uj)2,

for some between and . In other words, the discrete chain rule is accurate up to an error term of order .

Finally, let denote the discrete forward and backward difference operator in the time, i.e.,

 Dt±unj=∓un±1j−unjΔt.

The following identity is readily verified:

 (3.1) unjDt+unj=12Dt+(unj)2−Δt2(Dt+unj)2.

We propose the following fully-discrete (in space and time) finite difference scheme approximating the limiting solutions generated by the equation (1.1)-(1.2)

 (3.2) Dt+unj+D−hnj+12 =βΔxD+D−unj+γμ(Δx)Dt+D+D−unj,j∈Z,n∈N0, (3.3) u0j =1Δx∫xj+12xj−12u0(θ)dθ,j∈Z,

where are fixed parameters, and as . More specifically, we will either use or depending on the quest for the convergence of approximate solution towards a weak solution or the entropy solution, respectively.

###### Remark 3.1.

Here we used the notation to denote quantities that depend on and are bounded above by , where is a constant independent of . Likewise, we used the notation to denote quantities that depend on and are bounded above by , where is a constant independent of and .

The numerical flux corresponding to the flux function is given by

 hnj+12=kj+12^fnj+12,withkj+12=1Δx∫xj+1xjk(x)dx,

where is based on a two-point monotone numerical flux, i.e., non-decreasing with respect to the first argument and non-increasing with respect to the second argument, and consistent with the actual flux, i.e., . Moreover, in order to maintain monotonicity of the scheme (3.2) without the higher order terms (corresponds to ), the arguments of the numerical flux are transposed when the coefficient is negative. More specifically, we choose

 ^fnj+12=⎧⎪⎨⎪⎩^f(unj,unj+1),if kj+12≥0^f(unj+1,unj),if kj+12<0.

Summing up, the numerical flux is given by

 hnj+12=⎧⎪⎨⎪⎩kj+12^f(unj,unj+1),if kj+12≥0kj+12^f(unj+1,unj),if kj+12<0.

In particular, we focus on Engquist-Osher (EO) numerical flux given by

 (3.4) ^f(u,v)=12(f(u)+f(v))−12∫vu∣∣f′(s)∣∣ds.
###### Remark 3.2.

We have chosen to analyse the scheme (3.2)–(3.3) with EO flux because of its apparent simplicity. One can, however, adopt the method of proof developed in this paper and obtain similar results for other schemes (e.g., all monotone schemes).

To this end, observe that the EO flux given by (3.4) is Lipschitz continuous . In fact, for , it has continuous partial derivatives satisfying

 (3.5) f′−(v)=^fv(v,u)≤0≤^fu(v,u)=f′+(u),

using the conventional notations that and . It is also clear that serves as a Lipschitz constant for EO flux.

For a given initial data , we define the initial grid function by (3.3). Moreover, for the sequence , we associate the function defined by

 uΔx(x,t)=∑j∈Z,n≥0unj1Inj(x,t),

where denotes the characteristic function of the set . Similarly, for the coefficient approximated at each cell boundary, we associate the function defined by

 kΔx(x)=∑j∈Zkj+121Ij+12(x)

Note that and are discretized on grids that are staggered with respect to each other. This indeed results in a reduction in complexity, compared with the approach where two discretizations are aligned.

Throughout this paper, we use the notation to denote the functions associated with the sequence . For later use, recall that the discrete , and norms, and BV semi-norm for a lattice function are defined respectively as

 ∥uΔx(⋅,tn)∥ℓ2(R)=√Δx∑j∈Z∣∣unj∣∣2,|uΔx(⋅,tn)|BV=∑j∈Z∣∣unj+1−unj∣∣

For the sake of simplified notations, unless specified, we shall use the notation to denote the discrete norm.

## 4. A Priori Estimates

This section is devoted to the derivation of a priori estimates which turns out to be useful to prove “strong compactness” of the approximate solution . To begin with, following Coclite et al. [4], we assume that the approximate solutions generated by the scheme are uniformly bounded, i.e., . In other words, we assume that

1. [label=B.0,resume]

2. For almost every , , for some fixed constants ;

###### Remark 4.1.

It is worth mentioning that the above assumption on the approximate solutions is the manifestation of the specific structure of the flux function (depends explicitly on space variable). In fact, to obtain bound on the solution, one requires a priori bound on the solution (cf. (4.4)). This assumption can be toppled by replacing the “space dependent flux function” to a flux function which depends explicitly only on the solution. In such a scenario, one can use framework of the compensated compactness result [5], to reproduce all the results in this paper.

To proceed further, we first collect all the available estimates on the approximate solutions in the following lemma.

###### Lemma 4.1.

Let be a sequence of approximations generated by the scheme (3.2). Moreover, assume that the initial data lies in . Then the following estimate holds

 (4.1) 12Dt+∥un∥2+(γμ(Δx)2+βΔx22) Dt+∥D−un∥2+δΔx∥D−un∥2 +δΔx∥∥Dt+un∥∥2+δγΔxμ(Δx)∥∥Dt+D−un∥∥2≤C,

provided and satisfies the following CFL condition

 (4.2) max{2max{∥f∥∞,∥∥f′∥∥∞,∥k∥∞}+λ2,λ2(1+βΔx2γμ(Δx))}≤min(1−δ,β−δ),

with . Here and the constant is independent of .

In particular, the estimate (4.1) guarantees following space-time estimates:

 (4.3a) ∀n∈N,Δx∑j(unj)2 ≤C, (4.3b) Δx2Δt∑j∑n(D−unj)2 ≤C, (4.3c) Δx2Δt∑j∑n(Dt+unj)2 ≤C, (4.3d) ΔtΔx2μ(Δx)∑j∑n(Dt+D−unj)2 ≤C.
###### Remark 4.2.

In light of the CFL condition (4.2), we want to emphasize that if (required to prove convergence towards a weak solution, cf. Theorem 5.1), then we need to be kept fixed. On the other hand, if (required to prove convergence towards the entropy solution, cf. Theorem 6.1), then we need , with , to be kept fixed. To sum up, we need a stronger CFL condition to prove convergence of approximate solutions towards the unique entropy solution of (1.1). To this end, we mention that in the subsequent analysis the CFL condition (4.2) is always assumed to hold.

###### Proof.

To start with, we multiply the scheme (3.2) by and subsequently sum over . Then, using summation-by-parts formula and the identity (3.1), we obtain

 12Dt+∑jΔx(unj)2 −Δt2∑jΔx(Dt+unj)2+Δx∑junjD−hnj+12IΔx(f) =−βΔx∑jΔx∣∣D−uj∣∣2−γμ(Δx)Δx∑jD−unjDt+(D−unj).

Note that, the identity (3.1) also implies that

 μ(Δx)Δx∑jD−unjDt+(D−unj) =Δxμ(Δx)2∑jDt+(D−unj)2−Δxμ(Δx)Δt2∑j(Dt+D−unj)2.

Next, we move on to estimate the term . Using summation-by-parts formula, we obtain

 −IΔx(f) =−∑jΔxunjD+hnj−12=∑j(unj−unj−1)kj−12^fnj−12 =∑jkj−12[(unj−unj−1)^fnj−12−(F(unj)−F(unj−1))]+(F(unj)−F(unj−1)) =∑jkj−12[(unj−unj−1)^fnj−12−(F(unj)−F(unj−1))]Ej,n(f)−∑jF(unj)(kj+12−kj−12),

where is the primitive of , i.e., . A first order Taylor’s expansion together with the monotonicity of the numerical flux function gives us an estimate of . To see this, notice that

 Ej,n(f)=kj−12(unj−unj−1)(^fnj−12−f(unj−12)),

where lies in between and . To proceed further, we consider the following two cases:
Case 1: Assume that , then by the definition of numerical flux . If , then

 ^fnj−12≥^f(unj−12,unj)≥^f(unj−12,unj−12)=f(unj−12).

On the other hand, if , then

 ^fnj−12≤^f(unj−12,unj)≤^f(unj−12,unj−12)=f(unj−12).

Thus, in any case we have .
Case 2: Assume that , i.e., . If , then

 ^fnj−12≤^f(unj−12,unj−1)≤^f(unj−12,unj−12)=f(unj−12).

On the other hand, if , then

 ^fnj−12≥^f(unj−12,unj−1)≥^f(unj−12,unj−12)=f(unj−12).

Therefore, in this case also, we have . Having this in mind and making use of the Assumption 1, we conclude that

 (4.4) −IΔx(f)≤∣∣ ∣∣−∑j(kj−12−kj+12)F(unj)∣∣ ∣∣≤BV(k)maxj∣∣F(unj)∣∣≤C,

where is a constant independent of . Finally, combining all the above estimates, we obtain

 (4.5) 12Dt+∥un∥2−Δt2∥∥Dt+un∥∥2

Next, we multiply the scheme (3.2) by and sum over to obtain,

 (4.6) Δx2∑j(Dt+unj)2 +Δx2∑jD−hnj+12Dt+unj (4.7) =−β(Δx)3∑jD−unjDt+(D−unj)−γμ(Δx)Δx2∑j(Dt+D−unj)2.

Again, use of the identity (3.1) reveals that

 β(Δx)3∑jD−unjDt+(D−unj)=βΔx32∑jDt+(D−unj)2−βΔx