Discrete Shocks for Finite Difference Approximations to Scalar Conservation Laws.
Numerical simulations often provide strong evidences for the existence and stability of discrete shocks for certain finite difference schemes approximating conservation laws. This paper presents a framework for converting such numerical observations to mathematical proofs.  The framework is applicable to conservative schemes approximating stationary shocks of one-dimensional scalar conservation laws. The numerical flux function of the scheme is assumed to be twice differentiable but the scheme can be nonlinear and of any order of accuracy. To prove existence and stability, we show that it would suffice to verify some simple inequalities, which can usually be done using computers. As examples, we use the framework to give an unified proof of the existence of continuous discrete shock profiles for a modified first-order Lax--Friedrichs scheme and the second-order Lax--Wendroff scheme. We also show the existence and stability of discrete shocks for a third-order weighted essentially nonoscillatory (ENO) scheme.
Introduction
In this paper, we provide a general framework for proving the existence and stability of continuous
discrete shock profiles for conservative finite difference schemes which approximate
scalar conservation laws
Research supported by ARPA/GNR grant N00014-92-J-1890.
We consider schemes of conservative form:
with
Here \Delta) is the numerical flux of the scheme which satisfies g(u; u;
(consistency) and is twice continuously differentiable w.r.t. its arguments; u n
j is an approximation
to u(j \Deltax; n\Deltat)
is a constant integer such that (2p+1) is
the stencil width of the scheme. Schemes with such flux functions include the first order
Lax-Friedrichs scheme and some of its modified versions, the second order Lax-Wendroff
scheme and a class of high resolution weighted ENO schemes[7].
Let be two constants such that Eq. (1.1) with the initial data
admits an entropy satisfying shock given by
where s is the shock speed and by Rankine-Hugoniot condition,
In this paper, we will assume
Let us clarify the concepts we will use frequently in the paper:
Definition 1.1: For
1. sup
2. lim
3.
we call ' an approximate stationary discrete shock with parameter q. If, furthermore, ' also
satisfies
4.
we call it an exact stationary discrete shock for scheme (1.2) with parameter q.
Definition 1.2: If a function '(x)(x 2 R) is bounded, uniformly Lipschitz continuous in R
and for any q 2 [0; 1], is an exact stationary discrete shock for scheme
(1.2) with parameter q, then ' is called an continuous stationary discrete shock profile for
scheme (1.2).
Remark: An approximate discrete shock is not related to the scheme (1.2). However, we
will only be interested in approximate discrete shocks which are so accurate that condition
4 in Definition 1.1 is almost satisfied. In the following discussions, we will omit "stationary"
when referring to stationary discrete shocks. We will often refer to exact discrete shocks
plainly by "discrete shocks".
Existence and stability of discrete shocks are essential for the error analysis of difference
schemes approximating (1.1). It is well known that solutions to (1.1) generally contain
shocks and numerical schemes unavoidably commit O(1) error around the shocks. It is
important to understand whether this O(1) error will destroy the accuracy of the scheme in
smooth parts of the solution. For conservation laws whose solutions are sufficiently smooth
away from isolated shocks, Engquist and Yu [4] proved that, the O(1) error committed by a
finite difference scheme around shocks will not pollute the accuracy of the scheme in smooth
regions up to a certain time, provided that, (i) the scheme is linearly stable and (ii) the
scheme possesses stable discrete shocks.
This paper was motivated by the work of Liu and Yu [9]. Their approach was to linearize
the scheme around some constructed approximate discrete shocks. Existence and stability
of exact discrete shocks were then obtained by proving that this linearized scheme defines
a contractive mapping for small perturbations on the approximated discrete shock and the
orginal scheme behaves closely like the linearized scheme.
Our main observation is that, using computers, one can easily obtain approximate discrete
shocks and most importantly, these approximate discrete shocks can be made as accurate
as the machine limit allows. As we know, if a scheme possesses a discrete shock, it can be
often observed from numerical experiments that the scheme converges quickly to a numerical
discrete shock after a number of time iterations over an initial guess (this is often equivalently
stated as "the residue quickly settles down to machine zero"). When it is linearized around
such accurate approximate discrete shocks, the finite difference scheme can be expected
to behave very closely like the linearized scheme (e.g. in terms of contractiveness of the
induced mapping). The obvious advantage of using numerically computed approximate
discrete shocks is that it can be applied to almost all schemes and all conservation laws with
little efforts.
A brief review on discrete shocks for finite difference schemes is as follows: The existence
of a discrete shock was first studied by Jennings [6] for a monotone scheme by an L 1 -
contraction mapping and the Brower's fixed point theorem. For a first order system, Majda
and Ralston [10] used a center manifold theory and proved the existence of a discrete shock.
Yu [15] and Michaelson [11] followed the center manifold approach and showed the existence
and stability of a discrete shock for the Lax-Wendroff scheme and a third order scheme,
resp. In [14], Smylris and Yu studied the continuous dependence of the discrete shocks by
extending the functional space for finite difference schemes from L 1 (Z) to L 1 (R) and a
fixed point theorem. All existence theorems above require an artificial assumption on the
shock speed. In [8], Liu and Xin proved the existence and stability of a stationary discrete
shock for a modified Lax-Friedrichs scheme. For a first order system, Liu and Yu [9] showed
both the existence and stability of a discrete shock using a pointwise estimate and a fixed
point theorem as well as the continuous dependence of the discrete shock on the end states
Our paper is organized as follows:
First we prove a basic fixed point theorem in Section 2. Then we show in Section 3 how
the existence and stability problem can be formulated into a fixed point problem once an
approximate discrete shock is available. In Section 4, we derive sufficient conditions for the
scheme (1.2) to possess a single stable discrete shock and in Section 5, we derive sufficient
conditions for scheme (1.2) to possess a continuous discrete shock profile. In Section 6, we
discuss how to obtain approximate discrete shocks and how to verify the conditions derived
in earlier sections by computers. In Section 7, we apply our framework to give a unified proof
of the existence of continuous discrete shock profiles for a first order modified Lax-Friedrichs
scheme and the second order Lax-Wendroff scheme. We will also show the proof of existence
and stability of discrete shocks for a third order weighted ENO scheme. Some remarks will
be given in Section 8.
A Basic Fixed Point Theorem
In this section, we prove a basic fixed point theorem. First let us define a weighted l 2 norm.
Suppose that ff ? 1 and fi ? 1 are two constants. For any infinite dimensional vector
becomes the regular l 2 norm. We denote the corresponding space by
l 2
Because ff ? 1 and fi ? 1, it is easy to see that any vector v 2 l 2
lim
We denote a closed ball with radius r around a vector v 0 in l 2
Let F be a mapping from l 2
ff;fi to l 2
ff;fi and have the form
where
1. L is mapping which is linear, i.e. L[av
ff;fi , and contractive, i.e.
2. N is generally a nonlinear mapping and N[- where - is the null vector. Moreover,
v;w2Br (-)
3. E is a constant vector in l 2
independent of v.
We have the following basic fixed point theorem:
Theorem 2.1 fixed point theorem) If there exists oe ? 0 such that
Then the mapping F
1. is contractive in B oe (-) under the norm jj \Delta jj ff;fi ;
2. has a unique fixed point, i.e. there exists only one - (-) such that -
Moreover,
Proof: For any v; w 2 B oe (-), we have
Therefore
Moreover, for any v 2 B oe (-), since
From the condition (2.5), we obtain
Condition (2.5) implies that jjLjj ff;fi
ff;fi is strictly less than 1. Therefore (2.7) and (2.8)
together imply that the mapping F maps B oe (-) into itself and is a contractive mapping.
From Banach fixed point theorem, there exists a unique -
(-) such that F In
addition, we have
from which (2.6) follows. Q.E.D.
3 Formulation of a Fixed Point Problem
For a fixed q 2 [0; 1], assume -
' is an approximate stationary discrete shock with parameter
q and is given. By Definition 1:1, we know
lim
Assume that there exists an exact discrete shock for scheme (1.2) with parameter q
denoted as '. If we set
then condition 2 in Definition 1.1 and Eq. (3.3), resp. imply
lim
Condition 4 in Definition 1.1 gives
Here -
u means a vector sum of -
' and -
u, i.e. ( -
we have -
lim
Using (1.3) to write (3.5) in terms of -
v, we have
io
where we have used
lim
which is implied by (3.2), (3.6) and the consistency and continuity of the flux function g.
Let us define the right-hand-side (RHS) of (3.7) as a mapping in some subspace of the
space of infinite dimensional vectors, namely, let
with v be any vector satisfying (3.6), i.e. lim
0: Then Eq. (3.7) gives
which means that -
v is a fixed point of F .
If we reverse the above arguments, namely, assume that there exists an infinite dimensional
vector - v such that it satisfies (3.6) and is a fixed point of the mapping F defined in
(3.8), then it is easy to see that ' j -
is an exact stationary discrete shock for scheme
(1.2) with parameter q. So to prove the existence of a stationary discrete shock is equivalent
to prove the existence of a fixed point for the mapping F . We will restrict the space for the
search of fixed points of F to l 2
ff;fi in our study, where ff and fi are some suitable constants.
We can rewrite the mapping F in (3.8) in the form of (2.2) with
Here,
are the first order partial
derivatives of the numerical flux function g.
It is easy to see that the mapping L is linear and the mapping N is generally nonlinear
and satisfies N[-. In addition, E is a vector depending on -
but independent of v. The
linear mapping L is just the linearization of the mapping F around the approximate discrete
shock -
'. If -
' is accurate enough, we hope L becomes contractive under the weighted l 2 norm
To see the mapping N would satisfy (2.4), we write, for
any v; w
Z 1Z 1g 00
j+l djd-
where
and g 00
are the second order
derivatives of the flux function g. Notice that the summation on the RHS of (3.13) is a
double summation over k; Introducing the shifting operator (or a mapping)
which is defined as E k [v] vector v, we can rewrite (3.13) as
Z 1Z 1g 00
j+l djd-
If we apply the norm jj \Delta jj ff;fi on both sides of (3.15), due to the fact that v   is O(r), the mapping
would satisfy (2.4) provided that the second order derivatives of g are bounded. We will
give precise estimates of jjN jj r
ff;fi in the next two sections for different choices of approximate
discrete shocks -
' and slightly different forms of the mapping N . Then we derive sufficient
conditions on -
' and the first and second derivatives of g which will guarantee the existence
and stability of exact discrete shocks.
Our estimates will be based upon the following two bounds on the first two derivatives
of g:
The functions can be obtained analytically from the given flux function g.
Without loss of generality, we assume that both \Gamma 1 (r) and \Gamma 2 (r) are non-decreasing function
for r - 0.
4 A Single Discrete Shock Profile
For any fixed q 2 [0; 1], we estimate the upper bounds of jjLjj ff;fi and jjN jj r
assuming
that an approximate discrete shock -
' with parameter q is known. We then give a sufficient
condition which ensures the existence and stability of an exact discrete shock for scheme (1.2
with parameter q. The sense of stability will be made precise at the end of this section.
First, we estimate the upper bound of jjLjj ff;fi . Let us write the linear mapping L in the
form of matrix vector product:
thinking v as a column vector with the j th row entry being v j (j 2 Z). According to (3.10),
the infinite dimensional matrix A is given by
refers to the entry of A on i th row and j th column.
otherwise. We define D to be an infinite dimensional diagonal matrix with the i th diagonal
entry being
Use jj \Delta jj 2 and (\Delta; \Delta) 2 to stand for the norm and the inner product in l 2 . It is easy to see that,
for any v 2 l 2
ff;fi , we have jjvjj
(D \Gamma1 is the inverse of D), we
have
A T ~
A T ~
A T ~
Since ~
A T ~
A is symmetric, its l 2 norm is just its spectral radius, ae( ~
A T ~
A). Using Gerschgorin
Circle Theorem from matrix theory, we have
ae( ~
A T ~
A T ~
Notice that A is banded with bandwidth 2p + 1. Therefore ~
A T ~
A is also banded with band-width
not more than 4p + 1. We have the following upper bound for jjLjj
Lemma 4.1
~
A T ~
A
For later use, we define
~
A T ~
A
Here ~
with A and D given by (4.1) and (4.2).
Next we estimate the upper bound of jjN jj r
We start with a simple lemma on the norm of the shifting operator Z). Due to
the non-unitary weight in the norm jj \Delta jj ff;fi , the shifting operator is not unitary.
Lemma 4.2
Proof: Assume k ? 0, for any v 2 l 2
last equality implies
Similarly, for k ! 0, we have
and (4.5) follows by combining the above two inequalities. Q.E.D.
Recall (3.14), we have for any v; w
sup
k is defined in (3.1). Thus we have
j+l
is the function defined in (3.17). From (3.15) and Lemma 4.2, we have
k=\Gammap
Thus we obtain the following upper bound for jjN jj r
Lemma 4.3
are defined in (3.1) and (3.17), resp.
Combining Theorem 2.1 and Lemma 4.1 and 4.3, we obtain:
Theorem 4.1 (Existence and stability of a single discrete shock) If there exist ff ?
and an approximate discrete shock -
' with parameter q such that
oen
where ffi is defined in (4.4), then
1. the finite difference scheme (1.2) possesses an exact stationary discrete shock ' with
parameter q; Moreover, k j sup
2. ' is stable in B oe(-) in the following sense: For any v 2 B oe(-), lim
under maximum norm. Here L n [\Delta] means iterating the finite difference scheme (1.2) n
times using its argument as the initial vector.
Proof:
1. Based on the discussions in Section 3, the existence of an exact discrete shock for
scheme (1.2) is equivalent to the existence of a fixed point for the mapping F in (3.8).
This mapping can be put in the form of L[v] +N [v]
(3.12). By Lemma 4.1 and 4.3, condition (4.7) implies condition (2.5) in Theorem 2.1
(with an extra factor of 1
therefore the mapping (3.8) is contractive in B oe (-) and,
as a result of this, possesses a fixed point -
v. It is easy
to see that ' is an exact stationary discrete shock for scheme (1.2) with parameter q.
Moreover,
Due to the extra factor 1in (4.7), by the second conclusion in Theorem 2.1, namely,
(2.6), we have jj-vjj ff;fi - oe
2. For any v 2 B oe(-), we write
oe and the fact that the mapping F is contractive in B oe (-), we
know that, after applying the mapping F infinitely many times on -
v, the mapping
will converge to the fixed point - v. By the equivalence between the application of the
mapping F on iteration of the scheme (1.2) with initial vector
established in Section 3, we conclude that ' is stable in B oe(-).
5 Continuous Discrete Shock Profiles
In this section, we derive sufficient conditions for the existence of a continuous discrete shock
profile for the scheme (1.2).
For any fixed q 0 2 [0; 1], assume now the conditions in Theorem 4.1 are satisfied. Namely,
there exist ff ? and an approximate discrete shock -
' q0 with parameter q 0
such that condition (4.7) is true. Then by Theorem 4.1, there exists an exact discrete shock
Here we have added superscripts and subscripts q 0 to indicate the dependence on parameter
. The constants ff and fi, however, will be chosen to be independent of q 0 .
For proper choices of ff and fi, we are interested in finding the conditions on the approximate
discrete shock -
and the first two derivatives of the numerical flux function g
such that an exact discrete shock ' q for the scheme (1.2) is guaranteed to exist for any q
in a small neighborhood of q 0 , say [q
Once such
conditions are found and satisfied for a finite number of values of q 0 , e.g.
it becomes clear that for any q 2 [0; 1], there exists an exact discrete shock for scheme (1.2).
A continuous discrete shock profile is then obtained by properly arranging the family of exact
discrete shocks which are parameterized by q 2 [0; 1].
Let us take an approximate discrete shock with parameter q 2 [q
2M ] to be
Zg. It is easy to check that conditions 1,2,3 in Definition 1.1 are
satisfied. We have
sup
Define a mapping based on the approximate discrete shock -
' q in (5.1):
where v is any vector in l 2
ff;fi . We can rewrite the mapping F in the form of (2.2) with
It is easy to see L is a linear mapping and its norm can be estimated similarly as in Lemma 4.1,
namely, we have
A T
~
A q0
and similar to (4.4), we define
A T
~
A q 0
Here ~
is defined in (4.2); A q0 is given in (4.1) with -
' replaced by -
The mapping N satisfies N[- and its norm jjN jj r
ff;fi can be estimated as follows: For
any v; w same as (3.15), we have
Z 1Z 1g 00
j+l djd-
but with (different from (5.11))
Using (5.2) and (5.3), we have
sup
sup
Thus similar to Lemma 4.3, we have
Now we estimate jjE jj ff;fi . Since ' q0 is an exact discrete shock profile, we have '
which implies
We can write
In the third equality above, we have used the definition of -
' q in (5.1) and in the second
j has the bound
sup
Recall the function \Gamma 1 (\Delta) defined in (3.16), we get
Using the fact that jq \Gamma q
2M , we get a larger upper for jjE jj ff;fi which does not depend
on q,
We want to note that the factor jq \Gamma q 0 j on the RHS of (5) will enable us to prove uniform
Lipschitz continuity of a family of exact discrete shocks when they exist.
Based on the bounds found in (5.8), (5.11) and (5.13), we obtain the following sufficient
condition for the existence of a continuous discrete shock profile for scheme (1.2):
Theorem 5.1 (Existence of a continuous discrete shock profile) If there exist ff ?
(an integer) such that for each q
there exist oe q0 ? 0
and an approximate discrete shock profile -
' q0 for which the following inequality is true:
oe q0(
Here E q0 is given by
is given in are two functions defined in (3.16) and
(3.17). Then
1. for any q 2 [0; 1], there exists an exact discrete shock ' q for scheme (1.2);
2. if for any x 2 R, we define '(x)
are uniquely determined by
is a continuous discrete shock
profile for scheme (1.2).
Schematic proof:
1. Condition (5.14) clearly implies condition (4.7) in Theorem 4.1 for Therefore
there exists an exact discrete shock for scheme (1.2) with parameter q 0 . For any
], we can define an approximate discrete shock -
' q as in (5.1).
According to the estimates (5.8), (5.11) and (5.13), condition (5.14) implies (2.5) for
the mapping F in (5.4). By the same logic used in Theorem 4.1, we see that there
exists an exact discrete shock ' q for scheme (1.2) for any q 2 [q
condition (5.14) is true for all q
2. We only need to check that '(x) is bounded and uniformly Lipschitz continuous. For
all
are uniformly bounded by the definition of approximate
discrete shocks and the finiteness of M . Each ' q0 differs from -
' q0 by a vector whose
maximum can be bounded by oe q 0
, so does ' q differ from -
. Due to
choice of -
' q in (5.1), it is easy to see that ' q is uniformly bounded for q 2 [0; 1]. To
prove that '(x) is uniformly Lipschitz continuous, we first give two observations which
can be shown easily:
Observation #1: For any q 2 [0; 1], there exists q
M for some integer i between 0
and M , such that
sup
where the constant C are independent of q and q 0 . This is a result of the estimate
(5.12), the bound (2.6) in Theorem 2.1 and condition (5.14).
Observation #2:
This is due to the way we parameterize the family of discrete shocks, namely the
parameter q in Definition 1.1.
An easy generalization of Observation #1 is that for any q 1
or
for any j 2 Z. Let x
only need to consider the case which we have
Because each term on the RHS of the last equality is of the form of (5.16), the uniform
Lipschitz continuity of '(x) follows.
6 Algorithms for Computer Verification
In this section, we discuss how to use a computer to verify condition (4.7) in Theorem 4.1 to
prove the existence and stability of an exact discrete shock or condition (5.14) in Theorem 5.1
to prove the existence of a continuous discrete shock profile, for scheme (1.2).
6.1 Computing an approximate discrete shock
We start with providing a method of obtaining an accurate approximate discrete shock -
for any fixed q 2 [0; 1] using scheme (1.2). Let J be an integer. Set
where - j is given by (1.4) and f- are chosen such that
For example, we can take
We then apply the finite difference scheme (1.2) to the initial data u 0 repeatedly for sufficiently
many times. Note that when the scheme is applied to u n
need
values of u n
in order to compute values of u n+1
j for all jjj - J . We can
simply set u p. Although this makes
the scheme non-conservative in the bounded region (i.e. jjj - J), it actually does not make
an error much bigger than the machine accuracy if J is taken to be large enough. This is
because an exact discrete shock is generally believed to be converging to the two end states
exponentially fast. For the sake of rigorousness, we can modify the value of u n 0such that
to make the procedure conservative in the bounded region.
Here n 0 is assumed to the number of applications of the scheme on u 0 . To determine how
large n 0 should be, we can monitor sup
is defined in (3.12) with -
' replaced by u n ) to
see if it is small enough, say close to machine accuracy, for the purpose of our verification of
the conditions in Theorem 4.1 or Theorem 5.1. Finally, we can set the approximate discrete
shock to be
It is easy to check that -
' q satisfies the conditions 1 to 3 in Definition 1.1.
Remark: Theoretically the larger J is, the more accurate -
' q one can get. However, larger
J means longer computer time that it takes to verify the conditions. If
to 80 is believed to be good enough. If very small, J needs to be very
large and it may even exceed the computer power. In the latter case, the framework in this
paper may be improper.
6.2 Choosing the constants ff and fi
Once we have an approximate discrete shock -
' q , we can decide what to choose for ff and
fi. The criterion for this is to make the norm jjLjj ff;fi as much below 1 as possible, or in our
estimates, to maximize ffi in (4.4). The range of possible values of ff can be obtained by
studying the linearized scheme of (1.2) around u Similarly, fi can be obtained by studying
the linearized scheme of (1.2) around u \Gamma . See Smylris [13] for details. We can use one
approximate discrete shock or a few such approximate discrete shocks (corresponding to
different values of q in [0; 1]) to choose the constants ff and fi. In the latter case, we should
make the minimum of ffi (over different values of q) above 0 as far as possible. Note that the
matrix A in (4.4), which is given by (4.1) with -
' replaced by -
essentially finite due to
constancy of -
we have a finite row sums to take a maximum of
in (4.4).
6.3 Strategy for verification
We suggest the following strategy for verifying the condition (4.7) in Theorem 4.1. For a
given q 2 [0; 1],
1. Find the function \Gamma 2 (r). Make sure it is non-decreasing in r.
2. Compute an approximate discrete shock following the method described in subsection
6.1. Make it as accurate as possible should the condition (4.7) seems very demanding

3. Find the range of the constants ff and fi from the linear analysis of the scheme and
choose ff ? in (4.4) is positive and maximized.
4. Find the largest possible value for oe for which condition (4.7) is true. Mostly, we can
take
If the above four steps are through, we can conclude that scheme (1.2) does possess an exact
discrete shock with parameter q and it is stable in the sense stated in the second part of
Theorem 4.1.
To verify condition (5.14), we suggest the following steps:
1. Find the functions
2. For several q values, compute an approximate discrete shock -
' q0 for each q 0 . Find the
proper constants ff ? 1 and fi ? 1 such that minimum of the values of ffi q 0
based on
each -
positive and maximized. We then set M - oe \Gamma1
q0 and find oe q 0
such that the
RHS of (5.14) is maximized. Using this oe q0 and replacing - k q0 by
can obtain an estimate of the size of M by requiring the second term on the LHS of
(5.14) to be less than the RHS of (5.14). One can even replace M inside the functions
q0 as long as we eventually use an M - oe \Gamma1
q0 for all values of q 0
sampled (this is due to the monotonicity of \Gamma 1 (\Delta) and \Gamma 2 (\Delta)). We suggest one always
use a bigger M than necessary to attain a bigger margin of the RHS of (5.14) over the
LHS. Usually, one can take
3. For each q
we compute an approximate discrete shock -
parameter q 0 and check if for this -
q0 , there exists oe q0 ? 0 such that (5.14) is true.
If (5.14) is true for every q
conclude that the finite difference
scheme does have a continuous discrete shock profile.
7 Some Examples
In this section, we apply Theorem 5.1 to give an unified proof of the existence of a continuous
discrete shock profile for a modified Lax-Friedrichs scheme and the Lax-Wendroff scheme.
For a third order WENO scheme, we apply Theorem 4.1 to prove the existence and stability
of an exact discrete shock for some sample values of q in [0; 1].
As an example, we take the conservation law to be the Burgers' Eq.:
x
and take the end states \Gamma1. It is well known that an entropy satisfying
stationary shock exists for these two end states.
7.1 A modified Lax-Friedrichs scheme
The flux function for a modified Lax-Friedrichs scheme is
It is clear that the stencil width constant We take the upper bound of the first and
second derivatives of this flux function to be \Gamma 1

Appendix

A.1).
For 0:5, we use 41 points (i.e. to compute approximate discrete shocks and
choose to define the weighted l 2 norm. Roughly oe q 0
12ff maximizes the RHS of
and M can be estimated by 64ff 2
. We take 10000 and is able to verify condition
. Thus this modified Lax-Friedrichs scheme possesses a
continuous stationary discrete shock profile for Burgers' Eq. with end states
\Gamma1. The constant oe 0 , which represents the size of the stability region for the discrete
shocks (see conclusion 2 in Theorem 4.1), is approximately 7 \Theta 10 \Gamma3 . We plot values of the
LHS and RHS of (5.14) for 100 even spaced samples of q 0 in [0; 1] in Figure 1a. The discrete
shock profile is plotted in Figure 1b.
7.2 The Lax-Wendroff scheme
The flux function for the Lax-Wendroff scheme is
The stencil width constant

Appendix


A.2). For 0:5, we use 61 points (or to compute approximate discrete shocks.
For the weighted l 2 norm, we take Roughly oe q 0
24ff maximizes the RHS of
and M can be estimated by 96ff 2
. We take and is able to verify (5.14)
for
. Thus the Lax-Wendroff scheme possesses a continuous stationary
discrete shock profile for Burgers' Eq. with end states \Gamma1. The constant
for the size of the stability region, namely, oe 0 , is approximately 2:3 \Theta 10 \Gamma3 . We plot the LHS
and RHS of (5.14) for 100 even spaced samples of q 0 2 [0; 1] in Figure 1c. The discrete shock
profile is plotted in Figure 1d.
7.3 The third order WENO scheme
The WENO (weighted essentially non-oscillatory) schemes [7] are variations of the ENO
(esstially non-oscilltory) schemes [12]. They both achieve essentially non-oscillatory property
by favoring information from the smoother part of the stencil over that from the less smooth
(a)
(b) -20
(d)
Figure

1: (a) The modified Lax-Friedrichs scheme. Solid line: the value of the LHS of (5.14);
Dashed line: the value of the RHS of (5.14). (b) A discrete shock profile for the modified
Lax-Friedrichs scheme. (c) Same as (a) but for the Lax-Wendroff scheme. (d) Same as (b)
but for the Lax-Wendroff scheme.
or discontinuous part. However, the numerical flux function of the ENO schemes is at most
Lipschitz contiuous while the numerical flux function of the WENO schemes are infiinitely
smooth (if one takes ffl w appearing below to be nonzero).
The numerical flux function for the third order WENO scheme with global Lax-Friedrichs
flux splitting is
and
f \Sigma (z) =2 (f(z) \Sigma \Lambdaz)
Here, ffl w is a small constant to avoid the denominator to be zero and is taken as ffl
is a constant which is the maximum of jf 0 (u)j over all possible values of u. In our case, we
take   to be slightly above the maximum of the modulus of the two end states, Namely, we
set  1:1. The first and second order partial derivatives of the numerical flux g are shown
in the

Appendix

.
This scheme is third order accurate in space where the solution is monotone and smooth.
It degenerates to second order at smooth extrema. See [7] for details. If we use Euler Forward
in time, the scheme has the form
However this scheme is linearly unstable for any constant - ? 0. Abbreviate the RHS of
express the scheme with the third order Runge-Kutta scheme [12] in
time as
We abbreviate (7.6) as
We have two observations: (i) if f'
as in (3.8), if this mapping is contractive
under jj \Delta jj ff;fi for some ff ? 1 and fi ? 1, then the mapping derived from
contractive under the same norm.
The first observation is obvious. The second one is due to the fact that each stage in the
third order Runge-Kutta scheme (7.6) is a convex combination of u n and E[u] where u is u n
in stage one, u (1) in stage two and u (2) in the final stage. See [12] for details.
Therefore, in order to prove the existence and stability of exact discrete shocks for (7.6),
it suffices to prove the existence and stability of exact discrete shocks for (7.4). We have
attempted to apply Theorem 5.1 to prove the existence of a continuous discrete shock profile
for (7.4) and found that we need sample roughly 10 19 even spaced values of q 0 in [0; 1] which
is far beyond the computer power. Nevertheless, we are able to use Theorem 4.1 to prove the
existence and stability of exact discrete shocks for (7.4) for many sample values of q 0 2 [0; 1].
Our computer verification strongly indicates that a continuous discrete shock profile does
exist for this scheme.
Here are the details of the computer verification:
It is clear that the stencil width constant p equals 2. We take the upper bound for the
second derivatives of the numerical flux function to be
Its derivation is detailed in Appendix A.3. We have taken used 161 points (or
to compute the approximate discrete shocks. ff and fi are both taken to be 1:8. The
condition in Theorem 4.1 is verified for q This verification is done on
CRAY C-90 using double precision.
We plot the LHS and RHS of (4.7) for the 1000 even spaced samples of q 0 in Figure 2a. It
can be seen that the curve for the LHS is properly below the curve for the RHS. We believe
it is true for all q 0 2 [0; 1]. The discrete shock profile with 40 even spaced samples of q 0 is
plotted in Figure 2b.
The discrete shock profile appears to be monotone. In Figure 2cd, we show the profile
zoomed around between the grid points indexed from \Gamma10 to 0 and around
between points indexed from 0 to 10, resp. Notice that the profile contains very small
oscillations of magnitude around 10 \Gamma4 on both sides of the shock. However the profile looks
very smooth, which leads us to believe that a continuous discrete shock profile does exist
for this third order WENO scheme. We have yet been able to find a less stringent condition
than that in Theorem 5.1 in order to prove this by computer.
Concluding Remarks
We have provided sufficient conditions for a conservative scheme (1.2) to have a single discrete
shock (Theorem 4.1) or a continuous discrete shock profile (Theorem 5.1) for a scalar
conservation law in one dimension. These conditions can usually be verified by computers as
(a)
(b) -80 -60 -40 -20
-0.50.5(c)
(d)

Figure

2: The third order WENO scheme. (a) Solid line: the value of LHS of (4.7); Dashed
line: the value of the RHS of (4.7). (b) The discrete shock profile. (c) Zoom of (b) around
Zoom of (b) around
demonstrated in the last section. The key idea here was to linearize the scheme around accurate
numerical approximations of the discrete shocks and find suitable weighted l 2 norms
for this linear part to define a contractive mapping. If we can find sufficiently accurate
approximate discrete shocks, the orginal scheme behaves closely like the linearized scheme
around this approximate discrete shock in terms of contractiveness of the induced mapping.
Several generalizations or implications of Theorem 4.1 or Theorem 5.1 are immediate.
For example, we can find sufficent conditions, in the form of a single inequality, which assure
the existence of discrete shocks or a continuous discrete shock profile for a range of the
time-space ratio - and the end states Using the result in [4], we can generalize the
theorems to non-stationary shocks with sufficiently small shock speed. Pointwise convergence
rate estimates can also be obtained for schemes that possess stable discrete shocks and are
linearly stable, when used to approximate scalar conservation laws whose solution is smooth
except for some isolated shocks. However, we will not elaborate on such extensions.
9

Acknowledgments

We would like to thank Bj-orn Engquist and Stanley Osher for their support in this research.
The first author wants to thank Chi-Wang Shu for valuable suggestions.
Functions
Let
jzj-r
where f is the flux function of the conservation law (1.1). We derive the functions \Gamma 1 (r) and
for the numerical flux functions of the schemes discussed in Section 7.
For simplicity, we use g 0
k to denote the first order derivative of g(z
and use g 00
k;l to denote its second order derivative w.r.t z k and z l , where k;
A.1 The modified Lax-Friedrichs scheme
The numerical flux function is given in (7.1). Its first order derivatives are
Its second order derivatives are
Therefore, we have
A.2 Lax-Wendroff scheme
The numerical flux function is given in (7.2). Its first order derivatives are
Its second order derivatives are
1). Therefore, we have
A.3 The third order WENO scheme
The flux function for the third order WENO scheme is given in (7.3). To obtain \Gamma 1 (r) and
first find the first two derivatives of the flux function.
The derivatives of the function
are:
r 00
r 00
r 00
The derivatives of the function
a /
a
a
aa /)
a
a r 0
ab /)
Let
2. We define
Similarly we define by / \Sigma
aa / \Sigma
ab
bb the first and second order derivatives of / with
arguments, resp. same as / \Sigma . We drop 0 and 00 in the notation for derivatives of / for
simplicity.
The first order derivatives of the numerical flux function (7.3) are
a
a
2. The second order derivatives are
aa
ab
ab
ab
bb
ab
ab
bb
ab
aa
Other second order derivatives are know by symmetry, i.e. g 00
l;k .
We have the following simple observations:
rr 00
where r 0 can be r 0
a or r 0
can be r 00
aa or r 00
bb or r 00
ab and similarly, / 0 can be / 0
a or / 0
be / 00
aa or / 00
bb or / 00
ab .
Based on the above observations, we have



--R

Convergence of the Lax-Friedrichs scheme for isentropic gas dynamics III
Convergence of difference scheme with high resolution to conservation laws.
Courant and Friedrichs
Convergence of a finite difference scheme for a piece-wise smooth solution
Viscous limits for piecewise smooth solution to system of conservation laws
Gray Jennings Discrete shocks
Efficient implementation of Weighted ENO schemes

Construction and nonlinear stability of shocks for discrete conservation laws
Discrete shock profiles for systems of conservation laws
Discrete shocks for difference approximations to system of conservation laws
Efficient Implementation of essentially non-oscillatory shock-capturing schemes
Existence and Stability of Stationary profiles of the LW Scheme
On the existence and stability of a discrete shock profile for a finite difference scheme

--TR

--CTR
Hailiang Liu, The l1 global decay to discrete shocks for scalar monotone schemes, Mathematics of Computation, v.72 n.241, p.227-245, 01 January
conservation law;weighted ENO;discrete shock