\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