\documentclass[12pt,reqno]{amsart}
\newif\ifpdf
\ifx\pdfoutput\undefined
\pdffalse % we are not running PDFLaTeX
\else
\pdfoutput=1 % we are running PDFLaTeX
\pdftrue
\fi

\ifpdf
\usepackage[pdftex]{graphicx}
\else
\usepackage{graphicx}
\fi

\usepackage{psfig}
\usepackage{graphics}
\usepackage{verbatim}
\textwidth=6.2in
\textheight=8in
\advance\oddsidemargin by -.5 true in
\advance\evensidemargin by -.5 true in
\headheight=.5in
\numberwithin{equation}{section}
\newtheorem{thm}{Theorem}[section]

\newtheorem{lem}{Lemma}[section]
\newtheorem{prop}[thm]{Proposition}
\newtheorem{cor}[thm]{Corollary}

\theoremstyle{definition}
\newtheorem{mydef}{Definition}[section]
\newtheorem*{condition}{Conditions (H)}
\theoremstyle{remark}
\newtheorem{rem}{Remark}[section]
%\numberwithin{rem}{section}
%\numberwithin{thm}{section}
%\numberwithin{mydef}{section}
%\numberwithin{lem}{section}
%\numberwithin{cor}{section}
\renewcommand{\baselinestretch}{1.2}

\newcommand{\tst}{\textstyle}
\newcommand{\myf}[2]{\frac{\tst #1}{\tst{#2}}}
%\renewcommand{\tfrac}[2]{\frac{\tst#1}{\tst#2}}
\newcommand{\F}{\mathcal F}
\newcommand{\D}{\mathcal D}
\newcommand{\DD}{\mathbf D}
\newcommand{\K}{\tilde{K}}
\newcommand{\RH}{Rankine-Hugoniot }
\providecommand{\n}[1]{\left\lvert#1\right\rvert}
\providecommand{\bign}[1]{\bigl\lvert#1\bigr\rvert}
\providecommand{\Bign}[1]{\Bigl\lvert#1\Bigr\rvert}
\newcommand{\N}[1]{\left\lVert#1\right\rVert}
\newcommand{\BN}[1]{\Bigl\lVert#1\Bigr\rVert}
\newcommand{\bN}[1]{\bigl\lVert#1\bigr\rVert}
\newcommand{\p}[1]{\left(#1\right)}
\newcommand{\rtvert}[1]{\left.#1\right\rvert}
\newcommand{\bp}[1]{\bigl(#1\bigr)}% Used for parens smaller than
% the automatic ones given by \p
\newcommand{\Bp}[1]{\Bigl(#1\Bigr)}% Used for big parens
\newcommand{\bc}{\mathbf c}
\newcommand{\ba}{\mathbf a}
\newcommand{\bb}{\mathbf b}
\newcommand{\bs}{\mathbf s}
\newcommand{\ubs}{S_*}
\newcommand{\bt}{\mathbf t}
\newcommand{\bz}{\mathbf z}
\newcommand{\tq}{\tilde{q}}
\newcommand{\supnorm}[1]{\N{#1}_{\infty}}
\newcommand{\Bsupnorm}[1]{\BN{#1}_{\infty}}
\newcommand{\bsupnorm}[1]{\bN{#1}_{\infty}}
\newcommand{\onenorm}[1]{\N{#1}_{1}}
\newcommand{\bigonenorm}[1]{\bN{#1}_{1}}
\newcommand{\Bigonenorm}[1]{\BN{#1}_{1}}
%\def \ubar{\bar{u}}
\newcommand{\Ubar}{\bar{U}}
\newcommand{\g}{\gamma}
\newcommand{\la}{\lambda}
\newcommand{\R}{\mathbb R}
\newcommand{\Rn}{\R^n}
\newcommand{\br}[1]{\left[#1\right]}
\newcommand{\bigbr}[1]{\bigl[#1\bigr]}
\newcommand{\mybrace}[1]{\left\{#1\right\}}
\newcommand{\bigbrace}[1]{\bigl\{#1\bigr\}}
\newcommand{\angbr}[1]{\left\langle#1\right\rangle}
%\newcommand{\Su}[3]{\sum_{#1}^{#2}#3}
%\newcommand{\Chi}{\raise 2pt\hbox{$\chi$}}
\newcommand{\grad}{\nabla}
\newcommand{\tv}{{\operatornamewithlimits{TV}}\,}%
\newcommand{\supp}{{\operatornamewithlimits{supp}}\,}%
\newcommand{\eps}{\epsilon}%
\newcommand{\vareps}{\varepsilon}
%\newcommand{\v}{\varphi}
%partial derivative of #1 with respect to #2
\newcommand{\pd}[2]{\myf{\partial #1}{\partial #2}}
\newcommand{\pdtwo}[2]{\myf{\partial^2 #1}{\partial #2^2}}
\newcommand{\pdtwomix}[3]{\myf{\partial^2 #1}{\partial #2\partial #3}}
\newcommand{\talpha}{\tilde\alpha}
\newcommand{\tbeta}{\tilde\beta}
\newcommand{\tdelta}{\tilde\delta}
\newcommand{\tzeta}{\tilde\zeta}
\newcommand{\teta}{\tilde\eta}
\newcommand{\tmu}{\tilde\mu}
\newcommand{\tU}{\tilde{U}}
\newcommand{\teps}{\tilde\eps}
\newcommand{\tM}{\tilde{M}}
\newcommand{\tS}{\tilde{S}^*}
\newcommand{\utS}{\tilde{S}_*}
\newcommand{\Var}{{\operatornamewithlimits{Var}}}%
\begin{document}
\ifpdf
        \DeclareGraphicsExtensions{.pdf, .jpg, .tif}
        \else
        \DeclareGraphicsExtensions{.eps, .jpg}
        \fi
\title[Reacting Compressible Euler Equations]{
Global Entropy Solutions
to Exothermically Reacting, Compressible Euler Equations}

\author
{Gui-Qiang Chen}
\address[G.-Q.~Chen]{
Department of Mathematics\\
Northwestern University\\
Lunt Hall, 2033 Sheridan Road\\
Evanston, IL 60208-2730, USA}
\email[G.-Q.~Chen]
{gqchen@math.northwestern.edu
}

\author{David H. Wagner}
\address[D.~H.~Wagner]{
Department of Mathematics\\
University of Houston\\
Houston, TX 77204-3008, USA
}
\email[D.~H.~Wagner]{DWagner@uh.edu}

\date{April 20, 2002}
\begin{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.
\end{abstract}

\maketitle

\section{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:
\begin{subequations}\label{1.1}
\begin{align}
\rho_t + (\rho u)_x &= 0,\\
(\rho u)_t + \p{\rho u^2 + p} _x  &= 0,\\
(\rho E)_t + \left((\rho E + p)u\right)_x &=
q\rho Y\phi(T),\\
(\rho Y)_t + (\rho uY)_x &= - \rho Y\phi(T),\\
(\rho,\ u,\ E,\ Y)(x,0) &= (\rho_0,\ u_0,\ E_0,\ Y_0)(x).
\end{align}
\end{subequations}
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 system are
\begin{subequations}\label{1.2}
\begin{align}
&v = \frac{1}{\rho} \, (\text{specific volume)}, 
      \qquad E = e(v,S) + \frac{u^2}{2},\\
%%
%%Change
%&p = -\pd{e}{v}(v,S) > 0,\quad T = \pd{e}{S}(v,S) = \myf{p}{R\rho} > 0,\\
%To
%%
&p =-e_v(v,S) > 0,\qquad T = e_S(v,S) = \myf{p}{R\rho} > 0,\\
&p_v (v,S) < 0, \quad p_{vv}(v,S) > 0,
\end{align}
\end{subequations}
on any compact set in $v>0$.  We assume that the reaction rate
function $\phi$ is monotonically increasing and Lipschitz continuous.
In addition, inadmissible discontinuous solutions are eliminated by
requiring the following {\em entropy condition:}
\begin{equation}
\p{\rho S}_t + \p{\rho uS}_x \geq \myf{q \rho Y \phi(T)}{T}.
\end{equation}
The ideal gas assumption $p= R\rho T$ implies \cite{hyp:Whit} that the internal
energy can be written as a function of $T$ alone,
$e(v,S) = \tilde{e}(T)$.
Then $\tilde{e}'(T) = c_v(T)$ is the specific heat at constant volume. 
\begin{comment}
One can then compute
\begin{equation}\label{eq:1.2c}
    \begin{aligned}
    T &=\pd{e}{S}(v,S) = c_{v}(T)\pd{T}{S},\\
    p &= -\pd{e}{v}(v,S) = -c_{v}(T)\pd{T}{v}.
    \end{aligned}
\end{equation}
>From the ideal gas law, we calculate $dp = (R/v)dT - (RT/v^{2})dv$.
In combination with \eqref{eq:1.2c}, this implies
\begin{equation*}
    \begin{aligned}
        \pd{S}{T}(T,p) &= \myf{c_{v}(T)+R}{T},\\
        \pd{S}{p}(T,p) &= -\myf{R}{p}.\\
    \end{aligned}
\end{equation*}
\end{comment}
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=\g>1$ 
is the adiabatic exponent, a parameter that determines the equation 
of state. 

%%
%%Change
%This system of equations  (\ref{1.1}--\ref{1.2}) is useful for
%To
%%
This system of equations  (\ref{1.1})--(\ref{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.  
%%
%%Change
%This justifies the use of the Euler equations in (\ref{1.1}--\ref{1.2}), 
%To
%%
This justifies the use of the Euler equations in (\ref{1.1})--(\ref{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 \cite{M}, 
linearized stability analysis, as well as numerical and physical 
experiments, have shown that certain steady detonation waves are unstable
\cite{hyp:BMR,comb:Erp64,comb:FW,comb:LS,comb:OS}.  
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 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.~\ref{fig:1.1}.  
\begin{figure}[ht]\label{fig:1.1}
\hfil \psfig{figure=detpict.ps}\hfil
\caption{Computation of a pulsating detonation.}
\end{figure}
%\begin{figure}[ht]\label{fig:1.1}
%\ifpdf
%\includegraphics[viewport=33 1 458 188]{detpict2}
%\else
%\hfil \psfig{figure=detpict2.ps}\hfil
%\fi
%\caption{Computation of a pulsating detonation.}
%\end{figure}
The possibility of such oscillation clearly indicates that we are 
presented with both an interesting challenge and the possibility of 
discovering new and interesting mathematics.

%%
%%Change
% The system (\ref{1.1}-\ref{1.2}) is 
%To
%%
The system (\ref{1.1})-(\ref{1.2}) is 
a {\em hyperbolic system of balance laws}.
That is, it has the form
\begin{equation}\label{1.3}
U_t + F(U)_x = G(U),
\end{equation}
where $U(x, t)\in\Rn,\ F(U), G(U) \in \Rn$, and for each $U$, the
$n\times n$ matrix $DF(U)$ has a set of real eigenvectors which form a
basis of $\Rn$.  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 (\ref{1.3})
is called a system of {\em conservation} laws.  Such a system is {\em strictly}
hyperbolic if the $n$ eigenvalues of $DF(U)$  are real and distinct:
$\la_1 < \cdots < \la_n$.  The right eigenvectors 
$r_1, \cdots ,r_n$ correspond to the fields of simple waves admitted
by the system (\ref{1.3}).  The corresponding eigenvalues give the
characteristic speeds of propagation associated with these fields.
We say that the $j^{th}$ wave field is {\it genuinely nonlinear} if
$\grad \la_j(U) \cdot r_j(U) \neq 0$ for all $U$.  The $j^{th}$ wave
field is {\it linearly degenerate} if $\grad \la_j(U) \cdot 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 \cite{hyp:Gl65}.  The equations
of gas dynamics, that is, the first three equations of (\ref{1.1}),
together with (\ref{1.2}), omitting any terms containing $Y$, is not
completely genuinely nonlinear.  One may compute that
$$
\la_1 = u-c,\quad \la_2 = u,\quad \la_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 \cite{hyp:Liu77,hyp:Temp81}.

The system of conservation laws for (\ref{1.1}), that is, omitting
the reaction rate term $\rho Y\phi(T)$, has a fourth wave field which
is linearly degenerate with $\la_4 = \la_2 = u$.  Thus the system
(\ref{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:
\begin{subequations}\label{1.4}
\begin{align}
v_t - u_y &= 0,\\
u_t + p_y &= 0,\\
E_t + \p{pu}_y &=  q Y\phi(T),\\
Y_t &= - Y\phi(T) .
\end{align}
\end{subequations}
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 (\ref{1.4}).  Any existence theorem proved for this
system can be translated into an existence theorem for (\ref{1.1})
(cf. see \cite{hyp:Wg87}).

Several results have been obtained regarding the existence of entropy
solutions to hyperbolic systems of balance laws
\cite{hyp:DH,hyp:DC,hyp:LusTem,hyp:WY1,hyp:WY2}.
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 \cite{hyp:DH} to prevent resonance.
Some $2\times2$ physical systems were discussed in
\cite{hyp:DC,hyp:LusTem,hyp:WY1,hyp:WY2}.  
The papers \cite{hyp:DH,hyp:LusTem} had in view applications in
which the lower-order terms act in a way that reduces the spatial
total variation of the solution as time increases.  
In \cite{hyp:WY2}, the lower-order term has a coefficient $e^{-Kt}$,
so that the lower-order term decays uniformly without the need for
{\em a priori} estimates.  In \cite{hyp:WY1}, the first order terms
constitute a system of conservation laws for which solutions exist
with arbitrarily large initial data \cite{hyp:Ni68}.  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.
 
%%Change
%These remarks do not apply to the system (\ref{1.1}--\ref{1.2}), or
%To
%%
These remarks do not apply to the system (\ref{1.1})--(\ref{1.2}), or
(\ref{1.4}).  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 $\phi(T)$ never vanishes, so that
there is a positive minimum value $\Phi := \phi(T')>0$.  In a sense,
this is a very realistic condition.  Typically, $\phi(T)$ has the
Arrhenius form:
\begin{equation}\label{1.5}
\phi(T) = T^{\alpha}e^{-E/RT},
\end{equation}
which vanishes only at absolute zero temperature.  However, in a
typical unburned state, $\phi(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 \cite{hyp:Temp81}, we consider a one-parameter family of
functions $e(v,S,\eps)$, $\eps \geq 0$, which is $C^{5}$ and satisfies
\eqref{1.2}.   We assume that, 
when $\eps=0$, the equation of state is that of an isothermal gas:
\begin{equation}
    e(v,S,0) = -\ln(v) + \myf{S}{R}.
\end{equation}
%%Change
%For a polytropic gas, $\eps = \g-1$ and for $\eps > 0$,
%\begin{equation}
%   e(v,S,\eps) = \myf{1}{\eps}\p{\p{v\exp\p{-S/R}}^{-\eps}-1}
%\end{equation}
%To
%%
For a polytropic gas, $\eps = \g-1$, and for $\eps > 0$,
\begin{equation}
   e(v,S,\eps) = \myf{1}{\eps}\p{\p{v\exp\p{-S/R}}^{-\eps}-1}.
\end{equation}


One may easily check that this function is $C^{\infty}$ and that, 
as $\eps\rightarrow 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\^opital's rule to calculate 
\begin{equation}
  \pd{e}{\eps}(v,S,0) = \myf{1}{2}\p{-\ln(v) + \myf{S}{R}}^{2},
\end{equation}
and see that $\pd{e}{\eps}(v,S,\eps)$ is continuous at $\eps=0,\ v > 0$.

%%Change
%The value $\eps=0$ is mathematically special because at this value,
%To
%%
The value $\eps=0$ is mathematically special because, at this value,
the system \eqref{1.1} has a complete system of Riemann invariants:
\begin{equation}\label{eq:1.9}
  (r,s,S,Y) = (u - \ln(p), u + \ln(p),S,Y).
\end{equation}
Furthermore, all shock, rarefaction, and contact discontinuity curves
%%
%%Change
%in the $(r,s,S,Y)$-space are invariant under translation of the base point.
%To
%%
in $(r,s,S,Y)$-space are invariant under translation of the base point.
Following \cite{hyp:Temp81}, we use $(r,s,S,Y)$ as the coordinates for
our analysis in $\eps\geq 0$.  
Note that since $p=-\pd{e}{v}(v,S,\eps)$
and $e(v,S,\eps)$ is $C^{5}$,
the transformation between $(v,u,S)$ and $(r,s,S)$ is $C^{4}$.
Moreover, Temple \cite{hyp:Temp81} showed that this transformation
is a diffeomorphism.

Our principal result, Theorem \ref{thm:3.1}, is somewhat complex, but
a simple version may be stated as follows.
\begin{thm}\label{thm:1.1}
%%
%%Change
%Let $w-$ be a point in $rsS$-space, and let $\eps\in\br{0,1}$.
%To
%%
Let $(r_{-\infty},s_{-\infty},S_{-\infty})$ be a point in $rsS$-space, 
and let $\eps\in\br{0,1}$.
 Let $w_{0}(x)$ be given initial data for \eqref{1.4}, expressed in
$(r,s,S,Y)$-coordinates,
%%
%%Change
%and with $\lim_{x\to=\infty}w_0(x)=(w-,0)$,
%To
%%
and with $\lim_{x\to -\infty}w_0(x)=(r_{-\infty},s_{-\infty},S_{-\infty},0)$.
Then there is a function $C\p{\eps,\Phi,q,\supnorm{Y_0}}>0$ such that
$C\p{\eps,\Phi,q,\supnorm{Y_0}}\to \infty$ as $\eps \to 0$, and such
that, if
\begin{equation*}
Var_{rsSY}\p{w_0} < C\p{\eps,\Phi,q,\supnorm{Y_0}},
\end{equation*}
then there exists a global BV entropy solution to the Cauchy problem
%% delete "for"
%%
%%Change
%\eqref{1.1} with initial data determined by $w_0$.
%To
%%
\eqref{1.1}--\eqref{1.2} with initial data determined by $w_0$.

Moreover, for $\eps$ small,  $C\p{\eps,\Phi,q,\supnorm{Y_0}}$ has the form
\begin{equation*}
C\p{\eps,\Phi,q,\supnorm{Y_0}} = \myf{B}{\eps^{1/3}}\exp\p{-\myf{K q\supnorm{Y_0}}{\Phi}}.
\end{equation*}
for some constant $B$.
\end{thm}


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 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 \ref{thm:1.1} (and of the more sophisticated Theorem
\ref{thm:3.1}) may be outlined as follows.  In
Section \ref{sect:3}, we construct approximate solutions to the Cauchy
problem using a fractional--step method based on Glimm's scheme.  In
Section \ref{sect:4}, 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
\cite{hyp:Temp81}.  
In Section \ref{sect:5}, we use our bounds to prove that a subsequence
from our approximate solutions converges to an entropy solution of the
Cauchy problem.

\section{Approximate Solutions} \label{sect:3}
In this section, we construct approximate solutions to the Cauchy
problem for \eqref{1.4} and \eqref{1.2}. Using $(r,s,S,Y)$ as coordinates, 
the initial data are given by
\begin{equation}
(r, s, S, Y)|_{t=0} =(r_0, s_0, S_0, Y_0)(x).
\label{2.1}
\end{equation}
We may represent (\ref{1.4}), (\ref{1.2}), and (\ref{2.1}) in the
compact form:
\begin{equation}\label{2.2}
\begin{aligned}
U_t + F(U)_x &= G(U), \\
U|_{t=0}&= U_0(x).
\end{aligned}
\end{equation}

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.

\subsection{The Riemann Problem for the Non-reacting Gas}
The Riemann problem for the non-reacting gas is the following Cauchy
problem:
\begin{subequations}\label{2.3}
\begin{align}
&v_t - u_x = 0,\\
&u_t + p(v,S,\eps)_x = 0,\\
&\p{e(v,S,\eps) + u^2/2}_t + \bp{p(v,S,\eps)u}_x =  0,\\
&Y_t = 0 ,
\end{align}
\end{subequations}
with the initial data:
\begin{equation}
(r,\ s,\ S,\ Y)|_{t=0}=
 \begin{cases}
(r_L,\ s_L,\ S_L,\ Y_L), \qquad &x<0,\\
(r_R,\ s_R,\ S_R,\ Y_R), \qquad &x>0,
\end{cases}
\label{2.4}
\end{equation}
where $\eps\in \br{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 (\ref{2.3}d) is independent of (\ref{2.3}a,b,c).
Consequently, the solution of the Riemann problem (\ref{2.3})--(\ref{2.4}) 
decouples into the solution of two Riemann problems: 
one for non-reacting gas dynamics, (\ref{2.3}a,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 (\ref{2.3}
a,b,c) for polytropic gases can be found in \cite{hyp:CF,hyp:Smo}.
For the more general case, with the internal energy given by
$e(v,S,\eps)$, the solution of the Riemann problem was proved to exist
for small data in \cite{hyp:Lax57}, and for large data in
\cite{hyp:Smi,hyp:ChangHsiao79}.  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.

\begin{comment}
The system \eqref{1.4}, with $Y\geq 0$ and $\phi \equiv 0$, 
is also of the form \eqref{2.2}.  
In this case, the equation for $Y$ decouples from the rest of the system.
Therefore, we will simply review the results of \cite{hyp:Temp81}
for non-reacting gases.  We will see that the equation for $Y$
affects our analysis only through the reaction term.
\end{comment}

A backward, or 1-rarefaction wave, is a solution which is constant
on the 1-characteristics centered at a point $(x_{0},t_{0})$:
$\myf{x-t_{0}}{t-t_{0}}=\myf{dx}{dt}=\lambda_1\equiv -\sqrt{-p_v(v,S)}$, where
$\n{\lambda_{1}}=\rho c = c/v$ is the material, or Lagrangian, sound speed of
the gas.
In addition, the 1-Riemann invariants $S$ and $u-\int^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
$\lambda_3 = - \lambda_1$, and the 3-Riemann invariants are
$S$ and $u+\int^v \sqrt{-p_v(v,S)}dv$.

Forward and backward shock waves are simple jump discontinuities of $U$ along
a line $\myf{x-x_{0}}{t-t_{0}} = dx/dt = \sigma$, and satisfy the 
\RH conditions:
\begin{equation}\label{eq:2.5a}
\begin{aligned}
\sigma\br{v} = -\br{u},\\
\sigma\br{u} = \br{p},\\
\sigma\br{u^2/2+ e} = \br{pu}.\\
\end{aligned}
\end{equation}
Here $\br{f}$ denotes the jump in the quantity $f$ across the shock
wave. In addition to \eqref{eq:2.5a}, shock waves must satisfy the
{\em 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 \eqref{1.2}, Weyl
\cite{hyp:Weyl49} 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 $\lambda_{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,\eps)$
emanating from $w_L$ associated with each characteristic field such
that $w = H_{i}(z_{i},w_L,\eps)$ is connected to $w_L$ by an i--wave of
strength $z_i$.

%%
%%I have changes "parameterize" to "parametrize" below.
%%
We parametrize 1--shock curves by $z_{1}=r-r_L \leq 0$ and 3--shock curves
by $z_{3}=s-s_L \leq 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\geq 0$ and 3--rarefaction curves by
$z_{3}=s-s_L\geq 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}\neq 0$, these curves are at least $C^{3}$
\cite{hyp:Lax57,hyp:Temp81}.  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} \neq 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 $\bz = (z_{1},z_{2},z_{3})$ and 
\begin{equation*}
 H_{g}(\bz,w_L,\eps)= H_{3}\p{z_{3},
 H_{2}\p{z_{2},H_{1}\p{z_{1},w_L,\eps},\eps},\eps}.
\end{equation*}
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 $\bz \rightarrow H_{g}(\bz,w_L,\eps)$.  In
\cite{hyp:Temp81}, it is shown, using the implicit function theorem,
that for every $w=(r,s,S)\in \R^3$ there exists a neighborhood
$\Omega$ of $w$ such that the Riemann problem is solvable for $(w_L,\
w_R)\in\Omega\times\Omega$ and $\eps\in\br{0,1}$, and that one can
solve for the vector of signed wave strengths $\bz = (z_1,\ z_2,\
z_3)$ as a function of $(w_L,w_R,\eps)$,
\begin{equation*}
\bz = B_{g}(w_L,\ w_R,\ \eps),
\end{equation*}
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,\ \eps)$ the solution of the Riemann problem for
\eqref{2.3} with data $(w_L,w_R)\in \p{\Omega\times \br{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. 

\subsection{The Glimm Fractional--Step Scheme}\label{sect:2.2}
We will use a fractional--step scheme such as that described in
\cite{hyp:DH,hyp:LusTem}
based on the Glimm scheme \cite{hyp:Gl65}.
The construction of the fractional--step scheme for the inhomogeneous
system (\ref{1.3}) 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:
\begin{equation}
\Lambda = \max\limits_{1\le j\le 4}\n{\lambda_j(u)} \le \myf{l}{2h}
\label{2.15}
\end{equation}
is satisfied.

Partition $\R^+$ by the sequence $t_k=k h,\ k\in \mathbb Z^+$,
and partition $\R$ into cells with the $j$-th cell centered at
\begin{equation*}
x_j = jl, \qquad j=0,\pm1, \pm2, \cdots.
\end{equation*}
We begin by approximating the initial data $U_{0}(x)$ by a function
$U^{h}(x,0)$ which is constant for $x_{j-1}\leq 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\rightarrow 0$.  Choose a random
sequence $\chi_k,\ k =0,\ 1,\ 2,\ \cdots,$ 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:
\begin{alignat*}{2}%\label{2.16}
&U^h(x,t) = \mathcal{F}(t-kh, \mathcal{S}_0(t-k h,U^h(x, k h+0))),&
 \qquad &k h \le t< (k+1)h,\hfil \\
&U^h(x,k h+0) = U_j^k,&\quad j+k\mbox{ even,} \quad &(j-1)l<x<(j+1)l,\hfil\\
&U^{k+1}_j=\mathcal{R}^{k+1}_j\circ U^h(x,(k+1)h-).&&
\end{alignat*}
Here $\mathcal{F}(\tau,\cdot)$ is the fractional--step operator, which
advances the chemical reaction for $0<\tau\leq h$, 
$\mathcal{S}_0(\tau,\cdot)$
is the solution operator for the homogeneous system (\ref{2.3})
(i.e. $G(U)\equiv 0$) for $0 < \tau \leq h$, and $\mathcal{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 <k h$. Then we define
$U^h(x,k h+0)$ as follows.
\subsubsection*{Random Operator $\mathcal{R}_j^k$}
We define
\begin{align*}%\label{ 2.17}
&U_j^k = \mathcal{R}_j^k \circ U^h(x, k h-)\equiv U^h((j+\chi_k)l,
kh-),\\
&U^h(x,k h+0) \equiv U_j^k, \quad (j-1)l < x < (j+1)l,
\end{align*}
where $j+k$ is even and $\chi_k$ is the $k^{th}$ element of the random
sequence.
\subsubsection*{Homogeneous Solution Operator $S_0(t,\cdot)$}
In the strip $k h\le t <(k+1)h$,
we define
\begin{equation*}
U_0^h(\cdot,t) \equiv \mathcal{S}_0(t-k h, U^h(\cdot, k h+0)),
%\label{ 2.18}
\end{equation*}
where $\mathcal{S}_0(\tau,W)$ is the solution $U(\cdot,\tau)$ of the
following Cauchy problem:
\begin{equation}\label{eq:2.19}
\begin{split}
&U_{\tau} + F(U)_x = 0, \\
&U|_{\tau=0}= W(x).
\end{split}
\end{equation}

This Cauchy problem can be solved by
constructing the Riemann solutions of the Riemann problems 
(\ref{2.3})--(\ref{2.4}) with Riemann data:
\begin{equation*}
%\label{2.20}
  U_L =  U_j^k, \quad U_R = U_{j + 2}^k,
\qquad j= 0, \pm1, \pm2,\cdots.
\end{equation*}

This construction determines the unique solution of the Cauchy problem
(\ref{eq:2.19}), 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 (\ref{2.15}).

\subsubsection*{Fractional--Step Operator $\mathcal{F}(\tau,\cdot)$}
The fractional--step operator uses an approximation
to the solution of the initial value problem:
\begin{equation}\label{ivp}
\begin{split}
&\myf{dU}{d\tau} = G(U(\tau)),\\
&U|_{\tau=0}= U_0.
\end{split}
\end{equation}
This approximation must be
{\em  consistent} and {\em total variation stable}.  These
requirements are
defined as follows.
\begin{mydef}\label{def:2.1}
The explicit one-step approximation algorithm $\bar{U}(\tau)
= \F(\tau,U_0)$ is a {\em consistent first order} approximation to the solution
$U(\tau)$ of the initial value problem (\ref{ivp}) in the domain
$\DD\subset\R^n$ if $\F$ is Lipschitz continuous on
$\p{-\delta,\delta}\times \DD\to \R^n$ and there exists an increasing
function $\nu(\tau)>0$
such that $\nu(\tau)$ converges to $0$ as $\tau\to 0$
and  
\begin{equation}\label{eq:2.14b}
\N{\bar{U}(\tau)-U(\tau)} \leq \nu(\tau)\N{G\bp{U_{0}}}\n{\tau},
\end{equation}
for every $U_0\in \DD$.
\end{mydef}

\begin{rem} The requirement that $\nu(s)$ converges to $0$ as $\tau\to
0$ implies that $\F$ is ``{\em uniformly differentiable}'' with
respect to $\tau$, at $\tau=0$, on $\DD$.  The additional requirement in
(\ref{eq:2.14b}) regarding the factor $\N{G\bp{U_{0}}}$ is necessary
to obtain exponential decay in $Y$ as $t \to \infty$.
\end{rem}

\begin{mydef}\label{def:2.2}
The explicit one-step approximation algorithm $\bar{U}(\tau)=\F(\tau,U_0)$
is a {\em total variation stable} approximation to
the solution $U(\tau)$ of the initial value problem (\ref{ivp})
in the domain $\DD\subset\R^n$
if $\F,\ \pd{\F}{\tau},$ and $\pdtwo{\F}{\tau}$ are 
Lipschitz continuous on
$\p{-\delta,\delta} \times \DD \to \R^n$.
\end{mydef}

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.

\begin{mydef}\label{def:2.2a}
The explicit one-step approximation algorithm $\bar{U}(\tau)=\F(\tau,U_0)$
is a {\em conservative} approximation to the initial value problem
\begin{align}\label{eq:2.15a}
\myf{d\ }{d\tau}\begin{pmatrix}v(\tau)\\u(\tau)\\E(\tau)\\Y(\tau)\\
\end{pmatrix} 
&= \begin{pmatrix}0\\ 0\\ qY(\tau)\phi\bp{T(\tau)}\\-Y(\tau)\phi\bp{T(\tau)}\\
\end{pmatrix}, & 
\begin{pmatrix}v(0)\\u(0)\\E(0)\\Y(0)\\
\end{pmatrix} = 
\begin{pmatrix}v_{0}\\u_{0}\\E_{0}\\Y_{0}\\
\end{pmatrix}.
\end{align} 
if there is a function $f(\tau,U_{0}) \geq 0$ such that
$\F(\tau,U_{0})$ has the form:
\begin{equation}\label{eq:2.17a}
\F(\tau,\ v_{0},\ u_{0},\ E_{0},\ Y_{0}) = \p{v_{0},\ u_{0},\ E_{0}+q
f(\tau,U_{0}),
\ Y_{0}-f(\tau,U_{0})}.
\end{equation} 
\end{mydef}

\begin{lem}\label{lem:2.1}
If $\F(\tau,U_{0})$ is a conservative, consistent first order one-step
approximation algorithm with $f$ as in {\rm \eqref{eq:2.17a}}, then
$f(0,U_{0}) = 0$, and
\begin{equation*}
 \pd{f}{\tau}(0,v_{0},u_{0},T_{0},Y_{0}) = Y_{0}\phi(T_{0}).
\end{equation*} 
\end{lem}

\begin{rem}
 When $G$ is $C^2$, standard one-step approximation methods, such as
 Euler's method and Runge-Kutta methods, are {\em consistent first
 order, conservative,} and {\em 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 {\em consistent first order}, but will not meet our
 criteria for {\em 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 (\ref{1.1})
and the corresponding Navier--Stokes model
(cf. \cite{hyp:Wg89:1,hyp:Wg89:2,hyp:Chen92}),
it is customary to modify $\phi(T)$ to be zero for $T$ less than
an ``ignition temperature'' $T_i$.  The observation above indicates
that it may be inadvisable to make such a modification in a manner
that makes $\phi$ less than $C^2$---particularly in conjunction
with a numerical computation of time-dependent flow.
\end{rem}

Let $\mathcal{F}(\tau,U_0),\ 0 \leq \tau \leq h$, be a conservative,
consistent first order, total variation stable, one-step approximation
algorithm to the initial value problem (\ref{eq:2.15a}).  Given
$U^h_0(x,t)$, $k h\le t <(k+1)h$, we define $U^h(x,t) =
\mathcal{F}(t-kh, U^h_0(x,t))$.

We re-approximate $U^h(x,(k+1)h)$ using $\mathcal{R}^{k+1}_{j}$:
\begin{equation*}%\label{eq:2.24}
\begin{aligned}
U^{k+1}_{j}&=\mathcal{R}_{j}^{k+1}\circ U^h(x, (k+1)h-) \\
      &\equiv\mathcal{F}\bp{h; U_0^h\bp{(j+\chi_{k+1})l,\ (k+1)h-}},
    \quad jl < x < (j+2)l.
\end{aligned}
\end{equation*}

We will guarantee condition (\ref{2.15}) by showing that %deleted ","
for all $(x,t)$, $U^h(x,t)\in \D\times\br{0,1}$, where $\D$
is a compact set in $(r,s,S)$--space.
It follows that $\myf{h}{l}$ can be chosen so
that the Courant-Friedrichs-Levy condition (\ref{2.15}) is satisfied
for all time-steps.  Therefore, our approximate solutions for the
Cauchy problem of the system (\ref{1.1}--\ref{1.2}) are defined
unambiguously.

\section{Total Variation Stability}\label{sect:4}
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  $\ell_1$
norm, $\onenorm{\mathbf{v}}= \n{v_1}+\n{v_2}+\n{v_3}+M_{4}\n{v_4}$,
for a vector $\mathbf{v}\in \R^4$, and where $M_{4}>0$ is defined
later.  In addition, for a vector function $U(x)$ with components
$\p{r,s,S,Y}(x)$, we define the total variation using the weighted $\ell_1$ norm:
\begin{equation}\label{Vardef}
\Var_{rsSY}(U) = \onenorm{\p{\Var(r),
\Var(s),\Var(S),\Var(Y)}}.
\end{equation}
Other equivalent definitions are
possible; for example, one could use the Euclidean norm instead of the
$\ell_{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 $\myf{1}{h}$. Thus, when we sum these increases, we
obtain a finite, non-zero change in total variation.  However, this
same argument leaves us free to neglect terms of order
$h^2TV(U^h(\cdot,t))$ or higher, because in the limit as $h\rightarrow
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\times\br{0,1}$, where $\D$ 
will be defined later.  
Note that after the random choice step, the values of the solution
between $((j+\chi_k)l, kh)$ and $((j+1+\chi_k)l,(k\pm 1)h)$
depend only on the values of $U$ at these points.  We denote these
values by $\p{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$ {\em 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
\begin{itemize}
\item The Riemann problem for the non-reacting equations with data
$(U_L,U_R)$ is solvable for every ordered pair $(U_{gL},U_{gR}) \in
\D\times \D$.
\item The function $B_g(U_{gL},U_{gR},\eps)$, 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\times\D$.
\end{itemize}

\begin{lem}\label{lem:3.0.1}
Let
\begin{equation*}
\bz = (z_1,\ z_2,\ z_3,\ z_4) = B(U_L,\ U_R,\ \eps)
\end{equation*}
be the vector of signed wave strengths in the solution of the Riemann
problem with data $(U_L,\ U_R)$.  Let $\F(\tau,U)$ be a
consistent first order and total variation stable approximation to the
initial value problem {\rm\eqref{ivp}}. Let
\begin{equation}\label{eq:3.0.0}
\Gamma\bp{U_L,\ \bz,\ h}=B\bp{\F(h,U_L),\ \F(h,H(\bz,\ U_L)),\ \eps}.
\end{equation}
Then 
\begin{equation*}
\Gamma(U_L,\ \bz,\ h) = \bz + O(\onenorm{\bz})h.
\end{equation*}
\end{lem}

\begin{proof}
Note that $\Gamma(U_L,\bz,0)=\bz$ and $\Gamma(U_L,0,h) = 0$.
Thus we compute that
\begin{equation*}
\pd{\Gamma}{h}(U_L,\ 0,\ h) = \pdtwo{\Gamma}{h}(U_L,\ 0,\ h) =
0_{4\times 1}.
\end{equation*}
Let $L_1$ be a Lipschitz constant for $\pd{\Gamma}{h}$ with respect
to the $\ell_{1}$ norm.
Then
\begin{equation*}
\Bigonenorm{\pd{\Gamma}{h}(U_L,\ \bz,\ h)} \leq L_{1}\onenorm{\bz},
\end{equation*}
and
\begin{equation*}
\onenorm{\Gamma(U_L,\ \bz,\ h) -\bz}=
\Bigonenorm{\int_0^1 \pd{\Gamma}{h}(U_L,\ \bz,\ \theta h) h\ d\theta}
\leq L_1\onenorm{\bz}h.
\end{equation*}
\end{proof}

Lemma \ref{lem:3.0.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.

\begin{lem}\label{lem:3.2}
Let $L_2$ be a Lipschitz constant for $\pdtwo{\Gamma}{h}$.  Since
$\Gamma(U_L,0,h) = 0$, then
\begin{equation}\label{eq:3.00a}
\Bigonenorm{\Gamma(U_L,\ \bz,\ h) 
- \Bp{\bz + \pd{\Gamma}{h}(U_L,\ \bz,\ 0)h}} 
\leq L_2\onenorm{\bz}\myf{h^2}{2}.
\end{equation}
\end{lem}
\begin{proof}
By Taylor's Theorem,
\begin{equation*}
\begin{aligned}
\Gamma(U_L,\ \bz,\ h) -\bz
&=\pd{\Gamma}{h}\p{U_L,\ \bz,\ 0}h-\int_0^1\pdtwo{\Gamma}{h}\p{U_L,\
\bz,
\ \theta h}
(\theta-1)\ d\theta\ h^2.\\
\end{aligned}
\label{eq:3.00}
\end{equation*}
\end{proof}

Lemma \ref{lem: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 (\ref{eq:3.00a}) and Lemma \ref{lem:3.0.1}, we see that the 
term $\pd{\Gamma}{h}(U_L,\ \bz,\ 0)h$ is $O(\onenorm{\bz})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 $\supnorm{Y}$ decays exponentially as
$t\rightarrow \infty$.  Then we show that any increase in total
variation, due to the reaction step, is proportional to
$Y\onenorm{\bz}h$. This, together with the exponential decay of
$\supnorm{Y}$, will enable us to derive the conditions under which the
total increase of total variation is bounded.

\subsection{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 (\ref{1.4}),
this operator takes the form
\begin{equation*}\label{3.1.1}
\mathcal{F}\bp{\tau,(v_0^{h},\ u_0^{h},\ E_0^{h},\ Y_0^{h})} =
(v^{h},\ u^{h},\ E^{h},\ Y^{h}), \qquad 0 \le \tau < h,
\end{equation*}
where all quantities $v_0^h, v^h$, etc. are evaluated at $(x, kh+\tau)$,
and
\begin{equation} 
\begin{split}\label{3.1.2}
&v^{h}(x,k h +\tau) = v_0^{h}, \\
&u^{h}(x,k h +\tau) = u_0^{h},  \\
&E^{h}(x,k h+\tau) = E_0^{h}(x,\ k h+\tau)+q f(\tau,U_{0}^{h}(x,kh+\tau)),\\
&Y^{h}(x,k h + \tau) =  Y_0^{h}(x,\ k h+\tau)-f(\tau,U_{0}^{h}(x,kh+\tau)).
\end{split}
\end{equation} 
We need to estimate the change in $T$ due to $\mathcal{F}$.
By Lemma \ref{lem:3.2}, we only need to calculate
\begin{equation*}
\p{T^h-T_0^h}(x,k h +\tau) =\rtvert{\pd{T}{E}}_{U_0^h} 
          q f_{\tau}(0,U_0^h(x,kh+\tau))\tau =
          \myf{q}{e'(T)} f_{\tau}(0,U_{0}^{h}(x,kh+\tau))\tau.
\end{equation*}
By Lemma \ref{lem:2.1}, $f_\tau(0,U_0^h) = Y_0^h\phi(T_0^h)$.  Note
that this estimate is still valid for a non-polytropic 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.

\begin{lem}\label{lem:3.3}
Let $\mathcal{F}$ be a consistent, conservative one-step approximation
algorithm. Then, to first order in $\tau$,
\begin{equation}\label{3.8a}
\begin{aligned}
T^{h} &= T_0^{h}+
\myf{q}{c_v} Y_{0}^{h}\phi\bp{T_{0}^{h}}\tau,\\
Y^{h} &=  Y_0^{h}- Y_{0}^{h}\phi\bp{T_{0}^{h}}\tau,
\end{aligned}
\end{equation}
where $T^h,\ Y^h,\ T_0^h$, and $Y_0^h$ are all evaluated at $(x,k h +\tau)$.
\end{lem}

We have assumed that $T_0^{h}(x,t) \ge T'>0,\ 0 \le Y_0^{h}(x,t) \le
1$, and that $\phi$ is Lipschitz continuous, non-negative, and
increasing.  In particular,
\begin{equation*}
T > T' \implies 0 < \Phi \le \phi(T) \le C < \infty.
\end{equation*}
Note that the Arrhenius law (\ref{1.5}) satisfies these conditions.
Then
\begin{equation}\label{3.1.9}
Y_0^h\phi(T_0^h)\tau \geq Y_{0}^h\Phi \tau.
\end{equation}
Equation \eqref{3.8a}, or inequality \eqref{3.1.9}, implies 
the following lemma.

\begin{lem}\label{lem3.1} In the limit as $h \to 0$, the functions
$Y^{h}(x,t)$ and $T^{h}(x,t)$ satisfy
\begin{align*}
&0 \le Y^{h}(x,k h+\tau) \le Y_0^{h}(x,k h+\tau)e^{-\Phi \tau},
\qquad 0 \le \tau < h,\\
&T^{h}(x,k h+\tau) \ge T_0^{h}(x,k h+\tau) \ge c_0 >0.
\end{align*}
Furthermore, $Y^{h}(x,t) \leq Y^{h}(x,0)e^{-\Phi t}$ for all
$t \geq 0$.  
\end{lem}
Henceforth, for simplicity of exposition, we shall use this estimate
in the form given by the limit as $h\to 0$.


\subsection{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+\chi_k)l, k h)$, where $k
\in \mathbb N$ and $j \in \mathbb Z$ 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+\chi_k)l,k h)$ to mesh points
$((j+1+\chi_k)l, (k\pm1)h)$ (see Fig.~\ref{fig:3.1}).  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$.
\begin{figure}[ht]
\ifpdf
\includegraphics[viewport=0 0 330 80]{glimm}
\else
\hfil \psfig{figure=glimm.ps}\hfil
\fi
\caption{The mesh curves $J(k)$ and $J(k+1)$.}\label{fig:3.1}
\end{figure}

Let $J(k)$ be the unique mesh curve which connects the mesh points on
$t=k h$ 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 \le 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
\begin{align*}
&\Var_{rsSY}(r_0, s_0, S_0, Y_0)(\cdot) \le N < \infty, \\
&(r_0, s_0, S_0)(x)\in E \subset \D.
\end{align*}
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,\eps) =\bz= (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 $\eps=0$, we let $\alpha$ be the one--shock
in $\bz$ and $\mu$ be the 1--rarefaction wave in $\bz$ so that
\begin{equation*}
    \begin{aligned}
  \alpha&=\begin{cases}
  & 0,  \qquad \text{if }z_{1} \geq 0,\\
                 &\n{z_{1}}, \qquad \text{if }z_{1} < 0,
         \end{cases} \\
   \mu &= z_1 - \alpha.
    \end{aligned}  
\end{equation*}
Similarly, let
\begin{equation*}
 \begin{aligned}
     \beta &= \text{ 3-shock in }\bz,\\
 \eta &= \text{3-rarefaction in }\bz,\\
     \delta &= z_{2} = \text{contact wave in }\bz.\\
 \end{aligned}
\end{equation*}
Finally, for $\eps>0$, we define (see \cite{hyp:Temp81}, eq. (3.49)):
\begin{equation} \label{eq:3.11}
\begin{alignedat}{3}
&\text{For 1-waves, }\alpha = \Var_r^-(\alpha) + \eps\Var(\alpha),\
 &\delta_{\alpha} = \Var_S(\alpha),
 \ &\mu  = \Var_r^+(\mu) + \eps\Var(\mu),&\\
&\text{For 3-waves, }\beta  = \Var_s^-(\beta)  + \eps\Var(\beta),\ 
 &\delta_{\beta}  = \Var_S(\beta),
 \ &\eta = \Var_s^+(\eta)  + \eps\Var(\eta),&\\
&\text{For 2-waves, }\delta = S_+-S_-,&&&\\
&\text{For 4-waves, }\zeta=Y_+ - Y_-. &&&\\
\end{alignedat}
\end{equation}
Here $\Var(\alpha)$ 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 $\alpha$, and similarly for  the other waves.  
Note that since $S$ is monotone along the shock curve \cite{ hyp:Weyl49}, 
$\Var_{S}(\alpha)=\n{S_{+}-S_{-}}$.


We let $C$ denote a constant that depends only on $\D$ and $\phi$, 
and is independent of $\eps$ and the mesh length $h$.

For any mesh curve $J$, Temple defined \cite{hyp:Temp81} (p. 144):
\begin{equation*}%\label{3.13}
L_{T}(J) = \sum_J \mybrace{\p{\alpha_i-M_0\delta_{\alpha_i}}
+\p{\beta_i-M_0\delta_{\beta_i}} +M_0\n{\delta_i}}
+\eps\sum_J\mybrace{\mu_i+\eta_i} +V_0,
\end{equation*}
where $\alpha, \beta, \delta, \zeta, \delta_{\alpha}, \delta_{\beta},
\mu$, and $\eta$ 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 \ref{lemtemp:4.1} and \ref{lemtemp:4.2}
(see Lemmas 4.1 and 4.2 in  \cite{hyp:Temp81}).

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 \ref{lem:3.7}.  Let $U_{\pm\infty}$
denote the values of $U$ at $\pm \infty$, along a mesh curve. Let
$\alpha_{\infty}$, etc. denote the wave strengths (as defined in
\eqref{eq:3.11}) in the solution of the Riemann problem with data
$\bp{U_{-\infty},U_{\infty}}$.  Let
\begin{equation*}%\label{3.13}
 \begin{split}
L(J) = &\sum_J \mybrace{\p{\alpha_i-M_0\delta_{\alpha_i}}
    +\p{\beta_i-M_0\delta_{\beta_i}}+M_0\n{\delta_i}+M_4 \n{\zeta_i}
    +\eps (\mu_{i}+\eta_{i})}\\
&+ \alpha_{\infty}+\beta_{\infty}+\mu_{\infty}+\eta_{\infty}+
   \n{\delta_{\infty}} +M_{4}\n{\zeta_{\infty}}.
 \end{split}
\end{equation*}
The term $M_{4}\n{\zeta_{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
$\n{r_{\infty}-r_{-\infty}}+\n{s_{\infty}-s_{-\infty}}$.  However,
during the reaction step, $r_{\pm\infty}$ and $s_{\pm\infty}$ 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
$\bp{U_{-\infty},U_{\infty}}$ as a similar bound which works well with
our estimates on the reaction step.    Note that for any
mesh curve $J$, $L_{T}(J) \leq L(J)+V_0$.

Let $a_i$ and $b_j,\ 1\leq i,j\leq 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}$ {\em approaches} the wave
$b_{j}$ if either $i>j$ ($a_{i}$ is faster
than $b_{j}$), or $i=j\in\mybrace{1,3}$ and at least one of the waves
is a shock wave---either $a_{i}<0$ or $b_{j}<0$.

Since $\zeta$ (that is,
$Y$)--waves do not interact with any other waves, we can use the same
definition of $Q(J)$ that was given in \cite{hyp:Temp81} as follows.
Let $p_i$ denote arbitrary waves, $q_i$ arbitrary shock or rarefaction
waves, and $R_i$ arbitrary rarefaction waves.  We define 
\begin{equation*}
%Q(J)= M_1 Q_1 + M_2 Q_2 +M_3 Q_3,
Q(J)= M_1 \sum_{App}p_i\n{\delta_j} + M_2 \sum_{App}q_iR_j  +
M_3 \p{\sum_{App}\alpha_i\beta_j+\sum_{App}\alpha_{i}\alpha_{j}+
\sum_{App}\beta_{i}\beta_{j}} ,
\end{equation*}
%where
%$$
%Q_1=\sum_{App}p_i\n{\delta_j},
%\quad Q_2=\sum_{App}q_iR_j,
%\quad
%Q_3=\p{\sum_{App}\alpha_i\beta_j+
%\sum_{App}\alpha_{i}\alpha_{j}+
%\sum_{App}\beta_{i}\beta_{j}},
%$$
where the sums are taken over all pairs of approaching waves, and $p_i$
can be either to the left or right of $\n{\delta_j}$, etc.  We have
corrected an apparent typographical error in \cite{hyp:Temp81}, in which
the terms for shock waves of the same field were omitted.

Finally, we define
\begin{equation*}
F(J)= L(J) + \eps Q(J).
%\label{3.15}
\end{equation*}

In \cite{hyp:Temp81}, the following functionals were also defined 
in order to estimate $\Var_{rs}(J)$ for small $\eps$:
\begin{equation*}
\begin{aligned}
 L_{T0}(J) &= \sum_J\mybrace{\alpha_i + \beta_i} +V_0,\\
 F_{T0}(J) &= L_0(J) + \eps Q(J).
\end{aligned}
\end{equation*}
We modify $L_{T0}$ in the same
way that we modified $L_T$:
\begin{equation*}
L_0(J) = \sum_J\bigbrace{ \p{\alpha_i+\beta_i} + M_4\n{\zeta_i}+
\alpha_{\infty}+\beta_{\infty}+\mu_{\infty}+\eta_{\infty}+\n{\delta_{\infty}}
  +M_{4}q\n{\zeta_{\infty}}}.
\end{equation*}

Let $F_0 = L_0 + \eps Q$, $F_{T} = L_{T}+\eps Q$ and $F_{T0}= L_{T0}+\eps Q$.
Note that for any mesh curve $J$, $F_{T}(J) \leq F(J)+V_0$, 
and $F_{T0}(J) \leq F_{0}(J)+V_0$.

We will see in Lemmas \ref{lemtemp:3.1}--\ref{lem:3.13a} and \ref{lem:3.7}
 that all of the constants $M_0, M_2, M_3, M_4$ depend only on $\D$.


\subsection{Temple's Results} 
Before we prove the total variation stability of the approximate
solutions, we first review the following lemmas from Temple
\cite{hyp:Temp81} on which our analysis depends.  Lemmas
\ref{lemtemp:3.1} to \ref{lemtemp:3.5} below correspond to Lemmas 3.1
to 3.5 of \cite{hyp:Temp81}.

\begin{lem}\label{lemtemp:3.1}
For every compact convex open set $\D_{rs}$ in the $rs$-plane, there 
exist constants $M>0, \frac{1}{2}\le C_0<1, $ and $G>1$ such that, at $\eps=0$,
the following estimates hold across any interaction
$<w_L, w_M>+<w_M,w_R>\, \to\, <w_L, w_R>$ of states whose projections 
onto the $rs$-plane lie in $\D_{rs}$:
$$
\alpha'-\alpha_1-\alpha_2=A, \qquad
\beta'-\beta_1-\beta_2=B,
$$
where $A+B\le (C_0-1)\xi$ and $A=-\xi$ or $B=-\xi$.
Moreover,
\begin{alignat*}{2}
&|\delta'|-|\delta_1|-|\delta_2|
  +(\delta_{\alpha_1}+\delta_{\alpha_2}-\delta_{\alpha'})
  +(\delta_{\beta_1}+\delta_{\beta_2}-\delta_{\beta'})
  \le -M (A+B),\\
&|\delta'|-|\delta_1|-|\delta_2|\le G(D_2+D_3),\\
&\mu'-\mu_1-\mu_2\le G D_3,\qquad
  \eta'-\eta_1-\eta_2\le G D_3.
\end{alignat*}
\end{lem}
Here $D_2$ and $D_3$ are quadratic wave interaction terms which are dominated by the decrease in $Q$, as is shown in \cite{hyp:Temp81}.  Our concern is with the constant $M$.

\begin{lem}\label{lemtemp:3.2}
Let $\D=\D_{rs}\times [S_*, S^*]$ in $rsS$-space, there exists
$\eps_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:
\begin{equation*}
(\text{Change in strength at $\eps$})
\le (\text{Change in strength at $\eps=0$})
+G\eps Q.
\end{equation*}
\end{lem}

\begin{lem}\label{lemtemp:3.3}
There exists an $M>0$ depending only on $\D_{rs}$ and an $\eps_1>0$
such that, if $w_L,w_R \in \D$, the associated Riemann problem is
solvable for each $\eps \in \br{0,\eps_1}$, and the waves in these
Riemann problems satisfy the following estimates:
\begin{alignat*}{2}
&M\, \Var_r^-(\alpha) > \Var_{rsS}(\alpha),&\qquad&\Var_s^-(\alpha) < \alpha,\\
&M\, \Var_s^-(\beta) > \Var_{rsS}(\beta),&\qquad&\Var_r^-(\beta) < \beta,\\
&2\Var_r^+(\mu) > \Var_{rsS}(\mu),&\qquad &\Var_r^-(\mu)=0,\\
&2\Var_s^+(\eta) > \Var_{rsS}(\eta),&\qquad &\Var_s^-(\eta)=0,\\
&\Var_s(\mu) < \frac{1}{4}\mu,&\qquad &\Var_r(\eta) < \frac{1}{4}\eta.
\end{alignat*}
\end{lem}
Note that $r$ and $s$ are constant on contact waves so that
$\Var_{rsS}(\delta)=\n{\delta}$ and $\Var_{rsS}(\zeta)=\n{\zeta}$.

\begin{rem}
The above lemma imposes additional requirements on $M$, namely,
\begin{equation}\label{eq:3.13}
    M \geq 1 + \myf{\Var_{S}(\alpha)}{\Var_{r}(\alpha)}
    + \myf{\Var_{s}(\alpha)}{\Var_{r}(\alpha)},
\end{equation}
for all 1--shock curves in $\D$ and for $\eps\in\br{0,\eps_{1}}$. 
The constant $M$ must also satisfy a similar requirement with regard 
to 3--shock curves in $\D$.  Proposition 3.2 of \cite{hyp:Temp81} 
requires that, at $\eps=0$,
\begin{equation}\label{eq:3.14}
\begin{aligned}
 M &\geq \myf{\tilde{M}}{1-C_0},\\
 \tilde{M} &= 
2\max_{\eps=0}\p{\sup_{1-shocks}\n{\myf{dS}{dr}},\sup_{3-shocks}
\n{\myf{dS}{ds}}}.
\end{aligned}
\end{equation}
Here $\frac12\leq C_0<1$ is a constant defined in Lemma 3.1 and
Proposition 3.1 of \cite{hyp:Temp81}.  Note that we have corrected a
typographical error in \cite{hyp:Temp81} which defined $M =
(1-C_0)\tilde{M}$.

We require, in addition to the above, that
\begin{equation*}
M\geq  
\max_{0\leq\eps\leq\eps_1}\p{\sup_{1-shocks}\n{\myf{dS}{dr}},\sup_{3-shocks}
\n{\myf{dS}{ds}}}
\end{equation*}
holds for all 1-shock  and 3-shock curves in $\D$,  making $M$ larger as
needed. 
As a result, 
with $M_0 \leq \myf{1}{2M}$, $\alpha-M_0\delta_{\alpha}$ is monotone along
all 1--shock curves in $\D$ and $\beta-M_0\delta_{\beta}$ is monotone
along all 3--shock curves in $\D$, for $0\leq \eps\leq \eps_{1}$.
\end{rem}

\begin{lem}\label{lemtemp:3.4}
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 $\eps \in \br{0,1}$ with solution waves that
%%
%%Change
%satisfy the estimates of Lemma \ref{lemtemp:3.2}.
%To
%%Later we have made the same changes as well
%%
satisfy the estimates of Lemma {\rm \ref{lemtemp:3.2}}.
\end{lem}

Let $E$ be an arbitrary compact set in $rsS$-space.  
Let $\tilde{\D}=\tilde{\D}_{rs}\times\br{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 \ref{lemtemp:3.4}. 
Choose $\tilde{\eps}_{1}>0$ so that Lemmas \ref{lemtemp:3.1} to \ref{lemtemp:3.4} 
apply to $\tilde{\D}$.  Then Riemann problems $\angbr{w_{L},w_{R}}$ in
$\tilde{\D}$ are uniquely solvable if $\eps \leq \tilde{\eps}_{1}$,
or if $w_{L},\ w_{R}\in B_{C_{1}}(w)$ for some $w\in E$ and 
$0\leq\eps\leq 1$.  Let $\tilde{V}(w_L,w_R)$ denote the variation in the solution of
one of these Riemann problems at any time $t>0$.

\begin{lem}\label{lemtemp:3.5}
There exists a constant $K_{0}>1$
such that, for $(w_L, w_R) \in \tilde{\D}\times\tilde{\D}$,
\begin{equation*}
\tilde{V}(w_L,w_R)\leq K_{0}\onenorm{w_{L}-w_{R}}.
\end{equation*}
\end{lem}

Note that for the reacting equations, $\zeta$ or $Y$ waves are
decoupled from the acoustic and entropy waves so that Lemmas
{\ref{lemtemp:3.4}} and {\rm \ref{lemtemp:3.5}} generalize 
to our problem with
no change to the constants $C_{1}$ and $K_{0}$.

Lemmas \ref{lemtemp:4.1} and \ref{lemtemp:4.2} below correspond 
to Lemmas 4.1 and 4.2 of \cite{hyp:Temp81}.
\begin{lem}\label{lemtemp:4.1}
Let $\D=\D_{rs}\times\br{S_{*},S^{*}},\ \eps_1$, and $M$ satisfy the
conditions of Lemmas {\rm \ref{lemtemp:3.1}} 
to {\rm \ref{lemtemp:3.3}}.
Let $J_1$ be a mesh curve which evolves at $\eps \leq \eps_1$ from 
initial data $w_0(x)$ of variation $V_0$ through mesh curves J 
such that $w(U(J))\subset \D$.  Then the following
estimates hold for any $M_0\leq\myf{1}{2M}$ and 
%%
%Change
%any $J_2$ an immediate successor of $J_1$:
%To
%%
any $J_2$ which is an immediate successor of $J_1$:
\begin{alignat*}{3}
& (i) &&\Var_{rs}(J_2) \leq 20 L_{T0}(J_2);&&\\
& (ii)\quad &&\Var_{rsS}(J_2) \leq K L_T(J_2),&\mbox{ where }
  &K=\myf{20}{M_0^2},\ M_0 \leq \myf{1}{2M}.
\end{alignat*}
\end{lem}
Lemma \ref{lemtemp:4.1} easily generalizes to:
\begin{alignat*}{3}
& (i) &&\Var_{rs}(J_2) \leq 20 L_0(J_2);&&\\
& (ii)\quad &&\Var_{rsSY}(J_2) \leq K L(J_2),
\end{alignat*}
since $K>1$ and $\Var_Y(\zeta_i) = M_4\n{\zeta_i}$.  An examination of
the proof of this Lemma (Lemma 4.1 in \cite{hyp:Temp81}) 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 $\eps \in \br{0,1}$.  

\begin{lem}\label{lemtemp:4.2}
 Let $\tilde{\D}=\tilde{\D}_{rs}\times\bigbr{\utS,\tS}\subset E,
        \ \teps_1,\ \tM$, and $C_1$ satisfy the
conditions of Lemmas {\rm \ref{lemtemp:3.1}} to {\rm \ref{lemtemp:3.5}}.
Let $J_1$ be a mesh curve which evolves at $\eps \in \br{0,1}$
from initial data of total variation $V_0$, through mesh curves
$J<J_1$ for which $U(J) \subset B_{C_1}(w)$ for some $w\in 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 \leq \min\bp{\teps_1^2/2,\bp{2\tM}^{-1}} = \tM_0$:
\begin{equation*}
\Var_{rsS}(J_2) \leq K L_T(J_2),\qquad K = \myf{20}{M_0^2}.
\end{equation*}
\end{lem}
This lemma can similarly be generalized to:
\begin{equation*}
\Var_{rsSY}(J_2) \leq K L(J_2).
\end{equation*}

The main theorem in \cite{hyp:Temp81} (Theorem 4.1) is as follows:
\begin{lem}\label{thmtemp:4.1}
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)\subset E$ with $\Var(w_0(x))=V_0\leq N$ and
$Y\equiv 0$, if $\eps V_0 < C$, then $F_T(J)$ is non-increasing,
$F_T(J_0)<\infty$, and there exists a global weak solution to
{\rm (\ref{1.4}a,b,c)} with initial data $w_0(x)$.
\end{lem}

The proof of Lemma \ref{thmtemp:4.1} is in two different settings.
For small $\eps$, a compact convex set
$\tilde{\D}=\tilde{\D}_{rs}\times\br{\tilde{S}_{*},\tilde{S}^{*}}$ in
$rsS$--space is chosen, together with constants $\eps_{1},\
\tilde{M}_{0},\ K_{0}$, and $C_{1}$, such that $\tilde{\D}$ contains
all points within distance $C_{1}$ of the set $E$ (see
\cite{hyp:Temp81}, p. 142), and such that 
%%
%Change
%Lemmas \ref{lemtemp:3.1} to \ref{lemtemp:3.5}, and \ref{lemtemp:4.2} hold. 
%To
%%
Lemmas \ref{lemtemp:3.1}--\ref{lemtemp:3.5} and Lemma \ref{lemtemp:4.2} hold. 
Let $\D_{rs}$ be the
set of points in the $rs$--plane within a distance of $120K_{0}N$ of
$\tilde{\D}_{rs}$.  New values of $M$ and $C_{0}$ are chosen so that
%%
%Change
%Lemmas \ref{lemtemp:3.1} to \ref{lemtemp:3.3} and \ref{lemtemp:4.1}
%To
%%
Lemmas \ref{lemtemp:3.1}--\ref{lemtemp:3.3} and Lemma \ref{lemtemp:4.1}
are satisfied on $\D_{rs}$ for $\eps\in \br{0,\eps_{1}}$.  These
Lemmas require $M_{0}\le \myf{1}{2M}$ so that
\begin{equation*}
  M_{0} = \min\mybrace{\myf{1}{2M},\ \tilde{M}_{0}},\qquad
    K = \max\mybrace{\myf{20}{M_{0}^{2}},\ 1}
\end{equation*}
are defined.
Choose $S_{*}< \tilde{S}_{*}-6K_{0}KN$, and
$S^{*}>\tilde{S}^{*}+6K_{0}KN$, and define
$\D= \D_{rs}\times\br{S_{*},S^{*}}$.  
The mesh parameters $(l,h)$ are chosen so that the Courant-Friedrichs-Levy 
condition \eqref{2.15} is satisfied for values in $\D$. 
Later, a number $\eps_{0}\leq \eps_{1}$ is chosen. 
Then, for $0\le \eps \le \eps_{0}$, bounds on $F$ imply bounds 
on $\Var_{rsS}\bp{U^{h}(J)}$ which imply that $w(U^{h}(J))\subset \D$.  
The functional $F$ is proved to be finite and non-increasing for solutions which remain 
in $\D$.

For large $\eps$, that is, $\eps\in \br{0,1}$, 
smaller initial total variation is required, so that bounds
on $F$ ensure that $\Var_{rsS}\bp{U^{h}(J)} <C_{1}$.  This
implies that $w(U^{h}(J)) \subset B_{C_{1}}(\tilde{w})\subset \D$
for some $\tilde{w}\in E$ --- for example,
we could choose $\tilde{w} = \lim_{x\to-\infty}w_{0}(x)$.
Note that the ball $B_{C_1}(\tilde{w})$ should be defined in the $\ell_1$ norm.

We will need an estimate of the form
$F(J) \leq C\Var_{rsSY}(J)$. For this purpose, we extract the
following lemma from the proof of Lemma \ref{thmtemp:4.1}
\cite{hyp:Temp81}.

\begin{lem}\label{lem:3.13a}
Let $E$ be a compact set in $(r,s,S)$ space, and let $N\in(1,1/\eps]$ be given.
Suppose that the mesh curves $J$ and $J_0$ satisfy the conditions of 
Lemmas {\rm \ref{lemtemp:4.1}} or {\rm \ref{lemtemp:4.2}}.  
Then there are constants
$K_1(E,N)>1$ and $\eps_0(E,N)\in(0,1)$ such that, for initial
data $w_0(x)$ with $\Var_{rsS}(w_0) = V_0 <N$,  we have
\begin{alignat}{3}
F(J_0) &\leq 5 K_0 V_0,&\qquad&0<\eps\leq \eps_0,\\
F(J_0) &\leq K_1 V_0,&\qquad&\eps_0<\eps\leq 1.
\end{alignat}
Moreover, for approximate solutions $U^h(x,t)$, for which
$\Var_{rsSY}\bp{U^h(J)} \leq N$, we have
\begin{alignat}{3}
 F(J) &\leq 5 \Var_{rsSY}\bp{U^h(J)},&\qquad&0 
        <\eps\leq \eps_0,\label{eq:3.17a}\\
 F(J) &\leq \myf{K_1}{K_0}\Var_{rsSY}\bp{U^h(J)},
            &\qquad&\eps_0<\eps\leq 1.\label{eq:3.18a}
\end{alignat}
\end{lem}
\begin{proof}
We denote $L\bp{U^{h}(J)}$ and $\Var\bp{U^{h}(J)}$ by $L(J)$ and
$\Var(J)$.  Note that by \eqref{eq:3.11}, $p_{i}\leq
2\Var_{rsS}\p{p_{i}}$ for any wave $p_{i}$, $p_{i}\neq \zeta_{i}$,
while, by \eqref{Vardef}, $M_{4}\n{\zeta_{i}} = \Var(\zeta_{i})$.
Hence, for any mesh curve $J$,
\begin{equation}\label{A}
  \begin{gathered}
    \sum_{J,\ p_i\neq\zeta_i}\p{p_{i}+M_{4}\n{\zeta_{i}}} \leq \sum_{J}2\Var(p_{i}) = 2\Var_{rsSY}\p{J},\\
    L\bp{J}  \leq \sum_{J,\ p_i\neq\zeta_i}\p{p_{i}+M_4\n{\zeta_i}}+ 
\alpha_{\infty}+\beta_{\infty}+\mu_{\infty}+\eta_{\infty}+\n{\delta_{\infty}}
 +M_{4}\n{\zeta_{\infty}} \leq 4\Var_{rsSY}\bp{J}.
\end{gathered}
\end{equation}
Let $K = \myf{20}{M_{0}^{2}}>1$. Choose $G>1$ so that Lemmas
\ref{lemtemp:3.1} and \ref{lemtemp:3.2}
hold.  Let
\begin{equation}\label{C}
\begin{aligned}
  M_{1} &= 8G,\\
  M_{2} &= 8G + 2 M_{1}K(4K_{0}N)G,\\
  M_{3} &= 8G + 2 M_{1}K(4K_{0}N)G + 4 M_{2}K(4K_{0}N)G,\\
  K_{1} &= 4K_0 + 4M_3K_0^2.
\end{aligned}
\end{equation}
Lemma \ref{lemtemp:3.5} implies that 
$\Var_{rsSY}\bp{J_{0}}\leq K_{0}V_0$ if
$\eps\leq \eps_{0}$ or if $V_0 \leq C_{1}$.  Then, by \eqref{A},
\begin{equation}\label{E}
\begin{aligned}
 L(J_{0}) &\leq 4\Var\bp{J_{0}} \leq 4K_{0}V_0,\\
 L_{0}(J_{0}) &\leq 4\Var\bp{J_{0}} \leq 4K_{0}V_0.
\end{aligned}
\end{equation}
Furthermore, $Q(J_{0}) \leq M_{3}\p{2\Var(J_{0})}^{2} \leq 
4M_{3}K_{0}^{2}V_0^{2} \leq 4M_{3}K_{0}^{2}NV_0$.  Thus, for
$0 \leq \eps \leq 1$, since $\eps N \leq 1$,
\begin{equation}\label{G}
    F(J_{0}) \leq 4K_{0}V_0 + \eps 4M_{3}K_{0}^{2}NV_0\leq 
 K_{1}V_0.
\end{equation}
We now require $\eps_{0}\leq \myf{1}{K_{1}}$. Then 
$\eps\leq\eps_{0}\implies \eps Q(J_{0}) \leq 4\eps M_{3}K_{0}^{2}NV_0
\leq \frac{1}{2}\eps K_{1}V_0 \leq V_0$, and
\begin{equation}\label{H}
    \begin{aligned}
    F(J_{0}) &= L\bp{J_{0}} + \eps Q\bp{J_{0}} \leq 4 K_{0}V_0 + V_0 \leq 
    5K_{0}V_0,\\
    F_{0}(J_{0}) &= L_{0}\bp{J_{0}}+\eps Q\bp{J_{0}} \leq 5 K_{0}V_0,
    \end{aligned}
  \qquad\eps \leq \eps_{0}.  
\end{equation}
The proof of estimates \eqref{eq:3.17a} and \eqref{eq:3.18a} 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.
\end{proof}

Note that $\eps_{0}$ and $M_{i}$ depend only on $E$ and $N$. We now
continue with the derivation of the constant $C(E,N)$ from Lemma
\ref{thmtemp:4.1}, with modifications for small $\eps$ based on the
estimate $F(J_0) \leq 5 K_0 V_0$ .

For $\eps \leq \eps_0$, we define
\begin{equation}
\label{eq:C3}
C_0(E,N) =\myf{1}{20 M_3 K K_0}\min\p{\myf{1}{3 G -1},1-C_0}.
\end{equation}
Using Lemma \ref{lem:3.13a}, if $\eps V_0 < C_0(E,N)$, then
$\eps F(J_0) \leq \eps 5 K_0 V_0 <5 K_0 C_0(E,N)$ so that
\begin{equation}
\label{eq:newF}
\eps F(J_0)  < \myf{1}{4 M_3 K}\min\p{\myf{1}{3 G -1},1-C_0}.
\end{equation}
For $\eps_0 < \eps \leq 1$, we proceed as follows.
Let  $C_{2}= \myf{\eps_{0}C_{1}}{K_{1}K}$.  If $\eps\geq\eps_{0}$
and $\eps V_0 < C_{2}$, then $KK_{1}V_0 < C_{1}$.  Thus, by \eqref{G},
\begin{equation}\label{I}
KF(J_{0}) < C_{1},\qquad \eps \geq \eps_{0}.
\end{equation}
Since $C_{1}\leq 1$ and $K\geq 1$, $\myf{C_{1}}{K} \leq 1$. This, plus
$K_{0}N > 1$, implies that $F\bp{J_{0}} \leq 5 K_{0}N$ for
$\eps_{0}\leq \eps \leq 1$.  Let
\begin{equation}\label{eq:C(E,N)}
 C_3 = \myf{1}{4M_3KK_1}\min\p{\myf{1}{3G},\ 1-C_{0}},
\end{equation}
and let $C(E,N) = \min(C_2,C_3)$. Then, if $\eps V_0 < C(E,N)$, we have
$\eps F\bp{J_{0}}\leq \eps K_{1}V_0 \leq K_{1}C(E,N)$ so that
\eqref{eq:newF} holds for $\eps_0 < \eps \leq 1$.


\subsection{TV Stability of the Reaction Step} 
We now prove the $TV$ 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=-\infty$, 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.

\begin{lem}\label{lem:3.5a} There exists $T'(\eps)$ such that, 
for $\eps\in\br{0,\eps_{0}}$, $T(w)\ge T'(\eps)$ for all $w\in \D$, 
and, for $\eps \in\br{\eps_{0},1}$, $T(w)\ge T'(\eps)$ for all
$w\in B_{C_{1}}(\tilde{w})$.  Moreover, $T'(\eps)$ has a positive
minimum for $\eps\in\br{0,1}$.
\end{lem}
\begin{proof}
Note that $\ln(p)=\myf{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 $\eps$ such that, for all $w\in \D$,
$0<p'\leq p(w)\leq p''< \infty$.  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 $\R\times\p{0,\infty}\times\R$.

Since the map from $(u,v,S)$ to $(r,s,S)$ is a $C^{4}$ diffeomorphism
%%
%Change
%(see Section 2), and $T = \pd{}{S}e(v,S,\eps)$ is $C^{4}$, we have
%To
%%
(see Section 2), and $T = e_S(v,S,\eps)$ is $C^{4}$, we have
that $T$ is a continuous function on the compact set
$\D\times\br{0,1}$.  
Therefore, for each $\eps\in\br{0,1}$, $T$ has a minimum value 
$T'(\eps)$ on this set.  For $\eps\in\br{\eps_{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
$\eps=0$, $T = \myf{1}{R}$), the minimum value of $T'$ on $\br{0,1}$
is strictly positive.
\end{proof}

Note that $M_{0}$ was defined so that, if the strength of $\alpha$
varies with fixed left state, $\alpha -M_{0}\delta_{\alpha}$
is strictly increasing in $\alpha$ for $0\leq\eps\leq\eps_{1}$.
We require a similar property to hold for $\eps\in\br{0,1}$.

\begin{lem}\label{lem:3.60}
Let $\alpha$ and $\talpha$ be the strengths of two {\rm 1}--shock waves in $\D$, as
defined in \eqref{eq:3.11} for a common value of $\eps\in\br{0,1}$.
Let the left states of $\alpha$ and $\talpha$ be $U_{L}$ and
$\tilde{U}_{L}$, respectively.
Then there is a constant $C$ such that
\begin{equation}\label{eq:3.140}
\n{\alpha-M_{0}\delta_{\alpha} -\bp{\talpha
-M_{0}\delta_{\talpha}}} \leq \n{\alpha-\talpha} +
C\onenorm{U_{L}-\tilde{U}_{L}}\talpha.
\end{equation}
Similarly, if $\beta$ and $\tbeta$ are the strengths of
two {\rm 3}--shock waves in $\D$ with right states $U_{R}$ and $\tilde{U}_{R}$, 
respectively, then there is a constant $C$ such that
 \begin{equation}\label{eq:3.141}
\n{\beta-M_{0}\delta_{\beta} -\bp{\tbeta
-M_{0}\delta_{\tbeta}}} \leq \n{\beta-\tbeta} +
C\onenorm{U_{R}-\tilde{U}_{R}}\tbeta .
\end{equation}
\end{lem}

\begin{proof}
Let $\alpha_{0}$ be a 1--shock wave with the same left state as
$\alpha$ and such that $\alpha_{0}=\talpha$ as wave strengths.
Since $S$ is monotone along any shock curve \cite{hyp:Weyl49},
$\Var_{S}(\alpha)=\delta_{\alpha}$. 
We have that 
 \begin{equation}
    \begin{aligned}
 \n{\alpha-M_{0}\delta_{\alpha} -\bp{\talpha  -M_{0}\delta_{\talpha}}}
 &\leq \n{\alpha-M_{0}\delta_{\alpha} -
 \p{\alpha_{0}-M_{0}\delta_{\alpha_{0}}}}
 + M_{0}\n{\delta_{\alpha_{0}}-\delta_{\talpha}}.
\end{aligned}
\end{equation}
In accordance with Lemma \ref{lemtemp:4.2}, $M_{0}$ is chosen so that
$M_0 \leq\myf{\eps_{1}}{2}$.  Thus, for $\eps\leq M_{0}<\eps_{1}$,
since $\alpha$ and $\alpha_{0}$ have the same left state, and since
$\alpha-M_{0}\delta_{\alpha}$ is monotone in $\alpha$ with derivative
less than $1$,
\begin{equation*}
 \n{\alpha-M_{0}\delta_{\alpha} -
 \p{\alpha_{0}-M_{0}\delta_{\alpha_{0}}}} \leq
 \n{\alpha-\alpha_{0}}=\n{\alpha-\talpha}.
\end{equation*}
For $\eps \geq M_{0}$,
\begin{equation*}
\begin{split}
 &\n{\alpha-M_{0}\delta_{\alpha} -
\p{\alpha_{0}-M_{0}\delta_{\alpha_{0}}}}\\
&\ \qquad=\n{(1+\eps)\bp{\Var_{r}(\alpha)-\Var_{r}(\alpha_{0})} +
 \eps\bp{\Var_{s}(\alpha)-\Var_{s}(\alpha_{0})}
 + \p{\eps-M_{0}}\p{\delta_{\alpha}-\delta_{\alpha_{0}}}}.\\
\end{split}
\end{equation*}
Since $\alpha_{0}$ and $\alpha$ have the same left state, each of the
terms $\bp{\Var_{r}(\alpha)-\Var_{r}(\alpha_{0})}$,
$\bp{\Var_{s}(\alpha)-\Var_{s}(\talpha_{0})}$, and
$\p{\delta_{\alpha}-\delta_{\alpha_{0}}}$ have the same sign. Thus,
\begin{equation*}
     \begin{aligned}
 &\n{\alpha-M_{0}\delta_{\alpha}
 -\bp{\alpha_{0}-M_{0}\delta_{\alpha_{0}}}}\\
 &\qquad\leq \left\lvert(1+\eps)\p{\Var_{r}(\alpha)-\Var_{r}(\alpha_{0})}
 + \eps\p{\Var_{s}(\alpha)-\Var_{s}(\alpha_{0})}+
 \eps\p{\delta_{\alpha}-\delta_{\alpha_{0}}}\right\rvert\\
 &\qquad= \n{\alpha-\alpha_{0}}= \n{\alpha-\talpha}.
  \end{aligned}
\end{equation*}
Let
$G\bp{\tilde{U}_{L},U_{L},\talpha}=\delta_{\alpha_{0}}-\delta_{\talpha}.$
Observe that $G\bp{\tilde{U}_{L},U_{L},\talpha}=0$ if $\talpha=0$ or
$U_{L}=\tilde{U}_{L}$.  Since $G$ is $C^{3}$ for $\talpha\geq 0$, an
argument similar to the proof of Lemmas \ref{lem:3.0.1} and
\ref{lem:3.2} shows that there is a constant $C$ such that
$\bign{G\bp{\tilde{U}_{L},U_{L},\talpha}} \leq
C\talpha\bigonenorm{U_{L}-\tilde{U}_{L}}$.
\end{proof}


One consequence of Lemma \ref{lemtemp:3.3} is the following lemma.

\begin{lem}\label{lem:3.5}
Let $J_{\sigma}$ be a connected segment of a mesh curve $J$, and let
$U_{\sigma}$ be the restriction of an approximate solution $U^{h}$ to
$J_{\sigma}$.  Then there are positive constants $c$ and $C$ such that,
for all $U_{\sigma}$ with $w\bp{U_{\sigma}}\in\D$,
\begin{equation*}
c\ \Var_{rsSY}\bp{U_{\sigma}} \leq \sum_{J_{\sigma}}\p{\alpha +\beta+\mu
+\eta+\n{\delta}+\n{\zeta}} \leq C\ \Var_{rsSY}\bp{U_{\sigma}}.
\end{equation*}
\end{lem}

\begin{cor}\label{cor:3.6}
Let $J$ be any mesh curve, and let $U_{\sigma}$ be the restriction of
an approximate solution $U^{h}$ to $J$.  
Let
\begin{equation*}
L_{1}(J) = \sum_{J} \p{\alpha +\beta+\mu+\eta+\n{\delta}+M_{4}\n{\zeta}},
\end{equation*}
where the sum is understood to include the ``extra'' Riemann problem
for $\bp{U_{-\infty},U_{\infty}}$.  Then there are positive constants
$c$ and $C$ such that, for all $U_{\sigma}$ with $w\bp{U_{\sigma}}\in\D$,
\begin{equation*}
c\ \Var_{rsSY}\bp{U_{\sigma}}\leq L_{1}(J)\leq C\ \Var_{rsSY}\bp{U_{\sigma}}.
\end{equation*}
\end{cor}
Note that \eqref{eq:3.11} implies that, for $\eps \leq 1$,
\begin{equation}
L(J) \leq 2\Var_{rsSY}\bp{U_{\sigma}}.
\end{equation}
This estimate, together with the conclusions of Lemmas
\ref{lemtemp:4.1} and \ref{lemtemp:4.2}, implies that $L(J)$ is equivalent
to $\Var_{rsSY}\bp{U_{\sigma}}$, just as Corollary \ref{cor:3.6} 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
\ref{thmtemp:4.1} on $V_0$, namely $V_0 \leq N$ and $\eps V_0 < C(E,N)$;
these conditions are refined in Theorem \ref{thm:3.1}.  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 \ref{lemtemp:4.2} implies that $\Var_{rsS}(J) \leq K
L(J)$.  By definition, $L(J) \leq F(J)$.  Thus, we require $K F(J)
\leq N$ and $\eps K F(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.
\begin{lem}\label{lem:3.8}
In the $(r,s,S,Y)$-coordinates,
for the system {\rm \eqref{1.4}}, we have
  \begin{equation*}
        G(U) = q\bc(U)\cdot Y\phi(T),
        %\label{eq:3.0a}
  \end{equation*} 
where $\bc(U)$ is the vector
$\bp{-\myf{1}{c_{v}T},\myf{1}{c_{v}T},\myf{1}{T},-\myf{1}{q}}^T$.
If the gas is not polytropic, we let $c_{v}$ denote $\tilde{e}'(T)$.
\end{lem}
\begin{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,
\begin{equation}
  \myf{dr}{dY} = \myf{q}{c_{v}T}.
\end{equation}
Similar calculations hold for $r$ and $S$.
\end{proof}

We will find it convenient to denote
$\bp{-\myf{1}{c_{v}T},\myf{1}{c_{v}T},\myf{1}{T}}^\top$ by $\bc_{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'' $\tilde{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 $\tilde{J}$
to represent values before the reaction
step and $J$ to represent values after the reaction step.

\begin{lem}\label{lem:3.7} Let $\D$ be given as above.
Let $\F(\tau,U)$ be a consistent, total variation stable, and conservative 
approximation algorithm for the solution of the initial value 
problem {\rm (\ref{ivp})}.  
There is a constant $C$, which depends only on $\D$, such that 
if 
$U(J_k)$ has evolved at $\eps$ from
initial data of variation $V$ through mesh curves $J$ for which
$w(U(J))\subset \D$, then, for each $k$, to first order in
$h$,
\begin{equation*}%\label{eq:3.17}
    \begin{aligned}
L\bp{J_{k}} 
&\leq L(\tilde{J_k})+C hq\supnorm{Y_0}e^{-\Phi kh}L(\tilde{J_k}),\\
L_{0}\p{J_{k}} &\leq L_{0}\bp{\tilde{J}_{k}}+
C hq\supnorm{Y_0}e^{-\Phi kh}L(\tilde{J_k}).
\end{aligned}
\end{equation*}
\end{lem}

\begin{proof}
Note that  
\begin{equation}
 \begin{aligned}
 L\bp{J_{k}}-L\bp{\tilde{J_{k}}}  =& \sum_{-\infty<i< \infty}
 \bigl\{\alpha_{i}-M_{0}\delta_{\alpha_{i}}
 +\beta_{i}-M_{0}\delta_{\beta_{i}} +
 M_{0}\n{\delta_{i}}+\eps\p{\mu_{i}+\eta_{i}} +M_{4}\zeta_{i}
     \bigr.\\
 &\qquad\quad \bigl. -\bp{\talpha_{i}-M_{0}\delta_{\talpha_{i}} +
 \tbeta_{i}-M_{0}\delta_{\tbeta_{i}}+M_0\bign{\tdelta_{i}}
 +\eps\p{\tmu_{i}+\teta_{i}}+M_4\tzeta_{i}}\bigr\}\\
 &+\alpha_{\infty}+\beta_{\infty}+\mu_{\infty}+\eta_{\infty}
  +\n{\delta_{\infty}} + M_{4}\zeta_{\infty}\\
  &\qquad-\bp{\talpha_{\infty}+\tbeta_{\infty}+\tmu_{\infty}+\teta_{\infty}
  +\bign{\tdelta_{\infty}} + M_{4}\tzeta_{\infty}}.
\end{aligned}
\end{equation}
We denote the vector of the  end-terms by $\bz_{\infty}=B\bp{U_{-},U_{+}}$.
>From Lemma \ref{lem:3.60}, this sum is less than
\begin{equation}
    \begin{aligned}
  \sum_{-\infty<i\leq\infty}&\Bigl\{\n{\alpha_{i}-\talpha_i}
+ C\N{\F\bp{h,U_i}-U_i}\talpha +\bign{\beta_{i}-\tbeta_i}
    +C\N{\F\bp{h,U_{i+1}}-U_{i+1}}\tbeta\Bigr.\\
  &\qquad\qquad\qquad\qquad+ \Bigl.\Bign{\bign{\delta_{i}}-\bign{\tdelta_i}} +\n{\mu_{i}-\tmu_{i}}
  +\n{\eta_{i}-\teta_{i}}+M_{4}\n{\bign{\zeta_{i}}-\bign{\tzeta_{i}}}\Bigr\}\\
  &\leq \sum_{-\infty<i\leq \infty}\Bigl\{\onenorm{\Gamma\bp{U_{i},\bz_{i},h}
-\Gamma\bp{U_{i},\bz_{i},0}}
+C\N{\F\bp{h,U_i}-U_i}\bp{\talpha_{i}+\tbeta_{i-1}}\Bigr\}.
\end{aligned}
\end{equation}
In the limit as $h\to 0+$, we obtain
\begin{equation}
  L\bp{J_{k}}-L\bp{\tilde{J_{k}}}\leq \sum_{-\infty<i\leq \infty}
 \left\{ \onenorm{\pd{\Gamma}{h}\bp{U_{i},\bz_{i},0}}h
  + C q Y_{i}\phi\p{T_{i}}\bp{\talpha_{i}+\tbeta_{i-1}}h\right\}.
\end{equation}
Note that 
\begin{equation}\label{eq:3.23}
  \sum_{-\infty<i\leq \infty}C qY_{i}\phi\p{T_{i}}
  \bp{\talpha_{i}+\tbeta_{i-1}}h \leq
  C q\supnorm{Y}L_{1}\bp{\tilde{J}_{k}}h.
\end{equation}
We calculate (recall \eqref{eq:3.0.0})
\begin{equation}\label{3.36a}
\pd{\Gamma}{h}(U_i,\bz_{i},0) = \pd{B}{U_L}\pd{\F}{h}(0,U_i)
        + \pd{B}{U_R}\pd{\F}{h}\bp{0,H\p{U_i,\bz_{i}}}.
\end{equation}
For brevity, we note that $U_{i+1}=H\bp{U_{i},\bz_{i}}$.
By Lemmas \ref{lem:3.3} and \ref{lem:3.8},
\begin{equation*}
\pd{\F}{h}(0,U) = G(U)=q\bc(U) Y\phi(T).
\end{equation*}
We denote $\bc\bp{U_{i}}$ by $\bc_{i}$. Thus,
\begin{equation}\label{eq:3.48}
\begin{aligned}
&\pd{\Gamma}{h}(U_i,\bz_{i},0)
   = q\p{\pd{B}{U_L}\bc_{i} Y_i\phi(T_i)+\pd{B}{U_R}\bc_{i+1}
   Y_{i+1}\phi\bp{T_{i+1}}}\\
&= q\p{Y_i \p{\pd{B}{U_L}\phi(T_i)\bc_i + \pd{B}{U_R}
     \phi\bp{T_{i+1}}\bc_{i+1}} + \p{Y_{i+1} -
     Y_i}\pd{B}{U_R}\phi\bp{T_{i+1}}\bc_{i+1}}.\\
\end{aligned}
\end{equation}
As we noted in Section \ref{sect:3},  $Y$ is decoupled from $(r,s,S)$
in the solution of the non--reacting Riemann problem.  This
means that $\pd{B}{U_{L}}$ and $\pd{B}{U_R}$ are
block diagonal $4\times 4$ matrices with the upper left
$3\times 3$ block relating to non-reacting gas dynamics, {\em i.e.,} the
derivatives of wave strengths for acoustic waves and entropy waves,
with respect to $r,s$, and $S$.  The remaining $1\times 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 
$\partial B_{4}/\partial Y_{R}$ and $-1$ for $\partial B_{4}/\partial Y_{L}$.
We have that
\begin{equation}\label{eq:3.42}
\pd{B}{U_R}\cdot\bc_{i+1} = \pd{B}{\p{r,s,S}_R}\cdot \bc_{g(i+1)}
-\p{\begin{smallmatrix} 0\\0\\0\\1\\
\end{smallmatrix}}.
\end{equation}
Due to this block structure, $\pd{B}{\p{r,s,S}_R}$ has no $Y$-wave
component.  
Thus the $(r,s,S)$ components of \eqref{eq:3.48} have the form:
\begin{equation}\label{eq:3.48a}
\begin{split}
\pd{\Gamma_{g}}{h}(U_i,\bz_{i},0)
        &= qY_{i}\left(\pd{B_{g}}{\p{r,s,S}_{L}} \phi(T_{i})\bc_{gi}
+\pd{B_{g}}{\p{r,s,S}_{R}}\phi(T_{i+1})\bc_{g(i+1)}\right)\\
&\qquad+q\p{Y_{i+1} - Y_i}\pd{B_{g}}{\p{r,s,S}_{i+1}}
\phi(T_{i+1})\bc_{g(i+1)}.
\end{split}
\end{equation}

Let 
\begin{equation}
A(U_{i},\bz_{i}) = \pd{B_{g}}{\p{r,s,S}_{R}} \phi(T_{i})\bc_{gi}
+\pd{B_{g}}{\p{r,s,S}_R}\phi\bp{T_{i+1}}\bc_{g(i+1)}.
\end{equation}
Then
\begin{equation}\label{eq:3.15}
\pd{\Gamma_{g}}{h}(U_i,\bz_{i},0)= qY_i A(U_i,\bz_{i})
+ q\p{Y_{i+1} - Y_i}\pd{B_{g}}{\p{r,s,S}_R}\phi\bp{T_{i+1}}\bc_{g(i+1)}.
\end{equation}
Since $B$ is Lipschitz continuous, together with its first and second
derivative, the function $A(U_L,\bz)$ is Lipschitz continuous,
together with its first derivatives.
Furthermore, if $\bz_{gi}=\mathbf{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,
\begin{equation*}
A_{g}(U_{g}, \mathbf{0}) =
\phi(T)\pd{B_g}{\p{r,s,S}}(U_{g},U_{g})\bc_{g}=\mathbf{0}.
\end{equation*}
Thus, if $C$ is a Lipschitz constant for $\bz_{g}\to A_{g}(U_{g},\bz_{g})$ 
on $\D$ 
taken with respect to the $\ell_{1}$ norm of $\bz_{g}$, the first term in
\eqref{eq:3.15} can be estimated by
\begin{equation*}
\onenorm{Y_i A_{g}(U_{gi},\bz_{gi})q}
\leq C Y_i \onenorm{\bz_{gi}}q.
\end{equation*}

We now examine the last term of (\ref{eq:3.15}), because it is not
bounded by $\supnorm{Y}\sum_{i}\onenorm{\bz_{i}}$.  
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 (\ref{eq:3.15}) has the form
\begin{equation}\label{eq:3.31}
\p{Y_{i+1}-Y_i}
\pd{B_{g}}{\p{r,s,S}_{R}}\bc_{g(i+1)}
\phi(T_{i+1})q.
\end{equation}

The fourth component of equation (\ref{eq:3.48}) is the equation for
the strength of the $Y$-wave.  This equation is
\begin{equation*}
\pd{}{h}\p{Y_{i+1}-Y_i} = Y_{i}
\bp{\phi(T_{i})-\phi(T_{i+1})}-\p{Y_{i+1}-Y_{i}}\phi(T_{i+1}),
\end{equation*}
so that
\begin{equation*}
\pd{}{h}\n{Y_{i+1}-Y_{i}} \leq Y_{i} C\n{\bz_{g}}
-\n{Y_{i+1}-Y_i}\phi(T_{i+1}).
\end{equation*}

Thus, the reaction step produces possible increases in total
variation, which are bounded by 
\begin{equation*}
C Y_i \onenorm{\bz_{gi}}qh
 +\bigonenorm{\pd{B_{g}}{\p{r,s,S}_{R}}\bc_{g(i+1)}}
\n{Y_{i+1}-Y_i}\phi(T_{i+1})qh.
\end{equation*} 
The reaction step also produces a decrease in total variation for the
$Y$ component --- the fourth component of
$\pd{\Gamma}{h}(U_{i},\bz_{i},0)$ --- in the amount
$\n{Y_{i+1}-Y_i}\phi(T_{i+1})h$.  We now use the decrease in the $Y$
component to offset the ``bad'' increase in other components
proportional to $\n{Y_{i+1}-Y_i}$.  Since $\pd{B_{g}}{\p{r,s,S}_{R}}$
is Lipschitz continuous and $\D$ is bounded, there exists a finite
upper bound $M_4$ for
$\bigonenorm{\pd{B_{g}}{\p{r,s,S}_{R}}\bc_{g(i+1)}}q$ on $\D$.  Thus
the the effect of the term \eqref{eq:3.31} on the $(r,s,S)$ components
of $\pd{\Gamma}{h}(U_{i},\bz_{i},0)$ is bounded by
$M_4\n{\zeta}\phi(T_{i+1})h$, and this increase is offset by a
decrease in the term $M_4 \n{\zeta}$.

\begin{comment}
Note that $M_{0}\leq \myf{1}{2M}$ and $M\geq 1$, so that $M_{0}\leq
\myf12$.  Also note that $M_0$ has been chosen so that for
$\eps\in\br{0,1}$,
\begin{equation}\label{eq:3.98}
\myf{1}{2}< \myf{d}{d\alpha}\p{\alpha-M_0\delta_{\alpha}} \leq 1.
\end{equation}
\end{comment}

Thus, the change in $L$ is estimated as follows:
\begin{equation}
\begin{aligned}
L(J_k) 
&\leq L\bp{\tilde{J_k}} + C q\supnorm{Y}L_{1}\bp{\tilde{J}_{k}}h\\
 &+\sum_{-\infty<i\leq\infty}\left\{
 \bp{\talpha_{i}+ \tbeta_{i} + \bign{\tdelta_{i}}
     +\tilde{\mu_i}+\tilde{\eta_i}} Cq Y_{i} 
    + M_{4}\bign{\tzeta_{i}}\phi(T_{i+})
   -M_{4}\bign{\tzeta_{i}}\phi(T_{i+})\right\}h\\
&\leq L(\tilde{J_k}) + 2C\bsupnorm{\tilde{Y}}L_{1}\bp{\tilde{J_{k}}}qh.
 \end{aligned}
\end{equation}
Then, using the equivalence of $L_{1}(J),\ L(J)$, and $\Var_{rsSY}(J)$,
and redefining the constant $C$, we obtain
\begin{equation}\label{eq:3.49}
L(J_k) \leq L(\tilde{J_k})\p{1 + C\bsupnorm{\tilde{Y}}qh}.
\end{equation}
Similarly,
\begin{equation}\label{eq:3.49a}
 \begin{aligned}
    L_{0}(J_k) &\leq L_{0}(\tilde{J_k})+
    C\bsupnorm{\tilde{Y}}L(J_{k})qh,\\
     L_{1}(J_k) &\leq L_{1}(\tilde{J_k})+C\bsupnorm{\tilde{Y}}L(J_{k})qh.\\
\end{aligned}
 \end{equation}

Furthermore,
on the mesh curve $J(k)$, by Lemma \ref{lem3.1}, we have
\begin{align*}\label{eq:3.50}
L\p{J_{k}}
&\leq L\bp{\tilde{J_k}}+C q h\supnorm{Y_0}e^{-\Phi k h}
L\bp{\tilde{J_k}},\\
L_{0}\p{J_{k}}
&\leq L_{0}\bp{\tilde{J_k}}+C q h\supnorm{Y_0}e^{-\Phi k h}
L\bp{\tilde{J_k}}.\qedhere
\end{align*}
\end{proof}

Since $L_{0} \leq L\leq F,\ F_{0}\leq 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\supnorm{Y_0}$ by $\tq$.
Now we estimate the functional $F\p{J_{k}}$.  Recall that
\begin{equation*}
F(J_{k}) = L(J_{k}) + \eps Q(J_{k}).
\end{equation*}
We have that, to first order in $h$,
\begin{equation*}
    \begin{aligned}
    Q(J_{k})
  &\leq Q\bp{\tilde{J_{k}}} +
   M_{1}\sum_{App}\left\{\p{p_{i}-\tilde{p}_{i}}\n{\tdelta_{j}}
  +\tilde{p}_{i}\bp{\n{\delta_{j}}-\bign{\tilde{\delta}_{j}}}\right\}\\
  &\qquad+  M_{2}\sum_{App}\left\{ \p{q_{i}-\tilde{q}_{i}}\tilde{R}_{j}
        +\tq_{i}\bp{R_{j}-\tilde{R}_{j}}\right\}\\
  &\qquad + 
   M_{3}\sum_{i<j}\left\{ \bp{\alpha_{i}-\tilde{\alpha}_{i}
  +\beta_{i}-\tilde{\beta}_{i}}\p{\talpha_{j}+\tbeta_{j}}
  +\p{\talpha_{i}+\tbeta_{j}}\bp{\alpha_{j}-\talpha_{j}
   +\beta_{j}-\tbeta_{j}}\right\}\\
   &\leq Q\bp{\tilde{J}_{k}} +
   M_{3}\p{L_{1}\bp{J_{k}}-L_{1}\bp{\tilde{J}_{k}}}L_{1}(\tilde{J}_{k}),
    \end{aligned}
\end{equation*}
where we have used $M_{1}< M_{2}< M_{3}$ as defined in \cite{hyp:Temp81}.
Next, using \eqref{eq:3.49a} and the equivalence of $L(J_{k}),\
L_{1}(J_{k})$, and 
$\Var_{rsSY}(J_{k})$,  we have
\begin{equation*}
\begin{aligned}
    Q(J_{k}) &\leq Q\bp{\tilde{J}_{k}} +
    M_{3}C\bsupnorm{\tilde{Y}_{k}} L_{1}\bp{\tilde{J}_{k}}^{2}q h\\
    &\leq Q\bp{\tilde{J}_{k}} +  C\bsupnorm{\tilde{Y}_{k}}
    F\bp{\tilde{J}_{k}}^{2}q h.
\end{aligned}
\end{equation*}

\begin{lem}\label{lem:3.7.1}
Let $J_k$ be a mesh curve between $t=k h$ and $t=(k+1)h$. Then
\begin{equation}
F(J_k) 
\le F(\tilde{J_{k}})\left(1 + K_3  \tq e^{-\Phi k h}\bp{1
   + \eps F\bp{\tilde{J_{k}}}}h\right),
\label{3.22}
\end{equation} 
where $K_3$ is a constant depending only on $\D$ and $\phi$,
independent of $\eps\in \br{0,1}$ and the mesh lengths $l$ and $h$.
\end{lem}

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\to 0$, $F(J_{k})$ approaches a
solution to a nonlinear differential inequality,
\begin{equation*}
  \myf{dF}{dt} \leq K_{3}\tq e^{-\Phi t}\p{F+\eps F^{2}}.
\end{equation*}
Without the coefficient $e^{-\Phi 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) \leq
Z$.
Then, by (\ref{3.22}), we have
\begin{equation*}
F(J_k) \leq F(\tilde{J_{k}})\p{1 + K_3 \tq
e^{-\Phi k h}\p{1 + \eps Z}h}.
\end{equation*}
Let us abbreviate $e^{-\Phi k h}$ by $d^k$.  Then we have 
\begin{equation*}
F(J_k) \leq F(J_0)\prod_{j=0}^k \p{1+ K_3 h \tq
d^j\p{1+\eps Z}}.
\end{equation*}
Taking logarithms,
\begin{align*}
\ln\p{F(J_k)/F(J_0)}
&\leq\sum_{j=0}^k\ln\p{1+K_3 h\tq d^j\p{1+\eps Z}}\\
&\leq \sum_{j=0}^k K_3 h \tq d^j\p{1 + \eps Z}\\
  &\leq K_3 h\tq \p{1 + \eps Z}\myf{1}{(1-d)}.\\
\end{align*}
Thus we conclude
\begin{equation}\label{eq:3.25}
F(J_k) \leq F(J_0)\exp\p{\myf{K_3h\tq\p{1+\eps Z}}{1-e^{-\Phi h}}}.
\end{equation}
The function $f(h) = h/\p{1-e^{-\Phi h}}$ is increasing for $h>0$ and
tends to $1/\Phi$ as $h\to 0$.  Thus, for $h$ sufficiently small, we
obtain
\begin{equation}\label{eq:3.25.1}
       F(J_{k}) \leq F(J_0)\exp\p{\myf{2 K_3\tq\p{1+ \eps Z}}{\Phi}}.
\end{equation}
Let $K_4=2K_3$.  The estimate \eqref{eq:3.25.1} is valid as long as
$F(J_k)\leq Z$.  The condition required for this result is that
\begin{equation}\label{eq:3.22.1}
F(J_0)\leq \exp\p{-\myf{K_4 \tq\p{1 + \eps Z}}
{\Phi}} Z =g(Z).
\end{equation}
The value of $Z$ which maximizes $g(Z)$
is $Z = \myf{\Phi}{K_4\tq \eps }$. Thus, our least restrictive
condition on $F(J_0)$ is
\begin{equation}\label{eq:3.22.2}
F(J_0) \leq 
\exp\p{-1-\myf{K_4\tq}{\Phi}}\myf{\Phi}{K_4\eps \tq}
=\exp\p{-1-\myf{K_4 q\supnorm{Y_{0}}}{\Phi}}\myf{\Phi}{K_4\eps
q\supnorm{Y_0}}.
\end{equation}

We summarize these estimates with the following lemma.
\begin{lem}\label{lem:3.8.0} If $F(J_0)$ satisfies {\rm\eqref{eq:3.22.1}},
then, for all $k \geq 1$,
$F(J_k) \leq Z$.  In particular, if
$F(J_0)$ satisfies {\rm\eqref{eq:3.22.2}}, then, for all $k \geq 1$,
\begin{equation}
F(J_k) \leq Z = \myf{\Phi}{K_4 \eps q\supnorm{Y_0}}.
\end{equation}
Furthermore, if $F(J_{0})$ satisfies \eqref{eq:3.22.2}, then by
\eqref{eq:3.25},
\begin{equation}\label{eq:3.29}
  F(J_{k}) \leq F(J_{0})\exp\p{\myf{K_{4}q\supnorm{Y_{0}}}{\Phi}+1}.
\end{equation}
\end{lem}

Next, we estimate the amount that the solution ``drifts'' from its
original base point. We use $\tilde{w} =\lim_{x\to-\infty}w_{0}(x)$
as our base point.  Let $\D_{\eps} = \D$ for $0\leq \eps\leq\eps_{0}$, 
and $\D_{\eps}= B_{C_{1}}(\tilde{w})$ for $\eps_{0}< \eps \leq 1$.

\begin{lem}
Let $U^{h}(x,t)$ be an approximate solution to \eqref{2.2},
constructed with the fractional--step Glimm scheme, with
initial data $U_{0} \in BV(\R)$.  If $U^{h}(J) \in BV$ for every
mesh curve $J$, then
\begin{equation*}
U_{\pm\infty}=\lim_{x\to\pm\infty}U(x,t)
\end{equation*}
exists for all $t$ and satisfies
\begin{equation*}
\myf{dU_{\pm\infty}}{dt} = G\bp{U_{\pm\infty}}.
\end{equation*}
\end{lem}

For our application \eqref{1.4}, by Lemma \ref{lem:3.8},
$G(U) = Y\phi(T)\bp{\myf{-q}{c_{v}T},\myf{q}{c_{v}T},\myf{q}{T},-1}^T$.
Let $\bar{\phi} = \sup_{\D_{\eps}}\phi(T)$.

\begin{lem}\label{lem:3.13}
In the $(r,s,S)$-coordinates,
\begin{equation*}
    \begin{gathered}
    \n{r_{-\infty}(t)-r_{-\infty}(0)}+\n{s_{-\infty}(t)-s_{-\infty}(0)} \leq
        \int_{0}^{Y_{-\infty}(0)}\myf{2q}{c_{v}T}\ dY,\\
    \n{S_{-\infty}(t)-S_{-\infty}(0)} \leq
\int_{0}^{Y_{-\infty}(0)}\myf{q}{T}\ dY.\\
\end{gathered}
\end{equation*}
As $\eps \to 0$, $c_{v}\to \infty$.  Recall that, by Lemma
{\rm \ref{lem:3.5a}}, $T\geq T'(\eps)$.  Let $C_{4}$ be an upper bound for
$1/c_{v}$ on $\D_{\eps}$. Denote $Y_{-\infty}(0)$ by $y_{0}$. Then
\begin{equation*}
    \begin{gathered}
    \n{r_{-\infty}(t)-r_{-\infty}(0)}+\n{s_{-\infty}(t)-s_{-\infty}(0)} \leq
        C_{4}\myf{2y_{0}q}{T'},\\
         \n{S_{-\infty}(t)-S_{-\infty}(0)} \leq \myf{y_{0}q}{T'}.\\
\end{gathered}
\end{equation*}
\end{lem}

Applying Lemma \ref{lem:3.13}, we have
\begin{lem}\label{lem:3.14}
Let $U^{h}(x,t)$ be an approximate solution to \eqref{1.4}
constructed with the fractional--step Glimm scheme, with
initial data $U_{0} \in BV(\R)$ and $w(U_{0})(-\infty)= \tilde{w}$.
Then, for all $(x,t)\in J$,
\begin{equation*}
\begin{aligned}
 \onenorm{(r,s)^{h}(x,t) -(r,s)(\tilde{w})} &\leq
  2C_4 \myf{y_{0}q}{T'} + \Var_{rs}\bp{U^h(J_k)},\\
 \onenorm{(r,s,S)^{h}(x,t) -(r,s,S)(\tilde{w})}
&\leq \p{2C_4 +1}\myf{y_{0}q}{T'} + \Var_{rsS}\bp{U^h(J_k)}.
\end{aligned}
\end{equation*}
\end{lem}
\begin{comment}
For $0\leq \eps < \eps_0$ we need $ \onenorm{(r,s)^{h}(x,t)
-(r,s)(\tilde{w})} \leq 120 K_0 V$ to ensure that $(r,s)$ remains
within $\D_{rs}$, and we need $\onenorm{(r,s,S)^{h}(x,t)
-(r,s,S)(\tilde{w})}\leq 6K K_0 V$ to ensure that $(r,s,S)$ remains
within $\D$.  Furthermore, for $\eps_0 \leq \eps \leq 1$, we need
$$
\onenorm{(r,s,S)^{h}(x,t) -(r,s,S)(\tilde{w})} \leq C_1
$$
to ensure that $(r,s,S)$ remains within $B_{C_1}(\tilde{w})$.

Applying Lemmas \ref{lem:3.8.0}, \ref{lem:3.14}, \ref{lemtemp:4.1} and
\ref{lemtemp:4.2}, which implies $\Var_{rsSY}(J) \leq K F(J)$, we have:

\begin{lem}\label{lem:3.15}
Let $U_{0}(x)$ be given with $w(U_{0})(-\infty)=\tilde{w}$, and let
\begin{equation*}
\eps F(J_0) \leq 
 \exp\p{-1-\myf{K_4 q\supnorm{Y_{0}}}{\Phi}}\myf{\Phi}{K_4 q\supnorm{Y_0}}.
\end{equation*}
Then \eqref{eq:3.29} holds and, for $(x,t)$ on any mesh curve $J$,
\begin{equation*}
\begin{split}
&\onenorm{(r,s,S)^{h}(x,t) -(r,s,S)(\tilde{w})}\\
&\qquad\qquad\qquad    \leq  \p{2C_{4}+1}\myf{y_{0}q}{T'} +
          K F(J_{0})\exp\p{\myf{K_{4}q\supnorm{Y_{0}}}{\Phi}+1}.
\end{split}
\end{equation*}
\end{lem}
\end{comment}

\begin{thm}\label{thm:3.1}
Let $E$ be any compact set in $rsS$-space, let $\eps\in\br{0,1}$, and
let $N>1$ be any positive constant.  Let $\D$, $C_{1}$, and
$K_0$ be as determined in Lemma {\rm \ref{thmtemp:4.1}}. 
Let $C_0(E,N)$ and $C(E,N)$  be given by 
{\rm (\ref{eq:C3})} and {\rm (\ref{eq:C(E,N)})}.  
Let $K= 20/M_{0}^{2}$. Let $K_1$ and $\eps_0$ be as determined in Lemma
{\rm \ref{lem:3.13a}}.  Let $w_{0}(x)$ be given initial data for
\eqref{1.4}, expressed in $(r,s,S,Y)$-coordinates, and with
$(r,s,S,Y)$ values in $E\times\br{0,1}$.  Let $V_0=\Var_{rsSY}(w_{0})$.

If $0 \leq \eps < \eps_0$, and
\begin{align}
&  V_0 \leq
   \exp\p{-\myf{K_{4}q\supnorm{Y_{0}}}{\Phi}-1}
   \min\p{ \myf{C_0(E,N)}{\eps}, \myf{\Phi}{5\eps K_0  K_4 q\supnorm{Y_0}}},
   \label{eq:3.32} \\
&   V_0 \leq \exp\p{-\myf{K_{4}q\supnorm{Y_{0}}}{\Phi}-1}
\p{ N - \p{2C_{4}+1}\myf{y_{0}q}{5 K K_0 T'}},
   \label{eq:3.32b}\\ 
& V_0 \leq  \exp\p{-\myf{K_{4}q\supnorm{Y_{0}}}{\Phi}-1} 
%%
%Change
%        \p{ N - 2C_{4}\myf{y_{0}q}{100 K_{0}T'}} ,
%To
%%
\p{ N - C_{4}\myf{y_{0}q}{50 K_{0}T'}} ,
   \label{eq:3.32c}
\end{align}
then the fractional--step Glimm scheme constructs approximate solutions
$U^h(x,t)$ for which $w(U^h)(x,t)\in \D$ and with uniformly bounded spatial
total variation.

If $\eps_0 \leq \eps \leq 1$, and
\begin{align}
 &  V_0 \leq \exp\p{-\myf{K_{4}q\supnorm{Y_{0}}}{\Phi}-1}
    \min\p{\myf{6K_0 N}{K_1},\myf{C_3}{\eps}, \myf{\Phi}{\eps K_1 K_4 q\supnorm{Y_0}}},
    \label{eq:3.33}\\
   &    V_0 \leq \exp\p{-\myf{K_{4}q\supnorm{Y_{0}}}{\Phi}-1}\myf{1}{KK_1}
\p{C_1-\p{2C_{4}+1}\myf{y_{0}q}{T'}},
     \label{eq:3.33a}
\end{align}
then the fractional--step Glimm scheme constructs approximate solutions
$U^h(x,t)$ for which $w(U^h)(x,t)\in B_{C_1}(\tilde{w})$
and with uniformly bounded spatial total variation.
\end{thm}

\begin{proof}
    Our functionals $F$ and $F_{0}$ are nearly the same as Temple's.
    Our modifications consist of adding the terms $M_{4}\n{\zeta}$ and
    $\alpha_{\infty}+ \beta_{\infty}+\mu_{\infty}+\eta_{\infty}
    +\n{\delta_{\infty}}+M_{4}\n{\zeta_{\infty}}-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
    \ref{thmtemp:4.1} apply to our non-reacting step, as long as the
    conditions for that proof remain valid.  Those conditions are:
\begin{equation}
    \Var_{rsS}(J) \leq 6KK_0 N,\label{eq:3.44}\qquad\qquad
    \Var_{rs}(J) \leq 120 K_0 N,
\end{equation}
for $0 < \eps \leq \eps_0$, to ensure $(r,s,S)(x,t) \in \D$, and
\begin{equation}\label{eq:3.45}
    \Var_{rsS}(J) \leq C_1,
\end{equation}
for $\eps_0 < \eps \leq 1$, to ensure $(r,s,S)(x,t)\in B_{C_1}(w)$, and finally
\begin{equation}
\label{eq:3.46}
\eps F(J_0)  < \myf{1}{4 M_3 K}\min\p{\myf{1}{3 G -1},1-C_0},\qquad F(J) \leq 6 K_0 N,
\end{equation}
for $0< \eps \leq 1$ to ensure that $F(J)$ does not increase in a
non-reacting step.

By Lemmas \ref{lemtemp:4.1} and \ref{lemtemp:4.2} (and their
generalizations to our problem), we have
\begin{equation*}
 \begin{alignedat}{2} &\Var_{rsSY}(J) \leq K L(J) \leq K F(J),&\qquad
    &0\leq \eps\leq 1 ,\\ &\Var_{rs}(J) \leq 20 L_{0}(J) \leq 20
    F_{0}(J),&\qquad &0 \leq \eps \leq \eps_{0}.
\end{alignedat}
\end{equation*} 
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 \ref{lem:3.13a}, \ref{lem:3.8.0}, and
\ref{lem:3.14}, \eqref{eq:3.44} becomes \eqref{eq:3.32b} and
\eqref{eq:3.32c}, and \eqref{eq:3.45} becomes \eqref{eq:3.33a}.  For
$J_0$, the bounds \eqref{eq:3.46} are implied by Lemma \ref{lem:3.13a}
if $V_0\leq N$ and either $\eps V_0 \leq C_0(E,N)$ for $0<\eps \leq
\eps_0$, or $\eps V_0 \leq C(E,N)$ for $\eps_0< \eps \leq 1$. For
$J_k$, the bounds \eqref{eq:3.46} are implied for $0<\eps\leq \eps_0$
by \eqref{eq:3.32}. For $\eps_0<\eps \leq 1$, \eqref{eq:3.33} implies
$\eps V_0\leq C_3$ and $K_1 V_0 \leq 6K_0N$, or
\begin{equation*}
V_0 \leq \myf{6K_0 N}{8 K_0^2 M_3 N} \leq N . 
\end{equation*}
For $J> J_0$, using Lemma
\ref{lem:3.8.0}, we find that, if
\begin{equation}\label{eq:3.34}
  Z_{0} =  \min\p{6K_0 N,K_1 \myf{C_3}{\eps}},
\end{equation}
and
\begin{equation}\label{eq:3.34a}
   F\bp{J_{0}} \leq \exp\p{-\myf{K_{4}q\supnorm{Y_{0}}}{\Phi}-1}Z_{0},
\end{equation}
then, for all $k$, $F\bp{J_{k}} \leq Z_{0}$ so that 
the bounds \eqref{eq:3.46} are satisfied.  When we express 
(\ref{eq:3.34} and \ref{eq:3.34a}) in terms of $V_0$, using 
Lemma \ref{lem:3.13a}, we obtain \eqref{eq:3.32} and \eqref{eq:3.33}.
\begin{comment}    
Note that in the non-reacting step, the $Y$ equation is linear and
decoupled from the rest of the system.  Thus we do not need to show
that the total variation of $Y$ is small.  We need only show that
the total variation of $Y$ remains bounded---this is established by
any bound on $F$.
\end{comment}
\end{proof}


\section{Convergence to the Entropy Solution}\label{sect:5}
A scheme for constructing approximate solutions to a differential
equation is called {\em 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
\cite{hyp:LusTem} it is shown that the fractional--step scheme as
constructed in Section \ref{sect:2.2}, by using Euler's method to
approximate solutions to \eqref{ivp}, is consistent with the notion of
weak solutions for the nonlinear water-hammer problem.  It is also
shown in \cite{hyp:DH} that fractional--step schemes based on Glimm's
method and the first order approximation to \eqref{ivp} 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 \eqref{1.3}, as constructed in Section \ref{sect:2.2},
using any consistent approximation to \eqref{ivp},
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^{\infty}$ and converges
pointwise a.e. to a function $U(x,t)$, then $U(x,t)$ is an
entropy solution of the 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 (\ref{1.2}) and (\ref{1.4}).

Let $d\chi_k$ denote the uniform probability measure on
$\br{-1,1}$,
and let $d\chi$ denote the induced 
product probability measure for the random sample
$\{\chi_k\}_{k = 1}^{\infty}$ in the Cartesian product space
${\mathcal A}=\Pi_{k=1}^\infty\br{-1,1}$.

\begin{thm}\label{thm:4.1}
Suppose that
\begin{enumerate}
\item The sequence $U^{h_\ell}(x,t)$ is constructed by using the
Glimm fractional--step scheme, as described in
Section {\rm \ref{sect:2.2}}, with a fractional--step operator which is
consistent in the sense of Definition {\rm \ref{def:2.1}}, 
and with the random sample 
$\mybrace{\chi_k}_{k=1}^{\infty}$ chosen from $\mathcal{A}$.
\item The sequence $U^{h_\ell}(x,t)$ is uniformly bounded in
$L^{\infty}$ and converges pointwise a.e.  to the function $U(x,t)$.
\end{enumerate}
Then there exists a null set $\mathcal{N}\subset \mathcal{A}$ such that,
for $\mybrace{\chi_k}\in\mathcal{A}-\mathcal{N}$, 
the function $U(x,t)$ is an entropy solution of the Cauchy 
problem {\rm (\ref{1.3})} and {\rm (\ref{2.1})}. 
That is, for any convex entropy
pair $(\eta,q)$ with respect to $U$, the following inequality
\begin{equation}
\eta(U)_t +q(U)_x \le \grad\eta(U)G(U)
\label{eq:4.2}
\end{equation}
%%
%Change
%holds in the sense of distributions:
%To
%%
holds in the sense of distributions, that is,
\begin{equation}
\int_0^\infty\int_{-\infty}^\infty
\left(\eta(U)\phi_t + q(U)\phi_x+\grad\eta(U)G(U)\phi\, \right) dxdt
+\int_{-\infty}^{\infty}\eta\bp{U_{0}(x)}\phi(x,0)\ dx \ge 0,
\label{eq:4.3}
\end{equation}
where $\phi \in C_0^\infty((-\infty,\infty)\times[0,\infty))$
and $\phi(x,t)\ge 0$.
\end{thm}

\subsection*{Proof}
For any convex entropy pair $(\eta, q)$ with respect to $U$,
we define
\begin{equation}
\begin{split}
\mathcal{J}(\chi,h,\phi)
=- &\int_{0}^{\infty}\int_{-\infty}^{\infty}\left(\eta\bp{U^h}\phi_t +
q\bp{U^h}\phi_x +
                \grad \eta\bp{U^h}G\bp{U^h}\phi\right)dxdt,\\
   &-\int_{-\infty}^{\infty}\eta\bp{U_{0}(x)}\phi(x,0)\ dx,
\end{split}
\label{eq:4.4}
\end{equation}
where $\phi \in C_0^\infty ((-\infty,\infty)\times[0,\infty))$
and $\phi(x,t) \ge 0$.

Note that for $k\geq 0$,
\begin{equation}
\begin{split}
-\int_{k h}^{(k+1)h}\int_{-\infty}^{\infty} &\eta\bp{U^h}\phi_t\ dxdt
\\
= & -\int_{k h}^{(k+1)h}\int_{-\infty}^{\infty} \eta\bp{U_0^h}\phi_t\
dxdt
  -\int_{k h}^{(k+1) h}\int_{-\infty}^{\infty}
   \left(\eta\bp{U^h}
      -\eta\bp{U_0^h}\right)\phi_t\ dxdt.
\end{split}
\label{eq:4.6}
\end{equation}
Since $U_0^h$ is an entropy solution of the
conservation laws (\ref{eq:2.19}),
\begin{equation}
\begin{split}
&-\int_{k h}^{(k+1)h}\int_{-\infty}^\infty \eta\bp{U_0^h}\phi_t\ dxdt\\
\le &\int_{k h}^{(k+1)h}\int_{-\infty}^\infty q\bp{U_0^h}\phi_x\
dxdt\\&+
\int_{-\infty}^\infty\p{\eta\bp{U_{0}^h(x,k h+0)}\phi(x,k h)
              -\eta\bp{U_0^h(x,(k+1)h-)}\phi(x,(k+1)h)}\ dx.
\end{split}
\label{eq:4.7}
\end{equation}

Since $\F$ is consistent in the sense of Definition \ref{def:2.1},
there exists a function $\nu(s)$ such that
$\nu(s) \to 0$ as $s \to 0$, and
\begin{equation}
\N{\F\bp{s,U_{0}^{h}(x,t)}-U_{0}^{h}(x,t)-G\bp{U_{0}^{h}(x,t)}s} \le
\nu(s)\N{G\bp{U_{0}(x,t)}}\n{s}.
\label{eq:4.5}
\end{equation}



Equation (\ref{eq:4.5}) states that, as a function of $s$,
$\F\bp{s,U_{0}^{h}(x,t)}$
is differentiable at $s=0$, and has the partial derivative
$G\bp{U_{0}^{h}(x,t)}$.  Since $\eta$ is a differentiable
function of $U$, the chain rule states that, as a function of $s$,
$\eta\p{\F\bp{s,U_{0}^{h}(x,t)}}$ is differentiable
at $s=0$, and has partial derivative
$\grad\eta\bp{U_{0}^{h}(x,t)}G\bp{U_{0}^{h}(x,t)}$.  Therefore, there
is a function $e(s;x,t)$, which converges uniformly to $0$ as
$s\rightarrow 0$, such that
\begin{equation}
\eta\p{\F\bp{s,U_{0}^{h}(x,t)}} - \eta\bp{U_{0}^{h}(x,t)}
-\grad\eta\bp{U_{0}^{h}(x,t)}G\bp{U_{0}^{h}(x,t)}s 
=\eps(s;x,t)s.
\end{equation}
Thus, for $k h \leq t \leq (k+1)h$,
\begin{equation}
\begin{split}
\eta\bp{U^{h}(x,t)} =& \eta\bp{U_{0}^{h}(x,t)} \\
&+ \grad\eta\bp{U_{0}^{h}(x,t)}G\bp{U_{0}^{h}(x,t)}(t-k h)
+ \eps(t-k h;x,t)(t-k h).
\end{split}
\end{equation}

We can now compute
\begin{equation}\label{eq:4.8a}
\begin{split}
-&\int_{k h}^{(k+1)h}\int_{-\infty}^{\infty}
\left(\eta\bp{U^h}-\eta\bp{U_0^h}\right)\phi_t\ dxdt =
-\int_{k h}^{(k+1)h}\int_{-\infty}^{\infty}
\eps(t-k h; x,t)(t-k h)\phi_{t}\ dxdt\\
&\mspace{100mu}+\int_{-\infty}^{\infty}\int_{k h}^{(k+1) h}
\pd{}{t}\left(\grad\eta\bp{U_{0}^{h}(x,t)}G\bp{U_{0}^{h}(x,t)}(t-k h)\right)
\phi\ dxdt.
\end{split}
\end{equation}

Thus,
\begin{equation}\label{eq:4.8b}
\begin{split}
-\int_{k h}^{(k+1)h}&\int_{-\infty}^{\infty}
\p{\eta\bp{U^h}-\eta\bp{U_0^h}}\phi_t\ dxdt\\
=& -\int_{k h}^{(k+1)h}\int_{-\infty}^{\infty}
   \eps(t-k h; x,t)(t-k h)\phi_{t}\ dxdt\\
&+\int_{-\infty}^{\infty}\int_{k h}^{(k+1) h}
\left[\grad\eta\bp{U_{0}^{h}(x,t)}G\bp{U_{0}^{h}(x,t)}\right.\\
&\mspace{120mu} +\left.\pd{}{t}\p{\grad\eta\bp{U_{0}^{h}(x,t)}
G\bp{U_{0}^{h}(x,t)}}(t-k h)\right]\phi\ dxdt.
\end{split}
\end{equation}

Plugging (\ref{eq:4.7}) and (\ref{eq:4.8b}) into (\ref{eq:4.6}), we
obtain
\begin{equation*}
\begin{split}
-\int_{k h}^{(k+1)h}\int_{-\infty}^\infty &\eta\bp{U^h}\phi_t\ dxdt \\
\leq& \int_{k h}^{(k+1)h}\int_{-\infty}^\infty
 \br{q\bp{U_0^h}\phi_x+\grad\eta\bp{U_{0}^h}G\bp{U_0^h}\phi}\ dxdt\\
%%
%Change
%&+\int_{-\infty}^{\infty}\eta\bp{U_{0}^h(x,k h+0)}\phi(x,k h)
%-\eta\bp{U_{0}^h(x,(k+1)h-)} \phi(x,(k+1)h)\ dx\\
%To
%%
&+\int_{-\infty}^{\infty}\br{\eta\bp{U_{0}^h(x,k h+0)}\phi(x,k h)
  -\eta\bp{U_{0}^h(x,(k+1)h-)} \phi(x,(k+1)h)}\ dx\\
&+\int_{-\infty}^\infty\int_{k h}^{(k+1)h}
\pd{}{t}\br{\grad\eta\bp{U_{0}^{h}(x,t)}
G\bp{U_{0}^{h}(x,t)}}(t-k h)\phi(x,t)\ dxdt\\
&-\int_{k h}^{(k+1)h}\int_{-\infty}^{\infty}
 \eps(t-k h; x,t)(t-k h)\phi_{t}\ dxdt.
\end{split}
\label{eq:4.9}
\end{equation*}

Summing over $k$, we have
\begin{equation*}
\mathcal{J}(\chi,h,\phi)
\le \mathcal{I}(\chi,h,\phi)
  +\sum_{k=0}^\infty\mathcal{R}_j(\chi,h,\phi)
+\sum_{k=0}^\infty\mathcal{D}_j(\chi,h,\phi),
\label{eq:4.10}
\end{equation*}
where
\begin{align*}
\ &\mathcal{I}(\chi,h,\phi)
   = \sum_{k=0}^{\infty}\mathcal{I}_{k}(\chi,h,\phi),\\
\ &\mathcal{I}_{0}(\chi,h,\phi)
=\int_{-\infty}^\infty\p{\eta\bp{U^{h}(x,0)}-\eta\bp{U_{0}(x)}}
   \phi(x,k h)\ dx,\\
\ &\mathcal{I}_{k}(\chi,h,\phi)
=\int_{-\infty}^\infty\p{\eta\bp{U_{0}^h(x,k h+0)}
   -\eta\bp{U_{0}^h(x,k h-)}
  } \phi(x,k h)\ dx,\ k = 1,\cdots,
\end{align*}
\begin{align*}
\mathcal{R}_{k}(\chi,h,\phi)
=&\int_{k h}^{(k+1)h}\int_{-\infty}^\infty
        \p{q\bp{U_0^h}-q\bp{U^h}}\phi_x\ dxdt\\
&+\int_{k h}^{(k+1)h}\int_{-\infty}^\infty
\p{\grad\eta\bp{U_{0}^h}G\bp{U_0^h}-\grad\eta\bp{U^{h}}
                           G\bp{U^h}}\phi \, dxdt\\
&-\int_{k h}^{(k+1) h}\int_{-\infty}^{\infty}
  \eps(t-k h; x,t)(t-k h)\phi_{t}\ dxdt,
\end{align*}
and
\begin{equation*}
\mathcal{D}_{k}(\chi,h,\phi)
=\int_{k h}^{(k+1) h}\int_{-\infty}^\infty
 \pd{}{t}\Bp{\grad\eta\bp{U_{0}^{h}(x,t)}
G\bp{U_{0}^{h}(x,t)}}(t-k h)\phi(x,t)\ dxdt.
\end{equation*}
We now analyze each of these components for convergence to zero as
$h$ tends to zero.

\begin{lem}\label{lem:4.0}
There is a null subset ${\mathcal N}\subset {\mathcal A}$ and a
positive sequence $h_{m}$, which converges to zero, such
that, for $\chi \in {\mathcal A}-{\mathcal N}$, the functionals
$\mathcal{I}(\chi,h_m,\phi)$ converge weakly
to zero as measures on compact subsets of $\R\times\R^{+}$.
\end{lem}
\begin{proof}
Note that $\mathcal{I}_{0}\rightarrow 0$ since $U^{h}(x,0)$ was
chosen to converge to $U_{0}(x)$ in $L^{1}(\R)$.
Since
$$
\mathcal{I}_{k}(\chi,h,\phi)
=\sum_{j=-N}^{N}\int_{(j-1)l}^{(j+1)l}
 \p{\eta\bp{U_{0}^h((j+\chi_{k})l, k h-)}
  -\eta\bp{U_{0}^h(x,k h-)}}\phi(x,k h))dx,
$$
and
\begin{equation}
\n{\mathcal{I}_{k}(\chi,h,\phi)}
\le C TV\bp{U_{0}^h(x, k h-)}h\supnorm{\phi},
\label{eq:4.13}
\end{equation}
we have 
\begin{equation}
\n{\mathcal{I}(\chi,h,\phi)}
\le  C L d \supnorm{\phi},
\label{eq:4.14}
\end{equation}
where $L$ is an upper bound for $TV\bp{U_{0}^h(\cdot,(k-1)h)}$, and
$d$ is the diameter of $\supp \phi(x,t)$.  Thus $\phi \to
\mathcal{I}(\chi,h,\phi)$ is a distribution of order zero.  By the
Riesz representation theorem, $\mathcal{I}(\chi,h,\phi)$ corresponds
to a Radon measure \cite{an:Don69}.

Let $t_{k}$ denote $k h-$.  Notice that
\begin{equation}\label{eq:4.15}
\begin{split}
&\n{\int_{-1}^1\mathcal{I}_{k}(\chi,h,\phi)\ d\chi_{k}}\\
&=\n{\sum_{j=-N}^{N}\int_{-1}^{1}\int_{(j-1)l}^{(j+1)l}
 \p{\eta\bp{U_{0}^h((j+\chi_{k})l,t_{k})}
  -\eta\bp{U_{0}^h(x,t_{k})}}\phi(x,t_{k})\ dx d\chi_{k}}\\
&\le \n{\sum_{j=-N}^{N}\phi_j^k\int_{-1}^{1}\int_{(j-1)l}^{(j+1)l}
 \p{\eta\bp{U_{0}^h((j+\chi_{k})l,t_{k})}
  -\eta\bp{U_{0}^h(x,t_{k})}}\ dx d\chi_{k}}\\
&\quad +\supnorm{\grad\eta}\supnorm{\phi_x}l
   \sum_{j=-N}^{N}\int_{-1}^{1}\int_{(j-1)l}^{(j+1)l}
  \n{U_{0}^h((j+\chi_{k})l,t_{k})
   -U_{0}^h(x,t_{k})}\ dx d\chi_{k}\\
&\le C\supnorm{\grad\eta}\supnorm{\phi_x}l
\int_{-1}^{1}
  \int_{-l}^{l}\sum_{j=-N}^{N}\n{U_{0}^h((j+\chi_{k})l,t_{k})
                          -U_{0}^h(jl +x,t_{k})}\ dx d\chi_{k}\\
&\le  C\supnorm{\grad\eta}\supnorm{\phi_x}TV\bp{U_0^h(\cdot,t_{k})}h^2,
\end{split}
\end{equation}
where
$\phi_j^k=\phi(jl, t_{k})$.
\begin{comment}
\begin{equation}
\n{\int_{-1}^1\mathcal{I}_{k+1}(\chi,h,\phi)d\chi_{k+1}}
\le Ch^2\supnorm{\eta\bp{U^h}}\supnorm{\phi_x}
\le Ch^2\supnorm{\phi_x}.
\label{eq:4.15}
\end{equation}
\end{comment}
Therefore, for $0<k_1<k_2,$ we have
\begin{equation}
\aligned
\n{\int 
\mathcal{I}_{k_1}(\chi,h,\phi)\mathcal{I}_{k_2}(\chi,h,\phi)d\chi}
&=\n{\int\left(\int\mathcal{I}_{k_2}(\chi,h,\phi)d\chi_{k_2}\right)
\mathcal{I}_{k_1}(\chi,h,\phi)d\bar{\chi}_{k_2}}\\
&\le C h^3\supnorm{\phi}\supnorm{\phi_x},
\endaligned
\label{eq:4.16}
\end{equation}
where  $\bar{\chi}_{k_2}$ denotes the sequence obtained from
$\chi$ by
deleting the $k_2$-th element. Moreover, we have from (\ref{eq:4.13})
that
\begin{equation}
\int (\mathcal{I}_k(\chi,h,\phi))^2d\chi \le Ch^2\supnorm{\phi}^2.
\label{eq:4.17}
\end{equation}
Thus,
\begin{equation}
\aligned
\int(\mathcal{I}(\chi,h,\phi))^2d\chi
&=2 \sum\limits_{k_1<k_2}\int\mathcal{I}_{k_1}\mathcal{I}_{k_2}d\chi
 +\sum\limits_k\int\mathcal{I}_k^2d\chi \\
&\le C h \supnorm{\phi}\supnorm{\phi_x} + C h\supnorm{\phi}^2,
\endaligned
\label{eq:4.18}
\end{equation}
so that $\mathcal{I}(\chi,h,\phi)\to 0$ in $L^2(\mathcal{A},d\chi)$.  As a consequence, 
there exists a sequence $h_m\to 0$ such that $\mathcal{I}(\chi,h_m,\phi)\to 0$ pointwise a.e. in 
$(\mathcal{A},d\chi)$. In other words,
 for any $\phi\in C_0^\infty((-\infty,\infty)\times
[0,\infty))$,
there exists a sequence $h_m \to 0$ and a null subset ${\mathcal N}\subset {\mathcal A}$
such that, if $\chi \in {\mathcal A}-{\mathcal N}$, then
\begin{equation}
\mathcal{I}(\chi,h_m,\phi) \to 0, \quad\text{as } h_m \to 0.
\label{eq:4.19}
\end{equation}
Let $\{\phi_\ell\}_{\ell\in \Lambda}$ be a countable dense set in
$C_0^\infty((-\infty,\infty)\times [0,\infty))$ with respect to
the norm $\supnorm{\ }$.
Then, after refining $h_m$ above by a diagonal
argument and after taking the appropriate countable union of
null sets ${\mathcal N}_\ell$,
we can conclude that there exists a sequence
$h_m\to 0$ and a null subset ${\mathcal N}\subset {\mathcal A}$
such that, if $\chi\in {\mathcal A}-{\mathcal N}, \ell\in \Lambda$, then
\begin{equation}
\mathcal{I}(\chi,h_m,\phi_\ell)\to 0, \quad\text{as } m \to \infty.
\label{eq:4.20}
\end{equation}
Let $K$ be a compact subset of $\R \times {\R}^+$.
We see from (\ref{eq:4.14}) that
$\mathcal{I}(\chi,h_m,\phi)$ is a uniformly
bounded (i.e., Lipschitz continuous) sequence of linear
functionals on $C_0\br{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\br{K}\cap C^1\br{K}$.
This means that the functionals $\mathcal{I}(\chi,h_m,\phi)$
converge weakly
to zero as measures on $K$.
\end{proof}

Since, for $(k-1)h \leq t < k h$, $\n{U^{h}(x,t)-U_{0}^{h}(x,t)}
\le C (t-(k-1)h)\leq Ch$, we have
\begin{equation}
\sum_{k=1}^{\infty}\mathcal{R}_{k}(\chi,h,\phi) \le C
h\p{\onenorm{\phi_x}+
\onenorm{\phi} + \onenorm{\phi_{t}}}.
\label{eq:4.22}
\end{equation}

We now estimate 
\begin{equation*}
\mathcal{D}_{k}(\chi,h,\phi)
=\int_{k h}^{(k+1) h}\int_{-\infty}^\infty
 \pd{}{t}\Bp{\grad\eta\bp{U_{0}^{h}(x,t)}
G\bp{U_{0}^{h}(x,t)}}(t-k h)\phi(x,t)\ dxdt.
\end{equation*}
Let $H\bp{U_{0}^{h}}$ denote the derivative of $\grad\eta\cdot G$,
evaluated at $U_{0}^{h}$.

Then
\begin{equation*}
\aligned
\pd{}{t}\p{\grad\eta\bp{U_{0}^{h}}G\bp{U_{0}^{h}}}(x,t)
&=\widehat{H}\bp{U_{0}^{h}(x,t)}\cdot {U_{0}^{h}}_{t}(x,t)
=-\widehat{H}\bp{U_{0}^{h}(x,t)}\cdot F\bp{U_{0}^{h}(x,t)}_{x}\\
&=-\widehat{H}\bp{U_{0}^{h}(x,t)}\cdot
 \widehat{\grad F}(U_0^h)\bp{U_{0}^{h}(x,t)}_{x},
\endaligned
\end{equation*}
by the chain rule of differentiation.
Here $\widehat{H}(U^h(x,t))$ denotes Vol'pert's ``functional
superposition'' \cite{hyp:Vol67}, which is given by
\begin{equation}
\widehat{H}(U^h(x,t)) =
\int_0^1 H\p{\theta U^h(x,t-)+(1-\theta)U^h(x,t+0)}\ d\theta,
\end{equation}
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\in BV\br{\R\times{\R}^+;{\R}^n}$ and
$g \in C\br{{\R}^n;\R}$, $\widehat{H}\circ f$ is
well-defined almost everywhere with respect to one--dimensional
Hausdorff measure.

Thus, if $C_{1}$ is an 
$L_{\infty}$ matrix bound for
$\widehat{H}(U)\widehat{\grad F}(U)$,
\begin{equation}\label{eq:4.23}
\mathcal{D}_{k}(\chi,h,\phi) \leq C_{1}
\tv\bp{U_{0}^{h}}\supnorm{\phi}h^{2}.
\end{equation}
Let $C_{2}$ be the diameter of $\left\{t\mid \exists x \text{ such
that } 
\phi(x,t)\neq 0 \right\}$. Then
\begin{equation}
\sum_{k=1}^{\infty}\mathcal{D}_{k}(\chi,h,\phi) \leq
C_{1}C_{2}\supnorm{\phi} h.
\label{eq:4.24}
\end{equation}
\begin{comment}
Thus it follows
from (\ref{eq:4.22})-(\ref{eq:4.24}) that
\begin{equation}
\n{\sum_{k=1}^\infty \mathcal{R}_k(\chi,h,\phi)}
\leq CD h\p{\supnorm{\phi}+\supnorm{\phi_x}}.
\label{eq:4.25}
\end{equation}
\end{comment}
Therefore, we can conclude from (\ref{eq:4.22}), (\ref{eq:4.24}), and
Lemma \ref{lem:4.0}
that, if $\chi\in A$ and
$\phi\in C_0\br{K}\cup C^1\br{K}$, then
\begin{equation*}
\mathcal{J}(\chi, h_m,\phi)\to 0, \qquad\text{as } m\to \infty.
\end{equation*}

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\in(0,\infty)$. 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 $(\eta,q)$,
\begin{equation*}
\int_0^\infty\int_{-\infty}^\infty
  (\eta(U)\phi_t + q(U)\phi_x +\grad\eta(U)G(U)\phi)dxdt
+\int_{-\infty}^{\infty}\eta\bp{U(x,0)}\phi(x,0)dx \ge 0,
\end{equation*}
where $\phi\in C_0^\infty((-\infty,\infty)\times[0,\infty))$ and
$\phi(x,t)\ge 0$.

In particular, we choose
$(\eta(U), q(U)) =\pm\bp{U_j, F_j(U)}, 1\le j\le n,$
in (\ref{eq:4.2}) and conclude that $U$ satisfies the equations
\eqref{1.3} in the sense of distributions.\hfill\qedsymbol \smallskip


As a direct corollary to Theorem \ref{thm:4.1}, we obtain

\begin{thm}\label{thm:4.2} If $U(x,t)=\p{v,u,e+u^2/2,Y}$ is an entropy
solution of \eqref{1.4}, then 
%\begin{thm}\label{thm:4.2} Suppose that the function
%$U(x,t)=(v, u, e+\frac{u^2}{2}, Y)$ is the limit of a sequence of
%approximate solutions $U^{h_\ell}(x,t)$ as constructed in Section
%{\rm \ref{sect:2.2}} by using any consistent approximation algorithm to
%\eqref{ivp}.
%Then the function $U(x,t)$ is an entropy solution {\rm \eqref{eq:4.2}-\eqref{eq:4.3}}
%satisfying
\begin{gather}\label{4.26}
S_t \ge \myf{q\phi(T)Y}{T},\\
%%Change:
%\p{Y^2}_t + \phi(T)Y^2 \le 0,
%%To:
\p{\myf{Y^2}{2}}_t + \phi(T)Y^2 \le 0,\label{eq:4.26a}
\end{gather}
in the sense of distributions.
\end{thm}

%%Change:
%In fact, one can check that the entropy pairs $(\eta_*, q_*)$:
%\begin{equation*}
%\eta_*(U)= -S(v, e)+ \frac{1}{2}Y^2, \quad  q_*=0
%\end{equation*}
%is a convex entropy pair by a careful calculation, and (\ref{4.26})
%follows.
%%To:
In fact, one can check that both $\bp{-S(v,e),0}$ and $(Y^2/2,0)$
are convex entropy pairs by a careful calculation. Then  \eqref{4.26} and \eqref{eq:4.26a}
follow from \eqref{eq:4.2}.

%%Change:
%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 what is customary for homogeneous hyperbolic 
%systems of conservation laws.
%%To:
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 $\eta(U)$ be a convex ``entropy'' function (a.k.a. convex
extension) for a hyperbolic system of balance laws
\begin{equation*}
U_t + F(U)_x = G(U).
\end{equation*}
After multiplying by $\grad \eta(U)$, one obtains the following extra
balance
law for classical solutions:
\begin{equation}
\eta(U)_t + q(U)_x = \grad\eta(U)\cdot G(U).
\end{equation} For discontinuous weak solutions, the corresponding
entropy condition is \eqref{eq:4.2}.
As long as $G(U)$ is a locally integrable function of $(x,t)$, and
not a singular measure, 
the right-hand side of \eqref{eq:4.2} 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$ ($\eta(U)=-S$),
then the
inequality \eqref{eq:4.2} takes on more significance.


\section{Conclusions and Remarks}\label{sect:6}
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
\eqref{eq:3.32c} for $0< \eps\leq \eps_0$, or \eqref{eq:3.33a} for
$\eps_0 < \eps\leq 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\to \infty$ or as $x\to -\infty$--for
example, if $Y_0(x)$ has compact support or vanishes
at either $\infty$ or $-\infty$.

Secondly, we would like to be able to state a simple result in the
same form as that obtained in \cite{hyp:Ni68}, namely that
$\eps\Var\bp{U_0} \leq 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 \cite{hyp:Temp81} do not have this form.  Those results 
required:
\begin{equation}\label{eq:5.1}
V_0=\Var\bp{U_0} \leq \min\p{N,\myf{C(E,N)}{\eps}}.
\end{equation}
We recall from \cite{hyp:Temp81} that
\begin{align*}
C(E,N) &= \myf{1}{K_1 K}\min\p{\eps_0 C_1,\myf{1}{4M_3}
          \min\p{\myf{1}{3G},1-C_0}},\\
 M_3 &= G\p{8+64KK_0GN}\p{1+16KK_0GN},\\
 K_1 &= 8K_0^2 M_3 N.   
  \end{align*}
Thus $C(E,N) \approx
\myf{1}{N(A+BN)^4}$.  For $\eps$ small and $V_0$ large, we can let
$N=V_0$ and \eqref{eq:5.1} becomes:
\begin{equation*}
\eps V_0^2\p{A+B V_0}^4 \leq 1.
\end{equation*}

However, our refinement of the bounds on $F(J_0)$ in terms of $V_0$ in
Lemma \ref{lem:3.13a} improves this requirement somewhat for small
$\eps$.  This Lemma states that for $\eps \leq \eps_0$ $F(J_0) \leq 5
K_0 V_0\leq 5 K_0 N$.  
%%
%Change
%This estimate, plus \eqref{eq:newF} are the
%To
%%
This estimate and \eqref{eq:newF} are the
estimates required in \cite{hyp:Temp81} to prove that $F(J)$ is
non-increasing.
 
Thus, for $\eps$ small, we have $C(E,N)\approx \myf{1}{(A+B N)^2}$, so
that uniform total variation bounds apply to approximate solutions of
the non-reacting equations when
\begin{equation*}
\eps V_0\p{A+B V_0}^2 \leq 1,
\end{equation*}
or, ignoring lower-order terms,
\begin{equation*}
\eps B V_0^3 \leq 1,
\end{equation*}
\begin{comment}
We summarize the results of the previous sections with the following
theorem.

\begin{thm}\label{thm:5.1} There exists a constant $C_0>0$ such that,
when
\begin{equation}\label{eq:6.1}
(\gamma-1)TV(v_0,u_0,p_0,Y_0)(x) \le C_0 < \infty,
\end{equation} there exists a global BV entropy solution $(v,u,p,Y)$
for the Cauchy problem (2.1) satisfying  \eqref{eq:4.2} and (4.3).
\end{thm}
\end{comment}
\begin{rem}
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 $\phi(T) = 0$ for $T <
T_i$, which means that the reacting rate function is discontinuous.
In this case, $T_i$ is called an {\it 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 
\cite{hyp:Chen92,hyp:CHT}.
The existence of traveling wave solutions, and the ``ZND'', or
vanishing viscosity limit of such solutions, was discussed in
\cite{hyp:Wg89:1,hyp:Wg89:2}.

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.
\end{rem}

\begin{rem}
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
(\ref{1.4}), or (\ref{1.2}), and (\ref{2.1}). 
In this context we refer the reader to Liu
\cite{hyp:Liu79} for linear and nonlinear large-time behavior of
solutions of general systems of conservation laws, and to
Glimm-Lax \cite{hyp:GlLa}, DiPerna \cite{hyp:Di75}, Dafermos
\cite{hyp:Da89} 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 \cite{hyp:Chenfrid1,hyp:Chenfrid2}
for new analytical frameworks developed recently.
\end{rem}

\noindent
{\bf 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 \cite{hyp:SW}.

\renewcommand{\baselinestretch}{1.1}

\bibliography{hyperbolic,combustion,analysis}
\bibliographystyle{plain}

\end{document}


