\magnification = \magstephalf
\hsize=14truecm
\hoffset=1truecm
\parskip=5pt
\overfullrule=0pt
\nopagenumbers \pageno=181
\input amssym.def % The R for Real nunbers.
\input epsf
\font\eightrm=cmr8 \font\eighti=cmti8 \font\eightbf=cmbx8
\headline={\ifnum\pageno=1 \hfill\else%
{\tenrm\ifodd\pageno\rightheadline \else
\leftheadline\fi}\fi}
\def\rightheadline{\hfil TVD METHODS \hfil\folio}
\def\leftheadline{\folio\hfil Xuefeng Li \hfil}
\voffset=2\baselineskip
\vbox {\eightrm\noindent\baselineskip 9pt %
Differential Equations and Computational Simulations III\hfill\break
J. Graef, R. Shivaji, B. Soni \& J. Zhu (Editors)\hfill\break
Electronic Journal of Differential Equations, Conference~01, 1997, pp. 171--191.\hfill\break
ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu
\hfil\break ftp 147.26.103.110 or 129.120.3.113 (login: ftp)}
\footnote{}{\vbox{\hsize=10cm\eightrm\noindent\baselineskip 9pt %
1991 {\eighti Subject Classification:} 65C20, 65M12, 65M06.
\hfil\break
{\eighti Key words and phrases:} Conservation Laws, Godunov Method,
Entropy Condition, Convergence, High Accuracy.
\hfil\break
\copyright 1998 Southwest Texas State University and
University of North Texas. \hfil\break
Published November 12, 1998.} }
\bigskip\bigskip
\centerline{ENTROPY CONSISTENT, TVD METHODS WITH}
\centerline{HIGH ACCURACY FOR CONSERVATION LAWS}
\medskip
\centerline{Xuefeng Li}
\bigskip
{\eightrm\baselineskip=10pt \narrower
\centerline{\eightbf Abstract}
The Godunov method for conservation laws produces numerical solutions
that are total-variation diminishing (TVD) and converge to weak solutions
which satisfy the entropy condition (Entropy Consistency),
but the method is only first order accurate.
Many second and higher order accurate Godunov--type methods have been
developed by various researchers.
Although these high order methods perform very well numerically,
convergence and entropy-consistency has not been proven,
maybe due to the highly nonlinear approach.
In this paper, we develop a new class of Godunov--type methods that are
TVD, converge to weak solutions of conservation laws, and satisfy the
entropy condition. The error produced by these methods are theoretically
controllable by the choice the piecewise constant functions used in the
numerical approximation. Numerical experiments confirm that our methods
produce numerical solutions that are comparable to those produced by higher
order methods, while maintaining all the good characteristics of the
Godunov method.
\bigskip}
\def\half{{1\over 2}}
\bigbreak
\centerline{\bf 1. Introduction} \medskip\nobreak
A new class of numerical methods for nonlinear systems of hyperbolic
conservation laws,
$$\eqalignno{\partial_t u +\partial_x f(u)=0, & \quad
(x,t)\in (-\infty,+\infty)\times (t_0,+\infty),&(1.1a)\cr
u(x,t_0)=u_0(x), &\quad x\in(-\infty,+\infty), &(1.1b)\cr}$$
is presented in this paper. The solution vector $u$ has $m$ components
and is a function of two variables $x$ and $t$. That is,
$u=u(x,t)=(u_1,\ldots,u_m)^T$, $f(u)=(f_1(u),\ldots,f_m(u))^T$,
and matrix $\partial_uf$ has $m$ distinct real eigenvalues.
It is well known that a typical solution of (1.1) develops jump
discontinuities (called shocks) in finite time [Lax, 1973].
Thus solution at large
exists only in a weak sense. A weak solution of (1.1) is defined as a
function $u(x,t)$ that satisfies
$$\int_{t_0}^{+\infty}\!\!\int_{-\infty}^{+\infty}[(\partial_t\phi) u
+(\partial_x\phi) f(u)]\,dx\,dt+\int_{-\infty}^{+\infty}
\phi(x,t_0)u_0(x)\,dx=0, \eqno(1.2)$$
for all $\phi(x,t)\in C^{\infty}$ with compact support in the half plane
$(-\infty,+\infty)\times (t_0,+\infty)$.
Since weak solutions of (1.1) are not unique [Lax, 1973], additional conditions must be
used to identify the physically relevant solution, or the entropy
solution $u(x,t)$ (piecewise continuous)
that satisfies the entropy
conditions
$$\partial_t U(u)+\partial_x F(u)\leq 0,\quad
(x,t)\in (-\infty,+\infty)\times (t_0,+\infty),\eqno(1.3a)$$
in weak sense, or
$$-\int_{t_0}^{+\infty}\!\!\int_{-\infty}^{+\infty}[(\partial_t\phi) U(u)
+(\partial_x\phi) F(u)]\,dx\,dt-\int_{-\infty}^{+\infty}\phi(x,t_0)U(u_0(x))\,dx\leq 0,
\eqno(1.3b)$$
for all nonnegative $\phi(x,t)\in C^{\infty}$ with compact support in the half
plane $(-\infty,+\infty)\times (t_0,+\infty)$. In addition, across
a discontinuity of $u(x,t)$ with the speed of propagation being $S$,
there must be
$$F(u_r)-F(u_l)-S [U(u_r)-U(u_l)] \leq 0\,.\eqno(1.3c)$$
Here, $u_l$ and $u_r$ are the states of $u$ on the left and right
of the discontinuity of $u(x,t)$; $U(u)$ is the entropy
function of (1.1) and $U$ is convex, that is,
$$U(\half (p+q))\leq \half (U(p)+U(q))\,.\eqno(1.4)$$
$F(u)$ is the entropy flux of (1.1) and satisfies $\partial_u U\partial_u f=
\partial_u F$. The existence of entropy $U$ and entropy flux $F$ of (1.1) is
assumed here. For most conservation laws, this assumption is
indeed true [Harten, 1983].
In the case of scalar conservation laws, the entropy condition (1.3)
ensures the uniqueness of weak solution of (1.1) in the class of
piecewise continuous functions as stated in the following theorem.
\proclaim {Theorem 1.1 [Oleinik, 1957; Lax, 1973] (Uniqueness)}.
A weak piecewise continuous solution $u(x,t)$ of
a scalar conservation law in the form of (1.1) is unique
in the class of piecewise continuous functions (finitely
many discontinuities inside any compact set in
$x$--$t$ plane), if and only if
$u(x,t)$ satisfies (1.3).
Let $\{x_{j+\half}\}$ be a set of grid points on the $x-$axis. Denote
$[x_{j-\half},x_{j+\half}]$ by $I_j$. Define grid size
$h$ to be
$h=\sup_j|x_{j+\half}-x_{j-\half}|$. Let
$x_j=\half(x_{j-\half}+ x_{j+\half})$, $\Delta x_j=x_{j+\half}-x_{j-\half}$,
$\Delta x_{j+\half} =x_{j+1}-x_j$.
Let $\Delta t_n=\tau$ be the time
step, i.e., $t_{n+1}=t_n+\Delta t_n$ $=t_n+\tau$ ($n=0,1,\ldots$).
For ease of exposition, equal
spacing in $x$ is assumed from now on. That is, $\Delta
x_j=\Delta x_{j+\half}\equiv \Delta x\equiv h$.
Denote $\tau /h$ by
$\lambda$. If we define the averaged value of $u(x,t)$ over $I_j$ by
$u_j^n$, i.e.,
$$u_j^n={1\over h}\int_{I_j}u(x,t_n)dx,\eqno(1.5)$$
it can then be shown using Green's theorem over the rectangle $I_j\times
[t_n,t_{n+1}]$
that $\{u_j^n\}$ satisfy the following relations:
$$\eqalignno{u_j^{n+1}&=u_j^n-\lambda
[f_{j+\half}^n-f_{j-\half}^n],&(1.6a)\cr
f_{j+\half}^n&={1\over {\tau}}\int_{t_n}^{t_n+\tau}
f(u(x_{j+\half},t))dt.&(1.6b)\cr
}$$
In this paper, a new class of numerical methods is developed based on
(1.6). This new class of highly accurate numerical methods are
convergent, and the limits of their numerical solutions satisfy
the entropy condition of (1.3b). And their numerical accuracy
are comparable to those of high order methods.
A general discussion of numerical methods based on (1.6) is presented in
section 2. Approximations of functions using piecewise constant functions
are discussed in section 3. The formulation of the new methods
and further remarks are given in section 4.
Numerical implementation and tests are presented in section 5 in
comparison with results from Godunov method. It can be shown that the newly
developed first order methods
produce numerical solutions with sharp resolution
seen in those produced by high order
methods.
\bigbreak
\centerline{\bf 2. Godunov and Godunov Type methods} \medskip\nobreak
When a numerical solution to (1.1) is to be computed, it is
often the discrete
averaged values of $u(x,t)$ defined in (1.5) that are being calculated, using relations (1.6)
where $\{u_j^0\}$ are
initialized by the initial function $u_0(x)$ in (1.5) and $u(x,t_0)=u_0(x)$.
Many numerical methods are different only in the ways the integral
in (1.6b) is approximated. Evaluation of
the integral in (1.6b) often involves the evaluation of Riemann
problems, or generalized Riemann problems of (1.1). A Riemann problem
of (1.1) is defined as:
$$\eqalignno{\partial_t u +\partial_x f(u)=0,&\quad
(x,t)\in (-\infty,+\infty)\times (t_0,+\infty),\cr
u(x,t_0) =u_l,&\quad xx_0. \cr}$$
That is a conservation law with the initial function as a step function. The
solution of (2.1) is self--similar with respect to the point $(x_0,t_0)$ and is
denoted by $R((x-x_0)/(t-t_0);u_l,u_r)$. And a generalized Riemann problem of
(1.1) is defined as:
$$\eqalignno{\partial_t u +\partial_x f(u)=0,&\quad
(x,t)\in (-\infty,+\infty)\times (t_0,+\infty),\cr
u(x,t_0) =u_l(x),&\quad xx_0. \cr}$$
The solution of a generalized Riemann problem is denoted by
$G(x,t;x_0,t_0,u_l(\cdot),u_r(\cdot))$.
A difference scheme is said to be consistent with conservation law (1.1)
(called a conservative scheme) if it can be represented in the following
format
$$\eqalignno{v_j^{n+1}&=v_j^n-\lambda
[f_{j+\half}^n-f_{j-\half}^n], \cr
f_{j+\half}^n&=g(v_{j-l+1}^n,\ldots,v_{j+l}^n), &(2.3)\cr
g(u,\ldots,u)&\equiv f(u). \cr}$$
A difference scheme is said to be consistent with the entropy condition
(1.3a) if
$$\eqalignno{U_j^{n+1}&\leq U_j^n-\lambda
[F_{j+\half}^n-F_{j-\half}^n], \cr
U_j^n&=U(v_j^n),&(2.4)\cr
F_{j+\half}^n&=g(v_{j-l+1}^n,\ldots,v_{j+l}^n), \cr
g(u,\ldots,u)&\equiv F(u). \cr}$$
The lattice function $\{v_j^n\}$ is usually extended to continuous values
of $x, t$ by setting:
$$v_h(x,t)=v_j^n,\quad (x,t)\in I_j\times [t_n,t_{n+1}).\eqno(2.5)$$
A difference scheme is said to be total variation (TV) stable if
the total variation in $x$ of $v_h(x,t)$,
$$TV(v_h(\cdot,t))\equiv TV(v^n)=\sum_j|v_{j+1}^n-v_j^n|,\eqno(2.6)$$
is uniformly bounded in $t$ and $h$; here $n$ is the integer part of
$t/\tau$. A difference scheme is said to be total variation diminishing (TVD) if:
$$TV(v^{n+1})\leq TV(v^n).\eqno(2.7)$$
Here are two important theorems concerning the
above consistent difference schemes.
\proclaim {Theorem 2.1 [Lax, Wendroff, 1960; Harten, Lax,
Van Leer, 1983]}.
Suppose a difference scheme
is conservative and consistent with the entropy condition (1.3a)
in the form of (2.3). Let
$\{v_j^n\}$ be a solution of (2.3), with initial values $v_j^0=u_j^0$
as defined in (1.5). Let
$v_h(x,t)$ be the extended function of $\{v_j^n\}$ as defined in (2.5).
Suppose that for some sequence $h_k\rightarrow 0^+$,
$\tau_k/h_k=\lambda=$constant, the limit:
$$\lim_{h_k\rightarrow 0^+}v_{h_k}(x,t)=u(x,t),\eqno(2.8)$$
exists in the sense of bounded, $L^1_{loc}$ convergence. Then the limit
$u(x,t)$ satisfies the weak form (1.2) of the conservation law, and the
weak form (1.3b) of the entropy condition.\par
\proclaim {Theorem 2.2 [Harten, 1984]}. Suppose the difference scheme
(2.3) is conservative and TV stable. Then the scheme produces
numerical solution with a convergent subsequence
whose limit (in the sense of bounded, $L^1_{loc}$) is a weak solution of
(1.1).\par
In the case of scalar conservation laws, Theorem 2.1 ensures the uniqueness
of the limit solution when the entropy condition is satisfied, while Theorem
2.2 guarantees the existence of a limit solution if the difference scheme
is TV stable. Therefore, {\it a conservative difference scheme is convergent
and its limit function is the
unique solution of (1.1) if the scheme is TV stable and entropy consistent.
The goal of this paper is to develop a new class of difference methods that
are TV stable and entropy consistent. Furthermore, they produce numerical
solutions with accuracy comparable to those produced by high order
numerical methods.}
Godunov [1959], Van Leer (MUSCL) [1979], Colella (MUSCL) [1985], Colella and
Woodward (PPM) [1982], Harten, Engquist, Osher
and Chakravarthy (ENO) [1987]
developed numerical methods for solving (1.1) which
use relations (1.6) in the numerical approximations. In particular, Godunov
method uses piecewise constant functions (constant over each interval $I_j$ with
possible jumps at $\{x_{j+\half}\}$) to approximate $u(x,t_n)$; thus the
integral in (1.6b) can be evaluated exactly by solving Riemann problems of (1.1)
at $\{x_{j+\half}\}$. In the MUSCL (PPM) scheme, piecewise linear (quadratic)
functions (linear (quadratic) over each interval $I_j$ with possible jumps at
$\{x_{j+\half}\}$) are used to approximate $u(x,t_n)$; the integral in (1.6b) is
approximated using the trapezoidal rule. In the ENO schemes, piecewise
polynomial functions ($N$-th degree polynomial over each interval $I_j$ with
possible jumps at $\{x_{j+\half}\}$) are used to approximate $u(x,t_n)$; the integral in
(1.6b) is approximated using an appropriate $K$-point numerical quadrature. All
of those methods can be formulated in the following fashion:
$$\eqalignno{
v_j^{n+1}&=v_j^n-
\lambda
[f_{j+\half}^n-f_{j-\half}^n],&(2.9a)\cr
f_{j+\half}^n&=\sum_k\alpha_k f
(v_{j+\half}^{n,\gamma_k}),&(2.9b)\cr
\lambda&={{\tau}\over h}\leq {{\theta}\over
\Lambda},&(2.9c)\cr }$$
where $(2.9b)$ is an appropriate numerical quadrature approximating the integral
in $(1.6b)$; $\theta$ is a positive constant less than $1$, and it is called the
CFL number; $\alpha_k$ and $\gamma_k$ are coefficients of the numerical
quadrature; $\{v_{j+\half}^{n,\gamma_k}\}$ are values of $u$
needed in the quadrature; $\Lambda$ is the largest eigenvalue of the matrix
$\partial_u f$ in absolute value, which represents
the maximum speed of propagation of discontinuities of $u(x,t)$.
Equations listed in (2.9) represent a class of recursive algorithms
for solving (1.1). A major task of these algorithms is the evaluation
of $\{v_{j+\half}^{n,\gamma_k}\}$, values of $u(x,t)$
needed in the quadrature in (2.9b). Assuming $u(x,t_n)$ is approximated
by $u_j^n(x)$ over interval $I_j$. Then
$$v_{j+\half}^{n,\gamma_k}
=G(x_{j+\half},\gamma_k\tau;x_{j+\half},t_n,u_j^n(\cdot),u_{j+1}^n(\cdot)),\eqno(2.10)$$
as indicated in (2.2).
In the Godunov method, $\{u_j^n(x)\}$ are chosen to be constants. That is,
$u(x,t_n)$ is approximated using a piecewise constant function
$$P_n(x)=v_j^n={1\over {x_{j+\half}-x_{j-\half}}}
\int_{I_j}u(x,t_n)\,dx \quad \hbox{for } x\in I_j\,.$$
The piecewise constant function
$P_{n+1}(x)$ approximating $u(x,t_{n+1})$ is computed by averaging the exact
solutions of (1.1) over each $I_j$, using $P_n(x)$ as the initial function.
This exact solution consists of a sequence of Riemann solutions at
$\{x_{j+\half}\}$. Using Green's theorem over the rectangle $I_j\times
[t_n,t_{n+1}]$, one obtains
$$\eqalignno{\int_{I_j}P_n(x)dx&-\int_{I_j}u(x,t_{n+1})dx-\cr
\int_{t_n}^{t_n+\tau}f(u(x_{j+\half},t))dt
&+\int_{t_n}^{t_n+\tau}f(u(x_{j-\half},t))dt=0. &(2.11)\cr
}$$
Because the initial value $P_n(x)=v_j^n,\ x\in I_j$, is piecewise constant,
$$u(x_{j+\half},t)
=R((x-x_{j+\half})/(t-t_n);v_j^n,v_{j+1}^n)|_{x=x_{j+\half}}
=R(0;v_j^n,v_{j+1}^n)={\rm constant},\eqno(2.12)$$
provided $\tau\leq {{\theta}\over {\Lambda}}
h$, which ensures that shocks from $x_{j-\half}$ and
$x_{j+1+\half}$ will not reach $x_{j+\half}$. Therefore,
the approximated averaged value of $u(x,t_{n+1})$ over $I_j$ satisfies:
$$\eqalignno{v_j^{n+1}
&={1\over {\Delta x}}\int_{I_j}u(x,t_{n+1})\,dx\cr
&=v_j^n-\lambda[f(R(0;v_j^n,v_{j+1}^n))-f(R(0;v_{j-1}^n,v_j^n))],&(2.13)\cr
}$$
which indeed can be formulated in terms of (2.9):
$$\eqalignno{
v_j^{n+1}&=v_j^n-\lambda
[f_{j+\half}^n-f_{j-\half}^n],\cr
f_{j+\half}^n&=f
(R(0;v_j^n,v_{j+1}^n)),&(2.14)\cr
\lambda&={{\tau}\over h}\leq {{\theta}\over
{\Lambda}},\cr }$$
where the numerical quadrature used here is just the rectangle rule.
One concludes that $\{v_j^{n+1}\}$ is computed from $\{v_j^n\}$
through the use of two types of operations: Riemann Problem Solver (in (2.12))
and Integral Averaging (in (2.13)). The two processes are TV non--increasing
in the case of scalar conservation laws. Therefore, Godunov method is
TV stable (it is actually total variation diminishing).
For a convex entropy $U(u)$ and entropy
flux $F(u)$,
$${1\over {\Delta x}}\int_{I_j}U(u(x,t_{n+1}))dx\leq U(v_j^n)-\lambda
[F(R(0;v_j^n,v_{j+1}^n))-F(R(0;v_{j-1}^n,v_j^n))]\eqno(2.15)$$
because $u(x,t_{n+1})$ is the unique solution of (1.1) with
initial value $P_n(x)$ (Theorem 1.1). Due to Jensen's
inequality for convex functions,
$$\eqalignno{{1\over {\Delta x}}\int_{I_j}U(u(x,t_{n+1}))dx&\geq
U({1\over {\Delta x}}\int_{I_j}u(x,t_{n+1})dx)\cr
&=U(v_j^{n+1})=U_j^{n+1}.&(2.16)\cr
}$$
One thus derives that:
$$U_j^{n+1}\leq U_j^n-\lambda[F_{j+\half}^n-F_{j-\half}^n], \eqno(2.17)$$
which indicates that Godunov method is also entropy consistent.
The Godunov method enjoys the advantages of being TVD and entropy
consistent because $u_j^n(x)$ is constant, which also results in
the method being just first order accurate. On the other hand,
$u_j^n(x)$ is a polynomial in $x$, in methods like MUSCL, PPM, ENO,
thus these methods enjoy the advantage of being highly accurate
but lacking the stability property (such as TV stability)
and entropy consistency property.
TV stability and entropy consistency of those high order schemes are yet
to be further studied [Vila, 1989].
The goal of this paper is to develop a class of highly accurate methods that
use piecewise constant functions for approximations. Because they use
piecewise constant functions for approximations, they enjoy all advantages
the Godunov method does: TV stability and entropy consistency.
Let's first consider approximating a function $f(x)$ using
piecewise constant functions in the next section.
\bigbreak
\centerline{\bf 3. Approximation using Piecewise Constant Functions}\medskip\nobreak
The following theorem states results concerning approximating a function
using piecewise constant function over a partition of an interval.
\proclaim Theorem 3.1. Let $f(x)$ be a function over interval $[x_l,x_r]$.
Given an integer $N$, let
$x_l=x_10$.\cr
}$$
The solution to this problem consists of one rarefaction wave
traveling to the left, one shock wave traveling to the right,
and a contact discontinuity staying in between. The density
profile in the exact solution is monotone decreasing
where there is a small (weak) jump in density
at the contact discontinuity.
Numerical solution at time $T=2$ is reported. Density profiles are shown
in figures 5.1. The number of constants used in approximation is
indicated in the figures where when ONE constant is used for
approximation, it is the first order Godunov method.
\topinsert
\centerline{%
\vbox{\tabskip=0pt \offinterlineskip
\def\tablerule{\noalign{\hrule}}
\halign to 6 true in{
\strut#& % Col 1
\vrule#\tabskip=1em plus2em&\hfil#\hfil& % Col 2 and 3
\vrule#& \hfil#\hfil& % Col 4 and 5
\vrule#\quad& \hfil#\hfil\quad& % Col 6 and 7
\vrule#& \hfil#\%& % Col 8 and 9
\vrule#\tabskip=0pt\cr\tablerule % Col 10
&&\multispan7\hfil Table 5.2. Error Chart for Test Problem
\#2 ($dhx=0.05000$)\hfil&\cr\tablerule
&&\multispan7\hfil in ${\cal L}^{1}$ Norm\hfil&\cr\tablerule
&& \# of Constants used&&
\omit\hidewidth Norm\hidewidth&&
\omit\hidewidth Absolute Error\hidewidth&&
\omit\hidewidth Relative Error\hidewidth&\cr\tablerule
&&1&& 75.84542 && 1.77505 && 2.34035&\cr\tablerule
&&2&& 75.84542 && 1.02258 && 1.34825&\cr\tablerule
&&4&& 75.84542 && 0.71428 && 0.94176&\cr\tablerule
&&64&& 75.84542 && 0.67377 && 0.88834&\cr\tablerule
&&\sl 2nd Order Method
&&\sl 75.84542 &&\sl 1.29140 &&\sl 1.70276&\cr\tablerule
&&\multispan7\hfil in ${\cal L}^{2}$ Norm\hfil&\cr\tablerule
&&1&& 21.77542 && 1.02382 && 4.70173&\cr\tablerule
&&2&& 21.77542 && 0.74414 && 3.41734&\cr\tablerule
&&4&& 21.77542 && 0.64455 && 2.96001&\cr\tablerule
&&64&& 21.77542 && 0.61964 && 2.84561&\cr\tablerule
&&\sl 2nd Order Method
&&\sl 21.77542 &&\sl 0.81321 &&\sl 3.73455&\cr\tablerule
&&\multispan7\hfil in ${\cal L}^{\infty}$ Norm\hfil&\cr\tablerule
&&1&& 10.98673 && 3.35377 && 30.52561&\cr\tablerule
&&2&& 10.98673 && 3.09033 && 28.12782&\cr\tablerule
&&4&& 10.98673 && 3.01107 && 27.40644&\cr\tablerule
&&64&& 10.98673 && 2.93118 && 26.67929&\cr\tablerule
&&\sl 2nd Order Method
&&\sl 10.98673 &&\sl 3.91365 &&\sl 35.62160&\cr\tablerule
}}
}
%
\vbox{
\centerline{%
\epsfxsize=6.5 true cm
\epsfbox{htp1.eps}
\hfill
\epsfxsize=6.5 true cm
\epsfbox{htp2.eps}
}
\centerline{%
\epsfxsize=6.5 true cm
\epsfbox{htp3.eps}
\hfill
\epsfxsize=6.5 true cm
\epsfbox{htp4.eps}
}
\centerline{ Fig. 5.3. Pressure profiles of test \#2 (solid lines are the exact solutions)}
}
\endinsert
The second test problem [Harten, Engquist, Osher and Chakravarthy, 1987]
has initial values:
$$(q,p,\rho)=\cases{(0.698,3.528,0.445), &if\ $x<0$,\cr
(0,0.571,0.5), &if\ $x>0$.\cr
}$$
The solution to this problem consists of one rarefaction wave
traveling to the left, one shock wave traveling to the right,
and a contact discontinuity staying in between.
Different from the previous test problem, the density
profile of this problem has a built--up in the center part
of the solution. Therefore, this test problem provides
a different scenario to test all numerical methods.
Numerical solution at time $T=1.445$ is reported.
Density profiles are shown
in figures 5.2, and pressure profiles are shown in
figures 5.3. The number of constants used in approximation is
indicated in the figures where when ONE constant is used for
approximation, it is the first order Godunov method.
\topinsert
\centerline{%
\vbox{\tabskip=0pt \offinterlineskip
\def\tablerule{\noalign{\hrule}}
\halign to 6 true in{
\strut#& % Col 1
\vrule# \tabskip=1em plus2em
& \hfil#\hfil& % Col 2 and 3
\vrule#& \quad\hfil#\hfil\quad& % Col 4 and 5
\vrule#& \hfil#\hfil& % Col 6 and 7
\vrule#& \hfil#\hfil& % Col 8 and 9
\vrule#\tabskip=0pt\cr\tablerule % Col 10
&&\multispan7\hfil Table 5.3. Time Comparison Chart for Test Problem
\#1\hfil&\cr\tablerule
&&\multispan7\hfil Relative ${\cal L}^1$ error $\leq 1.0\%$\hfil&\cr\tablerule
&& \# of Constants used&&
\omit\hidewidth \# of grids used \hidewidth&&
\omit\hidewidth Ratio $N/\zeta$\hidewidth&&
\omit\hidewidth CPU Time\hidewidth&\cr\tablerule
&&1&& 500 && 1 && 14.0&\cr\tablerule
&&2&& 200 && 1 && 4.1&\cr\tablerule
&&4&& 120 && 1 && 1.7&\cr\tablerule
&&8&& 120 && 1 && 1.7&\cr\tablerule
&&16&& 120 && 16/14 && 2.2&\cr\tablerule
&&64&& 120 && 64/54 && 2.2&\cr\tablerule
&&\sl 2nd Order Method
&&\sl 120 &&\sl N/A &&\sl 2.0 &\cr\tablerule
}}
}
%
\medskip
\centerline{%
\vbox{\tabskip=0pt \offinterlineskip
\def\tablerule{\noalign{\hrule}}
\halign to 6 true in{
\strut#& % Col 1
\vrule# \tabskip=1em plus2em
& \hfil#\hfil& % Col 2 and 3
\vrule#& \quad\hfil#\hfil\quad& % Col 4 and 5
\vrule#& \hfil#\hfil& % Col 6 and 7
\vrule#& \hfil#\hfil& % Col 8 and 9
\vrule#\tabskip=0pt\cr\tablerule % Col 10
&&\multispan7\hfil Table 5.4. Time Comparison Chart for Test Problem
\#1\hfil&\cr\tablerule
&&\multispan7\hfil Relative ${\cal L}^2$ error $\leq 2.0\%$\hfil&\cr\tablerule
&& \# of Constants used&&
\omit\hidewidth \# of grids used \hidewidth&&
\omit\hidewidth Ratio $N/\zeta$\hidewidth&&
\omit\hidewidth CPU Time\hidewidth&\cr\tablerule
&&1&& 410 && 1 && 9.9&\cr\tablerule
&&2&& 140 && 1 && 2.0&\cr\tablerule
&&4&& 100 && 1 && 1.1&\cr\tablerule
&&8&& 100 && 1 && 1.2&\cr\tablerule
&&16&& 100 && 16/14 && 1.4&\cr\tablerule
&&64&& 80 && 64/54 && 0.9&\cr\tablerule
&&\sl 2nd Order Method
&&\sl 120 &&\sl N/A &&\sl 2.0 &\cr\tablerule
}}
}
%
\endinsert
The above profiles confirm that the newly developed methods indeed are able to
produce numerical solutions with better resolutions than those produced by a
first order method. The improvement is most significant in regions where
rarefaction waves exist. These profiles are comparable to those presented in
the paper by Harten, Engquist, Osher
and Chakravarthy (ENO) [1987].
Errors between the numerical solution and the exact solution at the
indicated times are computed in ${\cal L}^1$, ${\cal L}^2$ and
${\cal L}^{\infty}$ norms. They are reported in table 5.1 and 5.2 for the two
test problems. Similar
errors from a particular implementation of a second order Godunov--type
method [Li, 1994] are also reported in the tables for comparison purpose.
It can be seen that the current methods indeed produce numerical solutions
comparable to those produced by higher order methods which is also
apparent from the pictures in fig. 5.1--5.3. The results are much better than
those produced by a first order conservative method.
The above test problems were run on an IBM RS/6000 model 550
computer running AIX which is highly efficient in floating point number
computations. Table 5.3 records the CPU time used by the different
methods which are required to produce numerical solutions with
a relative ${\cal L}^1$ error less than or equal to $1.0\%$,
while Table 5.4 records the CPU time used by the different methods
which are required to produce numerical solutions with a relative
${\cal L}^2$ error less than or equal to $2.0\%$. One concludes from the
two tables that {\bf (a)}, to achieve a pre--determined accuracy,
the new methods are much more efficient than the first order Godunov
method in terms of CPU time and number of grid points used,
and they are comparable to a particular high order
accurate Godunov--type method, though the exact savings in CPU time
are computer--dependent and norm--dependent; {\bf (b)}, using
a pre--chosen number of grid points for numerical approximations,
the new methods produce numerical solutions
with much better accuracy than the one produced by the first
order Godunov method, and they are also comparable to the one
produced by a particular high order Godunov--type method.
A major problem
that hinders the efficiency of the new methods is the restriction
in $(4.12j)$, where a time step must be chosen to be smaller
than the one used in other Godunov--type methods as indicated
in $(2.9c)$. Otherwise the savings would have been much greater.
The author [Li, 1998] is working on an implicit version of those
new methods which will ease up the strict restriction
in $(4.12j)$, resulting in greater savings over other Godunov--type
methods. Implementation of the current new methods to conservation
laws in multiple space dimensions may also be investigated in the future.
\medskip
\noindent {\bf Acknowledgment}
The author thanks the referees for suggesting the inclusion of comparison
on the computing time among different methods which greatly
enhanced the quality of this manuscript.
\bigbreak
\centerline{\bf References } \medskip\nobreak
\item{[1]} Colella, P. and P. R. Woodward (1985). The
Piecewise-Parabolic Method(PPM) for Gas-Dynamics,
{\it J. Comp. Physics}, {54}, 174--201.\par
\item{[2]} Godunov, S. (1959). Finite Difference Method For Numerical
Computation of
Discontinuous Solutions of the Equations of Fluid Dynamics,
{\it Mat. Sbornik}, {47(89)}, Number 3, page 271.\par
\item{[3]} Harten, A. (1983). On the Symmetric Form of Systems of
Conservation Laws with Entropy,
{\it J. Comp. Physics}, {49}, 151--164.\par
\item{[4]} Harten, A., P. D. Lax and B. Van Leer (1983). On Upstream Differencing
and Godunov--type Schemes for Hyperbolic Conservation Laws,
{\it SIAM Review}, {25}, 35--61.\par
\item{[5]} Harten, A. (1984). On a Class of High Resolution Total--Variation--Stable
Finite--Difference Schemes,
{\it SIAM J. Numer. Anal.}, {21}, 1--23.\par
\item{[6]} Harten, A., B. Engquist, S. Osher and S. Chakravarthy (1987).
Uniformly High Order
Accurate Essentially Non-oscillatory Schemes,
III, {\it J. Comp. Physics}, {71}, 231--303.\par
\item{[7]} Lax, P. D. (1973). {\it Hyperbolic Systems of Conservation Laws and the
Mathematical Theory}
{\it of Shock Waves}, Society for Industrial and Applied
Mathematics.\par
\item{[8]} Lax, P. D., B. Wendroff (1960).
Systems of Conservation Laws, {\it Comm. Pure Appl. Math.}, {13}, 217--237.\par
\item{[9]} Li, X. (1994). A Numerical Method for Systems of Hyperbolic Conservation Laws with
Single Stencil Reconstructions, {\it Applied Mathematics and Computation},
{65}, 125--140.\par
\item{[10]} Li, X. (1998). Entropy Consistent, Implicit TVD Methods
with High Order Accuracy for Conservation Laws, {\it working paper}.\par
\item{[11]} Oleinik, O. A. (1957). On Discontinuous Solutions of Nonlinear
Differential Equations, Uspekhi Mat. Nauk, Vol. 12, pp.3--73.
English translation, Amer. Math. Soc. Trans., Ser. 2, No. 26, pp. 95--172.\par
\item{[12]} Sod, G. A. (1978). A Survey of Several
Finite Difference Methods for Systems of Nonlinear
Hyperbolic Conservation
Laws, {\it J. Comp. Physics}, {27}, 1--31.\par
\item{[13]} Van Leer, B. (1979). Towards
the Ultimate Conservative Differences Scheme, V. A Second
Order Sequel to
Godunov's Methods, {\it J. Comp. Physics}, {32}, 101--136.\par
\item{[14]} Vila, J. P. (1989). An Analysis of a Class of Second--Order
Accurate Godunov--type Schemes, {\it SIAM J. Numer. Anal.},
{26}, 830--853.
\bigskip
\noindent Xuefeng Li \hfil\break
Department of Mathematics and Computer Science, Loyola University\hfil\break
6363 St. Charles Avenue, New Orleans, LA 70118, USA.\hfil\break
E-mail address: Li@Loyno.edu
\bye