GLOBAL ENTROPY SOLUTIONS TO EXOTHERMICALLY REACTING, COMPRESSIBLE EULER EQUATIONS

GUI-QIANG CHEN AND DAVID H. WAGNER

Abstract. The global existence of entropy solutions is established for the compressible Euler equations for one-dimensional or plane-wave flow of an ideal gas, which undergoes a one-step exothermic chemical reaction under Arrhenius-type kinetics. We assume that the reaction rate is bounded away from zero and the total variation of the initial data is bounded by a parameter that grows arbitrarily large as the equation of state converges to that of an isothermal gas. The heat released by the reaction causes the spatial total variation of the solution to increase. However, the increase in total variation is proved to be bounded in t > 0 as a result of the uniform and exponential decay of the reactant to zero as t approaches infinity.

1. Introduction

We are concerned with the large-time existence of entropy solutions to the Cauchy problem for the equations of planar flow of an exothermically reacting ideal gas:

ρ_{t} + (ρu)_{x} = 0,  (1.1a)
(ρu)_{t} + (ρu^{2} + p)_{x} = 0,  (1.1b)
(ρE)_{t} + ((ρE + p)u)_{x} = qρY ϕ(T),  (1.1c)
(ρY)_{t} + (ρuY)_{x} = −ρY ϕ(T),  (1.1d)
(ρ, u, E, Y)(x, 0) = (ρ_{0}, u_{0}, E_{0}, Y_{0})(x).  (1.1e)

We assume, for simplicity, that the specific heats and molecular weights of the reactant and product gases are the same. Then the constitutive relations and conditions for this

Date: April 20, 2002.

system are

v = \frac{1}{ρ} (specific volume),  E = e(v, S) + \frac{u^{2}}{2} ,  (1.2a)
p = −e_{v}(v, S) > 0,  T = e_{S}(v, S) = \frac{p}{Rρ} > 0,  (1.2b)
p_{v}(v, S) < 0, p_{vv}(v, S) > 0,  (1.2c)

on any compact set in v > 0. We assume that the reaction rate function ϕ is monotonically increasing and Lipschitz continuous. In addition, inadmissible discontinuous solutions are eliminated by requiring the following entropy condition:

(ρS)_{t} + (ρuS)_{x} ≥ \frac{qρY ϕ(T)}{T}.  (1.3)

The ideal gas assumption p = RρT implies [34] that the internal energy can be written as a function of T alone, e(v, S) = ẽ(T). Then ẽ′(T) = c_{v}(T) is the specific heat at constant volume. We identify c_{v} + R = c_{p} as the specific heat at constant pressure.

When c_{v} is constant, then the gas is an ideal polytropic gas. For such a gas, c_{p} is also constant, and the ratio c_{p}/c_{v} = γ > 1 is the adiabatic exponent, a parameter that determines the equation of state.

This system of equations (1.1)–(1.2) is useful for studying the behavior of plane detonation waves. In a detonation wave, the effects of pressure gradients, which support the shock wave, and the conversion of chemical energy to mechanical energy are far greater than the diffusive effects such as viscosity, heat conduction, and diffusion of chemical species. This justifies the use of the Euler equations in (1.1)–(1.2), rather than the Navier-Stokes equations, in this context. The shock wave solutions one observes in this model are jump discontinuities. This is a very good representation of the shock waves one observes experimentally, which have a width of several molecular mean free paths. The reaction zone of a detonation wave, by way of contrast, has a width which is generally hundreds of mean free paths.

Our interest in this system is partly stimulated by an interest in new and different types of behavior exhibited by solutions of this system. Whereas non-reacting shock waves are known, under reasonable assumptions, to be stable [22], linearized stability analysis, as well as numerical and physical experiments, have shown that certain steady detonation waves are unstable [1, 13, 14, 18, 24]. One particular kind of instability that takes place within the context of one space dimension produces pulsating detonation waves. In certain parameter regimes, steady planar detonation waves are unstable and

REACTING COMPRESSIBLE EULER EQUATIONS 3

evolve into oscillating waves. These oscillating waves generate a steady stream of waves which propagate behind the wave. For example, a numerical calculation of such an evolution, performed by the second author, is presented in Fig. 1. The possibility of such

Figure 1. Computation of a pulsating detonation.

oscillation clearly indicates that we are presented with both an interesting challenge and the possibility of discovering new and interesting mathematics.

The system (1.1)-(1.2) is a hyperbolic system of balance laws. That is, it has the form

U_{t} + F (U)_{x} = G(U),  (1.4)

where U(x, t) ∈ ℝ^{n}, F (U), G(U) ∈ ℝ^{n}, and for each U, the n × n matrix DF (U) has a set of real eigenvectors which form a basis of ℝ^{n}. One of the principal features of the theory of quasilinear hyperbolic systems is the formation of shock waves. This shock formation makes difficult the establishment of theorems regarding the existence and/or the uniqueness of solutions. Classical, smooth solutions to the Cauchy problem will not usually exist for all t > 0. However, weak, discontinuous entropy solutions to this problem, under reasonable conditions, do exist for all t > 0, as we shall show in this paper.

In the case where G(U) is identically zero, the system (1.4) is called a system of conservation laws. Such a system is strictly hyperbolic if the n eigenvalues of DF (U) are real and distinct: λ_{1} < ⋯ < λ_{n}. The right eigenvectors r_{1}, ⋯, r_{n} correspond to the fields of simple waves admitted by the system (1.4). The corresponding eigenvalues

4 GUI-QIANG CHEN AND DAVID H. WAGNER

give the characteristic speeds of propagation associated with these fields. We say that the j^{th} wave field is genuinely nonlinear if ∇λ_{j}(U) · r_{j}(U) ≠ 0 for all U. The j^{th} wave field is linearly degenerate if ∇λ_{j}(U) · r_{j}(U) = 0 for all U.

The large-time existence of entropy solutions of strictly hyperbolic and genuinely nonlinear systems of conservation laws, with initial data of small total variation, was proved in [15]. The equations of gas dynamics, that is, the first three equations of (1.1), together with (1.2), omitting any terms containing Y, is not completely genuinely nonlinear. One may compute that

λ_{1} = u − c, λ_{2} = u, λ_{3} = u + c,

where c = v\sqrt{−p_{v}} is the speed of sound. The first and third wave fields propagate acoustic waves and shock waves and are genuinely nonlinear. The second field is linearly degenerate and propagates contact discontinuities. The large-time existence of entropy solutions of the Cauchy problem for gas dynamics is proved in [19, 28].

The system of conservation laws for (1.1), that is, omitting the reaction rate term ρY ϕ(T), has a fourth wave field which is linearly degenerate with λ_{4} = λ_{2} = u. Thus the system (1.1) is not strictly hyperbolic. However, it is hyperbolic, and the lack of strict hyperbolicity does not affect the existence theory for this system. In fact one may rewrite this system in Lagrangian coordinates:

v_{t} − u_{y} = 0,  (1.5a)
u_{t} + p_{y} = 0,  (1.5b)
E_{t} + (pu)_{y} = qY ϕ(T),  (1.5c)
Y_{t} = −Y ϕ(T).  (1.5d)

Because this choice of coordinates reduces the fourth equation of the system to an ordinary differential equation, which is coupled with the rest of the system only through the temperature function T, we will work primarily with the system (1.5). Any existence theorem proved for this system can be translated into an existence theorem for (1.1) (cf. see [30]).

Several results have been obtained regarding the existence of entropy solutions to hyperbolic systems of balance laws [9, 10, 21, 35, 36]. However, the equations of reacting flow considered here do not satisfy the hypotheses for these results. Our problem presents a double eigenvalue, which is ruled out in [9] to prevent resonance. Some 2 × 2 physical systems were discussed in [10, 21, 35, 36]. The papers [9, 21] had in view applications in

REACTING COMPRESSIBLE EULER EQUATIONS 5

which the lower-order terms act in a way that reduces the spatial total variation of the solution as time increases. In [36], the lower-order term has a coefficient e^{−Kt}, so that the lower-order term decays uniformly without the need for a priori estimates. In [35], the first order terms constitute a system of conservation laws for which solutions exist with arbitrarily large initial data [23]. Such a system is rather unusual, and no decay of the lower-order terms is required in order to obtain the large-time existence of entropy solutions.

These remarks do not apply to the system (1.1)–(1.2), or (1.5). For this system, the exothermic reaction can increase the total variation in a number of ways. For example, in the formation of a detonation wave, a chemical reaction behind a shock wave can increase the strength of that shock wave. More subtle phenomena are also possible. In a nearly constant, unreacted state, a very small variation in temperature can cause the gas in one region to react prior to the gas in nearby regions, resulting in a large increase in total variation. Moreover, the hot spot created by such an event would generate waves, some of which would be shock or rarefaction waves. These waves could propagate away from the hot spot before the remaining reactant ignites.

The theorem that we present in this paper is only a first step in dealing with these difficulties. We assume that the initial data are such that the reaction rate function ϕ(T) never vanishes, so that there is a positive minimum value Φ := ϕ(T′) > 0. In a sense, this is a very realistic condition. Typically, ϕ(T) has the Arrhenius form:

ϕ(T) = T^{α}e^{−E/RT} ,  (1.6)

which vanishes only at absolute zero temperature. However, in a typical unburned state, ϕ(T) is very small. We make this assumption in order to obtain the uniform decay of the reactant to zero. Thus, although the total variation of the solution may very well increase while the reaction is active, the reaction must eventually die out. Consequently, the increase in total variation can be estimated rigorously.

Following [28], we consider a one-parameter family of functions e(v, S, ϵ), ϵ ≥ 0, which is C^{5} and satisfies (1.2). We assume that, when ϵ = 0, the equation of state is that of an isothermal gas:

e(v, S, 0) = − ln(v) + \frac{S}{R}.  (1.7)

For a polytropic gas, ϵ = γ − 1, and for ϵ > 0,

e(v, S, ϵ) = \frac{1}{ϵ} ((v exp (−S/R))^{−ϵ} − 1) .  (1.8)

6 GUI-QIANG CHEN AND DAVID H. WAGNER

One may easily check that this function is C^{∞} and that, as ϵ → 0+, all partial derivatives converge uniformly on any compact set in v > 0 to the corresponding derivatives of e(v, S, 0). In particular, one may use L’Hôpital’s rule to calculate

\frac{∂e}{∂ϵ}(v, S, 0) = \frac{1}{2} (− ln(v) + \frac{S}{R})^{2} ,  (1.9)

and see that \frac{∂e}{∂ϵ}(v, S, ϵ) is continuous at ϵ = 0, v > 0.

The value ϵ = 0 is mathematically special because, at this value, the system (1.1) has a complete system of Riemann invariants:

(r, s, S, Y) = (u − ln(p), u + ln(p), S, Y).  (1.10)

Furthermore, all shock, rarefaction, and contact discontinuity curves in (r, s, S, Y)-space are invariant under translation of the base point. Following [28], we use (r, s, S, Y) as the coordinates for our analysis in ϵ ≥ 0. Note that since p = −\frac{∂e}{∂v}(v, S, ϵ) and e(v, S, ϵ) is C^{5}, the transformation between (v, u, S) and (r, s, S) is C^{4}. Moreover, Temple [28] showed that this transformation is a diffeomorphism.

Our principal result, Theorem 3.2, is somewhat complex, but a simple version may be stated as follows.

Theorem 1.1. Let (r_{−∞}, s_{−∞}, S_{−∞}) be a point in rsS-space, and let ϵ ∈ [0, 1]. Let w_{0}(x) be given initial data for (1.5), expressed in (r, s, S, Y)-coordinates, and with lim_{x→−∞} w_{0}(x) = (r_{−∞}, s_{−∞}, S_{−∞}, 0). Then there is a function C (ϵ, Φ, q, ∥Y_{0}∥_{∞}) > 0 such that C (ϵ, Φ, q, ∥Y_{0}∥_{∞}) → ∞ as ϵ → 0, and such that, if

V ar_{rsSY} (w_{0}) < C (ϵ, Φ, q, ∥Y_{0}∥_{∞}) ,

then there exists a global BV entropy solution to the Cauchy problem (1.1)–(1.2) with initial data determined by w_{0}.

Moreover, for ϵ small, C (ϵ, Φ, q, ∥Y_{0}∥_{∞}) has the form

C (ϵ, Φ, q, ∥Y_{0}∥_{∞}) = \frac{B}{ϵ^{1/3}} exp (−\frac{Kq ∥Y_{0}∥_{∞}}{Φ} ) .

for some constant B.

There is an interesting common thread connecting our results with previous ones concerning balance laws. While earlier results had in view lower-order terms that exerted a damping effect, or otherwise reduced the total variation, our result requires the decay

REACTING COMPRESSIBLE EULER EQUATIONS 7

of the lower-order term, even though the total variation may increase in the process. Thus, in either case, decay of some kind seems essential.

The proof of Theorem 1.1 (and of the more sophisticated Theorem 3.2) may be outlined as follows. In Section 2, we construct approximate solutions to the Cauchy problem using a fractional–step method based on Glimm’s scheme. In Section 3, we establish bounds on the spatial total variation of our approximate solutions. The proof of these bounds involves use of Temple’s results for non-reacting gas dynamics [28]. In Section 4, we use our bounds to prove that a subsequence from our approximate solutions converges to an entropy solution of the Cauchy problem.

2. Approximate Solutions

In this section, we construct approximate solutions to the Cauchy problem for (1.5) and (1.2). Using (r, s, S, Y) as coordinates, the initial data are given by

(r, s, S, Y)|_{t=0} = (r_{0}, s_{0}, S_{0}, Y_{0})(x).  (2.1)

We may represent (1.5), (1.2), and (2.1) in the compact form:

\begin{aligned}U_{t} + F (U)_{x} & = G(U), \\ U|_{t=0} & = U_{0}(x). \end{aligned}  (2.2)

The approximate solutions are constructed by using the Glimm scheme, combined with a fractional–step method which incorporates the reaction rate term. We begin by discussing the Riemann solutions for the homogeneous (non-reacting) system.

2.1. The Riemann Problem for the Non-reacting Gas. The Riemann problem for the non-reacting gas is the following Cauchy problem:

v_{t} − u_{x} = 0,  (2.3a)
u_{t} + p(v, S, ϵ)_{x} = 0,  (2.3b)
(e(v, S, ϵ) + u^{2}/2)_{t} + (p(v, S, ϵ)u)_{x} = 0,  (2.3c)
Y_{t} = 0,  (2.3d)

with the initial data:

(r, s, S, Y)|_{t=0} = \begin{cases}(r_{L}, s_{L}, S_{L}, Y_{L}), & x < 0, \\ (r_{R}, s_{R}, S_{R}, Y_{R}), & x > 0, \end{cases}  (2.4)

8 GUI-QIANG CHEN AND DAVID H. WAGNER

where ϵ ∈ [0, 1], and (r_{L}, s_{L}, S_{L}, Y_{L}) and (r_{R}, s_{R}, S_{R}, Y_{R}) are constant states. Note that (2.3d) is independent of (2.3a,b,c). Consequently, the solution of the Riemann problem (2.3)–(2.4) decouples into the solution of two Riemann problems: one for nonreacting gas dynamics, (2.3a,b,c), and a second trivial Riemann problem for Y. For the remainder of this section, we will omit any discussion of Y.

An explicit solution to the Riemann problem for the system (2.3 a,b,c) for polytropic gases can be found in [7, 27]. For the more general case, with the internal energy given by e(v, S, ϵ), the solution of the Riemann problem was proved to exist for small data in [17], and for large data in [26, 2]. The solution to the Riemann problem consists of forward and backward rarefaction and shock waves, together with contact discontinuities of speed zero. These waves are separated by regions in which the solution is constant.

A backward, or 1-rarefaction wave, is a solution which is constant on the 1-characteristics

centered at a point (x_{0}, t_{0}): \frac{x − t_{0}}{t − t_{0}} = \frac{dx}{dt} = λ_{1} ≡ −\sqrt{−p_{v}(v, S)}, where |λ_{1}| = ρc = c/v

is the material, or Lagrangian, sound speed of the gas. In addition, the 1-Riemann invariants S and u − ∫^{v} \sqrt{−p_{v}(v, S)}dv are constant within the wave.

A forward, or 3-rarefaction wave, is similar to the 1-wave except that λ_{3} = −λ_{1}, and the 3-Riemann invariants are S and u + ∫^{v} \sqrt{−p_{v}(v, S)}dv.

Forward and backward shock waves are simple jump discontinuities of U along a line \frac{x − x_{0}}{t − t_{0}} = dx/dt = σ, and satisfy the Rankine-Hugoniot conditions:

\begin{matrix}σ [v] = − [u] , \\ σ [u] = [p] , \\ σ [u^{2}/2 + e] = [pu] . \end{matrix}  (2.5)

Here [f] denotes the jump in the quantity f across the shock wave. In addition to (2.5), shock waves must satisfy the entropy condition, which states that the specific entropy S increases as t increases, or equivalently, as the gas crosses the shock wave. Under the conditions stated in (1.2), Weyl [33] showed that, for a fixed state U_{0} on the upstream side of a shock wave, the change in entropy across the wave is monotone increasing with the strength of the wave.

Contact discontinuities are simple jump discontinuities in S along a straight line of speed λ_{2} = 0. The quantities p and u are constant across contact discontinuities.

For given w_{L}, there are curves of states w = H_{i}(z_{i}, w_{L}, ϵ) emanating from w_{L} associated with each characteristic field such that w = H_{i}(z_{i}, w_{L}, ϵ) is connected to w_{L} by an i–wave of strength z_{i}.

REACTING COMPRESSIBLE EULER EQUATIONS 9

We parametrize 1–shock curves by z_{1} = r−r_{L} ≤ 0 and 3–shock curves by z_{3} = s−s_{L} ≤ 0, where r_{L} and s_{L} are the values of the coordinates r and s on the left side of the shock wave. We parametrize 1–rarefaction curves by z_{1} = r − r_{L} ≥ 0 and 3–rarefaction curves by z_{3} = s − s_{L} ≥ 0. For given w_{L}, the 1–rarefaction and 1–shock curve based at w_{L} meet at w_{L} with C^{2} contact to form a single curve parametrized by z_{1}. For z_{1} ≠ 0, these curves are at least C^{3} [17, 28]. Similarly, the 3–shock and 3–rarefaction curves meet at w_{L} to form a single curve parametrized by z_{3}, which is C^{2} at w_{L} and at least C^{3} for z_{3} ≠ 0.

There is a line of states in (r, s, S)-space which are connected to w_{L} by a contact discontinuity. The equations for this line are p = p_{L}, u = u_{L}, or simply r = r_{L}, s = s_{L}. We parametrize this line by z_{2} = S − S_{L}.

Let z = (z_{1}, z_{2}, z_{3}) and

H_{g}(z, w_{L}, ϵ) = H_{3} (z_{3}, H_{2} (z_{2}, H_{1} (z_{1}, w_{L}, ϵ) , ϵ) , ϵ) .

Henceforth, we use a subscript g to denote a function or vector relating to non-reacting gas dynamics. Then the Riemann problem with data (w_{L}, w_{R}) has a solution if and only if w_{R} is in the range of the map z → H_{g}(z, w_{L}, ϵ). In [28], it is shown, using the implicit function theorem, that for every w = (r, s, S) ∈ ℝ^{3} there exists a neighborhood Ω of w such that the Riemann problem is solvable for (w_{L}, w_{R}) ∈ Ω × Ω and ϵ ∈ [0, 1], and that one can solve for the vector of signed wave strengths z = (z_{1}, z_{2}, z_{3}) as a function of (w_{L}, w_{R}, ϵ),

z = B_{g}(w_{L}, w_{R}, ϵ),

where B_{g} is C^{2} with locally Lipschitz second derivatives. The solution is unique (among self-similar solutions) if this map is one–to–one.

To this solution of the Riemann problem for the non-reacting gas, we now add the decoupled solution of the linear Riemann problem for Y. We let H_{4}(z_{4}, r_{L}, s_{L}, S_{L}, Y_{L}) = (r_{L}, s_{L}, S_{L}, Y_{L} + z_{4}), and we let B(w_{L}, w_{R}, ϵ) the solution of the Riemann problem for (2.3) with data (w_{L}, w_{R}) ∈ (Ω × [0, 1])^{2}. Clearly, B will also be C^{2} with locally Lipschitz second derivatives.

We will make use of B and its regularity properties in Section 3, where we will derive conditions for the total variation stability of the fractional–step approximation scheme. This result will depend only on the regularity of B and the size of the data measured in total variation and in the uniform norm. Specifically, the result does not depend on the choice of parameters (z_{1}, z_{2}, z_{3}, z_{4}) used to describe H and B, as long as the regularity properties of B are maintained.

2.2. The Glimm Fractional–Step Scheme. We will use a fractional–step scheme such as that described in [9, 21] based on the Glimm scheme [15]. The construction of the fractional–step scheme for the inhomogeneous system (1.4) is as follows.

Choose mesh lengths h > 0 and l > 0 in the t and x directions, respectively, such that the Courant-Friedrichs-Levy condition:

Λ = max_{1≤j≤4} |λ_{j}(u)| ≤ \frac{l}{2h}  (2.6)

is satisfied.

Partition ℝ^{+} by the sequence t_{k} = kh, k ∈ ℤ^{+}, and partition ℝ into cells with the j-th cell centered at

x_{j} = jl, j = 0, ±1, ±2, · · · .

We begin by approximating the initial data U_{0}(x) by a function U^{h}(x, 0) which is constant for x_{j−1} ≤ x < x_{j+1} for j even, and which converges to U_{0}(x) pointwise a.e., and in L^{1} on all bounded intervals as h → 0. Choose a random sequence χ_{k}, k = 0, 1, 2, ⋯, from the uniform probability distribution on the interval (−1, 1). Then our approximation scheme and approximate solutions can be written in the following abstract form:

U^{h}(x, t) = F(t − kh, S_{0}(t − kh, U^{h}(x, kh + 0))),  kh ≤ t < (k + 1)h,
U^{h}(x, kh + 0) = U_{j}^{k},  j + k even, (j − 1)l < x < (j + 1)l,

U_{j}^{k+1} = R_{j}^{k+1} ∘ U^{h}(x, (k + 1)h−).

Here F(τ, ·) is the fractional–step operator, which advances the chemical reaction for 0 < τ ≤ h, S_{0}(τ, ·) is the solution operator for the homogeneous system (2.3) (i.e. G(U) ≡ 0) for 0 < τ ≤ h, and R_{j}^{k} is the random choice operator determined by the random sequence. The detailed description of these operators is given below.

Assume that U^{h}(x, t) is defined for t < kh. Then we define U^{h}(x, kh + 0) as follows.

Random Operator R_{j}^{k}. We define

U_{j}^{k} = R_{j}^{k} ∘ U^{h}(x, kh−) ≡ U^{h}((j + χ_{k})l, kh−),
U^{h}(x, kh + 0) ≡ U_{j}^{k}, (j − 1)l < x < (j + 1)l,

where j + k is even and χ_{k} is the k^{th} element of the random sequence.

REACTING COMPRESSIBLE EULER EQUATIONS 11

Homogeneous Solution Operator S_{0}(t, ·). In the strip kh ≤ t < (k + 1)h, we define

U_{0}^{h}(·, t) ≡ S_{0}(t − kh, U^{h}(·, kh + 0)),

where S_{0}(τ, W) is the solution U(·, τ) of the following Cauchy problem:

\begin{matrix}U_{τ} + F (U)_{x} = 0, \\ U|_{τ=0} = W (x). \end{matrix}  (2.7)

This Cauchy problem can be solved by constructing the Riemann solutions of the Riemann problems (2.3)–(2.4) with Riemann data:

U_{L} = U_{j}^{k}, U_{R} = U_{j+2}^{k}, j = 0, ±1, ±2, · · · .

This construction determines the unique solution of the Cauchy problem (2.7), based on the existence and uniqueness result for the Riemann problem of the system. This solution is valid as long as the waves of the different Riemann problems do not interact. The non-interaction of these waves is guaranteed by the Courant-Friedrichs-Levy condition (2.6).

Fractional–Step Operator F(τ, ·). The fractional–step operator uses an approximation to the solution of the initial value problem:

\begin{matrix}\frac{dU}{dτ} = G(U(τ)), \\ U|_{τ=0} = U_{0}. \end{matrix}  (2.8)

This approximation must be consistent and total variation stable. These requirements are defined as follows.

Definition 2.1. The explicit one-step approximation algorithm Ū(τ) = F(τ, U_{0}) is a consistent first order approximation to the solution U(τ) of the initial value problem (2.8) in the domain D ⊂ ℝ^{n} if F is Lipschitz continuous on (−δ, δ) × D → ℝ^{n} and there exists an increasing function ν(τ) > 0 such that ν(τ) converges to 0 as τ → 0 and

∥Ū(τ) − U(τ)∥ ≤ ν(τ) ∥G(U_{0})∥ |τ| ,  (2.9)

for every U_{0} ∈ D.

Remark 2.1. The requirement that ν(s) converges to 0 as τ → 0 implies that F is “uniformly differentiable” with respect to τ, at τ = 0, on D. The additional requirement in (2.9) regarding the factor ∥G(U_{0})∥ is necessary to obtain exponential decay in Y as t → ∞.

12 GUI-QIANG CHEN AND DAVID H. WAGNER

Definition 2.2. The explicit one-step approximation algorithm Ū(τ) = F(τ, U_{0}) is a total variation stable approximation to the solution U(τ) of the initial value problem (2.8) in the domain D ⊂ ℝ^{n} if F, \frac{∂F}{∂τ} , and \frac{∂^{2}F}{∂τ^{2}} are Lipschitz continuous on (−δ, δ)×D → ℝ^{n}.

In addition, we note one more desirable property for a one-step approximation algorithm. This property is specific to chemical reaction systems; it states that the algorithm conserves total energy just as the chemical reaction does.

Definition 2.3. The explicit one-step approximation algorithm Ū(τ) = F(τ, U_{0}) is a conservative approximation to the initial value problem

\frac{d}{dτ} (\begin{matrix}v(τ) \\ u(τ) \\ E(τ) \\ Y (τ) \end{matrix}) = (\begin{matrix}0 \\ 0 \\ qY (τ)ϕ(T (τ)) \\ −Y (τ)ϕ(T (τ)) \end{matrix}) ,  (\begin{matrix}v(0) \\ u(0) \\ E(0) \\ Y (0) \end{matrix}) = (\begin{matrix}v_{0} \\ u_{0} \\ E_{0} \\ Y_{0} \end{matrix}) .  (2.10)

if there is a function f(τ, U_{0}) ≥ 0 such that F(τ, U_{0}) has the form:

F(τ, v_{0}, u_{0}, E_{0}, Y_{0}) = (v_{0}, u_{0}, E_{0} + qf(τ, U_{0}), Y_{0} − f(τ, U_{0})) .  (2.11)

Lemma 2.1. If F(τ, U_{0}) is a conservative, consistent first order one-step approximation algorithm with f as in (2.11), then f(0, U_{0}) = 0, and

\frac{∂f}{∂τ} (0, v_{0}, u_{0}, T_{0}, Y_{0}) = Y_{0}ϕ(T_{0}).

Remark 2.2. When G is C^{2}, standard one-step approximation methods, such as Euler’s method and Runge-Kutta methods, are consistent first order, conservative, and total variation stable in this sense. However, when G is merely Lipschitz continuous, methods of order two or higher will generally not be C^{2} with respect to s. Such methods will be consistent first order, but will not meet our criteria for total variation stability. Although our definition appears to require an unusual degree of regularity in F, this regularity appears to be necessary for our proof of stability of the approximation scheme with respect to the total variation, as we will make clear in Section 3. Yet we have no proof that this degree of regularity is actually necessary for total variation stability. Furthermore, the regularity of F alone does not seem to be sufficient for global total variation stability because, as we have already remarked, some type of decay seems necessary.

In the study of traveling wave solutions of (1.1) and the corresponding Navier–Stokes model (cf. [32, 31, 3]), it is customary to modify ϕ(T) to be zero for T less than an

REACTING COMPRESSIBLE EULER EQUATIONS 13

“ignition temperature” T_{i}. The observation above indicates that it may be inadvisable to make such a modification in a manner that makes ϕ less than C^{2}—particularly in conjunction with a numerical computation of time-dependent flow.

Let F(τ, U_{0}), 0 ≤ τ ≤ h, be a conservative, consistent first order, total variation stable, one-step approximation algorithm to the initial value problem (2.10). Given U_{0}^{h}(x, t), kh ≤ t < (k + 1)h, we define U^{h}(x, t) = F(t − kh, U_{0}^{h}(x, t)).

We re-approximate U^{h}(x, (k + 1)h) using R_{j}^{k+1}:

\begin{aligned}U_{j}^{k+1} & = R_{j}^{k+1} ∘ U^{h}(x, (k + 1)h−) \\ & ≡ F(h; U_{0}^{h}((j + χ_{k+1})l, (k + 1)h−)), jl < x < (j + 2)l. \end{aligned}

We will guarantee condition (2.6) by showing that for all (x, t), U^{h}(x, t) ∈ D × [0, 1], where D is a compact set in (r, s, S)–space. It follows that \frac{h}{l} can be chosen so that the Courant-Friedrichs-Levy condition (2.6) is satisfied for all time-steps. Therefore, our approximate solutions for the Cauchy problem of the system (1.1–1.2) are defined unambiguously.

3. Total Variation Stability

In this section, we estimate the approximate solutions U^{h}(x, t) in the total variation norm and prove that the spatial total variation of the approximate solutions is uniformly bounded with respect to the mesh length h. We will measure the total variation of approximate solutions, using the sum of the absolute values of the strengths of waves in the solution of each Riemann problem occurring in the non-reacting step of the approximation scheme. In discussing this sum and its various terms, it is convenient to use a weighted ℓ_{1} norm, ∥v∥_{1} = |v_{1}| + |v_{2}| + |v_{3}| + M_{4} |v_{4}|, for a vector v ∈ ℝ^{4}, and where M_{4} > 0 is defined later. In addition, for a vector function U(x) with components (r, s, S, Y) (x), we define the total variation using the weighted ℓ_{1} norm:

Var_{rsSY} (U) = ∥(Var(r), Var(s), Var(S), Var(Y))∥_{1} .  (3.1)

Other equivalent definitions are possible; for example, one could use the Euclidean norm instead of the ℓ_{1} norm.

We will find that for each time step, the fractional step causes the total variation of the approximate solution to increase by an amount that is of order h. In any fixed time interval, the number of time-steps is of order \frac{1}{h}. Thus, when we sum these increases, we obtain a finite, non-zero change in total variation. However, this same argument leaves

14 GUI-QIANG CHEN AND DAVID H. WAGNER

us free to neglect terms of order h^{2}T V (U^{h}(·, t)) or higher, because in the limit as h → 0 the sum of these tends to 0 on any time interval.

We formalize this argument with the following lemma, which is applicable to any system of hyperbolic balance laws. In all that follows, we assume that the approximate solutions U^{h} and U_{0}^{h} remain in a compact domain D×[0, 1], where D will be defined later. Note that after the random choice step, the values of the solution between ((j +χ_{k})l, kh) and ((j + 1 + χ_{k})l, (k ± 1)h) depend only on the values of U at these points. We denote these values by (U_{L}, U_{R}). The key to our total variation estimate for the reaction step is our estimate of the change in the solution of the Riemann problem with data (U_{L}, U_{R}). Note that a bound for the total variation of U prior to the random choice step can be achieved by applying our methods to estimate the effect of the reacting step on the individual waves between (U_{L}, U_{R}). We will choose D such that

• The Riemann problem for the non-reacting equations with data (U_{L}, U_{R}) is solvable for every ordered pair (U_{gL}, U_{gR}) ∈ D × D.

• The function B_{g}(U_{gL}, U_{gR}, ϵ), which gives the wave strengths of the Riemann solution for the data (U_{gL}, U_{gR}), is Lipschitz continuous, together with its first and second derivatives, in D × D.

Lemma 3.1. Let

z = (z_{1}, z_{2}, z_{3}, z_{4}) = B(U_{L}, U_{R}, ϵ)

be the vector of signed wave strengths in the solution of the Riemann problem with data (U_{L}, U_{R}). Let F(τ, U) be a consistent first order and total variation stable approximation to the initial value problem (2.8). Let

Γ(U_{L}, z, h) = B(F(h, U_{L}), F(h, H(z, U_{L})), ϵ).  (3.2)

Then

Γ(U_{L}, z, h) = z + O(∥z∥_{1})h.

Proof. Note that Γ(U_{L}, z, 0) = z and Γ(U_{L}, 0, h) = 0. Thus we compute that

\frac{∂Γ}{∂h}(U_{L}, 0, h) = \frac{∂^{2}Γ}{∂h^{2}}(U_{L}, 0, h) = 0_{4×1}.

Let L_{1} be a Lipschitz constant for \frac{∂Γ}{∂h} with respect to the ℓ_{1} norm. Then

∥\frac{∂Γ}{∂h}(U_{L}, z, h)∥_{1} ≤ L_{1} ∥z∥_{1} ,

REACTING COMPRESSIBLE EULER EQUATIONS 15

and

∥Γ(U_{L}, z, h) − z∥_{1} = ∥∫_{0}^{1} \frac{∂Γ}{∂h}(U_{L}, z, θh)h dθ∥_{1} ≤ L_{1} ∥z∥_{1} h.

□

Lemma 3.1 implies that the fractional step increases the total variation of the approximate solution at no more than an exponential rate. For a semi-linear problem, such a rate of increase might be acceptable. However, the existence theory for entropy solutions of quasi–linear hyperbolic systems of conservation laws, which we shall apply, requires that the total variation remains bounded by a certain constant. In order to prove that the total variation remains bounded, we need more detailed estimates.

Lemma 3.2. Let L_{2} be a Lipschitz constant for \frac{∂^{2}Γ}{∂h^{2}}. Since Γ(U_{L}, 0, h) = 0, then

∥Γ(U_{L}, z, h) − (z + \frac{∂Γ}{∂h}(U_{L}, z, 0)h)∥_{1} ≤ L_{2} ∥z∥_{1} \frac{h^{2}}{2} .  (3.3)

Proof. By Taylor’s Theorem,

Γ(U_{L}, z, h) − z = \frac{∂Γ}{∂h} (U_{L}, z, 0) h − ∫_{0}^{1} \frac{∂^{2}Γ}{∂h^{2}} (U_{L}, z, θh) (θ − 1) dθ h^{2}.
□

Lemma 3.2 shows that we can estimate the increase in total variation for the reacting step by calculating first derivatives of the solution operator for the Riemann problem.

Comparing (3.3) and Lemma 3.1, we see that the term \frac{∂Γ}{∂h}(U_{L}, z, 0)h is O(∥z∥_{1})h. In order to prove our theorem, we need to show that the sum of all increases in total variation which are caused by the reaction, is bounded by a certain constant. In the next subsection, we show that ∥Y ∥_{∞} decays exponentially as t → ∞. Then we show that any increase in total variation, due to the reaction step, is proportional to Y ∥z∥_{1} h. This, together with the exponential decay of ∥Y ∥_{∞}, will enable us to derive the conditions under which the total increase of total variation is bounded.

3.1. Estimates on the Reacting Step. We first analyze the properties of a consistent first order conservative fractional–step operator for the chemical reaction. As it acts upon the conserved densities (v, u, E, Y) of (1.5), this operator takes the form

F(τ, (v_{0}^{h}, u_{0}^{h}, E_{0}^{h}, Y_{0}^{h})) = (v^{h}, u^{h}, E^{h}, Y^{h}), 0 ≤ τ < h,

where all quantities v_{0}^{h}, v^{h}, etc. are evaluated at (x, kh + τ), and

\begin{matrix}v^{h}(x, kh + τ) = v_{0}^{h}, \\ u^{h}(x, kh + τ) = u_{0}^{h}, \\ E^{h}(x, kh + τ) = E_{0}^{h}(x, kh + τ) + qf(τ, U_{0}^{h}(x, kh + τ)), \\ Y^{h}(x, kh + τ) = Y_{0}^{h}(x, kh + τ) − f(τ, U_{0}^{h}(x, kh + τ)). \end{matrix}  (3.4)

We need to estimate the change in T due to F. By Lemma 3.2, we only need to calculate

(T^{h} − T_{0}^{h}) (x, kh + τ) = \frac{∂T}{∂E} |_{U_{0}^{h}} qf_{τ}(0, U_{0}^{h}(x, kh + τ))τ = \frac{q}{e′(T)}f_{τ}(0, U_{0}^{h}(x, kh + τ))τ.

By Lemma 2.1, f_{τ}(0, U_{0}^{h}) = Y_{0}^{h}ϕ(T_{0}^{h}). Note that this estimate is still valid for a nonpolytropic gas even though c_{v} is not constant but varies with T.

A similar calculation applies to Y^{h}. Thus we have proved the following lemma.

Lemma 3.3. Let F be a consistent, conservative one-step approximation algorithm. Then, to first order in τ,

\begin{aligned}T^{h} & = T_{0}^{h} + \frac{q}{c_{v}} Y_{0}^{h}ϕ(T_{0}^{h})τ, \\ Y^{h} & = Y_{0}^{h} − Y_{0}^{h}ϕ(T_{0}^{h})τ, \end{aligned}  (3.5)

where T^{h}, Y^{h}, T_{0}^{h}, and Y_{0}^{h} are all evaluated at (x, kh + τ).

We have assumed that T_{0}^{h}(x, t) ≥ T′ > 0, 0 ≤ Y_{0}^{h}(x, t) ≤ 1, and that ϕ is Lipschitz continuous, non-negative, and increasing. In particular,

T > T′ ⟹ 0 < Φ ≤ ϕ(T) ≤ C < ∞.

Note that the Arrhenius law (1.6) satisfies these conditions. Then

Y_{0}^{h}ϕ(T_{0}^{h})τ ≥ Y_{0}^{h}Φτ.  (3.6)

Equation (3.5), or inequality (3.6), implies the following lemma.

Lemma 3.4. In the limit as h → 0, the functions Y^{h}(x, t) and T^{h}(x, t) satisfy

0 ≤ Y^{h}(x, kh + τ) ≤ Y_{0}^{h}(x, kh + τ)e^{−Φτ}, 0 ≤ τ < h,
T^{h}(x, kh + τ) ≥ T_{0}^{h}(x, kh + τ) ≥ c_{0} > 0.

Furthermore, Y^{h}(x, t) ≤ Y^{h}(x, 0)e^{−Φt} for all t ≥ 0.

Henceforth, for simplicity of exposition, we shall use this estimate in the form given by the limit as h → 0.

REACTING COMPRESSIBLE EULER EQUATIONS 17

3.2. Glimm Functional for the Fractional–Step Scheme. Our proof of the BV stability of the fractional–step scheme is based on Glimm’s method. Following that method, we define a functional on the restriction of the approximate solutions U^{h} to certain “mesh curves” J. We define a mesh point to be a point (x, t) = ((j + χ_{k})l, kh), where k ∈ ℕ and j ∈ ℤ such that j + k is even. A mesh curve J is a piecewise linear curve in the (x, t)-space, which successively connects mesh points ((j +χ_{k})l, kh) to mesh points ((j + 1 + χ_{k})l, (k ± 1)h) (see Fig. 2). We define a partial order on the set of mesh curves by stating that larger curves lie toward larger time. We call J_{2} an immediate successor of J_{1} if J_{2} connects the same mesh points as J_{1}, except for one mesh point, and if J_{2} > J_{1}.

Figure 2. The mesh curves J(k) and J(k + 1).

Let J(k) be the unique mesh curve which connects the mesh points on t = kh to the mesh points on t = (k + 1)h. Note that J(k) crosses all the waves in the Riemann solutions of U_{0}^{h}(x, t) in the strip kh ≤ t < (k + 1)h.

We will define a functional F on the set of mesh curves. The coefficients in F will depend on a set D containing the (r, s, S) values of all approximate solutions. We will assume that the initial data satisfy

Var_{rsSY} (r_{0}, s_{0}, S_{0}, Y_{0})(·) ≤ N < ∞,
(r_{0}, s_{0}, S_{0})(x) ∈ E ⊂ D.

Given values of U^{h} on adjacent mesh points at t_{k} = kh: U_{L} = U^{h}(x_{j}, t_{k}) and U_{R} = U^{h}(x_{j+2}, t_{k}), with B(U_{L}, U_{R}, ϵ) = z = (z_{1}, z_{2}, z_{3}, z_{4}), we define the wave strengths of the approximate solution U^{h}(x, t) between (x_{j}, t_{k}) and (x_{j+2}, t_{k}) as follows. For any wave p, let Var_{r}^{−}(p) denote the decreasing variation in r across p. Then, at ϵ = 0, we let α be the one–shock in z and μ be the 1–rarefaction wave in z so that

\begin{aligned}α & = \{ \begin{matrix}0,  if z_{1} ≥ 0, \\ |z_{1}| ,  if z_{1} < 0, \end{matrix} \\ μ & = z_{1} − α. \end{aligned}

Similarly, let

β = 3-shock in z,
η = 3-rarefaction in z,
δ = z_{2} = contact wave in z.

Finally, for ϵ > 0, we define (see [28], eq. (3.49)):

For 1-waves, α = Var_{r}^{−}(α) + ϵVar(α), δ_{α} = Var_{S}(α), μ = Var_{r}^{+}(μ) + ϵVar(μ),
For 3-waves, β = Var_{s}^{−}(β) + ϵVar(β), δ_{β} = Var_{S}(β), η = Var_{s}^{+}(η) + ϵVar(η),

For 2-waves, δ = S_{+} − S_{−},

For 4-waves, ζ = Y_{+} − Y_{−}.

(3.7)

Here Var(α) is taken to mean the total variation of (r, s, S) along the 1–shock curve, between the left and right states of a 1–shock wave α, and similarly for the other waves. Note that since S is monotone along the shock curve [33], Var_{S}(α) = |S_{+} − S_{−}|.

We let C denote a constant that depends only on D and ϕ, and is independent of ϵ and the mesh length h.

For any mesh curve J, Temple defined [28] (p. 144):

L_{T} (J) = ∑_{J} {(α_{i} − M_{0}δ_{α_{i}}) + (β_{i} − M_{0}δ_{β_{i}}) + M_{0} |δ_{i}|} + ϵ ∑_{J} {μ_{i} + η_{i}} + V_{0},

where α, β, δ, ζ, δ_{α}, δ_{β}, μ, and η are as defined above, for each shock wave or contact discontinuity crossing J. The constant V_{0} is the total variation of (r, s, S) along the initial data w_{0}(x). The constant M_{0} is defined in Lemmas 3.10 and 3.11 (see Lemmas 4.1 and 4.2 in [28]).

We make the following modifications of L_{T}. Let M_{4} be a constant, the value of which will be specified in the proof of Lemma 3.18. Let U_{±∞} denote the values of U at ±∞, along a mesh curve. Let α_{∞}, etc. denote the wave strengths (as defined in (3.7)) in the solution of the Riemann problem with data (U_{−∞}, U_{∞}). Let

\begin{aligned}L(J) = & ∑_{J} \{(α_{i} − M_{0}δ_{α_{i}}) + (β_{i} − M_{0}δ_{β_{i}}) + M_{0} |δ_{i}| + M_{4} |ζ_{i}| + ϵ(μ_{i} + η_{i})\} \\ & + α_{∞} + β_{∞} + μ_{∞} + η_{∞} + |δ_{∞}| + M_{4} |ζ_{∞}| . \end{aligned}

The term M_{4} |ζ_{i}| is added as a measurement of Var_{Y} (J), in order to make L(J) equivalent to Var_{rsSY} (J). In L_{T}, V_{0} was used as an upper estimate for |r_{∞} − r_{−∞}| + |s_{∞} − s_{−∞}|. However, during the reaction step, r_{±∞} and s_{±∞} change, and this change must be estimated. For this reason, instead of V_{0}, we use the wave strengths in the Riemann problem with data (U_{−∞}, U_{∞}) as a similar bound which works well with our estimates on the reaction step. Note that for any mesh curve J, L_{T} (J) ≤ L(J) + V_{0}.

Let a_{i} and b_{j}, 1 ≤ i, j ≤ 3, be two waves in an approximate solution U^{h}, which cross a mesh curve J, with a_{i} to the left of b_{j} on J. We say that the wave a_{i} approaches the wave b_{j} if either i > j (a_{i} is faster than b_{j}), or i = j ∈ {1, 3} and at least one of the waves is a shock wave—either a_{i} < 0 or b_{j} < 0.

Since ζ (that is, Y)–waves do not interact with any other waves, we can use the same definition of Q(J) that was given in [28] as follows. Let p_{i} denote arbitrary waves, q_{i} arbitrary shock or rarefaction waves, and R_{i} arbitrary rarefaction waves. We define

Q(J) = M_{1} ∑_{App} p_{i} |δ_{j}| + M_{2} ∑_{App} q_{i}R_{j} + M_{3} (∑_{App} α_{i}β_{j} + ∑_{App} α_{i}α_{j} + ∑_{App} β_{i}β_{j}) ,

where the sums are taken over all pairs of approaching waves, and p_{i} can be either to the left or right of |δ_{j}|, etc. We have corrected an apparent typographical error in [28], in which the terms for shock waves of the same field were omitted.

Finally, we define

F (J) = L(J) + ϵQ(J).

In [28], the following functionals were also defined in order to estimate Var_{rs}(J) for small ϵ:

\begin{aligned}L_{T 0}(J) & = ∑_{J} \{α_{i} + β_{i}\} + V_{0}, \\ F_{T 0}(J) & = L_{0}(J) + ϵQ(J). \end{aligned}

We modify L_{T 0} in the same way that we modified L_{T}:

L_{0}(J) = ∑_{J}{(α_{i} + β_{i}) + M_{4} |ζ_{i}| + α_{∞} + β_{∞} + μ_{∞} + η_{∞} + |δ_{∞}| + M_{4}q |ζ_{∞}|}.

20 GUI-QIANG CHEN AND DAVID H. WAGNER

Let F_{0} = L_{0} + ϵQ, F_{T} = L_{T} + ϵQ and F_{T 0} = L_{T 0} + ϵQ. Note that for any mesh curve J, F_{T} (J) ≤ F (J) + V_{0}, and F_{T 0}(J) ≤ F_{0}(J) + V_{0}.

We will see in Lemmas 3.5–3.13 and 3.18 that all of the constants M_{0}, M_{2}, M_{3}, M_{4} depend only on D.

3.3. Temple’s Results. Before we prove the total variation stability of the approximate solutions, we first review the following lemmas from Temple [28] on which our analysis depends. Lemmas 3.5 to 3.9 below correspond to Lemmas 3.1 to 3.5 of [28].

Lemma 3.5. For every compact convex open set D_{rs} in the rs-plane, there exist constants M > 0, \frac{1}{2} ≤ C_{0} < 1, and G > 1 such that, at ϵ = 0, the following estimates hold across any interaction < w_{L}, w_{M} > + < w_{M}, w_{R} > → < w_{L}, w_{R} > of states whose projections onto the rs-plane lie in D_{rs}:

α′ − α_{1} − α_{2} = A, β′ − β_{1} − β_{2} = B,

where A + B ≤ (C_{0} − 1)ξ and A = −ξ or B = −ξ. Moreover,

|δ′| − |δ_{1}| − |δ_{2}| + (δ_{α_{1}} + δ_{α_{2}} − δ_{α′}) + (δ_{β_{1}} + δ_{β_{2}} − δ_{β′}) ≤ −M(A + B),
|δ′| − |δ_{1}| − |δ_{2}| ≤ G(D_{2} + D_{3}),
μ′ − μ_{1} − μ_{2} ≤ GD_{3}, η′ − η_{1} − η_{2} ≤ GD_{3}.

Here D_{2} and D_{3} are quadratic wave interaction terms which are dominated by the decrease in Q, as is shown in [28]. Our concern is with the constant M.

Lemma 3.6. Let D = D_{rs} × [S_{∗}, S^{∗}] in rsS-space, there exists ϵ_{1} > 0 and G > 1 such that interactions are defined for every w_{L}, w_{M}, and w_{R} in D and such that, for each wave field, the following estimate holds in any interaction:

(Change in strength at ϵ) ≤ (Change in strength at ϵ = 0) + GϵQ.

Lemma 3.7. There exists an M > 0 depending only on D_{rs} and an ϵ_{1} > 0 such that, if w_{L}, w_{R} ∈ D, the associated Riemann problem is solvable for each ϵ ∈ [0, ϵ_{1}], and the waves in these Riemann problems satisfy the following estimates:

M Var_{r}^{−}(α) > Var_{rsS}(α),  Var_{s}^{−}(α) < α,
M Var_{s}^{−}(β) > Var_{rsS}(β),  Var_{r}^{−}(β) < β,
2Var_{r}^{+}(μ) > Var_{rsS}(μ),  Var_{r}^{−}(μ) = 0,
2Var_{s}^{+}(η) > Var_{rsS}(η),  Var_{s}^{−}(η) = 0,
Var_{s}(μ) < \frac{1}{4}μ,  Var_{r}(η) < \frac{1}{4}η.

Note that r and s are constant on contact waves so that Var_{rsS}(δ) = |δ| and Var_{rsS}(ζ) = |ζ|.

Remark 3.1. The above lemma imposes additional requirements on M, namely,

M ≥ 1 + \frac{Var_{S}(α)}{Var_{r}(α)} + \frac{Var_{s}(α)}{Var_{r}(α)},  (3.8)

for all 1–shock curves in D and for ϵ ∈ [0, ϵ_{1}]. The constant M must also satisfy a similar requirement with regard to 3–shock curves in D. Proposition 3.2 of [28] requires that, at ϵ = 0,

M ≥ \frac{M̃}{1 − C_{0}} ,
M̃ = 2 max_{ϵ=0} ( sup_{1−shocks} |\frac{dS}{dr} | , sup_{3−shocks} |\frac{dS}{ds} |) .

(3.9)

Here \frac{1}{2} ≤ C_{0} < 1 is a constant defined in Lemma 3.1 and Proposition 3.1 of [28]. Note that we have corrected a typographical error in [28] which defined M = (1 − C_{0})M̃.

We require, in addition to the above, that

M ≥ max_{0≤ϵ≤ϵ_{1}} ( sup_{1−shocks} |\frac{dS}{dr} | , sup_{3−shocks} |\frac{dS}{ds} |)

holds for all 1-shock and 3-shock curves in D, making M larger as needed. As a result, with M_{0} ≤ \frac{1}{2M}, α − M_{0}δ_{α} is monotone along all 1–shock curves in D and β − M_{0}δ_{β} is monotone along all 3–shock curves in D, for 0 ≤ ϵ ≤ ϵ_{1}.

Lemma 3.8. For every compact set E in rsS–space, there exists a constant 0 < C_{1} < 1 such that, for every B_{C_{1}}(w) with w in E (B_{C_{1}}(w) = ball of radius C_{1} with center w), interaction problems in B_{C_{1}}(w) are solvable for each ϵ ∈ [0, 1] with solution waves that satisfy the estimates of Lemma 3.6.

22 GUI-QIANG CHEN AND DAVID H. WAGNER

Let E be an arbitrary compact set in rsS-space. Let D̃ = D̃_{rs} × [S_{∗}, S^{∗}] be a compact set in rsS-space that contains the points within a distance C_{1} of E, where C_{1} is given by Lemma 3.8. Choose ϵ̃_{1} > 0 so that Lemmas 3.5 to 3.8 apply to D̃. Then Riemann problems ⟨w_{L}, w_{R}⟩ in D̃ are uniquely solvable if ϵ ≤ ϵ̃_{1}, or if w_{L}, w_{R} ∈ B_{C_{1}}(w) for some w ∈ E and 0 ≤ ϵ ≤ 1. Let Ṽ (w_{L}, w_{R}) denote the variation in the solution of one of these Riemann problems at any time t > 0.

Lemma 3.9. There exists a constant K_{0} > 1 such that, for (w_{L}, w_{R}) ∈ D̃ × D̃,

Ṽ (w_{L}, w_{R}) ≤ K_{0} ∥w_{L} − w_{R}∥_{1} .

Note that for the reacting equations, ζ or Y waves are decoupled from the acoustic and entropy waves so that Lemmas 3.8 and 3.9 generalize to our problem with no change to the constants C_{1} and K_{0}.

Lemmas 3.10 and 3.11 below correspond to Lemmas 4.1 and 4.2 of [28].

Lemma 3.10. Let D = D_{rs} × [S_{∗}, S^{∗}] , ϵ_{1}, and M satisfy the conditions of Lemmas 3.5 to 3.7. Let J_{1} be a mesh curve which evolves at ϵ ≤ ϵ_{1} from initial data w_{0}(x) of variation V_{0} through mesh curves J such that w(U(J)) ⊂ D. Then the following estimates hold for any M_{0} ≤ \frac{1}{2M} and any J_{2} which is an immediate successor of J_{1}:

(i) Var_{rs}(J_{2}) ≤ 20L_{T 0}(J_{2});
(ii) Var_{rsS}(J_{2}) ≤ KL_{T} (J_{2}), where K = \frac{20}{M_{0}^{2}} , M_{0} ≤ \frac{1}{2M} .

Lemma 3.10 easily generalizes to:

(i) Var_{rs}(J_{2}) ≤ 20L_{0}(J_{2});
(ii) Var_{rsSY} (J_{2}) ≤ KL(J_{2}),

since K > 1 and Var_{Y} (ζ_{i}) = M_{4} |ζ_{i}|. An examination of the proof of this Lemma (Lemma 4.1 in [28]) reveals that our substitution of the “waves at infinity” for V_{0} in L has no effect on the proof.

The next lemma is similar to the previous one, except that it concerns ϵ ∈ [0, 1].

Lemma 3.11. Let D̃ = D̃_{rs} × [S̃_{∗}, S̃^{∗}] ⊂ E, ϵ̃_{1}, M̃, and C_{1} satisfy the conditions of Lemmas 3.5 to 3.9. Let J_{1} be a mesh curve which evolves at ϵ ∈ [0, 1] from initial data of total variation V_{0}, through mesh curves J < J_{1} for which U(J) ⊂ B_{C_{1}}(w) for some w ∈ E. Then the following estimate holds for any J_{2} which is an immediate successor of J_{1}, so long as M_{0} satisfies M_{0} ≤ min(ϵ̃_{1}^{2}/2, (2M̃)^{−1}) = M̃_{0}:

Var_{rsS}(J_{2}) ≤ KL_{T} (J_{2}), K = \frac{20}{M_{0}^{2}} .

This lemma can similarly be generalized to:

Var_{rsSY} (J_{2}) ≤ KL(J_{2}).

The main theorem in [28] (Theorem 4.1) is as follows:

Lemma 3.12. Let E be any compact set in rsS-space, and let N > 1 be any positive constant. Then there exists a constant C = C(E, N) such that, for initial data w_{0}(x) ⊂ E with Var(w_{0}(x)) = V_{0} ≤ N and Y ≡ 0, if ϵV_{0} < C, then F_{T} (J) is non-increasing, F_{T} (J_{0}) < ∞, and there exists a global weak solution to (1.5a,b,c) with initial data w_{0}(x).

The proof of Lemma 3.12 is in two different settings. For small ϵ, a compact convex set D̃ = D̃_{rs} × [S̃_{∗}, S̃^{∗}] in rsS–space is chosen, together with constants ϵ_{1}, M̃_{0}, K_{0}, and C_{1}, such that D̃ contains all points within distance C_{1} of the set E (see [28], p. 142), and such that Lemmas 3.5–3.9 and Lemma 3.11 hold. Let D_{rs} be the set of points in the rs–plane within a distance of 120K_{0}N of D̃_{rs}. New values of M and C_{0} are chosen so that Lemmas 3.5–3.7 and Lemma 3.10 are satisfied on D_{rs} for ϵ ∈ [0, ϵ_{1}]. These Lemmas require M_{0} ≤ \frac{1}{2M} so that

M_{0} = min \{\frac{1}{2M} , M̃_{0}\} , K = max \{\frac{20}{M_{0}^{2}} , 1\}

are defined. Choose S_{∗} < S̃_{∗} − 6K_{0}KN, and S^{∗} > S̃^{∗} + 6K_{0}KN, and define D = D_{rs} × [S_{∗}, S^{∗}]. The mesh parameters (l, h) are chosen so that the Courant-Friedrichs-Levy condition (2.6) is satisfied for values in D. Later, a number ϵ_{0} ≤ ϵ_{1} is chosen. Then, for 0 ≤ ϵ ≤ ϵ_{0}, bounds on F imply bounds on Var_{rsS}(U^{h}(J)) which imply that w(U^{h}(J)) ⊂ D. The functional F is proved to be finite and non-increasing for solutions which remain in D.

For large ϵ, that is, ϵ ∈ [0, 1], smaller initial total variation is required, so that bounds on F ensure that Var_{rsS}(U^{h}(J)) < C_{1}. This implies that w(U^{h}(J)) ⊂ B_{C_{1}}(w̃) ⊂ D for some w̃ ∈ E — for example, we could choose w̃ = lim_{x→−∞} w_{0}(x). Note that the ball B_{C_{1}}(w̃) should be defined in the ℓ_{1} norm.

We will need an estimate of the form F (J) ≤ CVar_{rsSY} (J). For this purpose, we extract the following lemma from the proof of Lemma 3.12 [28].

24 GUI-QIANG CHEN AND DAVID H. WAGNER

Lemma 3.13. Let E be a compact set in (r, s, S) space, and let N ∈ (1, 1/ϵ] be given. Suppose that the mesh curves J and J_{0} satisfy the conditions of Lemmas 3.10 or 3.11. Then there are constants K_{1}(E, N) > 1 and ϵ_{0}(E, N) ∈ (0, 1) such that, for initial data w_{0}(x) with Var_{rsS}(w_{0}) = V_{0} < N, we have

F (J_{0}) ≤ 5K_{0}V_{0},  0 < ϵ ≤ ϵ_{0},  (3.10)
F (J_{0}) ≤ K_{1}V_{0},  ϵ_{0} < ϵ ≤ 1.  (3.11)

Moreover, for approximate solutions U^{h}(x, t), for which Var_{rsSY} (U^{h}(J)) ≤ N, we have

F (J) ≤ 5Var_{rsSY} (U^{h}(J)),  0 < ϵ ≤ ϵ_{0},  (3.12)
F (J) ≤ \frac{K_{1}}{K_{0}} Var_{rsSY} (U^{h}(J)),  ϵ_{0} < ϵ ≤ 1.  (3.13)

Proof. We denote L(U^{h}(J)) and Var(U^{h}(J)) by L(J) and Var(J). Note that by (3.7), p_{i} ≤ 2Var_{rsS} (p_{i}) for any wave p_{i}, p_{i} ≠ ζ_{i}, while, by (3.1), M_{4} |ζ_{i}| = Var(ζ_{i}). Hence, for any mesh curve J,

∑_{J, p_{i}≠ζ_{i}} (p_{i} + M_{4} |ζ_{i}|) ≤ ∑_{J} 2Var(p_{i}) = 2Var_{rsSY} (J) ,

L(J) ≤ ∑_{J, p_{i}≠ζ_{i}} (p_{i} + M_{4} |ζ_{i}|) + α_{∞} + β_{∞} + μ_{∞} + η_{∞} + |δ_{∞}| + M_{4} |ζ_{∞}| ≤ 4Var_{rsSY} (J).

(3.14)

Let K = \frac{20}{M_{0}^{2}} > 1. Choose G > 1 so that Lemmas 3.5 and 3.6 hold. Let

\begin{aligned}M_{1} & = 8G, \\ M_{2} & = 8G + 2M_{1}K(4K_{0}N)G, \\ M_{3} & = 8G + 2M_{1}K(4K_{0}N)G + 4M_{2}K(4K_{0}N)G, \\ K_{1} & = 4K_{0} + 4M_{3}K_{0}^{2}. \end{aligned}  (3.15)

Lemma 3.9 implies that Var_{rsSY} (J_{0}) ≤ K_{0}V_{0} if ϵ ≤ ϵ_{0} or if V_{0} ≤ C_{1}. Then, by (3.14),

\begin{aligned}L(J_{0}) & ≤ 4Var(J_{0}) ≤ 4K_{0}V_{0}, \\ L_{0}(J_{0}) & ≤ 4Var(J_{0}) ≤ 4K_{0}V_{0}. \end{aligned}  (3.16)

Furthermore, Q(J_{0}) ≤ M_{3} (2Var(J_{0}))^{2} ≤ 4M_{3}K_{0}^{2}V_{0}^{2} ≤ 4M_{3}K_{0}^{2}NV_{0}. Thus, for 0 ≤ ϵ ≤ 1, since ϵN ≤ 1,

F (J_{0}) ≤ 4K_{0}V_{0} + ϵ4M_{3}K_{0}^{2}NV_{0} ≤ K_{1}V_{0}.  (3.17)

We now require ϵ_{0} ≤ \frac{1}{K_{1}} . Then ϵ ≤ ϵ_{0} ⟹ ϵQ(J_{0}) ≤ 4ϵM_{3}K_{0}^{2}NV_{0} ≤ \frac{1}{2}ϵK_{1}V_{0} ≤ V_{0},

and

\begin{aligned}F (J_{0}) & = L(J_{0}) + ϵQ(J_{0}) ≤ 4K_{0}V_{0} + V_{0} ≤ 5K_{0}V_{0}, \\ F_{0}(J_{0}) & = L_{0}(J_{0}) + ϵQ(J_{0}) ≤ 5K_{0}V_{0}, \end{aligned} ϵ ≤ ϵ_{0}.  (3.18)

The proof of estimates (3.12) and (3.13) is nearly identical, except that the constant K_{0} is not needed since in this case we do not need to estimate total variation of Riemann solutions from initial data. □

Note that ϵ_{0} and M_{i} depend only on E and N. We now continue with the derivation of the constant C(E, N) from Lemma 3.12, with modifications for small ϵ based on the estimate F (J_{0}) ≤ 5K_{0}V_{0} .

For ϵ ≤ ϵ_{0}, we define

C_{0}(E, N) = \frac{1}{20M_{3}KK_{0}} min (\frac{1}{3G − 1}, 1 − C_{0}) .  (3.19)

Using Lemma 3.13, if ϵV_{0} < C_{0}(E, N), then ϵF (J_{0}) ≤ ϵ5K_{0}V_{0} < 5K_{0}C_{0}(E, N) so that

ϵF (J_{0}) < \frac{1}{4M_{3}K} min (\frac{1}{3G − 1}, 1 − C_{0}) .  (3.20)

For ϵ_{0} < ϵ ≤ 1, we proceed as follows. Let C_{2} = \frac{ϵ_{0}C_{1}}{K_{1}K}. If ϵ ≥ ϵ_{0} and ϵV_{0} < C_{2}, then KK_{1}V_{0} < C_{1}. Thus, by (3.17),

KF (J_{0}) < C_{1},  ϵ ≥ ϵ_{0}.  (3.21)

Since C_{1} ≤ 1 and K ≥ 1, \frac{C_{1}}{K} ≤ 1. This, plus K_{0}N > 1, implies that F (J_{0}) ≤ 5K_{0}N for ϵ_{0} ≤ ϵ ≤ 1. Let

C_{3} = \frac{1}{4M_{3}KK_{1}} min (\frac{1}{3G}, 1 − C_{0}) ,  (3.22)

and let C(E, N) = min(C_{2}, C_{3}). Then, if ϵV_{0} < C(E, N), we have ϵF (J_{0}) ≤ ϵK_{1}V_{0} ≤ K_{1}C(E, N) so that (3.20) holds for ϵ_{0} < ϵ ≤ 1.

3.4. TV Stability of the Reaction Step. We now prove the T V stability of the approximate solutions during the reaction step. Our total variation bounds imply bounds on the length of U^{h}(J), but we must also deal with the “drift” of the solution due to the source term G(U). We solve this problem by deriving conditions, under which the sum of the drift at one point, say, x = −∞, plus the total variation, remains less than Temple’s total variation bounds.

As we have already noted, we require a positive lower bound T′ for T in order to assure the uniform decay of the reactant.

Lemma 3.14. There exists T′(ϵ) such that, for ϵ ∈ [0, ϵ_{0}], T (w) ≥ T′(ϵ) for all w ∈ D, and, for ϵ ∈ [ϵ_{0}, 1], T (w) ≥ T′(ϵ) for all w ∈ B_{C_{1}}(w̃). Moreover, T′(ϵ) has a positive minimum for ϵ ∈ [0, 1].

Proof. Note that ln(p) = \frac{s − r}{2}, so that bounds on s − r imply bounds on ln(p). Such bounds imply that there exist constants p′ and p′′ independent of ϵ such that, for all w ∈ D, 0 < p′ ≤ p(w) ≤ p′′ < ∞. Similarly, bounds on s + r = 2u imply bounds on the velocity u. Thus, a compact set in (r, s, S)-space, considered in the (u, p, S)-coordinates, is a compact set in ℝ × (0, ∞) × ℝ.

Since the map from (u, v, S) to (r, s, S) is a C^{4} diffeomorphism (see Section 2), and T = e_{S}(v, S, ϵ) is C^{4}, we have that T is a continuous function on the compact set D × [0, 1]. Therefore, for each ϵ ∈ [0, 1], T has a minimum value T′(ϵ) on this set. For ϵ ∈ [ϵ_{0}, 1], we can take this minimum value on the smaller set B_{C_{1}}(w_{0}). Since all values of T on D are strictly positive (note, for ϵ = 0, T = \frac{1}{R}), the minimum value of T′ on [0, 1] is strictly positive. □

Note that M_{0} was defined so that, if the strength of α varies with fixed left state, α − M_{0}δ_{α} is strictly increasing in α for 0 ≤ ϵ ≤ ϵ_{1}. We require a similar property to hold for ϵ ∈ [0, 1].

Lemma 3.15. Let α and α̃ be the strengths of two 1–shock waves in D, as defined in (3.7) for a common value of ϵ ∈ [0, 1]. Let the left states of α and α̃ be U_{L} and Ũ_{L}, respectively. Then there is a constant C such that

|α − M_{0}δ_{α} − (α̃ − M_{0}δα̃)| ≤ |α − α̃| + C ∥U_{L} − Ũ_{L}∥_{1} α̃.  (3.23)

Similarly, if β and β̃ are the strengths of two 3–shock waves in D with right states U_{R} and Ũ_{R}, respectively, then there is a constant C such that

|β − M_{0}δ_{β} − (β̃ − M_{0}δβ̃)| ≤ |β − β̃| + C ∥U_{R} − Ũ_{R}∥_{1} β̃.  (3.24)

Proof. Let α_{0} be a 1–shock wave with the same left state as α and such that α_{0} = α̃ as wave strengths. Since S is monotone along any shock curve [33], Var_{S}(α) = δ_{α}. We have that

|α − M_{0}δ_{α} − (α̃ − M_{0}δα̃)| ≤ |α − M_{0}δ_{α} − (α_{0} − M_{0}δ_{α_{0}})| + M_{0} |δ_{α_{0}} − δα̃| .  (3.25)

REACTING COMPRESSIBLE EULER EQUATIONS 27

In accordance with Lemma 3.11, M_{0} is chosen so that M_{0} ≤ \frac{ϵ_{1}}{2} . Thus, for ϵ ≤ M_{0} < ϵ_{1}, since α and α_{0} have the same left state, and since α − M_{0}δ_{α} is monotone in α with derivative less than 1,

|α − M_{0}δ_{α} − (α_{0} − M_{0}δ_{α_{0}})| ≤ |α − α_{0}| = |α − α̃| .

For ϵ ≥ M_{0},

|α − M_{0}δ_{α} − (α_{0} − M_{0}δ_{α_{0}})| = |(1 + ϵ)(Var_{r}(α) − Var_{r}(α_{0})) + ϵ(Var_{s}(α) − Var_{s}(α_{0})) + (ϵ − M_{0}) (δ_{α} − δ_{α_{0}})| .

Since α_{0} and α have the same left state, each of the terms (Var_{r}(α) − Var_{r}(α_{0})), (Var_{s}(α) − Var_{s}(α̃_{0})), and (δ_{α} − δ_{α_{0}}) have the same sign. Thus,

|α − M_{0}δ_{α} − (α_{0} − M_{0}δ_{α_{0}})|
≤ |(1 + ϵ) (Var_{r}(α) − Var_{r}(α_{0})) + ϵ (Var_{s}(α) − Var_{s}(α_{0})) + ϵ (δ_{α} − δ_{α_{0}})|
= |α − α_{0}| = |α − α̃| .

Let G(Ũ_{L}, U_{L}, α̃) = δ_{α_{0}} − δα̃. Observe that G(Ũ_{L}, U_{L}, α̃) = 0 if α̃ = 0 or U_{L} = Ũ_{L}. Since G is C^{3} for α̃ ≥ 0, an argument similar to the proof of Lemmas 3.1 and 3.2 shows that there is a constant C such that |G(Ũ_{L}, U_{L}, α̃)| ≤ Cα̃∥U_{L} − Ũ_{L}∥_{1}. □

One consequence of Lemma 3.7 is the following lemma.

Lemma 3.16. Let J_{σ} be a connected segment of a mesh curve J, and let U_{σ} be the restriction of an approximate solution U^{h} to J_{σ}. Then there are positive constants c and C such that, for all U_{σ} with w(U_{σ}) ∈ D,

c Var_{rsSY} (U_{σ}) ≤ ∑_{J_{σ}} (α + β + μ + η + |δ| + |ζ|) ≤ C Var_{rsSY} (U_{σ}).

Corollary 3.1. Let J be any mesh curve, and let U_{σ} be the restriction of an approximate solution U^{h} to J. Let

L_{1}(J) = ∑_{J} (α + β + μ + η + |δ| + M_{4} |ζ|) ,

where the sum is understood to include the “extra” Riemann problem for (U_{−∞}, U_{∞}). Then there are positive constants c and C such that, for all U_{σ} with w(U_{σ}) ∈ D,

c Var_{rsSY} (U_{σ}) ≤ L_{1}(J) ≤ C Var_{rsSY} (U_{σ}).

Note that (3.7) implies that, for ϵ ≤ 1,

L(J) ≤ 2Var_{rsSY} (U_{σ}).  (3.26)

This estimate, together with the conclusions of Lemmas 3.10 and 3.11, implies that L(J) is equivalent to Var_{rsSY} (U_{σ}), just as Corollary 3.1 does for L_{1}(J).

We now make a more detailed estimate of the increase in total variation for the solution of the Riemann problem as a result of the reaction step. To be specific, we estimate the increase in F (J). Our goal is to derive conditions on F (J_{0}) that ensure that Var_{rsS}(J) satisfies the conditions imposed by Lemma 3.12 on V_{0}, namely V_{0} ≤ N and ϵV_{0} < C(E, N); these conditions are refined in Theorem 3.2. We will assume that these conditions are satisfied for all predecessors of J. Under these conditions, F does not increase in the non-reacting step. Lemma 3.11 implies that Var_{rsS}(J) ≤ KL(J). By definition, L(J) ≤ F (J). Thus, we require KF (J) ≤ N and ϵKF (J) < C(E, N). These bounds are implicitly assumed in the following lemmas.

In order to be consistent with the estimates in Section 3.3, we must work in the (r, s, S, Y)-coordinates.

Lemma 3.17. In the (r, s, S, Y)-coordinates, for the system (1.5), we have

G(U) = qc(U) · Y ϕ(T),

where c(U) is the vector (−\frac{1}{c_{v}T} , \frac{1}{c_{v}T} , \frac{1}{T} , −\frac{1}{q})^{T}. If the gas is not polytropic, we let c_{v} denote ẽ′(T).

Proof. We have that r = u − ln(p) = u − ln(R) − ln(T) + ln(v). In the reaction step de(T)/dY = c_{v}(T)dT/dY = −q, and (v, u) remain unchanged. Thus,

\frac{dr}{dY} = \frac{q}{c_{v}T} .  (3.27)

Similar calculations hold for r and S. □

We will find it convenient to denote (−\frac{1}{c_{v}T} , \frac{1}{c_{v}T} , \frac{1}{T} )^{⊤} by c_{g}(U).

In order to discuss the effect of the exothermic reaction on the functionals L and Q, it is convenient to identify a new “mesh curve” J̃ which, as a curve, is the same as a given mesh curve J, but upon which the value of U differs from the value of U on J by a single reaction step along all of J. We take J̃ to represent values before the reaction step and J to represent values after the reaction step.

REACTING COMPRESSIBLE EULER EQUATIONS 29

Lemma 3.18. Let D be given as above. Let F(τ, U) be a consistent, total variation stable, and conservative approximation algorithm for the solution of the initial value problem (2.8). There is a constant C, which depends only on D, such that if U(J_{k}) has evolved at ϵ from initial data of variation V through mesh curves J for which w(U(J)) ⊂ D, then, for each k, to first order in h,

\begin{aligned}L(J_{k}) & ≤ L(J̃_{k}) + Chq ∥Y_{0}∥_{∞} e^{−Φkh}L(J̃_{k}), \\ L_{0} (J_{k}) & ≤ L_{0}(J̃_{k}) + Chq ∥Y_{0}∥_{∞} e^{−Φkh}L(J̃_{k}). \end{aligned}

Proof. Note that

L(J_{k}) − L(J̃_{k}) = ∑_{−∞<i<∞}{α_{i} − M_{0}δ_{α_{i}} + β_{i} − M_{0}δ_{β_{i}} + M_{0} |δ_{i}| + ϵ (μ_{i} + η_{i}) + M_{4}ζ_{i}
−(α̃_{i} − M_{0}δα̃_{i} + β̃_{i} − M_{0}δβ̃_{i} + M_{0}|δ̃_{i}| + ϵ (μ̃_{i} + η̃_{i}) + M_{4}ζ̃_{i})}
+ α_{∞} + β_{∞} + μ_{∞} + η_{∞} + |δ_{∞}| + M_{4}ζ_{∞}
− (α̃_{∞} + β̃_{∞} + μ̃_{∞} + η̃_{∞} + |δ̃_{∞}| + M_{4}ζ̃_{∞}).

(3.28)

We denote the vector of the end-terms by z_{∞} = B(U_{−}, U_{+}). ¿From Lemma 3.15, this sum is less than

∑_{−∞<i≤∞}{|α_{i} − α̃_{i}| + C ∥F(h, U_{i}) − U_{i}∥ α̃ + |β_{i} − β̃_{i}| + C ∥F(h, U_{i+1}) − U_{i+1}∥ β̃
+ ||δ_{i}| − |δ̃_{i}|| + |μ_{i} − μ̃_{i}| + |η_{i} − η̃_{i}| + M_{4} ||ζ_{i}| − |ζ̃_{i}||}

≤ ∑_{−∞<i≤∞}{∥Γ(U_{i}, z_{i}, h) − Γ(U_{i}, z_{i}, 0)∥_{1} + C ∥F(h, U_{i}) − U_{i}∥ (α̃_{i} + β̃_{i−1})}.

(3.29)

In the limit as h → 0+, we obtain

L(J_{k}) − L(J̃_{k}) ≤ ∑_{−∞<i≤∞} \{∥\frac{∂Γ}{∂h}(U_{i}, z_{i}, 0)∥_{1} h + CqY_{i}ϕ (T_{i}) (α̃_{i} + β̃_{i−1})h\} . (3.30)

Note that

∑_{−∞<i≤∞} CqY_{i}ϕ (T_{i}) (α̃_{i} + β̃_{i−1})h ≤ Cq ∥Y ∥_{∞} L_{1}(J̃_{k})h.  (3.31)

We calculate (recall (3.2))

\frac{∂Γ}{∂h}(U_{i}, z_{i}, 0) = \frac{∂B}{∂U_{L}} \frac{∂F}{∂h} (0, U_{i}) + \frac{∂B}{∂U_{R}} \frac{∂F}{∂h} (0, H (U_{i}, z_{i})).  (3.32)

For brevity, we note that U_{i+1} = H(U_{i}, z_{i}). By Lemmas 3.3 and 3.17,

\frac{∂F}{∂h} (0, U) = G(U) = qc(U)Y ϕ(T).

We denote c(U_{i}) by c_{i}. Thus,

\begin{matrix}\frac{∂Γ}{∂h}(U_{i}, z_{i}, 0) = q (\frac{∂B}{∂U_{L}} c_{i}Y_{i}ϕ(T_{i}) + \frac{∂B}{∂U_{R}} c_{i+1}Y_{i+1}ϕ(T_{i+1})) \\ = q (Y_{i} (\frac{∂B}{∂U_{L}} ϕ(T_{i})c_{i} + \frac{∂B}{∂U_{R}} ϕ(T_{i+1})c_{i+1}) + (Y_{i+1} − Y_{i}) \frac{∂B}{∂U_{R}} ϕ(T_{i+1})c_{i+1}) . \end{matrix} (3.33)

As we noted in Section 2, Y is decoupled from (r, s, S) in the solution of the non–reacting Riemann problem. This means that \frac{∂B}{∂U_{L}} and \frac{∂B}{∂U_{R}} are block diagonal 4 × 4 matrices with the upper left 3×3 block relating to non-reacting gas dynamics, i.e., the derivatives of wave strengths for acoustic waves and entropy waves, with respect to r, s, and S. The remaining 1 × 1 block contains the derivative of the wave strength of the Y -contact with respect to Y —the value of this derivative is 1 for ∂B_{4}/∂Y_{R} and −1 for ∂B_{4}/∂Y_{L}. We have that

\frac{∂B}{∂U_{R}} · c_{i+1} = \frac{∂B}{∂ (r, s, S)_{R}} · c_{g(i+1)} − ( \begin{matrix}0 \\ 0 \\ 0 \\ 1 \end{matrix} ) .  (3.34)

Due to this block structure, \frac{∂B}{∂ (r, s, S)_{R}} has no Y -wave component. Thus the (r, s, S) components of (3.33) have the form:

\begin{aligned}\frac{∂Γ_{g}}{∂h} (U_{i}, z_{i}, 0) & = qY_{i} (\frac{∂B_{g}}{∂ (r, s, S)_{L}} ϕ(T_{i})c_{gi} + \frac{∂B_{g}}{∂ (r, s, S)_{R}} ϕ(T_{i+1})c_{g(i+1)}) \\ & + q (Y_{i+1} − Y_{i}) \frac{∂B_{g}}{∂ (r, s, S)_{i+1}} ϕ(T_{i+1})c_{g(i+1)}. \end{aligned}  (3.35)

Let

A(U_{i}, z_{i}) = \frac{∂B_{g}}{∂ (r, s, S)_{R}} ϕ(T_{i})c_{gi} + \frac{∂B_{g}}{∂ (r, s, S)_{R}} ϕ(T_{i+1})c_{g(i+1)}.  (3.36)

Then

\frac{∂Γ_{g}}{∂h} (U_{i}, z_{i}, 0) = qY_{i}A(U_{i}, z_{i}) + q (Y_{i+1} − Y_{i}) \frac{∂B_{g}}{∂ (r, s, S)_{R}} ϕ(T_{i+1})c_{g(i+1)}.  (3.37)

Since B is Lipschitz continuous, together with its first and second derivative, the function A(U_{L}, z) is Lipschitz continuous, together with its first derivatives. Furthermore, if z_{gi} = 0, then U_{g(i+1)} = U_{gi} and, in particular, T_{i+1} = T_{i}. Since B_{g}(U_{g}, U_{g}) is the vector of wave strengths for a Riemann problem with equal states,

A_{g}(U_{g}, 0) = ϕ(T)\frac{∂B_{g}}{∂ (r, s, S)}(U_{g}, U_{g})c_{g} = 0.

Thus, if C is a Lipschitz constant for z_{g} → A_{g}(U_{g}, z_{g}) on D taken with respect to the ℓ_{1} norm of z_{g}, the first term in (3.37) can be estimated by

∥Y_{i}A_{g}(U_{gi}, z_{gi})q∥_{1} ≤ CY_{i} ∥z_{gi}∥_{1} q.

We now examine the last term of (3.37), because it is not bounded by ∥Y ∥_{∞} ∑_{i} ∥z_{i}∥_{1}. In fact, it is the only such term in our estimate of the increase in total variation resulting from the reaction step. This term of (3.37) has the form

(Y_{i+1} − Y_{i}) \frac{∂B_{g}}{∂ (r, s, S)_{R}} c_{g(i+1)}ϕ(T_{i+1})q.  (3.38)

The fourth component of equation (3.33) is the equation for the strength of the Y - wave. This equation is

\frac{∂}{∂h} (Y_{i+1} − Y_{i}) = Y_{i}(ϕ(T_{i}) − ϕ(T_{i+1})) − (Y_{i+1} − Y_{i}) ϕ(T_{i+1}),

so that

\frac{∂}{∂h} |Y_{i+1} − Y_{i}| ≤ Y_{i}C |z_{g}| − |Y_{i+1} − Y_{i}| ϕ(T_{i+1}).

Thus, the reaction step produces possible increases in total variation, which are bounded by

CY_{i} ∥z_{gi}∥_{1} qh + ∥\frac{∂B_{g}}{∂ (r, s, S)_{R}} c_{g(i+1)}∥_{1} |Y_{i+1} − Y_{i}| ϕ(T_{i+1})qh.

The reaction step also produces a decrease in total variation for the Y component — the fourth component of \frac{∂Γ}{∂h}(U_{i}, z_{i}, 0) — in the amount |Y_{i+1} − Y_{i}| ϕ(T_{i+1})h. We now use the decrease in the Y component to offset the “bad” increase in other components proportional to |Y_{i+1} − Y_{i}|. Since \frac{∂B_{g}}{∂ (r, s, S)_{R}} is Lipschitz continuous and D is bounded, there

exists a finite upper bound M_{4} for ∥\frac{∂B_{g}}{∂ (r, s, S)_{R}} c_{g(i+1)}∥_{1}q on D. Thus the the effect of the

term (3.38) on the (r, s, S) components of \frac{∂Γ}{∂h}(U_{i}, z_{i}, 0) is bounded by M_{4} |ζ| ϕ(T_{i+1})h, and this increase is offset by a decrease in the term M_{4} |ζ|.

Thus, the change in L is estimated as follows:

L(J_{k}) ≤ L(J̃_{k}) + Cq ∥Y ∥_{∞} L_{1}(J̃_{k})h

+ ∑_{−∞<i≤∞} {(α̃_{i} + β̃_{i} + |δ̃_{i}| + μ̃_{i} + η̃_{i})CqY_{i} + M_{4}|ζ̃_{i}|ϕ(T_{i+}) − M_{4}|ζ̃_{i}|ϕ(T_{i+})} h

≤ L(J̃_{k}) + 2C∥Ỹ ∥_{∞}L_{1}(J̃_{k})qh.

(3.39)

Then, using the equivalence of L_{1}(J), L(J), and Var_{rsSY} (J), and redefining the constant C, we obtain

L(J_{k}) ≤ L(J̃_{k}) (1 + C∥Ỹ ∥_{∞}qh) .  (3.40)

Similarly,

\begin{aligned}L_{0}(J_{k}) & ≤ L_{0}(J̃_{k}) + C∥Ỹ ∥_{∞}L(J_{k})qh, \\ L_{1}(J_{k}) & ≤ L_{1}(J̃_{k}) + C∥Ỹ ∥_{∞}L(J_{k})qh. \end{aligned}  (3.41)

Furthermore, on the mesh curve J(k), by Lemma 3.4, we have

L (J_{k}) ≤ L(J̃_{k}) + Cqh ∥Y_{0}∥_{∞} e^{−Φkh}L(J̃_{k}),
L_{0} (J_{k}) ≤ L_{0}(J̃_{k}) + Cqh ∥Y_{0}∥_{∞} e^{−Φkh}L(J̃_{k}).  □

Since L_{0} ≤ L ≤ F, F_{0} ≤ F, and L is equivalent to L_{1}, it is clear that L_{1}, L_{0}, and F_{0} remain bounded as long as F remains bounded.

Henceforth, we will denote q ∥Y_{0}∥_{∞} by q̃. Now we estimate the functional F (J_{k}). Recall that

F (J_{k}) = L(J_{k}) + ϵQ(J_{k}).

We have that, to first order in h,

Q(J_{k}) ≤ Q(J̃_{k}) + M_{1} ∑_{App} {(p_{i} − p̃_{i}) |δ̃_{j}| + p̃_{i}(|δ_{j}| − |δ̃_{j}|)}

+ M_{2} ∑_{App} {(q_{i} − q̃_{i}) R̃_{j} + q̃_{i}(R_{j} − R̃_{j})}

+ M_{3} ∑_{i<j} {(α_{i} − α̃_{i} + β_{i} − β̃_{i}) (α̃_{j} + β̃_{j}) + (α̃_{i} + β̃_{j}) (α_{j} − α̃_{j} + β_{j} − β̃_{j})}

≤ Q(J̃_{k}) + M_{3} (L_{1}(J_{k}) − L_{1}(J̃_{k})) L_{1}(J̃_{k}),

where we have used M_{1} < M_{2} < M_{3} as defined in [28]. Next, using (3.41) and the equivalence of L(J_{k}), L_{1}(J_{k}), and Var_{rsSY} (J_{k}), we have

\begin{aligned}Q(J_{k}) & ≤ Q(J̃_{k}) + M_{3}C∥Ỹ_{k}∥_{∞}L_{1}(J̃_{k})^{2}qh \\ & ≤ Q(J̃_{k}) + C∥Ỹ_{k}∥_{∞}F (J̃_{k})^{2}qh. \end{aligned}

Lemma 3.19. Let J_{k} be a mesh curve between t = kh and t = (k + 1)h. Then

F (J_{k}) ≤ F (J̃_{k}) (1 + K_{3}q̃e^{−Φkh}(1 + ϵF (J̃_{k}))h) ,  (3.42)

where K_{3} is a constant depending only on D and ϕ, independent of ϵ ∈ [0, 1] and the mesh lengths l and h.

We now derive the conditions which ensure that F (J_{k}) remains bounded independent of k. The rate of increase in F (J_{k}) is nonlinear: in the limit as h → 0, F (J_{k}) approaches a solution to a nonlinear differential inequality,

\frac{dF}{dt} ≤ K_{3}q̃e^{−Φt} (F + ϵF^{2}) .

Without the coefficient e^{−Φt}, there would be no means of obtaining a uniform bound on F. However, it is easy to show that solutions of this differential inequality, with sufficiently small initial values, have bounded solutions. We show that a similar estimate holds for F independent of h.

We suppose, first of all, that a bound exists, namely, F (J_{k}) ≤ Z. Then, by (3.42), we have

F (J_{k}) ≤ F (J̃_{k}) (1 + K_{3}q̃e^{−Φkh} (1 + ϵZ) h) .

Let us abbreviate e^{−Φkh} by d^{k}. Then we have

F (J_{k}) ≤ F (J_{0}) ∏_{j=0}^{k} (1 + K_{3}hq̃d^{j} (1 + ϵZ)) .

Taking logarithms,

ln (F (J_{k})/F (J_{0})) ≤ ∑_{j=0}^{k} ln (1 + K_{3}hq̃d^{j} (1 + ϵZ))

≤ ∑_{j=0}^{k} K_{3}hq̃d^{j} (1 + ϵZ)

≤ K_{3}hq̃(1 + ϵZ) \frac{1}{(1 − d)}.

Thus we conclude

F (J_{k}) ≤ F (J_{0}) exp (\frac{K_{3}hq̃(1 + ϵZ)}{1 − e^{−Φh}}) .  (3.43)

The function f(h) = h/ (1 − e^{−Φh}) is increasing for h > 0 and tends to 1/Φ as h → 0. Thus, for h sufficiently small, we obtain

F (J_{k}) ≤ F (J_{0}) exp (\frac{2K_{3}q̃(1 + ϵZ)}{Φ}) .  (3.44)

Let K_{4} = 2K_{3}. The estimate (3.44) is valid as long as F (J_{k}) ≤ Z. The condition required for this result is that

F (J_{0}) ≤ exp (−\frac{K_{4}q̃(1 + ϵZ)}{Φ}) Z = g(Z).  (3.45)

The value of Z which maximizes g(Z) is Z = \frac{Φ}{K_{4}q̃ϵ}. Thus, our least restrictive condition on F (J_{0}) is

F (J_{0}) ≤ exp (−1 − \frac{K_{4}q̃}{Φ}) \frac{Φ}{K_{4}ϵq̃} = exp (−1 − \frac{K_{4}q ∥Y_{0}∥_{∞}}{Φ} ) \frac{Φ}{K_{4}ϵq ∥Y_{0}∥_{∞}} . (3.46)

We summarize these estimates with the following lemma.

Lemma 3.20. If F (J_{0}) satisfies (3.45), then, for all k ≥ 1, F (J_{k}) ≤ Z. In particular, if F (J_{0}) satisfies (3.46), then, for all k ≥ 1,

F (J_{k}) ≤ Z = \frac{Φ}{K_{4}ϵq ∥Y_{0}∥_{∞}} .  (3.47)

Furthermore, if F (J_{0}) satisfies (3.46), then by (3.43),

F (J_{k}) ≤ F (J_{0}) exp (\frac{K_{4}q ∥Y_{0}∥_{∞}}{Φ} + 1) .  (3.48)

REACTING COMPRESSIBLE EULER EQUATIONS 35

Next, we estimate the amount that the solution “drifts” from its original base point. We use w̃ = lim_{x→−∞} w_{0}(x) as our base point. Let D_{ϵ} = D for 0 ≤ ϵ ≤ ϵ_{0}, and D_{ϵ} = B_{C_{1}}(w̃) for ϵ_{0} < ϵ ≤ 1.

Lemma 3.21. Let U^{h}(x, t) be an approximate solution to (2.2), constructed with the fractional–step Glimm scheme, with initial data U_{0} ∈ BV (ℝ). If U^{h}(J) ∈ BV for every mesh curve J, then

U_{±∞} = lim_{x→±∞} U(x, t)

exists for all t and satisfies

\frac{dU_{±∞}}{dt} = G(U_{±∞}).

For our application (1.5), by Lemma 3.17, G(U) = Y ϕ(T)(\frac{−q}{c_{v}T} , \frac{q}{c_{v}T} , \frac{q}{T} , −1)^{T}. Let ϕ̄ = sup_{D_{ϵ}} ϕ(T).

Lemma 3.22. In the (r, s, S)-coordinates,

\begin{matrix}|r_{−∞}(t) − r_{−∞}(0)| + |s_{−∞}(t) − s_{−∞}(0)| ≤ ∫_{0}^{Y_{−∞}(0)} \frac{2q}{c_{v}T} dY, \\ |S_{−∞}(t) − S_{−∞}(0)| ≤ ∫_{0}^{Y_{−∞}(0)} \frac{q}{T} dY. \end{matrix}

As ϵ → 0, c_{v} → ∞. Recall that, by Lemma 3.14, T ≥ T′(ϵ). Let C_{4} be an upper bound for 1/c_{v} on D_{ϵ}. Denote Y_{−∞}(0) by y_{0}. Then

\begin{matrix}|r_{−∞}(t) − r_{−∞}(0)| + |s_{−∞}(t) − s_{−∞}(0)| ≤ C_{4} \frac{2y_{0}q}{T′} , \\ |S_{−∞}(t) − S_{−∞}(0)| ≤ \frac{y_{0}q}{T′} . \end{matrix}

Applying Lemma 3.22, we have

Lemma 3.23. Let U^{h}(x, t) be an approximate solution to (1.5) constructed with the fractional–step Glimm scheme, with initial data U_{0} ∈ BV (ℝ) and w(U_{0})(−∞) = w̃. Then, for all (x, t) ∈ J,

\begin{aligned}∥(r, s)^{h}(x, t) − (r, s)(w̃)∥_{1} & ≤ 2C_{4} \frac{y_{0}q}{T′} + Var_{rs}(U^{h}(J_{k})), \\ ∥(r, s, S)^{h}(x, t) − (r, s, S)(w̃)∥_{1} & ≤ (2C_{4} + 1) \frac{y_{0}q}{T′} + Var_{rsS}(U^{h}(J_{k})). \end{aligned}

36 GUI-QIANG CHEN AND DAVID H. WAGNER

Theorem 3.2. Let E be any compact set in rsS-space, let ϵ ∈ [0, 1], and let N > 1 be any positive constant. Let D, C_{1}, and K_{0} be as determined in Lemma 3.12. Let C_{0}(E, N) and C(E, N) be given by (3.19) and (3.22). Let K = 20/M_{0}^{2}. Let K_{1} and ϵ_{0} be as determined in Lemma 3.13. Let w_{0}(x) be given initial data for (1.5), expressed in (r, s, S, Y)-coordinates, and with (r, s, S, Y) values in E × [0, 1]. Let V_{0} = Var_{rsSY} (w_{0}).

If 0 ≤ ϵ < ϵ_{0}, and

V_{0} ≤ exp (−\frac{K_{4}q ∥Y_{0}∥_{∞}}{Φ} − 1) min (\frac{C_{0}(E, N)}{ϵ}, \frac{Φ}{5ϵK_{0}K_{4}q ∥Y_{0}∥_{∞}} ) ,  (3.49)
V_{0} ≤ exp (−\frac{K_{4}q ∥Y_{0}∥_{∞}}{Φ} − 1) (N − (2C_{4} + 1) \frac{y_{0}q}{5KK_{0}T′} ) ,  (3.50)
V_{0} ≤ exp (−\frac{K_{4}q ∥Y_{0}∥_{∞}}{Φ} − 1) (N − C_{4} \frac{y_{0}q}{50K_{0}T′} ) ,  (3.51)

then the fractional–step Glimm scheme constructs approximate solutions U^{h}(x, t) for which w(U^{h})(x, t) ∈ D and with uniformly bounded spatial total variation.

If ϵ_{0} ≤ ϵ ≤ 1, and

V_{0} ≤ exp (−\frac{K_{4}q ∥Y_{0}∥_{∞}}{Φ} − 1) min (\frac{6K_{0}N}{K_{1}} , \frac{C_{3}}{ϵ} , \frac{Φ}{ϵK_{1}K_{4}q ∥Y_{0}∥_{∞}} ) ,  (3.52)
V_{0} ≤ exp (−\frac{K_{4}q ∥Y_{0}∥_{∞}}{Φ} − 1) \frac{1}{KK_{1}} (C_{1} − (2C_{4} + 1) \frac{y_{0}q}{T′} ) ,  (3.53)

then the fractional–step Glimm scheme constructs approximate solutions U^{h}(x, t) for which w(U^{h})(x, t) ∈ B_{C_{1}}(w̃) and with uniformly bounded spatial total variation.

Proof. Our functionals F and F_{0} are nearly the same as Temple’s. Our modifications consist of adding the terms M_{4} |ζ| and α_{∞} + β_{∞} + μ_{∞} + η_{∞} + |δ_{∞}| + M_{4} |ζ_{∞}| − V to L. We note that none of these terms change during the non-reacting step. Therefore, estimates on F (J_{2}) − F (J_{1}) in the proof of Lemma 3.12 apply to our non-reacting step, as long as the conditions for that proof remain valid. Those conditions are:

Var_{rsS}(J) ≤ 6KK_{0}N,  Var_{rs}(J) ≤ 120K_{0}N,  (3.54)

for 0 < ϵ ≤ ϵ_{0}, to ensure (r, s, S)(x, t) ∈ D, and

Var_{rsS}(J) ≤ C_{1},  (3.55)

for ϵ_{0} < ϵ ≤ 1, to ensure (r, s, S)(x, t) ∈ B_{C_{1}}(w), and finally

ϵF (J_{0}) < \frac{1}{4M_{3}K} min (\frac{1}{3G − 1}, 1 − C_{0}) ,  F (J) ≤ 6K_{0}N,  (3.56)

for 0 < ϵ ≤ 1 to ensure that F (J) does not increase in a non-reacting step. By Lemmas 3.10 and 3.11 (and their generalizations to our problem), we have

\begin{matrix}Var_{rsSY} (J) ≤ KL(J) ≤ KF (J), & 0 ≤ ϵ ≤ 1, \\ Var_{rs}(J) ≤ 20L_{0}(J) ≤ 20F_{0}(J), & 0 ≤ ϵ ≤ ϵ_{0}. \end{matrix}

As we have noted, the drift due to chemical reactions requires additional conditions in order to ensure that the solution stays in D. Thus, using Lemmas 3.13, 3.20, and 3.23, (3.54) becomes (3.50) and (3.51), and (3.55) becomes (3.53). For J_{0}, the bounds (3.56) are implied by Lemma 3.13 if V_{0} ≤ N and either ϵV_{0} ≤ C_{0}(E, N) for 0 < ϵ ≤ ϵ_{0}, or ϵV_{0} ≤ C(E, N) for ϵ_{0} < ϵ ≤ 1. For J_{k}, the bounds (3.56) are implied for 0 < ϵ ≤ ϵ_{0} by (3.49). For ϵ_{0} < ϵ ≤ 1, (3.52) implies ϵV_{0} ≤ C_{3} and K_{1}V_{0} ≤ 6K_{0}N, or

V_{0} ≤ \frac{6K_{0}N}{8K_{0}^{2}M_{3}N} ≤ N.

For J > J_{0}, using Lemma 3.20, we find that, if

Z_{0} = min (6K_{0}N, K_{1} \frac{C_{3}}{ϵ} ) ,  (3.57)

and

F (J_{0}) ≤ exp (−\frac{K_{4}q ∥Y_{0}∥_{∞}}{Φ} − 1) Z_{0},  (3.58)

then, for all k, F (J_{k}) ≤ Z_{0} so that the bounds (3.56) are satisfied. When we express (3.57 and 3.58) in terms of V_{0}, using Lemma 3.13, we obtain (3.49) and (3.52). □

4. Convergence to the Entropy Solution

A scheme for constructing approximate solutions to a differential equation is called consistent if the convergence of the scheme (in a suitable sense) implies that the limit of the approximation is actually a solution of the differential equation. In [21] it is shown that the fractional–step scheme as constructed in Section 2.2, by using Euler’s method to approximate solutions to (2.8), is consistent with the notion of weak solutions for the nonlinear water-hammer problem. It is also shown in [9] that fractional–step schemes based on Glimm’s method and the first order approximation to (2.8) are consistent with the notion of entropy solutions. In this section we will show that a fractional–step Glimm scheme, for general hyperbolic systems of balance laws (1.4), as constructed in Section 2.2, using any consistent approximation to (2.8), is consistent with the notion of entropy solutions. That is, if a sequence U^{h_{n}}(x, t) of approximate solutions is bounded in L^{∞} and converges pointwise a.e. to a function U(x, t), then U(x, t) is an entropy solution of the

38 GUI-QIANG CHEN AND DAVID H. WAGNER

system of balance laws, which satisfies the entropy inequality. As a corollary, we obtain the existence of a global entropy solution to the exothermically reacting, compressible Euler equations (1.2) and (1.5).

Let dχ_{k} denote the uniform probability measure on [−1, 1], and let dχ denote the induced product probability measure for the random sample {χ_{k}}_{k=1}^{∞} in the Cartesian product space A = Π_{k=1}^{∞} [−1, 1].

Theorem 4.1. Suppose that

(1) The sequence U^{h_{ℓ}}(x, t) is constructed by using the Glimm fractional–step scheme, as described in Section 2.2, with a fractional–step operator which is consistent in the sense of Definition 2.1, and with the random sample {χ_{k}}_{k=1}^{∞} chosen from A.

(2) The sequence U^{h_{ℓ}}(x, t) is uniformly bounded in L^{∞} and converges pointwise a.e.

to the function U(x, t). Then there exists a null set N ⊂ A such that, for {χ_{k}} ∈ A − N, the function U(x, t) is an entropy solution of the Cauchy problem (1.4) and (2.1). That is, for any convex entropy pair (η, q) with respect to U, the following inequality

η(U)_{t} + q(U)_{x} ≤ ∇η(U)G(U)  (4.1)

holds in the sense of distributions, that is,

∫_{0}^{∞} ∫_{−∞}^{∞} (η(U)ϕ_{t} + q(U)ϕ_{x} + ∇η(U)G(U)ϕ) dxdt + ∫_{−∞}^{∞} η(U_{0}(x))ϕ(x, 0) dx ≥ 0, (4.2)

where ϕ ∈ C_{0}^{∞}((−∞, ∞) × [0, ∞)) and ϕ(x, t) ≥ 0.

Proof. For any convex entropy pair (η, q) with respect to U, we define

\begin{aligned}J (χ, h, ϕ) = − & ∫_{0}^{∞} ∫_{−∞}^{∞} (η(U^{h})ϕ_{t} + q(U^{h})ϕ_{x} + ∇η(U^{h})G(U^{h})ϕ) dxdt, \\ & − ∫_{−∞}^{∞} η(U_{0}(x))ϕ(x, 0) dx, \end{aligned}  (4.3)

where ϕ ∈ C_{0}^{∞}((−∞, ∞) × [0, ∞)) and ϕ(x, t) ≥ 0.

Note that for k ≥ 0,

\begin{matrix}− ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} & η(U^{h})ϕ_{t} dxdt \\ = & − ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} η(U_{0}^{h})ϕ_{t} dxdt − ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} (η(U^{h}) − η(U_{0}^{h})) ϕ_{t} dxdt. \end{matrix} ; (4.4)

Since U_{0}^{h} is an entropy solution of the conservation laws (2.7),

\begin{aligned}& − ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} η(U_{0}^{h})ϕ_{t} dxdt \\ ≤ & ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} q(U_{0}^{h})ϕ_{x} dxdt \\ & + ∫_{−∞}^{∞} (η(U_{0}^{h}(x, kh + 0))ϕ(x, kh) − η(U_{0}^{h}(x, (k + 1)h−))ϕ(x, (k + 1)h)) dx. \end{aligned} (4.5)

Since F is consistent in the sense of Definition 2.1, there exists a function ν(s) such that ν(s) → 0 as s → 0, and

∥F(s, U_{0}^{h}(x, t)) − U_{0}^{h}(x, t) − G(U_{0}^{h}(x, t))s∥ ≤ ν(s) ∥G(U_{0}(x, t))∥ |s| .  (4.6)

Equation (4.6) states that, as a function of s, F(s, U_{0}^{h}(x, t)) is differentiable at s = 0, and has the partial derivative G(U_{0}^{h}(x, t)). Since η is a differentiable function of U, the chain rule states that, as a function of s, η (F(s, U_{0}^{h}(x, t))) is differentiable at s = 0, and has partial derivative ∇η(U_{0}^{h}(x, t))G(U_{0}^{h}(x, t)). Therefore, there is a function e(s; x, t), which converges uniformly to 0 as s → 0, such that

η (F(s, U_{0}^{h}(x, t))) − η(U_{0}^{h}(x, t)) − ∇η(U_{0}^{h}(x, t))G(U_{0}^{h}(x, t))s = ϵ(s; x, t)s.  (4.7)

Thus, for kh ≤ t ≤ (k + 1)h,

\begin{matrix}η(U^{h}(x, t)) = & η(U_{0}^{h}(x, t)) \\ + ∇η(U_{0}^{h}(x, t))G(U_{0}^{h}(x, t))(t − kh) + ϵ(t − kh; x, t)(t − kh). \end{matrix} (4.8)

We can now compute

\begin{aligned}− & ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} (η(U^{h}) − η(U_{0}^{h})) ϕ_{t} dxdt = − ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} ϵ(t − kh; x, t)(t − kh)ϕ_{t} dxdt \\ & + ∫_{−∞}^{∞} ∫_{kh}^{(k+1)h} \frac{∂}{∂t} (∇η(U_{0}^{h}(x, t))G(U_{0}^{h}(x, t))(t − kh)) ϕ dxdt. \end{aligned} ; (4.9)

Thus,

\begin{aligned}− ∫_{kh}^{(k+1)h} & ∫_{−∞}^{∞} (η(U^{h}) − η(U_{0}^{h})) ϕ_{t} dxdt \\ = & − ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} ϵ(t − kh; x, t)(t − kh)ϕ_{t} dxdt \\ & + ∫_{−∞}^{∞} ∫_{kh}^{(k+1)h} [∇η(U_{0}^{h}(x, t))G(U_{0}^{h}(x, t)) \\ & + \frac{∂}{∂t} (∇η(U_{0}^{h}(x, t))G(U_{0}^{h}(x, t))) (t − kh)] ϕ dxdt. \end{aligned} (4.10)

Plugging (4.5) and (4.10) into (4.4), we obtain

− ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞}η(U^{h})ϕ_{t} dxdt
≤ ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} [q(U_{0}^{h})ϕ_{x} + ∇η(U_{0}^{h})G(U_{0}^{h})ϕ] dxdt

+ ∫_{−∞}^{∞} [η(U_{0}^{h}(x, kh + 0))ϕ(x, kh) − η(U_{0}^{h}(x, (k + 1)h−))ϕ(x, (k + 1)h)] dx

+ ∫_{−∞}^{∞} ∫_{kh}^{(k+1)h} \frac{∂}{∂t} [∇η(U_{0}^{h}(x, t))G(U_{0}^{h}(x, t))] (t − kh)ϕ(x, t) dxdt

− ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} ϵ(t − kh; x, t)(t − kh)ϕ_{t} dxdt.

Summing over k, we have

J (χ, h, ϕ) ≤ I(χ, h, ϕ) + ∑_{k=0}^{∞} R_{j}(χ, h, ϕ) + ∑_{k=0}^{∞} D_{j}(χ, h, ϕ),

where

I(χ, h, ϕ) = ∑_{k=0}^{∞} I_{k}(χ, h, ϕ),

I_{0}(χ, h, ϕ) = ∫_{−∞}^{∞} (η(U^{h}(x, 0)) − η(U_{0}(x))) ϕ(x, kh) dx,

I_{k}(χ, h, ϕ) = ∫_{−∞}^{∞} (η(U_{0}^{h}(x, kh + 0)) − η(U_{0}^{h}(x, kh−))) ϕ(x, kh) dx, k = 1, ⋯,

R_{k}(χ, h, ϕ) = ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} (q(U_{0}^{h}) − q(U^{h})) ϕ_{x} dxdt

+ ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} (∇η(U_{0}^{h})G(U_{0}^{h}) − ∇η(U^{h})G(U^{h})) ϕ dxdt

− ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} ϵ(t − kh; x, t)(t − kh)ϕ_{t} dxdt,

and

D_{k}(χ, h, ϕ) = ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} \frac{∂}{∂t}(∇η(U_{0}^{h}(x, t))G(U_{0}^{h}(x, t)))(t − kh)ϕ(x, t) dxdt.

We now analyze each of these components for convergence to zero as h tends to zero.

Lemma 4.1. There is a null subset N ⊂ A and a positive sequence h_{m}, which converges to zero, such that, for χ ∈ A − N, the functionals I(χ, h_{m}, ϕ) converge weakly to zero as measures on compact subsets of ℝ × ℝ^{+}.

Proof. Note that I_{0} → 0 since U^{h}(x, 0) was chosen to converge to U_{0}(x) in L^{1}(ℝ). Since

I_{k}(χ, h, ϕ) = ∑_{j=−N}^{N} ∫_{(j−1)l}^{(j+1)l} (η(U_{0}^{h}((j + χ_{k})l, kh−)) − η(U_{0}^{h}(x, kh−))) ϕ(x, kh))dx,

and

|I_{k}(χ, h, ϕ)| ≤ CT V (U_{0}^{h}(x, kh−))h ∥ϕ∥_{∞} ,  (4.11)

we have

|I(χ, h, ϕ)| ≤ CLd ∥ϕ∥_{∞} ,  (4.12)

where L is an upper bound for T V (U_{0}^{h}(·, (k−1)h)), and d is the diameter of supp ϕ(x, t). Thus ϕ → I(χ, h, ϕ) is a distribution of order zero. By the Riesz representation theorem, I(χ, h, ϕ) corresponds to a Radon measure [12].

Let t_{k} denote kh−. Notice that

\begin{matrix}|∫_{−1}^{1} I_{k}(χ, h, ϕ) dχ_{k}| \\ = | ∑_{j=−N}^{N} ∫_{−1}^{1} ∫_{(j−1)l}^{(j+1)l} (η(U_{0}^{h}((j + χ_{k})l, t_{k})) − η(U_{0}^{h}(x, t_{k}))) ϕ(x, t_{k}) dxdχ_{k}| \\ ≤ | ∑_{j=−N}^{N} ϕ_{j}^{k} ∫_{−1}^{1} ∫_{(j−1)l}^{(j+1)l} (η(U_{0}^{h}((j + χ_{k})l, t_{k})) − η(U_{0}^{h}(x, t_{k}))) dxdχ_{k}| \\ + ∥∇η∥_{∞} ∥ϕ_{x}∥_{∞} l ∑_{j=−N}^{N} ∫_{−1}^{1} ∫_{(j−1)l}^{(j+1)l} |U_{0}^{h}((j + χ_{k})l, t_{k}) − U_{0}^{h}(x, t_{k})| dxdχ_{k} \\ ≤ C ∥∇η∥_{∞} ∥ϕ_{x}∥_{∞} l ∫_{−1}^{1} ∫_{−l}^{l} ∑_{j=−N}^{N} |U_{0}^{h}((j + χ_{k})l, t_{k}) − U_{0}^{h}(jl + x, t_{k})| dxdχ_{k} \\ ≤ C ∥∇η∥_{∞} ∥ϕ_{x}∥_{∞} T V (U_{0}^{h}(·, t_{k}))h^{2}, \end{matrix} (4.13)

where ϕ_{j}^{k} = ϕ(jl, t_{k}). Therefore, for 0 < k_{1} < k_{2}, we have

\begin{aligned}|∫ I_{k_{1}}(χ, h, ϕ)I_{k_{2}}(χ, h, ϕ)dχ| & = |∫ (∫ I_{k_{2}}(χ, h, ϕ)dχ_{k_{2}}) I_{k_{1}}(χ, h, ϕ)dχ̄_{k_{2}} | \\ & ≤ Ch^{3} ∥ϕ∥_{∞} ∥ϕ_{x}∥_{∞} , \end{aligned} (4.14)

where χ̄_{k_{2}} denotes the sequence obtained from χ by deleting the k_{2}-th element. Moreover, we have from (4.11) that

∫ (I_{k}(χ, h, ϕ))^{2}dχ ≤ Ch^{2} ∥ϕ∥_{∞}^{2} .  (4.15)

Thus,

\begin{aligned}∫ (I(χ, h, ϕ))^{2}dχ & = 2 ∑_{k_{1}<k_{2}} ∫ I_{k_{1}}I_{k_{2}}dχ + ∑_{k} ∫ I_{k}^{2}dχ \\ & ≤ Ch ∥ϕ∥_{∞} ∥ϕ_{x}∥_{∞} + Ch ∥ϕ∥_{∞}^{2} , \end{aligned}  (4.16)

so that I(χ, h, ϕ) → 0 in L^{2}(A, dχ). As a consequence, there exists a sequence h_{m} → 0 such that I(χ, h_{m}, ϕ) → 0 pointwise a.e. in (A, dχ). In other words, for any ϕ ∈ C_{0}^{∞}((−∞, ∞) × [0, ∞)), there exists a sequence h_{m} → 0 and a null subset N ⊂ A such that, if χ ∈ A − N, then

I(χ, h_{m}, ϕ) → 0, as h_{m} → 0.  (4.17)

REACTING COMPRESSIBLE EULER EQUATIONS 43

Let {ϕ_{ℓ}}_{ℓ∈Λ} be a countable dense set in C_{0}^{∞}((−∞, ∞) × [0, ∞)) with respect to the norm ∥ ∥_{∞}. Then, after refining h_{m} above by a diagonal argument and after taking the appropriate countable union of null sets N_{ℓ}, we can conclude that there exists a sequence h_{m} → 0 and a null subset N ⊂ A such that, if χ ∈ A − N , ℓ ∈ Λ, then

I(χ, h_{m}, ϕ_{ℓ}) → 0, as m → ∞.  (4.18)

Let K be a compact subset of ℝ×ℝ^{+}. We see from (4.12) that I(χ, h_{m}, ϕ) is a uniformly bounded (i.e., Lipschitz continuous) sequence of linear functionals on C_{0} [K]. We have seen that this sequence converges pointwise to zero on a dense subset of this space. Hence the sequence must converge pointwise to zero on all of C_{0} [K] ∩ C^{1} [K]. This means that the functionals I(χ, h_{m}, ϕ) converge weakly to zero as measures on K. □

Since, for (k − 1)h ≤ t < kh, |U^{h}(x, t) − U_{0}^{h}(x, t)| ≤ C(t − (k − 1)h) ≤ Ch, we have

∑_{k=1}^{∞} R_{k}(χ, h, ϕ) ≤ Ch (∥ϕ_{x}∥_{1} + ∥ϕ∥_{1} + ∥ϕ_{t}∥_{1}) .  (4.19)

We now estimate

D_{k}(χ, h, ϕ) = ∫_{kh}^{(k+1)h} ∫_{−∞}^{∞} \frac{∂}{∂t}(∇η(U_{0}^{h}(x, t))G(U_{0}^{h}(x, t)))(t − kh)ϕ(x, t) dxdt.

Let H(U_{0}^{h}) denote the derivative of ∇η · G, evaluated at U_{0}^{h}.

Then

\begin{aligned}\frac{∂}{∂t} (∇η(U_{0}^{h})G(U_{0}^{h})) (x, t) & = Ĥ(U_{0}^{h}(x, t)) · U_{0}^{h}_{t}(x, t) = −Ĥ(U_{0}^{h}(x, t)) · F (U_{0}^{h}(x, t))_{x} \\ & = −Ĥ(U_{0}^{h}(x, t)) · ˆ∇F(U_{0}^{h})(U_{0}^{h}(x, t))_{x}, \end{aligned}

by the chain rule of differentiation. Here Ĥ(U^{h}(x, t)) denotes Vol’pert’s “functional superposition” [29], which is given by

Ĥ(U^{h}(x, t)) = ∫_{0}^{1} H (θU^{h}(x, t−) + (1 − θ)U^{h}(x, t + 0)) dθ,  (4.20)

at all points (x, t) at which U^{h} is either approximately continuous or has a well-defined jump discontinuity with a normal direction not parallel to the x-axis. For arbitrary functions f ∈ BV [ℝ × ℝ^{+}; ℝ^{n}] and g ∈ C [ℝ^{n}; ℝ], Ĥ∘f is well-defined almost everywhere with respect to one–dimensional Hausdorff measure.

Thus, if C_{1} is an L_{∞} matrix bound for Ĥ(U)ˆ∇F(U),

D_{k}(χ, h, ϕ) ≤ C_{1}TV (U_{0}^{h}) ∥ϕ∥_{∞} h^{2}.  (4.21)

Let C_{2} be the diameter of {t | ∃x such that ϕ(x, t) ≠ 0}. Then

∑_{k=1}^{∞} D_{k}(χ, h, ϕ) ≤ C_{1}C_{2} ∥ϕ∥_{∞} h.  (4.22)

Therefore, we can conclude from (4.19), (4.22), and Lemma 4.1 that, if χ ∈ A and ϕ ∈ C_{0} [K] ∪ C^{1} [K], then

J (χ, h_{m}, ϕ) → 0, as m → ∞.

Furthermore, we conclude from Theorems 3.1 and 3.2 that the approximate solution sequence {U^{h}(x, t)} converges strongly in L^{1} to a bounded variation function U(x, t) for any rational number t ∈ (0, ∞). Observing that U^{h}(x, t) has a finite speed of propagation and using the uniform estimate on the total variation in x of U^{h}, we conclude that {U^{h}(x, t)} is L^{1}-Lipschitz continuous in t, from which follows the convergence of U^{h} to the function U for all t.

Hence, for any convex entropy pair (η, q),

∫_{0}^{∞} ∫_{−∞}^{∞}(η(U)ϕ_{t} + q(U)ϕ_{x} + ∇η(U)G(U)ϕ)dxdt + ∫_{−∞}^{∞} η(U(x, 0))ϕ(x, 0)dx ≥ 0,

where ϕ ∈ C_{0}^{∞}((−∞, ∞) × [0, ∞)) and ϕ(x, t) ≥ 0.

In particular, we choose (η(U), q(U)) = ±(U_{j}, F_{j}(U)), 1 ≤ j ≤ n, in (4.1) and conclude that U satisfies the equations (1.4) in the sense of distributions. □

As a direct corollary to Theorem 4.1, we obtain

Theorem 4.2. If U(x, t) = (v, u, e + u^{2}/2, Y) is an entropy solution of (1.5), then

S_{t} ≥ \frac{qϕ(T)Y}{T} ,  (4.23)
(\frac{Y^{2}}{2} )_{t} + ϕ(T)Y^{2} ≤ 0,  (4.24)

in the sense of distributions.

In fact, one can check that both (−S(v, e), 0) and (Y^{2}/2, 0) are convex entropy pairs by a careful calculation. Then (4.23) and (4.24) follow from (4.1).

Note that due to the lower-order term, which represents the reaction rate in our case, the form of the entropy condition is somewhat different from that which is customary for homogeneous hyperbolic systems of conservation laws. Let η(U) be a convex “entropy” function (a.k.a. convex extension) for a hyperbolic system of balance laws

U_{t} + F (U)_{x} = G(U).

REACTING COMPRESSIBLE EULER EQUATIONS 45

After multiplying by ∇η(U), one obtains the following extra balance law for classical solutions:

η(U)_{t} + q(U)_{x} = ∇η(U) · G(U).  (4.25)

For discontinuous weak solutions, the corresponding entropy condition is (4.1). As long as G(U) is a locally integrable function of (x, t), and not a singular measure, the righthand side of (4.1) has no bearing on the admissibility of shock waves. However, if one is interested in the proper expression of the increase of physical entropy S (η(U) = −S), then the inequality (4.1) takes on more significance.

5. Conclusions and Remarks

We conclude with several observations. First, we note that the lower-order term representing the chemical reaction causes the solution to drift, and this drift, in turn, imposes a new requirement, namely (3.51) for 0 < ϵ ≤ ϵ_{0}, or (3.53) for ϵ_{0} < ϵ ≤ 1, in order to keep the solution values inside a given range. However, these conditions become redundant if the lower-order term tends to zero as x → ∞ or as x → −∞–for example, if Y_{0}(x) has compact support or vanishes at either ∞ or −∞.

Secondly, we would like to be able to state a simple result in the same form as that obtained in [23], namely that ϵVar(U_{0}) ≤ C guarantees the existence of solutions for the Cauchy problem. Unfortunately, we are not yet able to do so. The basic reason for this is that the results of [28] do not have this form. Those results required:

V_{0} = Var(U_{0}) ≤ min (N, \frac{C(E, N)}{ϵ}) .  (5.1)

We recall from [28] that

C(E, N) = \frac{1}{K_{1}K} min (ϵ_{0}C_{1}, \frac{1}{4M_{3}} min (\frac{1}{3G}, 1 − C_{0})) ,
M_{3} = G (8 + 64KK_{0}GN) (1 + 16KK_{0}GN) ,
K_{1} = 8K_{0}^{2}M_{3}N.

Thus C(E, N) ≈ \frac{1}{N(A + BN)^{4}} . For ϵ small and V_{0} large, we can let N = V_{0} and (5.1) becomes:

ϵV_{0}^{2} (A + BV_{0})^{4} ≤ 1.

However, our refinement of the bounds on F (J_{0}) in terms of V_{0} in Lemma 3.13 improves this requirement somewhat for small ϵ. This Lemma states that for ϵ ≤ ϵ_{0} F (J_{0}) ≤

46 GUI-QIANG CHEN AND DAVID H. WAGNER

5K_{0}V_{0} ≤ 5K_{0}N. This estimate and (3.20) are the estimates required in [28] to prove that F (J) is non-increasing.

Thus, for ϵ small, we have C(E, N) ≈ \frac{1}{(A + BN)^{2}} , so that uniform total variation bounds apply to approximate solutions of the non-reacting equations when

ϵV_{0} (A + BV_{0})^{2} ≤ 1,

or, ignoring lower-order terms,

ϵBV_{0}^{3} ≤ 1,

Remark 5.1. One would like to think that a more detailed analysis would improve these results to solutions such as steady strong detonation waves, steady weak deflagration waves, and perturbations of such waves. Such steady waves can only exist when there is some T_{i} > 0 such that ϕ(T) = 0 for T < T_{i}, which means that the reacting rate function is discontinuous. In this case, T_{i} is called an ignition temperature. The results we present here do not establish a theory for this class of solutions. However, an existence theory for such solutions with large initial data has been established for compressible Navier-Stokes models of combustion in [3, 6]. The existence of traveling wave solutions, and the “ZND”, or vanishing viscosity limit of such solutions, was discussed in [32, 31].

As an alternative to the ignition temperature assumption, one can study steady detonation or deflagration wave solutions to an initial–boundary value problem, where fresh reactant is supplied through a boundary value condition. However, the techniques developed in this paper are not directly applicable to this problem, because our method depends on the uniform decay of the reactant. It may be possible, however, that this difficulty can be overcome by analyzing the increase in total variation of the solution due to finite quantities of reactant, and observing that a given finite mass of reactant still decays to zero as it flows through the reaction zone—thus the damage done to a variation estimate is limited. Yet there is still a continuing increase in total variation, and the only apparent way to offset this increase is through decreases in total variation resulting from shock-rarefaction interactions. It would be interesting to make such an extension to the compressible Euler equations with such discontinuous rate function.

Remark 5.2. In this paper the large-time behavior of the solution is not explored. It would be interesting to investigate the asymptotic behavior of the generalized solutions to the Cauchy problem (1.5), or (1.2), and (2.1). In this context we refer the reader to Liu [20] for linear and nonlinear large-time behavior of solutions of general systems of conservation laws, and to Glimm-Lax [16], DiPerna [11], Dafermos [8] and references cited therein for decay results of solutions to systems with two conservation laws in the context of the Glimm scheme. We also refer to Chen-Frid [5, 4] for new analytical frameworks developed recently.

Acknowledgments. The research was initiated when the authors visited IMA at Minnesota with funding provided by the National Science Foundation. Gui-Qiang Chen’s research was supported in part by the National Science Foundation grant DMS-9971793, INT-9987378, and INT-9726215. David Wagner’s research was supported in part by by the Texas Coordinating Board for Higher Education, Advanced Research Program, and by the Institute for Mathematics and its Applications with funds provided by the National Science Foundation. The calculation depicted in Fig. 1 was performed by David Wagner using a fractional–step method based on the Sanders-Weiser scheme [25].

References

[1] Anne Bourlioux, Andrew Majda, and Victor Roytburd. Theoretical and numerical structure for unstable one–dimensional detonations. SIAM J. Appl. Math., 51:303–343, 1991.

[2] Tung Chang and Ling Hsiao. Riemann problem for one dimensional adiabatic flow without convexity. Acta Math. Sinica, 22:719–32, 1979.

[3] Gui-Qiang Chen. Global solutions to the compressible Navier-Stokes equations for a reacting mixture. SIAM J. Math. Anal., 23:609–634, 1992.

[4] Gui-Qiang Chen and Hermano Frid. Decay of entropy solutions of nonlinear conservation laws. Arch. Rat. Mech. Anal., 146:95–127, 1999.

[5] Gui-Qiang Chen and Hermano Frid. Uniqueness and asymptotic stability of riemann solutions for the compressible Euler equations. Transactions of the American Mathematical Society, 353:1103– 1117, 2001.

[6] Gui-Qiang Chen, David Hoff, and Konstantina Trivisa. Well–posedness and asymptotic behavior of solutions to the reacting, compressible Navier–Stokes equations with large discontinuous initial data. Preprint, January 2001.

[7] Richard Courant and Kurt O. Friedrichs. Supersonic Flow and Shock Waves, volume 12 of Applied Mathematical Sciences. Wiley–Interscience, New York, 1948.

[8] Constantine M. Dafermos. Generalized characteristics in hyperbolic systems of conservation laws. Arch. Rat. Mech. Anal., 107:127–155, 1989.

[9] Constantine M. Dafermos and Ling Hsiao. Hyberbolic systems of balance laws with inhomegeneity and dissipation. Indiana Univ. Math. J., 31:471–491, 1982.

[10] X. Ding, Gui-Qiang Chen, and P. Luo. Convergence of the fractional step Lax-Friedrichs scheme and Godunov scheme for the isentropic gas dynamics. Comm. Math. Phys., 121:63–84, 1989.

[11] Ronald J. DiPerna. Decay and asymptotic behavior of solutions to nonlinear hyperbolic systems of conservation laws. Indiana Univ. Math. J., 24:1047–1071, 1975.

[12] William F. Donoghue, Jr. Distributions and Fourier Transforms, volume 32 of Pure and Applied Mathematics. Academic Press, New York and London, 1969.

[13] J.J. Erpenbeck. Stability of idealized one-reaction detonations. Phys. Fluids, 7:684–696, 1964.

[14] W. Fickett and W.W. Wood. Flow calculation for pulsating one-dimensional detonation. Phys. Fluids, 9:903–916, 1966.

[15] James Glimm. Solutions in the large for nonlinear hyperbolic systems of equations. Comm. Pure Appl. Math., 18:697–715, 1965.

[16] James Glimm and Peter Lax. Decay of solutions of systems of nonlinear hyperbolic conservation laws. Mem. Amer. Math. Soc., 101, 1970.

[17] Peter D. Lax. Hyperbolic systems of conservations laws II. Comm. Pure Appl. Math., 101:537–566, 1957.

[18] H.I. Lee and D. Scott Stewart. Calculation of linear detonation instablility: one dimensional instability of plane detonation. J. Fluid Mech., 216:103–132, 1990.

[19] Tai-Ping Liu. Solutions in the large for the equations of nonisentropic gas dynamics. Indiana U. Math. J., 26:147–177, 1977.

[20] Tai-Ping Liu. Quasilinear hyperbolic systems. Comm. Math. Phys., 68:141–172, 1979.

[21] M. Luskin and J.B. Temple. The existence of a global weak solution to the nonlinear waterhammer problem. Comm. Pure Appl. Math., 34:697–735, 1982.

[22] Andrew Majda. The Stability of Multi-Dimensional Shock Fronts—A New Problem For Linear Hyberbolic Equations, volume 275 of AMS Memoirs. American Mathematical Society, Providence, RI, 1983.

[23] Takaaki Nishida. Global solution for an initial boundary value problem of a quasilinear hyperbolic system. Proc. Japan Acad., 44:642–646, 1968.

[24] A.K. Oppenheim and R.I. Soloukhin. Experiments in gasdynamics of explosions. Ann. Rev. Fluid Mech., 5:31–58, 1973.

[25] Richard Sanders and Alan Weiser. A high resolution staggered mesh approach for nonlinear hyperbolic systems of conservation laws. Jour. Comp. Phys., 101:314–329, 1992.

[26] R. Smith. The Riemann problem in gas dynamics. Trans. Amer. Math. Soc., 249:1–50, 1979.

[27] Joel A. Smoller. Shock Waves and Reaction-Diffusion Equations, volume 258 of Grundlehren der mathematischen Wissenschaften. Springer-Verlag, 1982.

[28] J. Blake Temple. Solutions in the large for the nonlinear hyperbolic conservation laws of gas dynamics. Jour. Diff. Eqs., 41:96–161, 1981.

[29] A. I. Vol’pert. The spaces BV and quasilinear equations. Math. USSR–Sbornik, 2:225–267, 1967.

[30] David H. Wagner. Equivalence of the Euler and Lagrangian equations of gas dynamics for weak solutions. Jour. Diff Eq., 68:118–136, 1987.

[31] David H. Wagner. Detonation waves and deflagration waves in the one dimensional ZND model for high Mach number combustion. In W. Brent Lindquist, editor, Current Progress in Hyberbolic Systems: Riemann Problems and Computations, volume 100 of Contemporary Mathematics, pages 277–285, Providence, 1989. American Mathematical Society.

[32] David H. Wagner. The existence and behavior of viscous structure for plane detonation waves. SIAM J. Math. Anal., 20:1035–1054, 1989.

[33] Hermann Weyl. Shock waves in arbitrary fluids. Comm. Pure Appl. Math., 2:102–122, 1949.

[34] G. B. Whitham. Linear and Nonlinear Waves. Pure & Applied Mathematics. Wiley-Interscience, 1974.

[35] L.-A. Ying and C.-H. Wang. Global solutions of the Cauchy problem for a nonhomogeneous quasilinear hyperbolic system. Comm. Pure Appl. Math., 33:579–597, 1980.

[36] L.-A. Ying and C.-H. Wang. Solution in the large for nonhomogeneous quasilinear hyperbolic systems of equations. Journal of Mathematical Analysis and Applications, 78:440–454, 1980.

(G.-Q. Chen) Department of Mathematics, Northwestern University, Lunt Hall, 2033 Sheridan Road, Evanston, IL 60208-2730, USA

E-mail address, G.-Q. Chen: gqchen@math.northwestern.edu

(D. H. Wagner) Department of Mathematics, University of Houston, Houston, TX 77204-3008, USA

E-mail address, D. H. Wagner: DWagner@uh.edu