Journal of Computational Physics 160, 241–282 (2000)
doi:10.1006/jcph.2000.6459, available online at http://www.idealibrary.com on

New High-Resolution Central Schemes
for Nonlinear Conservation Laws andConvection–Diffusion Equations
Alexander Kurganov∗ and Eitan Tadmor†∗ Department of Mathematics, University of Michigan, Ann Arbor, Michigan 48109;
and †Department of Mathematics, UCLA, Los Angeles, California 90095
E-mail: ∗ kurganov@math.lsa.umich.edu, †tadmor@math.ucla.edu
Received April 8, 1999; revised December 8, 1999

Central schemes may serve as universal finite-difference methods for solving nonlinear convection–diffusion equations in the sense that they are not tied to the specificeigenstructure of the problem, and hence can be implemented in a straightforwardmanner as black-box solvers for general conservation laws and related equations governing the spontaneous evolution of large gradient phenomena. The first-order Lax–Friedrichs scheme (P. D. Lax, 1954) is the forerunner for such central schemes. The
central Nessyahu–Tadmor (NT) scheme (H. Nessyahu and E. Tadmor, 1990) offers
higher resolution while retaining the simplicity of the Riemann-solver-free approach.
The numerical viscosity present in these central schemes is of order O((1x)2r /1t).
In the convective regime where 1t ∼ 1x, the improved resolution of the NT scheme
and its generalizations is achieved by lowering the amount of numerical viscositywith increasing r . At the same time, this family of central schemes suffers fromexcessive numerical viscosity when a sufficiently small time step is enforced, e.g.,
due to the presence of degenerate diffusion terms.
In this paper we introduce a new family of central schemes which retain the simplicity of being independent of the eigenstructure of the problem, yet which enjoy
a much smaller numerical viscosity (of the corresponding order O(1x)2r −1 )). In
particular, our new central schemes maintain their high-resolution independent ofO(1/1t), and letting 1t ↓ 0, they admit a particularly simple semi-discrete formulation. The main idea behind the construction of these central schemes is the use of
more precise information of the local propagation speeds. Beyond these CFL related
speeds, no characteristic information is required. As a second ingredient in theirconstruction, these central schemes realize the (nonsmooth part of the) approximate
solution in terms of its cell averages integrated over the Riemann fans of varying
size.
The semi-discrete central scheme is then extended to multidimensional problems,with or without degenerate diffusive terms. Fully discrete versions are obtained with241
0021-9991/00 $35.00
c 2000 by Academic Press
Copyright °All rights of reproduction in any form reserved.

242

KURGANOV AND TADMOR

Runge–Kutta solvers. We prove that a scalar version of our high-resolution central
scheme is nonoscillatory in the sense of satisfying the total-variation diminishingproperty in the one-dimensional case and the maximum principle in two-space dimensions. We conclude with a series of numerical examples, considering convex
and nonconvex problems with and without degenerate diffusion, and scalar and systems of equations in one- and two-space dimensions. Time evolution is carried out
by the third- and fourth-order explicit embedded integration Runge–Kutta methods
recently proposed by A. Medovikov (1998). These numerical studies demonstrate
the remarkable resolution of our new family of central scheme. °c 2000 Academic Press
Key Words: hyperbolic conservation laws; multidimensional systems; degenerate
diffusion; central difference schemes; non-oscillatory time differencing.

CONTENTS
1. Introduction.
2. Central schemes—A brief overview.
3. The fully discrete scheme—One-dimensional setup.
4. The reduction to semi-discrete formulation. 4.1. One-dimensional hyperbolic conservation laws. 4.2. Onedimensional convection-diffusion equations. 4.3. Multidimensional extensions.
5. From semi-discrete back to fully discrete—The general setup.
6. Numerical examples. 6.1. One-dimensional scalar linear hyperbolic equation. 6.2. One-dimensional
scalar hyperbolic conservation laws. 6.3. One-dimensional systems of hyperbolic conservation laws.
6.4. One-dimensional convection–diffusion equations. 6.5. Two-dimensional problems.

1. INTRODUCTION

During the past few decades there has been an enormous amount of activity related to theconstruction of approximate solutions for nonlinear conservation laws,∂∂u(x, t) +f (u(x, t)) = 0,
∂t
∂x

(1.1)

and for the closely related convection–diffusion equations,∂∂∂u(x, t) +f (u(x, t)) =Q[u(x, t), u x (x, t)].
∂t
∂x
∂x

(1.2)

Here, u(x, t) = (u 1 (x, t), . . . , u N (x, t)) is an N -vector of conserved quantities, f (u) is anonlinear convection flux, and Q(u, u x ) is a dissipation flux satisfying the (weak) parabolicity condition ∇s Q(u, s) ≥ 0 ∀u, s. In the general multidimensional case u is an N -vector in
the d-spatial variables x = (x1 , . . . , xd ), with the corresponding fluxes f (u) = ( f 1 , . . . , f d )and Q(u, ∇x u) = (Q 1 , . . . , Q d ).
These equations are of great practical importance since they govern a variety of physicalphenomena that appear in fluid mechanics, astrophysics, groundwater flow, meterology,
semiconductors, and reactive flows. Convection–diffusion equations (1.2) also arise in twophase flow in oil reservoirs, non-Newtonian flows, front propagation, traffic flow, financial
modeling, and several other areas.
In this work we present new second-order central difference approximations to (1.1)
and (1.2). These new schemes can be viewed as modifications of the Nessyahu–Tadmor(NT) scheme [38]. Our schemes enjoy the major advantages of the central schemes over

HIGH-RESOLUTION CENTRAL SCHEMES243

the upwind ones: first, no Riemann solvers are involved, and second—as a result of being
Riemann solver free—their realization and generalization for complicated multidimensional
systems (1.1) and (1.2) are considerably simpler than in the upwind case. At the same time,
the new schemes have a smaller amount of numerical viscosity than the original NT scheme,
and unlike other central schemes, they can be written and efficiently integrated in their semidiscrete form.
We would like to emphasize the importance of semi-discrete formulations for solving “real,” practical problems associated with multidimensional systems (1.1) and (1.2).
Semi-discrete schemes are especially effective when they combine high-resolution, nonoscillatory spatial discretization with high-order, large stepsize ODE solvers for their time
evolution.
The advantage of our semi-discrete scheme is clearly demonstrated later, in Figs. 6.21 and
6.22, where it is compared with the fully discrete NT solution of the degenerate convection–diffusion equation.
µu t + f (u)x =uxp

1 + u 2x

¶
.

(1.3)

x

This model, recently proposed in [28], describes high-gradient phenomena with possible
discontinuous subshock solutions; consult [12, 28] for details. When the NT scheme is used
to resolve these discontinuities, the computed subshocks are smeared due to the larger numerical dissipation which is accumulated ove the small time steps enforced by the restricted
CFL stability condition, 1t ∼ (1x)2 . This situation—of excessive numerical dissipation (oforder O((1x)2r /1t))—is typical for fully discrete central schemes with time steps muchsmaller than the convective CFL limitation. Alternatively, our new central scheme will accumulate less dissipation (of order O(1x)2r −1 ) and hence can be efficiently used with time
steps as small as required.
This paper is organized as follows. In Section 2 we provide a brief description of thecentral differencing approach for hyperbolic conservation laws.
In Section 3 we introduce our new fully discrete second-order central scheme, whichis constructed for systems of one-dimensional hyperbolic conservation laws. The limiting
case, 1t ↓ 0, brings us to the semi-discrete version presented in Subsection 4.1. Here we
prove that our second-order semi-discrete central scheme satisfies the scalar total-variation
diminishing (TVD) property; consult Theorem 4.1 below. In Subsections 4.2 and 4.3 oursemi-discrete scheme is extended, respectively, to one-dimensional convection–diffusion
equations and to multidimensional hyperbolic and (degenerate) parabolic problems.
In Section 5 we return to the fully discrete framework, discussing time discretizationfor our semi-discrete central scheme. Specifically, we use efficient Runge–Kutta ODE
solvers to integrate the semi-discrete schemes outlined earlier in Section 4. We retain theoverall Riemann-free simplicity without giving up high resolution. Here we prove thatthe resulting second-order fully discrete central scheme satisfies the scalar the maximum
principle; consult Theorem 5.1.
We end in Section 6 by presenting a number of numerical results. These results areconvincing illustrations that our new central schemes provide high resolution at a lowest cost, when applied both to hyperbolic systems of conservation laws and to a varietyof convection–diffusion models. These numerical results confirm an essential aspect of
the current approach—retaining high-resolution without the costly (approximate) Riemann

244

KURGANOV AND TADMOR

solvers, characteritic decompositions, etc. This aspect in the context of high-resolution
schemes was introduced in the Nessyahu–Tadmor scheme [38] and was extended to twodimensional problems in [18]. A nonstaggered and hence less dissipative version was presented in [17]. The relaxation scheme introduced in [19] is closely related to these staggered
central schemes; in fact, they coicide in the relaxation limit ε ↓ 0. The choice for a relaxation matrix A in [19] provides us a family of high-resolution schemes; we note in passingthat the special scalar choice A = ρ(∂ f (u)/∂u)I is an O(ε) perturbation of the central
scheme discussed in this paper. Other componentwise approaches were presented in [34].
The CUSP scheme presented in [16] is a semi-discrete scheme which avoids characteristic
decompositions. And more recently, Liu and Osher [35] introduced a semi-discrete scheme
based on a pointwise formulation of ENO which retains high resolution without Riemannsolvers.
2. CENTRAL SCHEMES—A BRIEF OVERVIEW

Central schemes offer universal finite-difference methods for solving hyperbolic conservation laws, in the sense that they are not tied up to the specific eigenstructure of the
problem and hence can be implemented in a straightforward manner as a black-box solver
for general systems (1.1). In particular, they do not involve characteristic decomposition of
the flux f . In fact, even computation of the Jacobian of f can be avoided; in the particularcase of the second-order NT scheme, for example, numerical derivatives of the flux in (2.6)
below can be implemented componentwise—consult [18, 36].
In 1954 Lax and Friedrichs [10, 29] introduced the first-order stable central scheme, thecelebrated Lax–Friedrichs (LxF) scheme:
=
u n+1
j

u nj+1 + u nj−1
2

−

¡
¢¤
λ£ ¡ n ¢f u j+1 − f u nj−1 .
2

(2.1)

Here, λ := 1t/1x is the fixed mesh ratio, and u nj is an approximate value of u(x = x j , t =t n ) at the grid point (x j := j1x, t n := n1t). Compared with the canonical first-order upwind
scheme of Godunov [11], the central LxF scheme has the advantage of simplicity, sinceno (approximate) Riemann solvers, e.g., [42], are involved in its construction. The main
disadvantage of the LxF scheme, however, lies in its large numerical dissipation, which
prevents sharp resolution of shock discontinuities and rarefaction tips.
A natural high-order extension of the LxF scheme—the NT scheme—was presented in1990 in [38]. The main idea of this generalization is replacing the first-order piecewise
constant solution which is behind the original LxF scheme with van Leer’s MUSCL-typepiecewise-linear second-order approximation, e.g., [31]. This is then combined with a LxF
solver—an alternative to the upwind solvers, which avoids the time-consuming resolutionof Riemann fans by staggered (x, t)-integration. Thus, the NT scheme retains the simplicityof the Riemann-free LxF framework while gaining high resolution, which eliminates the
disadvantage of excessive first-order dissipation.
Here is a brief readers’ digest based on the representation of the LxF and NT schemes
as Godunov-type schemes. We follow [36, Section 2]. To this end, we utilize the sliding
average of u(·, t),
½
¾
Z
1x
1
u(ξ, t) d ξ ,Ix := ξ : |ξ − x| ≤,ū(x, t) :=
|Ix |
2
Ix

245

HIGH-RESOLUTION CENTRAL SCHEMES

so that the integration of (1.1) over the rectangle I x × [t, t + 1t] yields an equivalent
reformulation of the conservation law, (1.1),
 t+1t

t+1t
¶¶
¶¶
µ µ
Z
Z µ µ
1x
1x
1 
, τ dτ −, τ dτ .
f u x+
f u x−ū(x, t + 1t) = ū(x, t) −
1x
2
2
τ =t

τ =t

(2.2)
We begin by assuming that we have already computed an approximation to the solution attime level t = t n —a piecewise linear approximation ũ(x, t n ) ≈ u(x, t n ) of the formũ(x, t n ) :=X£¤ū nj + (u x )nj (x − x j ) 1[x j−1/2 ,x j+1/2 ] ,x j ± 1/2 := x j ±

j

1x
.
2

(2.3)

R
Here, {ū nj } are the computed cell averages, ū nj ≈ ū(x j , t n ) = Ix u(ξ, t n ) dξ/1x, and {(u x )nj }jare approximations to the exact derivatives, u x (x j , t n ). These approximate derivatives arereconstructed from the computed cell averages. The nonoscillatory behavior of the central
schemes hinges on the appropriate choice of approximate derivatives, and there is a library ofrecipes for such nonoscillatory reconstructions. A total-variation (TV) stability, a maximum
principle, or a weaker nonoscillatory property of this piecewise-linear approximation (e.g.,decreasing the number of extrema [36, Section 4]) can be satisfied for a wide variety ofsuch scalar reconstructions proposed and discussed in [4, 13, 14, 27, 32, 36, 38, 41]. For
example, a scalar TVD reconstruction in (2.3) is obtained via the ubiquitous minmod limiter[13, 31, 41],¶µ n
u j − u nj−1 u nj+1 − u nj
,
,
(2.4)
(u x )nj = minmod
1x1xwith minmod(a, b) := 12 [sgn(a) + sgn(b)] · min(|a|, |b|). This is a particular case of a oneparameter family of limiters outlined in (5.2) below.
We then proceed to solve Eq. (2.2) subject to the piecewise-linear initial data (2.3) depicted
in Fig. 2.1.
The piecewise-linear interpolant, ũ(x, t n ), may be discontinuous at points {x j+1/2 }. Yetfor sufficiently small 1t, the solution of problem (2.2)–(2.3) will remain smooth around x j
for t ≤ t n + 1t =: t n+1 , due to the finite speed of propagation. Hence, if we take Ix to be thestaggered grid cell, [x j , x j+1 ] (see Fig. 2.1), we can compute ū(x, t) on the RHS of (2.2)
exactly, and the flux integrals there can be approximated by the midpoint rule. This resultsin the NT scheme [38]ū n+1
j+1/2 =ū nj + ū nj+1
2+¢
£ ¡ n+1/2 ¢¡ n+1/2 ¢¤1x ¡,(u x )nj − (u x )nj+1 − λ f u j+1 − f u j
8n+1/2

where the midpoint values, u j(2.5), are predicted by Taylor expansion,
n+1/2

uj

= ū nj −1t( f x )nj .
2

(2.6)

Thus, we have computed an approximate solution at the next time level t = t n+1 , a solution
which is realized by its (staggered) cell averages, ū n+1
j+1/2 .
Extensions of the second order central NT scheme (2.5), (2.6) to higher-order central
schemes can be found in [4, 32, 36]. Multidimensional extensions were introduced in246

KURGANOV AND TADMOR

FIG. 2.1. Central differencing approach—staggered integration over the local Riemann fan.

[2, 3, 18]. We also would like to mention the corresponding central schemes for incompressible flows in [24–26, 33] and applications to various systems in, e.g., [1, 7, 43]. Thus,
the use of higher-order reconstructions enables us to decrease the numerical dissipation
present in central schemes, and achieve a higher resolution of shocks, rarefactions, andother spontaneous evolution of large gradient phenomena.
Remarks. 1. Characteristic vs componentwise approach. A key advantage of central
schemes is their simplicity–one avoids here the intricate and time-consuming characteristic decompositions based on (approximate) Riemann solvers, which are necessary inhigh-resolution upwind formulations. For systems of conservation laws, the numericalderivatives (u x )nj can be implemented by componentwise extension of the scalar recipe fornonoscillatory limiters. Similarly, the predicted values in (2.6) are based on approximate
derivatives of the flux, ( f x )nj . These values can be computed in terms of the exact Jacobian,∂f(ū nj )(u x )nj . Alternatively, we can even avoid the use of the computationally expensive∂u(and sometimes inaccessible) exact Jacobian ∂∂uf . Instead, the approximate flux derivatives,( f x )nj , are computed in a componentwise manner based on the neighboring discrete valuesof f (ū(x j−1 , t n )), f (ū(x j , t n )) and f (ū(x j+1 , t n )). It was pointed out in [36, 18] that thisJacobian-free version of the central scheme does not deteriorate its high resolution.
2. Cell averages vs point values. Note that here one realizes the approximate solution
by its cell averages, ū n+1
j+1/2 . In general, when dealing with first- and second-order schemes,
n+1
the cell averages, ū n+1
j+1/2 , can be identified with the corresponding point values, u j+1/2 ,modulo a negligible second-order term. We therefore from now on omit the bar notation.
(Consult [4, 36], for example, for this distinction with higher-order central schemes).
3. Second- vs first-order. In the particular case of (u x )nj ≡ 0, the second-order NT
scheme is reduced to the staggered form of the first-order LxF scheme. The nonstaggeredversion of a second-order central scheme can be found in [17].

HIGH-RESOLUTION CENTRAL SCHEMES

247

The second-order NT scheme and its extensions owe their superior resolution to the lower
amount of numerical dissipation—considerably lower than in the first-order LxF scheme.
The dissipation present in these central scheme has an amplitude of order O((1x)2r /1t);
unfortunately, this does not circumvent the difficulties with small time steps which arise, e.g.,with convection–diffusion equations (1.2). Consider, for example, the degenerate parabolic
equation (1.3). Stability necessitates small time steps, 1t ∼ (1x)2 , and the influence of thenumerical dissipation, accumulated over the many steps of the NT scheme, can be clearlyseen in the smeared subshock computed in Fig. 6.22.
One possible way to overcome this difficulty is to use a semi-discrete formulation: whena semi-discrete scheme is coupled with an appropriate ODE solver, one ends up with small
numerical viscosity proportional to the vanishing size of the time step 1t. But in this context,
the central LxF scheme, NT scheme, and their extensions are of limited use, since theseschemes do not admit a semi-discrete form. To make our point, consider the LxF scheme(2.1) in its viscous form,− u nj
u n+1
j
1t

+

¡
¢
¡
¢
f u nj+1 − f u nj−1
21x

=

¢ ¡
¢¤
1 £¡ nu j+1 − u nj − u nj − u nj−1 . (2.7)
21t

Passing to the limit 1t → 0 (while leaving 1x to be fixed), we get the semi-discrete divergence on the left of (2.7), u̇ j (t) + { f (u j+1 (t)) − f (u j−1 (t))}/21x, which is balanced withan increasing amount of dissipation on the right ∼u x x (1x)2 /1t ↑ ∞, as we refine the time
step 1t ↓ 0. In the degenerate viscous case, for example, the CFL restriction 1t ¿ 1x is
responsible for the excessive smearing in the LxF scheme. The second-order NT scheme
has a considerably smaller numerical viscosity, with amplitude of order O((1x)4 /1t) awayfrom extrema cells.1 Nevertheless, the central NT scheme and its higher-order generalizations do not admit any semi-discrete versions, and hence are inappropriate for small time
step computations or steady-state calculations as t ↑ ∞.
This brings us to the new class of central schemes introduced in this paper. These newcentral schemes have smaller numerical dissipation and are the first fully discrete Godunovtype central schemes that admit a semi-discrete form.

3. THE FULLY DISCRETE SCHEME—ONE-DIMENSIONAL SETUPThe NT scheme is based on averaging over the nonsmooth Riemann fans using spatial
cells of the fixed width, 1x. The main idea in the construction of our new central schemes
is to use more precise information about the local speed of wave propagation, in order toaverage the nonsmooth parts of the computed solution over smaller cells of variable size oforder O(1t). We proceed as follows.
Assume that we have already computed the piecewise-linear solution at time level t n ,based on the cell averages u nj , and have reconstructed approximate derivatives (u x )nj in (2.3).
We now turn to evolve it in time. To begin with, we estimate the local speed of propagation at
the cell boundaries, x j+1/2 : the upper bound (disregarding the direction of the propagation)
1

The first two terms on the right of the NT scheme (2.5) yield for smooth u’sū nj − 2ū nj+1/2 + ū nj+1
21t

+

(1x)4
1x
s(u x )nj − (u x )nj+1 d ∼ 1t u x x x x .
81t

248

KURGANOV AND TADMOR

FIG. 3.2. Modified central differencing.

is denoted by a nj+1/2 and is given by2
¶
∂f
(u) ,
max
ρ
∂u
u∈C(u −j+1/2 ,u +j+1/2 )µa nj+1/2 =(3.1)

where u +j+1/2 := u nj+1 − 1x
(u x )nj+1 and u −j+1/2 := u nj + 1x
(u x )nj are the correspondent left22−n
and right intermediate values of ũ(x, t ) at x j+1/2 , and C(u j+1/2 , u +j+1/2 ) is a curve in phasespace connecting u −j+1/2 and u +j+1/2 via the Riemann fan.
Remark. In most practical applications, these local maximal speeds can be easily evaluated. For example, in the genuinely nonlinear or linearly degenerate case one finds that
(3.1) reduces to
½ µ
¶ µ
¶¾
∂f ¡ − ¢∂f ¡ + ¢.
(3.2)
u j+1/2 , ρ
u j+1/2
a nj+1/2 := max ρ∂u∂u
In fact, the maximal local speeds are related to the already calculated CFL number. We
emphasize that these local speeds are the only additional information required to modify
the NT scheme.
Our new scheme is constructed in two steps. First, we proceed along the lines of the
NT scheme. The NT scheme is based on averaging over the staggered control volumes[x j , x j+1 ] × [t n , t n+1 ] of fixed spatial width 1x. Instead, we now use narrower control volumes, where at each time step we integrate over the intervals [x nj+1/2,l , x nj+1/2,r ] × [t n , t n+1 ];see Fig. 3.2. Due to the finite speed of propagation, the points x nj+1/2,l := x j+1/2 − a nj+1/2, 1t
and x nj+1/2,r := x j+1/2 + a nj+1/2 1t separate between smooth and nonsmooth regions, and2Let λi (A) be the eigenvalues of A; then we use ρ(A) := maxi |λi (A)| to denote its spectral radius.

249

HIGH-RESOLUTION CENTRAL SCHEMES

hence the nonsmooth parts of the solution are contained inside these narrower control
volumes of spatial width 2a nj+1/2 1t.
We proceed with the exact evaluation of the new cell averages at t n+1 . Let 1x j+1/2 :=nx j+1/2,r − x nj+1/2,l denote the width of the Riemann fan which originates at x j+1/2 . Exactcomputation of the spatial integrals yields
x nj+1/2,rZ

1
1x j+1/2u(ξ, t n+1 ) dξx nj+1/2,lx nj+1/2,r

=

=

Z

1

ũ(ξ, t n ) dξ −1x j+1/2

x nj+1/2,lu nj + u nj+1
2
−

+

1
2a nj+1/2 1t

1
1x j+1/2

Zt n+1¢¢¡ ¡¢¢¤£ ¡ ¡ ndτf u x j+1/2,r , τ − f u x nj+1/2,l , τtn¢
1x − a nj+1/2 1t ¡(u x )nj − (u x )nj+1
4
Zt n+1¢¢¡ ¡¢¢¤£ ¡ ¡ ndτ.
f u x j+1/2,r , τ − f u x nj+1/2,l , τ

(3.3)

tn

Similarly, let 1x j := x nj+1/2,l − x nj−1/2,r = 1x − 1t (a nj−1/2 + a nj+1/2 ) denote the width of
strip around x j which is free of the neighboring Riemann fans. Then exact integration yields
x nj+1/2,l

1
1x j

Z

u(ξ, t n+1 ) dξx nj−1/2,rx nj+1/2,l

1
=
1x j

Z

x nj−1/2,r1ũ(ξ, t n ) dξ −1x j

Zt n+1¢¢¡ ¡¢¢¤£ ¡ ¡ ndτf u x j+1/2,l , τ dτ − f u x nj−1/2,r , τtn¢
1t ¡ na j−1/2 − a nj+1/2 (u x )nj
4
Zt n+1¢¢¡ ¡¢¢¤£ ¡ ¡ n1−dτ.
f u x j+1/2,l , τ − f u x nj−1/2,r , τ
1x j= u nj +

(3.4)

tn

Using the midpoint rule to approximate the flux integrals on the RHS of (3.3) and (3.4), we
conclude with the new cell averages at t = t n+1 ,w n+1
j+1/2 =

u nj + u nj+1
2−

1
2a nj+1/2

1t= u nj +
w n+1
j
2

−

+

1x − a nj+1/2 1t ¡(u x )nj − (u x )nj+1¢

4£ ¡ n+1/2 ¢¡ n+1/2 ¢¤f u j+1/2,r − f u j+1/2,l ,¡ n¢a j−1/2 − a nj+1/2 (u x )nj£ ¡¢¡ n+1/2 ¢¤λ¡¢ f u n+1/2
j+1/2,l − f u j−1/2,r .
1 − λ a nj−1/2 + a nj+1/2

(3.5)

250

KURGANOV AND TADMOR

Here, the midpoint values are obtained from the corresponding Taylor expansions,¢
1t ¡ nn+1/2
f u j+1/2,l x ,u j+1/2,l := u nj+1/2,l −
2

n+1/2

u j+1/2,r := u nj+1/2,r −¢1t ¡ nf u j+1/2,r x ,2µ
u nj+1/2,l := u nj + 1x(u x )nj¶
1n− λa j+1/2 ,2µu nj+1/2,r := u nj+1 − 1x(u x )nj+1

(3.6)
¶
1
n
− λa j+1/2 .
2

Again, if we want to, we can avoid the computation of the Jacobian of f while using acomponentwise evaluation of f x on the right of (3.6).
At this stage, we realize the solution at time level t = t n+1 in terms of the approximate cell
n+1
averages, w n+1
j+1/2 , w j . These averages spread over a nonuniform grid which is oversampled
by twice the number of the original cells at t = t n . In the second and final step of theconstruction of our scheme, we convert these nonuniform averages back into the original
grid we started with at t = t n , along the lines of the conversion recipe outlined in [17]. As
a by-product of this conversion, we avoid the staggered form of the original NT scheme(2.5)–(2.6).
To obtain the cell averages over the original grid of the uniform, nonstaggered cells[x j−1/2 , x j+1/2 ], we consider the piecewise-linear reconstruction over the nonuniform cells
at t = t n+1 , and following [17], we project its averages back onto the original uniform grid.
Note that we do not need to reconstruct the average of the smooth portion of the solution,
wn+1
j , as it will be averaged out (consult Fig. 3.2), and hence the required piecewise-linear
approximation takes the formw̃(x, t n+1 ) :=

X ©£

¡
¢¤
n+1
n
n
w n+1
j+1/2 + (u x ) j+1/2 x − x j+1/2 1[x j+1/2,l ,x j+1/2,r ]

j

ª
n
n
+ wn+1
j 1[x j−1/2,r ,x j+1/2,l ] .

(3.7)

Here, the exact spatial derivatives, u x (x j+1/2 , t n+1 ), are approximated by(u x )n+1
j+1/2 =

!
Ã
n+1
n+1
w n+1
wn+1
2
j+1 − w j+1/2
j+1/2 − w j¡¡¢,¢ . (3.8)· minmod1x
1 + λ a nj+1/2 − a nj+3/2 1 + λ a nj+1/2 − a nj−1/2

Finally, the desired cell averages, u n+1
j , are obtained by averaging the approximate solutionin (3.7). Our fully discrete second-order central scheme then recasts into the final form

=
u n+1
j

1
1x

xZj+1/2£¡ n¢¤ n+1nw̃(ξ, t n+1 ) dξ = λa nj−1/2 w n+1
wj
j−1/2 + 1 − λ a j−1/2 + a j+1/2

x j−1/2+ λa nj+1/2 wn+1
j+1/2 +

¡ n
¢2
¢2
¤
1x £¡ nn+1
λa j−1/2 (u x )n+1
j−1/2 − λa j+1/2 (u x ) j+1/2 ,
2

(3.9)

n+1
specified in (3.5) are expressed in terms
where the intermediate values of w n+1
j+1/2 and w j
n+1/2
n+1/2
n
of the local speed, a j±1/2 , the midvalues, u j+1/2,l , u j+1/2,r , and the reconstructed slopes,
n+1(u x ) j±1/2 , given in (3.2), (3.6), and (3.8), respectively.

251

HIGH-RESOLUTION CENTRAL SCHEMES

Remarks. 1. Central differencing. The approach taken here can be still viewed as central differencing in the sense that the Riemann fans are inside the domain of averaging.
Consequently, since no (approximate) Riemann solvers are involved, we retain one of themain advantages of the central schemes—simplicity; no Jacobians and charactristic decompositions are needed. At the same time, treating smooth and nonsmooth regions separately,
we gain smaller numerical viscosity (independent of O(1/1t)). In particular, the resulting central scheme (3.9), (3.5) admits a semi-discrete form which is discussed in the next
section.
2. Nonoscillatory properties. The exact entropy evolution operator associated with thescalar equation satisfies the TVD property, ku(·, t)kBV ≤ ku(·, 0)kBV . The various ingredients in the construction of our central scheme retain this TVD property—the nonoscillatory
reconstruction (with appropriate choice of approximate derivatives), exact evolution, andcell averaging. Thus, the only ingredient that is potentially oscillatory enters when we use
the midpoint quadrature rule for temporal integration of the fluxes, yet this does not seemto violate the overall TVD property of our fully discrete central scheme; see, e.g., the TVD
proof of the original NT scheme in [38].
In the particular semi-discrete case (discussed in Section 4 below), the midpoint rule is“exact,” and the TVD of the semi-discrete version of our scheme follows. A direct proof forthe semi-discrete scalar TVD is outlined in Theorem 4.1 below. Moreover, when this semidiscrete scheme is coupled with appropriate Runge–Kutta solvers, we arrive at fully discrete,second-order, central TVD schemes; consult Section 5 below. A maximum principle forthese schemes in two-space dimensions is outlined in Theorem 5.1.
3. Nonstaggered reconstruction. The piecewise linear reconstruction, (3.7), is necessary in order to ensure second-order accuracy, since simple averaging (without reconstruction) over [x j−1/2 , x j+1/2 ] reduces the order of the resulting scheme to first-order accuracy;see [17].
4. First-order version. We conclude by commenting on the first-order version of ourscheme. To this end, we set the slopes, both (u x )nj and (u x )n+1
j+1/2 , to be zero. Then the
staggered cell averages in (3.5) are reduced tow n+1
j+1/2 =

u nj + u nj+1
2−

1
2a nj+1/2£ ¡ n ¢
¡ ¢¤f u j+1 − f u nj ,= u nj .
wn+1
jInserting these values into (3.9) yields the first-order scheme, which takes the viscous
form
u n+1
= u nj −j¡¢¤ 1 £ n¡¢¡¢¤λ£ ¡ n ¢f u j+1 − f u nj−1 +λa j+1/2 u nj+1 − u nj − λa nj−1/2 u nj − u nj−1 ,2
2(3.10)

where a nj+1/2 are the maximal local speeds. This scheme was originally attributed to Rusanov[30]; it is a special case of the family of first-order Godunov-type scheme introduced in [15](based on a symmetric approximate Rieman solver) and it coincides with the so-called local
LxF scheme in [45]. It should be noted, however, that although this scheme is similar to the
LxF scheme, (2.7), its numerical viscosity coefficient, Q nj+1/2 = λa nj+1/2 , is always smaller
than the corresponding LxF one, Q LxF
j+1/2 ≡ 1, thanks to the CFL condition. In regions witha small local speed of propagation (e.g., near the sonic points), the numerical viscosity
present in (3.10) is in fact considerably smaller, 1ta nj+1/2 ¿ 1x.

252

KURGANOV AND TADMOR

4. THE REDUCTION TO SEMI-DISCRETE FORMULATION

4.1. One-Dimensional Hyperbolic Conservation Laws
We consider the fully discrete second-order central scheme (3.9), (3.5), expressed in
n+1
n+1
terms of the cell averages, wn+1
j±1/2 , w j , and the approximate derivates, (u x ) j±1/2 . Observethat except the original averages, u nj , which participates in wnj , all the terms on the right of(3.9), (3.5) are proportional to 1t (or λ). Rearranging the divided differences accordingly
while separating the vanishing terms proportional to λ (as 1t ↓ 0) we find− u nju n+1
j1tn(3.9) a j−1/2
↓

=

µ
wn+1
j−1/2 +¶
a nj−1/2 + a nj+1/2
a nj+1/2 n+1
1
1 n
−
w n+1
w j+1/2 −u + O(λ)
+
j
1x
1x
1x
1t j

1x
½ n
¢ 1
¡
¢
1 £ ¡ n+1/2 ¢ ¡ n+1/2 ¢¤(3.5) a j−1/2 ¡ n
u j−1 + u nj + a nj−1/2 (u x )nj−1 −(u x )nj −f u j−1/2,r − f u j−1/2,l
↓
= 21x
4
21x
¢
¡ n
¢
a j−1/2 + a nj+1/2 n 1 ¡ n
1 £ ¡ n+1/2 ¢−u j + a j−1/2 − a nj+1/2 (u x )nj −f u j+1/2,l1x2
1x¡ n+1/2 ¢¤ a nj+1/2 ¡ n¢ 1¡¢− f u j−1/2,r +u j + u nj+1 + a nj+1/2 (u x )nj − (u x )nj+1
21x
4¾¡ n+1/2 ¢¤1 £ ¡ n+1/2 ¢f u j+1/2,r − f u j+1/2,l+ O(λ)−
21x(£¡ ¡ n+1/2 ¢¡ n+1/2 ¢¢ ¡ ¡ n+1/2 ¢¡ n+1/2 ¢¢¤1− f u j+1/2,r + f u j+1/2,l − f u j−1/2,r + f u j−1/2,l
=
21x
·µ
¶ µ
¶¸
a nj+1/2
1x
1x+
u nj+1 −(u x )nj+1 − u nj +(u x )nj
1x
2
2
·µ
¶
µ
¶¸ )
a nj−1/2
1x
1x
n
x
n
n
uj −
(u x ) j − u j−1 +(u x ) j−1 . + O(λ).
−
1x22Note that as 1t → 0, the midvalues on the right approach (consult (3.6))
1x(u x ) j+1 (t) =: u +j+1/2 (t),
2
1x
n+1/2
(u x ) j (t) =: u −j+1/2 (t),u j+1/2,l → u j (t) +2
n+1/2

u j+1/2,r → u j+1 (t) −(4.1)

where (u x ) j (t) are the numerical derivatives reconstructed from the computed cell averages,u j (t). Thus, letting 1t ↓ 0, the resulting semi-discrete central scheme can be written in its
compact form¢¡¢¢ ¡ ¡
¢ ¡ ¡
¢¢
¡ ¡ +
f u j+1/2 (t) + f u −j+1/2 (t) − f u +j−1/2 (t) + f u −j−1/2 (t)du j (t) = −
dt
21x
£ +
¤
£
¤ª
1 ©
−
a j+1/2 (t) u j+1/2 (t) − u j+1/2 (t) − a j−1/2 (t) u +j−1/2 (t) − u −j−1/2 (t) .
+
21x
(4.2)

253

HIGH-RESOLUTION CENTRAL SCHEMES

Recall that a j+1/2 (t) is the maximal local speed, e.g., in the generic case one may take
½ µ
¶ µ
¶¾
¢
¢
∂f ¡ +∂f ¡ −uu(t) , ρ(t).
a j+1/2 (t) := max ρ∂u j+1/2∂u j+1/2
Remarks.
tive form,

1. Conservation form. The second-order scheme, (4.2), admits the conservaH j+1/2 (t) − H j−1/2 (t)du j (t) = −
,
dt
1x

(4.3)

with the numerical flux
¢
¡
¢
¡
¤
f u +j+1/2 (t) + f u −j+1/2 (t)a j+1/2 (t) £ +−u j+1/2 (t) − u −j+1/2 (t) . (4.4)
H j+1/2 (t) :=22Here, the intermediate values u ±j+1/2 are given by
u +j+1/2 := u j+1 (t) −1x(u x ) j+1 (t),2u −j+1/2 := u j (t) +1x(u x ) j (t).
2

(4.5)

One verifies that H j+1/2 (t) ≡ H (u j−1 (t), u j (t), u j+1 (t), u j+2 (t)) is a numerical flux consistent with Eq. (1.1), i.e., H (v, v, v, v) = f (v). In fact, with the minmod limiter, (2.4),
the corresponding approximate derivatives, (u x ) j (t), vanish at extrema values,sgn(u j+1 (t) − u j (t)) + sgn(u j (t) − u j−1 (t)) = 0 ⇒ (u x ) j (t) = 0,

(4.6)

and hence the corresponding numerical flux satisfies the essentially three-point consistency
H (·, v, v, ·) = f (v).

(4.7)

2. Numerical viscosity. The second expression on the right of (4.2) accounts for thenumerical viscosity of the scheme. Taylor’s expansion shows that for smooth u’s, this amountof numerical viscoity is of order ∼(1x)3 (a(u)u x x x )x /8. This O(1x)3 term, uniformly
bounded w.r.t. 1/1t, should be contrasted with the corresponding numerical viscosity termsof order O((1x)2 /1t) in the first-order LxF scheme (indicated earlier in (2.7)) and of orderO((1x)4 /1t) in the second-order NT scheme.
3. Simplicity. We again would like to emphasize the simplicity of the second-ordersemi-discrete central scheme, (4.2), so that it does not require any information about the
eigenstructure of the underlying problem beyond the CFL-related speeds, a j+1/2 (t). The
computation of the numerical derivatives, (u x ) j (t), is carried out componentwise; no specificknowledge of characteristic decomposition based on (approximate) Riemann solvers is
required.
4. First-order reduction. If we reset all the numerical derivatives, (u x ) j (t) = 0, then
(4.2) is reduced to the first-order semi-discrete central scheme corresponding to the 1t ↓0-limit of the Rusanov scheme, (3.10),
f (u j+1 (t)) − f (u j−1 )du j (t) = −
dt
21x
¤
1 £
+
a j+1/2 (t)(u j+1 (t) − u j (t)) − a j−1/2 (t)(u j (t) − u j−1 (t)) . (4.8)
21x
We conclude this section with the proof of the one-dimensional scalar TVD property forour new semi-discrete scheme.

254

KURGANOV AND TADMOR

THEOREM 4.1 (TVD of Semi-discrete Central Scheme). Consider the scalar semidiscrete central scheme (4.2), with intermediate values u ±j±1/2 in (4.5) based on approximate
derivatives (u x ) j (t) satisfying (4.6), e.g., the family of minmod-like limiters (with θ = 1 corresponding to (2.4))¶µ n
u j − u nj−1 u nj+1 − u nj−1 u nj+1 − u nj,,θ,1 ≤ θ ≤ 2. (4.9)(u x )nj := minmod θ
1x
21x
1x
Then the following TVD property holds:X|u j+1 (t) − u j (t)| ≤ ku(·, 0)kBV .
ku(·, t)kBV :=jRemark. Notice that in the scalar case the local propagation speeds are given by
a j+1/2 (t) :=maxu∈[u −j+1/2 (t),u +j+1/2 (t)]| f 0 (u)|,(4.10)and in the special case of convex f , this is further simplified:
a j+1/2 (t) := max{| f 0 (u −j+1/2 (t))|, | f 0 (u +j+1/2 (t))|}.

(4.11)

Proof. The second-order flux in (4.2), H j+1/2 (t), can be viewed as a generalized MUSCL
flux [40],¢¡H j+1/2 (t) = H Rus u +j+1/2 (t), u −j+1/2 (t) ,expressed in terms of the first-order E-flux H Rus = H Rus (u ` , u r ), associated with the firstorder Rusanov scheme (4.8),
H Rus (u ` , u r ) :=f (u ` ) + f (u r ) a`r−(u r − u ` ),22a`r = max | f 0 (u)|.
u∈[u ` ,u r ]According to [47], the TVD property of such scalar, semi-discrete generalized MUSCL
schemes is guaranteed if (consult [47, Example 2.4])¯¯¯ (u x ) j+1/2 ¯¯ ≤ 2.
¯¯¯ 1u(4.12)

j±1/2

This is clearly fulfilled by the choice of approximate derivatives in (4.9). (We note in passing
the necessity of the clipping phenomenon, (4.6), enforced by (4.12).) j
4.2. One-Dimensional Convection–Diffusion Equations
Consider the convection–diffusion equation (1.2). If the dissipation flux, Q(u, u x ), is anonlinear function, then Eq. (1.2) can be a strongly degenerate parabolic equation which
admits nonsmooth solutions. To solve it numerically is a highly challenging problem. In
this context, the operator splitting technique was used in, e.g., [6, 8, 12, 20, 22, 23], yet thisapproach suffers the familiar limitations of splitting, e.g., limited accuracy, etc.

255

HIGH-RESOLUTION CENTRAL SCHEMES

Our second-order semi-discrete scheme, (4.3)–(4.4), can be applied to Eq. (1.2) in astraightforward manner, since we can treat the hyperbolic and the parabolic parts of (1.2)
simultaneously. This results in the following conservative scheme:u̇ j (t) = −H j+1/2 (t) − H j−1/2 (t) P j+1/2 (t) − P j−1/2 (t)
+
.
1x
1x

(4.13)

Here, H j+1/2 (t) is our numerical convection flux, (4.4), and P j+1/2 (t) is a reasonable approximation to the diffusion flux, e.g., the simplest central difference approximation
· µ
¶
µ
¶¸
1
u j+1 (t) − u j (t)u j+1 (t) − u j (t)Q u j (t),+ Q u j+1 (t),. (4.14)
P j+1/2 (t) =21x1x4.3. Multidimensional Extensions
Our second-order semi-discrete schemes, (4.2) and (4.13), (4.4), (4.14), can be extended
to both multidimensional hyperbolic and parabolic problems. Without loss of generality, let
us consider the two-dimensional convection–diffusion equationu t + f (u)x + g(u) y = Q x (u, u x , u y )x + Q y (u, u x , u y ) y ,(4.15)

where Q x ≡ Q y ≡ 0 corresponds to the two-dimensional hyperbolic conservation law.
We use a uniform spatial grid, (x j , yk ) = ( j1x, k1y). Suppose that we have computed
the solution at some time level t and have reconstructed the two-dimensional, non-oscillatory
piecewise-linear polynomial approximationX[u j,k (t)+(u x ) j,k (t)(x − x j )+(u y ) j,k (t)(y − yk )]1[x j−1/2 ,x j+1/2 ] × [yk−1/2 ,yk+1/2 ] .
u(x, y, t) ≈j,k

, yk±1/2 := yk ± 1y
; (u x ) j,k (t) and (u y ) j,k (t) are numerical derivaHere, x j±1/2 := x j ± 1x
2
2
tives, which approximate the exact ones, u x (x j , yk , t) and u y (x j , yk , t), respectively. Witha proper choice of numerical derivatives, the reconstruction of piecewise polynomial approximation is nonoscillatory. For example, using the minmod limiter, (2.4), guarantees
the nonoscillatory property in the sense of satisfying a (local) scalar maximum principle;
consult [18, Thm. 1].
The 2D extension of the scheme (4.13), (4.4), (4.14) can be written in the conservative
form
y

y

x
x
− H j−1/2,k
H j+1/2,k
H j,k+1/2 − H j,k−1/2du j,k (t) = −
−
dt
1x
1y

+

x
x
P j+1/2,k− P j−1/2,k

1x
y

y

y

+

y

P j,k+1/2 − P j,k−1/2
1y.

(4.16)xx≡ H j+1/2,k(t) and H j,k+1/2 ≡ H j,k+1/2 (t) are x- and y-numerical convection
Here, H j+1/2,kfluxes, respectively (viewed as a generalization of the one-dimensional flux constructedabove in (4.4)),
¢
¡
¢
¡
¤
a xj+1/2,k (t) £ +f u +j+1/2,k (t) + f u −j+1/2,k (t)x
H j+1/2,k (t) :=−u j+1/2,k (t) − u −j+1/2,k (t) ,2
2(4.17)
¢
¡ −
¢
¡ +
y
¤
g u j,k+1/2 (t) + g u j,k+1/2 (t)a j,k+1/2 (t) £ +y−u j,k+1/2 (t) − u −j,k+1/2 (t) ,H j,k+1/2 (t) :=
2
2

256

KURGANOV AND TADMOR

which are expressed in terms of the intermediate values
1x(u x ) j+1/2±1/2,k (t),21y(u y ) j,k+1/2±1/2 (t),u ±j,k+1/2 (t) := u j,k+1 (t) ∓2

u ±j+1/2,k (t) := u j+1,k (t) ∓(4.18)

yand the local speeds, a xj+1/2,k (t) and a j,k+1/2 (t), are computed, e.g., by
µ
¶
µ
¶
¢
¢
∂f ¡ ±∂g ¡ ±
ya xj+1/2,k (t) := max ρ
a j,k+1/2 (t) := max ρ
u j+1/2,k (t) ,u j,k+1/2 (t) .
±
±
∂u
∂u
(4.19)yyxxSimilarly, P j+1/2,k≡ P j+1/2,k(t) and P j,k+1/2 ≡ P j,k+1/2 (t) are the corresponding x- andy-numerical diffusion fluxes, given by· µ¶1
u j+1,k − u j,kxx
P j+1/2,k :=Q u j,k ,, (u y ) j,k
2
1x
µ
¶¸
u j+1,k − u j,k, (u y ) j+1,k ,+ Q x u j+1,k ,1x(4.20)· µ¶1
u j,k+1 − u j,kyy
P j,k+1/2 :=Q u j,k , (u x ) j,k ,
2
1y
¶¸
µ
u j,k+1 − u j,k.
+ Q y u j,k+1 , (u x ) j,k+1 ,1y5. FROM SEMI-DISCRETE BACK TO FULLY DISCRETE—THE GENERAL SETUP

The two-dimensional semi-discrete central scheme (4.2) forms a system of nonlinear
ODEs, the so-called “method of lines” for the discrete unknowns {u j,k (t)}. To integrate
in time, one must introduce a variable time step, 1t n , stepping forward from time levelt n to t n+1 := t n + 1t n . We start by considering the simplest scenario of first-order timedifferencing. The nonoscillatory behavior of the forward Euler scheme is summarized in
the maximum principle stated in Theorem 5.1 below. To retain the overall high accuracy ofthe spatial differencing, however, higher-order stable time discretizations are required. To
this end, the forward Euler time differencing can be used as a building block for higherorder Runge–Kutta and multi-level ODE solvers. In particular, second- and third-order ODE
solvers can be constructed by convex combinations of the simple forward Euler differencing,retaining the overall maximum principle. Thus, we conclude in Corollaries 5.1 and 5.2
below with a fully discrete second-order central scheme satisfying the two-dimensional
scalar maximum principle.
We begin with
THEOREM 5.1 (Maximum Principle). Consider the two-dimensional central scheme¢ 1t n ¡ y¢
1t n ¡ xyn
x(t n ) −H j+1/2,k (t n ) − H j−1/2,k
H j,k+1/2 (t n ) − H j,k−1/2 (t n ) .
u n+1
j,k = u j,k −
1x
1y
(5.1)
Here, H x (t) and H y (t) are the numerical fluxes given in (4.17)–(4.19); let their numericalderivatives be determined by one of the following one-parameter family of minmod-like257

HIGH-RESOLUTION CENTRAL SCHEMES

limiters:µ(u x )nj,k := minmodµ(u y )nj,k := minmod¶
u nj,k − u nj−1,k u nj+1,k − u nj−1,k u nj+1,k − u nj,k
,
,θ
,
θ
1x
21x
1x
¶
u nj,k − u nj,k−1 u j,k+1 − u j,k−1 u nj,k+1 − u nj,k
,
,θ
,
θ
1y
21y
1y

1 ≤ θ ≤ 2,
(5.2)
1 ≤ θ ≤ 2.

Assume the following CFL condition holds:
µ n
¶
1t
1t n
1
max
max| f 0 (u)|,max|g 0 (u)| ≤ .
1x u
1y u
8

(5.3)

Then the resulting fully discrete central scheme satisfies the maximum principle
max{u j,k (t n+1 )} ≤ max{u j,k (t n )}.
j,k

(5.4)

j,k

Proof. With λn := 1t n /1x and µn := 1t n /1y denoting the x- and y-mesh ratios, theforward Euler time discretization of our scheme takes the explicit form
¢
¡
¢
¡
¢
¡
¢
λn £ ¡ +f u j+1/2,k (t n ) + f u −j+1/2,k (t n ) − f u +j−1/2,k (t n ) − f u −j−1/2,k (t n )
2
©
ª
©
x
− a j+1/2,k (t n ) u +j+1/2,k (t n ) − u −j+1/2,k (t n ) + a xj−1/2,k (t n ) u +j−1/2,k (t n )n
u n+1
j,k = u j,k −ª¤ µn £ ¡ +
¢
¡
¢
¡
¢
g u j,k+1/2 (t n ) + g u −j,k+1/2 (t n ) − g u +j,k−1/2 (t n )− u −j−1/2,k (t n ) −
2
¢
©
ª
¡ −
y
n
− g u j,k−1/2 (t ) − a j,k+1/2 (t n ) u +j,k+1/2 (t n ) − u −j,k+1/2 (t n )©
ª¤y+ a j,k−1/2 (t n ) u +j,k−1/2 (t n ) − u −j,k−1/2 (t n ) .
To simplify notations, we use the standard abbreviations¡¢¡¢
1xj,k f := f u −j+1/2,k (t n ) − f u +j−1/2,k (t n ) ,1xj±1/2,k u := u +j±1/2,k (t n )−u −j±1/2,k (t n ),with the similar notations for y-differences, e.g., 1 j,k g := g(u −j,k+1/2 (t n ))−g(u +j,k−1/2 (t n )),etc. Then our scheme can be rewritten as follows (where all the quantities on the right are
taken at time level t = t n ):
y

u n+1
j,k =

u −j+1/2,k + u +j−1/2,k + u −j,k+1/2 + u +j,k−1/2
4

−

·
¢
λn 1xj+1/2,k f ¡ +u j+1/2,k − u −j+1/2,kx
2 1 j+1/2,k u¢ 1xj−1/2,k f ¡ +¢1xj,k f ¡ −u j+1/2,k − u +j−1/2,k + x
u j−1/2,k − u −j−1/2,k+2 x
1 j,k u1 j−1/2,k u
¸

¡
¢
¡
¢
− a xj+1/2,k u +j+1/2,k − u −j+1/2,k + a xj−1/2,k u +j−1/2,k − u −j−1/2,k
−

· y
y
¢
¢
1 j,k g ¡ −µn 1 j,k+1/2 g ¡ +−+2uu j,k+1/2 − u +j,k−1/2−uyyj,k+1/2
j,k+1/2
2 1 j,k+1/2 u
1 j,k uy¢¡¢
1 j,k−1/2 g ¡ +yu j,k−1/2 − u −j,k−1/2 − a j,k+1/2 u +j,k+1/2 − u −j,k+1/2
y
1 j,k−1/2 u¸¡¢
y+ a j,k−1/2 u +j,k−1/2 − u −j,k−1/2 .
+258

KURGANOV AND TADMOR

Rearranging the terms, we find that u n+1
j,k is given by the following linear combination of
the intermediate values, u ±j±1/2,k and u ±j,k±1/2 :
u n+1
j,k =

½
¾
½
¾
1xj+1/2,k f
1xj−1/2,k fλnλna xj+1/2,k − x
u +j+1/2,k +a xj−1/2,k + x
u −j−1/2,k2
1 j+1/2,k u2
1 j−1/2,k u·¸¾½
1xj,k f1 λn 1xj+1/2,k fx+−au −j+1/2,k+−2
j+1/2,k4
2 1xj+1/2,k u
1xj,k u½·¸¾
1xj,k f1 λn 1xj+1/2,k f
x
+
−
2
−
+
a
u +j−1/2,kj−1/2,k4
2 1xj+1/2,k u
1xj,k u½½¾¾
yy
1 j,k+1/2 g +1 j,k−1/2 g −
µn
µn
y
y
+
a j,k+1/2 − ya j,k−1/2 + y
u j,k+1/2 +u j,k−1/2
2
2
1 j,k+1/2 u
1 j,k−1/2 u
½
· y
¸¾
y
1 j,k g
1 µn 1 j,k+1/2 gy++− a j,k+1/2 − 2 yu −j,k+1/2
y4
2 1 j,k+1/2 u
1 j,k u· y¸¾½
y
1 j,k g
1 µn 1 j,k−1/2 gy−+au +j,k−1/2 .
+
−
2
(5.5)
y
y
j,k−1/2
4
2 1 j,k−1/2 u
1 j,k uNote that all the coefficients in (5.5) are positive due to our CFL assumption, (5.3). Thismeans that the linear combination on the RHS of (5.5) is a convex combination and hence±±
the value of u n+1
j,k does not exceed the values of u j±1/2,k and u j,k±1/2 . And since ourchoice of minmod-like approximate derivatives in (5.2) guarantees that these intermediate values, u ± , satisfy a local maximum principle w.r.t. the original averages, u n , e.g., [18],max j,k {u ±j±1/2,k (t), u ±j,k±1/2 (t)} ≤ max j,k {u j,k (t)}, the result (5.4) then follows. j
The forward Euler scheme is limited to first-order accuracy. It can be used, however, asa building block for higher-order schemes based on Runge–Kutta (RK) or multi-level timedifferencing. Shu and Osher [44, 45] have identified a whole family of such schemes, basedon convex combinations of forward Euler steps.
To this end, we let C[w] denote our spatial recipe (4.17)–(4.19) for central differencing
a grid function w = {w j,k },·C[w] := −xx(w) − H j−1/2,k(w)
H j+1/2,k

1x

y

+

y

H j,k+1/2 (w) − H j,k−1/2 (w)

¸

1y

.

(5.6)

Expressed in terms of the forward Euler solver, w + 1tC[w], we consider the one-parameterfamily of RK schemesu (1) = u n + 1t n C[u n ]¡£ ¤¢u (`+1) = η` u n + (1 − η` ) u (`) + 1t n C u (`) ,` = 1, 2, . . . , s − 1,

(5.7)

u n+1 := u (s) .
In Table 5.1 we quote the preferred second- and third-order choices of [45]. We state
COROLLARY 5.1 (Maximum Principle for Runge–Kutta Time Differencing). Assume
that the CFL condition (5.3) holds. Then the fully discrete central scheme (4.17)–(4.19),(5.2),(5.7) with ηm specified in Table 5.1 satisfies the maximum principle (5.4).

259

HIGH-RESOLUTION CENTRAL SCHEMESTABLE 5.1
Runge–Kutta Methodsη1η2Second order time differencing
1
Two-step modified Euler (s = 2)2—

Third order time differencing
3
Three-step method (s = 3)4

1
3

A similar two-parameter family of multi-level methods was identified in [44]. They takethe particularly simple form
u n+1 = η(u n + c0 1t n C[u n ]) + (1 − η)(u n−s + cs 1t n C[u n−s ]),(5.8)

with positive coefficients given in Table 5.2. We state
COROLLARY 5.2 (Maximum Principle for Multi-level Time Differencing). Assume that
the CFL condition
µ n
¶
1t
1t n
1
0
0
max| f (u)|,max|g (u)| ≤ min
max
(5.9)
u
u
c
k
1x
1y
8ck
holds. Then the fully discrete multi-level central scheme (4.17)–(4.19), (5.2), (5.8), with ηand c’s specified in Table 5.2, satisfies the maximum principle (5.4).
We close by noting that Corollaries 5.1 and 5.2 extend the maximum principle for thesecond-order fully discrete two-dimensional scheme introduced in [18, Thm. 1].
6. NUMERICAL EXAMPLES

We conclude this paper with a number of numerical examples. In all the numericalresults presented below we have used the θ -dependent family of limiters corresponding
to (4.9). (These are in general the less dissipative limiter than the original minmod, (2.4),
corresponding to (4.9) with θ = 1). The spatial derivative of u(x, t) is approximated byµ¶ū(x + a) − ū(x) ū(x + a) − ū(x − b) ū(x) − ū(x − b),,θ,u x (x, t) ≈ minmod θ
a
a+b
b
(6.1)
TABLE 5.2
Multi-level Methods
η

c0

cs

Second-order time differencing
34-level method (s = 2)2
4

05-level method (s = 3)

3
2

0

Third-order time differencing
16
5-level method (s = 3)27

3

6-level method (s = 4)

2

12
11
10
7
30
17

7-level method (s = 5)

8
9

25
32
108
125

5
3

260

KURGANOV AND TADMOR

where a and b are appropriate grid scales, and the multivariable minmod function is definedby
if x j > 0 ∀ j,min j {x j },minmod(x1 , x2 , . . .) = max j {x j },if x j < 0 ∀ j,0,otherwise.
The parameter θ ∈ [1, 2] has been chosen in the optimal way in every example. Note thatθ = 2 corresponds to the least dissipative limiter (no new local extrema are introduced),
whereas θ = 1 ensures a nonoscillatory nature of the approximate solution in the sensethat there is no increase of the total-variation. The reconstruction depicted in Fig. 3.2, forexample, does not increase in the variation at the interface x j−1/2 but it is an oscillatory one
since new extrema is introduced at x j+1/2 . In the scalar examples below, θ = 2 has provided
a satisfactory results, but consult the nonconvex Buckley–Levertt equation in Fig. 6.7 as acounterexample; for systems the optimal values of θ vary between 1.1 and 1.5.
Another point we would like to stress concerns the time discretization of our semi-discrete
central schemes. In general, the Runge–Kutta time differencing is preferable over the multilevel differencing, since the former enables a straightforward use of variable time steps.
We note that for the standard explicit RK methods the time step can be very small due totheir strict stability restriction. There are two different approaches to increasing efficiency
at this point. First, one can use implicit or explicit–implicit ODE solvers. These methods areunconditionally stable, but they require inverting nonlinear operators (in the general case ofa nonlinear diffusion), which is a computationally expensive and analytically complicatedprocedure.
In all the numerical examples shown below, we preferred the second approach—to solve
systems of ODEs by means of the explicit embedded integration third-order RK method
recently introduced by Medovikov [37] (his original code, DUMKA3, was used). This highorder differencing produces accurate results, and its larger stability domains (in comparison
with the standard RK methods) allow us to use larger time steps; the explicit form retains
simplicity, and the embedded formulas permit an efficient stepsize control. In practice thesemethods preserve all the advantages of explicit methods and work as fast as implicit methods(see [37] for details).
Remark. Below, we abbreviate by FD1 and SD1 the Rusanov first-order fully andsemi-discrete schemes. We also use FD2 and SD2 notation for our second-order fully andsemi-discrete schemes. As previously, LxF and NT stand for the Lax–Friedrichs and the
Nessyahu–Tadmor schemes.
6.1. One-Dimensional Scalar Linear Hyperbolic Equation
EXAMPLE 1 (Linear Steady Shocks). First, consider the simplest linear equation withthe zero flux, f (u) = 0,
u t = 0,subject to the discontinuous initial data,½1,u(x, 0) =0,−0.5 < x < 0.5,
otherwise.

(6.2)

(6.3)

HIGH-RESOLUTION CENTRAL SCHEMES

261

FIG. 6.3. Problem (6.2)–(6.3). N = 100, T = 2; first-order methods.

Notice that due to their dissipative nature, neither the LxF scheme nor the NT scheme
provides a reasonable approximation to this problem. But using the Rusanov scheme andour second-order fully discrete scheme, which are, in some sense, the least dissipative
central schemes, we achieve the perfect resolution of the discontinuities (the same result isobtained by the corresponding semi-discrete schemes).
In Figs. 6.3 and 6.4 the solutions computed by our schemes are compared with the(staggered) LxF and NT schemes. It also can be easily checked analytically that both the
FD1 scheme and the FD2 schemes solve problem (6.2)–(6.3) exactly.

FIG. 6.4. Problem (6.2)–(6.3). N = 100, T = 2; second-order methods.

262

KURGANOV AND TADMORTABLE 6.1
Initial Value Problem (6.4), L1 - and L∞ -Norms of the Errors
L ∞ -error

L 1 -error
N

NT

Rate

FD2

Rate

NT

Rate

FD2

Rate

40
80
160
320
640
1280

2.920e-03
4.583e-04
1.115e-04
2.360e-05
5.273e-06
1.249e-06

—
2.67
2.04
2.24
2.16
2.08

8.716e-03
1.876e-03
3.892e-04
7.943e-05
1.659e-05
3.430e-06

—
2.22
2.27
2.29
2.26
2.27

3.151e-03
9.963e-04
3.704e-04
1.263e-04
4.463e-05
1.690e-05

—
1.66
1.43
1.55
1.50
1.40

7.818e-03
2.598e-03
9.262e-04
2.881e-04
1.028e-04
3.593e-05

—
1.59
1.49
1.68
1.49
1.52

EXAMPLE 2 (Accuracy Test). Let us consider the linear equation subject to periodic
initial datau t + u x = 0,u(x, 0) = sin x.

(6.4)

This problem admits the global smooth solution that was computed at time T = 1 with thevarying number of grid points, N .
In Table 6.1 we compare the accuracy of our second-order fully discrete scheme, FD2,
with the accuracy of the NT scheme. These results show that for the FD2 scheme the absolute
error is larger, but the rate of convergence is slightly higher than for the NT scheme.

6.2. One-Dimensional Scalar Hyperbolic Conservation Laws
EXAMPLE 3 (Burgers’ Equation: Pre- and Post-shock Solutions). In this example weapproximate solutions to the inviscid Burgers’ equation,
µ 2¶
u
= 0.
(6.5)
ut +
2 x
Let us start with the case of smooth periodic initial data, e.g.,u(x, 0) = 0.5 + sin x.

(6.6)

The well-known solution of (6.5)–(6.6) develops a shock discontinuity at the critical time
Tc = 1. Table 6.2 shows the L 1 - and L ∞ -norms of the errors at the pre-shock time T = 0.5
TABLE 6.2
Initial Value Problem (6.5)–(6.6), L1 - and L∞ -Norms of the Errors at T = 0.5
L ∞ -error

L 1 -error
N

NT

Rate

FD2

Rate

NT

Rate

FD2

Rate

40
80
160
320
640
1280

1.011e-02
2.116e-03
4.705e-04
1.095e-04
2.517e-05
5.926e-06

—
2.26
2.17
2.10
2.12
2.09

9.101e-03
1.843e-03
4.272e-04
9.334e-05
2.163e-05
4.867e-06

—
2.30
2.11
2.19
2.11
2.15

6.782e-03
2.951e-03
9.918e-04
3.727e-04
1.248e-04
4.433e-05

—
1.20
1.57
1.41
1.58
1.49

6.283e-03
2.333e-03
7.481e-04
2.603e-04
9.508e-05
3.132e-05

—
1.43
1.64
1.52
1.45
1.60

HIGH-RESOLUTION CENTRAL SCHEMES

263

FIG. 6.5. Problem (6.5)–(6.6). N = 100, T = 2; the FD2 scheme.

when the solution is still infinitely smooth. Unlike the linear case (Example 2), both the absolute errors and the convergence rates of the FD2 scheme are better than the corresponding
errors and rates of the NT scheme. This indicates a certain advantage of our fully discrete
second-order scheme over the NT scheme while applied to nonlinear problems.
In Figs. 6.5 and 6.6 we present the approximate solutions at the post-shock time T = 2,
when the shock is well developed. Second-order behavior is confirmed by the measuring
Lip´-errors, [39], which are recorded in Table 6.3. Again, the solution obtained by the FD2
scheme is slightly more accurate than the solution computed by the NT scheme.

FIG. 6.6. Problem (6.5)–(6.6). N = 100, T = 2; the NT scheme.

264

KURGANOV AND TADMORTABLE 6.3
Initial Value Problem (6.5)–(6.6), Lip0 -Norm of the Errors at T = 2
N

NT

Rate

FD2

Rate

40
80
160
320
640
1280

1.233e-02
1.904e-03
4.809e-04
1.253e-04
3.415e-05
1.000e-05

—
2.70
1.99
1.94
1.88
1.77

7.779e-03
1.982e-03
4.905e-04
1.222e-04
3.047e-05
7.558e-06

—
1.97
2.01
2.01
2.00
2.01

EXAMPLE 4 (Nonconvex Flux). In this example we show results of applying our fully
discrete second-order scheme to the following Riemann problem:(· 2¸(u − 1)(u 2 − 4)2,
if x < 0,= 0,u(x, 0) =
(6.7)
ut +
−2,
if x > 0.
4x
The solutions to this initial value problem are depicted at time T = 1.2. Figure 6.7 demonstrates the clear advantage of our new FD2 scheme over the NT scheme; in particular, thelatter seems to give a wrong solution (even after the grid refinement.3 We note, however,that when a more restrictive minmod limiter was used (corresponding to θ = 1 in (6.1)), theNT solution did converge to the entropy solution at the expense of additional smoothing on
the edges of the Riemann fan, which can be noticed in Fig. 6.8.
6.3. One-Dimensional Systems of Hyperbolic Conservation Laws
EXAMPLE 5 (Euler Equations of Gas Dynamics). Let us consider the one-dimensional
Euler System


 
m
µ
¶
ρ
∂ 
ρ 2
∂ 

2
p = (γ − 1) · E − u ,m + ρu + p  = 0,
∂t
∂x
2
E
u(E + p)where ρ, u, m = ρu, p, and E are the density, velocity, momentum, pressure, and totalenergy, respectively. Here, the conserved quantities are uE = (ρ, m, E)T , and the flux vector
function is f (Eu ) = (m, ρu 2 + p, u(E + p))T . We solve this system subject to Riemann
initial data,½
uE (x, 0) =uE L ,
uE R ,

x < 0,
x > 0.

We apply our scalar-designed schemes to this problem in a straightforward manner. We now
prefer the alternative approximation to the flux derivatives, needed in (3.6). The minmod
limiter, (6.1), is employed directly on the corresponding values of f (Eu ) to avoid an expensive
computation of the Jacobian, ∂∂uf .
3
A similar failure of convergence towards the entropy solution by upwind approximation of nonconvex equations
was reported in [5].

HIGH-RESOLUTION CENTRAL SCHEMES265

FIG. 6.7. Riemann IVP (6.7). N = 100.

We compute the solution to two different Riemann problems:• The first Riemann problem was proposed by Sod [46]. The initial data are givenbyuE L = (1, 0, 2.5)T ,uE R = (0.125, 0, 0.25)T .

The approximations to the density, velocity, and pressure obtained by the FD2 scheme arepresented in Figs. 6.9–6.14.

FIG. 6.8. Riemann IVP (6.7). N = 100.

266

KURGANOV AND TADMOR

FIG. 6.9. Sod problem—density. N = 200, T = 0.1644.

FIG. 6.10. Sod problem—density. N = 400, T = 0.1644.

FIG. 6.11. Sod problem—velocity. N = 200, T = 0.1644.

HIGH-RESOLUTION CENTRAL SCHEMES

FIG. 6.12. Sod problem—velocity. N = 400, T = 0.1644.

FIG. 6.13. Sod problem—pressure. N = 200, T = 0.1644.

FIG. 6.14. Sod problem—pressure. N = 400, T = 0.1644.

267

268

KURGANOV AND TADMOR

FIG. 6.15. Lax problem—density. N = 200, T = 0.16.

• The second Riemann problem was proposed by Lax [29]. The initial values are
given byuE L = (0.445, 0.311, 8.928)T ,uE R = (0.5, 0, 1.4275)T .

The results computed by the FD2 scheme are shown in Figs. 6.15–6.20.
Our numerical results for this system are comparable with the results obtained in[38]. We would like to stress again that as in the case of the original NT scheme thecharacteristic decomposition is not required; i.e., our new schemes still can be appliedcomponentwise.
Remark. In all the above 1D hyperbolic examples we have presented the numericalresults obtained by our fully discrete scheme. We also tested the corresponding

FIG. 6.16. Lax problem—density. N = 400, T = 0.16.

HIGH-RESOLUTION CENTRAL SCHEMES

FIG. 6.17. Lax problem—velocity. N = 200, T = 0.16.

FIG. 6.18. Lax problem—velocity. N = 400, T = 0.16.

FIG. 6.19. Lax problem—pressure. N = 200, T = 0.16.

269

270

KURGANOV AND TADMOR

FIG. 6.20. Lax problem—pressure. N = 400, T = 0.16.

semi-discrete scheme on the same examples. The results are very similar, but the semidiscrete scheme is slightly more dissipative than the fully discrete one.
As we have already mentioned, the main advantage of the semi-discrete approach can beseen while we apply our scheme to (degenerate) parabolic convection–diffusion equations.
Below, we show several examples of such problems.
6.4. One-Dimensional Convection–Diffusion Equations
EXAMPLE 6 (Burgers-Type Equation with Saturating Dissipation). We begin with theconvection–diffusion equation with bounded dissipation flux proposed in [28]. Consider
Eq. (1.3) with f (u) = u 2 subject to the Riemann initial data,½u(x, 0) =

1.2,
−1.2,

x < 0,
x > 0.

(6.8)

It was proved in [12] that the solution to this initial value problem contains a subshock
located at x = 0. This is why solving (1.3), (6.8) numerically is quite challenging problem.
Our second-order semi-discrete scheme, SD2, provides a very good resolution of the
discontinuity (Fig. 6.21) while the fully discrete NT scheme fails to resolve it (Fig. 6.22;
see also numerical results in [28]). The SD2 scheme was tested on all the examples from[12]. The numerical results are highly satisfactory, and using this semi-discrete approachno operator splitting is needed (consult [12] for details).
EXAMPLE 7 (Buckley–Leverett Equation). Next, let us consider the scalar convection–diffusion Buckley–Leverett equation (1.2) with Q(u, s) = ν(u)s,u t + f (u)x = ε(ν(u)u x )x ,εν(u) ≥ 0.

(6.9)

This is a prototype model for oil reservoir simulations (two-phase flow). Typically, ν(u)vanishes at some values of u, and (6.9) is a degenerate parabolic equation. Usually, theoperator splitting technique is used (see [6, 8, 20, 22, 23]) to solve it numerically, but thelimitations of such an approach are well known.

HIGH-RESOLUTION CENTRAL SCHEMES

271

FIG. 6.21. Sharp resolution by the SD2 scheme. N = 400, T = 1.5.

We take ε to be 0.01, f (u) to have an s-shaped form,u2,
f (u) = 2
u + (1 − u)2and ν(u) to vanish at u = 0, 1,
ν(u) = 4u(1 − u).

FIG. 6.22. Smeared discontinuity by the NT scheme. N = 400, T = 1.5.

(6.10)(6.11)

272

KURGANOV AND TADMOR

FIG. 6.23. Initial-boundary value problem (6.9)–(6.12). T = 0.2.

The initial function isu(x, 0) =(1 − 3x,0,0 ≤ x ≤ 13 ,1
< x ≤ 1,3

(6.12)

the boundary value of u(0, t) = 1 is kept fixed.
The numerical solution computed by the SD2 scheme for different numbers of grid
points is presented in Fig. 6.23. No exact solution to problem (6.9)–(6.12) is available, but
if compared with the numerical solutions reported in [20], our solutions seem to converge
to the correct entropy solution.
EXAMPLE 8 (Gravitational Effects). We now consider the Buckley–Leverett equation,(6.9), with the same ε = 0.01, the same diffusion coefficient, (6.11), and the flux function
f (u) including gravitational effects:
f (u) =

u2
(1 − 5(1 − u)2 ).
u 2 + (1 − u)2(6.13)

This equation is more complicated than the previous one since we should handle the flux
(6.13) where f 0 (u) changes sign. The numerical solutions to this equation and to Eqs. (6.9)–(6.11) subject to the Riemann initial data,(u(x, 0) =0,

0 ≤ x < 1 − √12 ,1,

1 − √12 ≤ x ≤ 1,

are shown in Fig. 6.24. It can be observed that our semi-discrete scheme provides the samehigh quality of numerical solutions for both of these problems.

HIGH-RESOLUTION CENTRAL SCHEMES

273

FIG. 6.24. Riemann problem for the Buckley–Leverett equation with and without gravitation. T = 0.2.

EXAMPLE 9 (Glacier Growth Model). In this example we consider a one-dimensionalmodel for glacier growth (see [9, 21]). Let a glacier of height h(x, t) rest upon a flat
mountain. Its evolution is described by the nonhomogenious convection–diffusion equationh t + f (h)x = ε(ν(h)h x )x + S(x, t, h).

(6.14)Let ε = 0.01. The typical flux and diffusion coefficient aref (h) =

h + 3h 6,
4

We first look at the Riemann problem with(
1,h(x, 0) =0,ν(h) = 3h 6 .

(6.15)

x < 0,
x > 0,

that describes an outlet into a valley disregarding seasonal variations. To complete this simple
model we use the source S(x, t, h) = S0 (x) if h(x, t) > 0, and S(x, t, h) = max{S0 (x), 0}
if h(x, t) = 0, where0,x < −0.4,

1−0.4 ≤ x ≤ −0.2,
S0 (x) = 2 (x + 0.4),
− 1 x,x > −0.2.
2

The numerical simulations for different number of grid points and at different times arepresented in Figs. 6.25–6.27. These solutions, obtained by our SD2 scheme, seem to be
more accurate than the solutions obtained by the operator splitting method in [21].

274

KURGANOV AND TADMOR

FIG. 6.25. Moving glacier at T = 1.

FIG. 6.26. Moving glacier at T = 2.

FIG. 6.27. Moving glacier at T = 4.

HIGH-RESOLUTION CENTRAL SCHEMES

275

FIG. 6.28. Growing glacier at T = 7.5.

Second, we look at the growth of a new glacier (h(x, 0) ≡ 0). Let the glacier be restricted
to the interval [−5, 5]. Then an appropriate source term isS(x, t) =½
0,−0.01x + 0.05 sin(2πt),x ≤ −5,x > −5.

Here the second, trigonometric term models seasonal variations.
Figures 6.28–6.33 illustrate the glacier growth computed by the Rusanov first-order andby our second-order semi-discrete schemes, SD1 and SD2, with different numbers of grid

FIG. 6.29. Growing glacier at T = 10.5.

276

KURGANOV AND TADMOR

FIG. 6.30. Growing glacier at T = 15.

points. The SD2 scheme provides more accurate resolution of the upstream front, but admits
some oscillations on the glacier downstream. The amplitude of these oscillations remains
small but does not diminish with the grid refinement. Moreover, they tend to propagate up
to the top of the glacier as N increases. At the same time, applying the SD1 scheme with alarge number of grid points gives a very accurate, nonoscillatory solution, comparable withthe one reported in [21].
EXAMPLE 10 (Hyperbolic–Parabolic Equation). We conclude this subsection with an
example of strongly degenerate parabolic (or, hyperbolic–parabolic) convection–diffusion

FIG. 6.31. Growing glacier at T = 22.5.

HIGH-RESOLUTION CENTRAL SCHEMES

277

FIG. 6.32. Growing glacier at T = 30.

equation. Consider Eq. (6.9) with ε = 0.1, f (u) = u 2 , and½
ν(u) =

0,1,|u| ≤ 0.25,|u| > 0.25.

(6.16)

This ν(u) is a discontinuous function, and the equation is therefore of hyperbolic naturewhen u ∈ [−0.25, 0.25] and parabolic elsewhere.

FIG. 6.33. Growing glacier at T = 37.5.

278

KURGANOV AND TADMOR

FIG. 6.34. Initial value problem (6.9), (6.16), (6.17). T = 0.7.

We apply our SD2 scheme to Eqs. (6.9), (6.16) subject to the initial data,

1,



u(x, 0) = −1,



0,

− √12 − 0.4 < x < − √12 + 0.4,
√1 − 0.4 < x < √1 + 0.4,
22

(6.17)

otherwise.

The results for two different number of grid points is shown in Fig. 6.34. We would liketo point out the high resolution of discontinuities and the accurate transition between thehyperbolic and parabolic regions.
6.5. Two-Dimensional Problems
EXAMPLE 11 (Two-Dimensional Burgers-Type Equations). Consider the two-dimensional extension of the equation from Example 10,
u t + (u 2 )x + (u 2 ) y = ε(ν(u)u x )x + ε(ν(u)u y ) y .

(6.18)The numerical results obtained by our SD2 scheme are presented in Figs. 6.35 and 6.36.
In these examples ε = 0.1; ν(u) is given by (6.16) (strongly degenerate parabolic problem)

FIG. 6.35. Pure hyperbolic problem—solution at time T = 0.5 on a 60 ∗ 60 grid.

HIGH-RESOLUTION CENTRAL SCHEMES

279

FIG. 6.36. Degenerate parabolic problem—solution at time T = 0.5 on a 60 ∗ 60 grid.

or ν(u) = 0 (hyperbolic problem). The initial data are equal to −1 and 1 inside two circles ofradius 0.4 centered at (0.5, 0.5) and (−0.5, −0.5), respectively, and zero elsewhere inside thesquare [−1.5, 1.5] × [−1.5, 1.5]. As in the one-dimensional examples, our scheme performs
well in both hyperbolic and hyperbolic–parabolic cases even with relatively small numberof grid points.
EXAMPLE 12 (Two-Dimensional Buckley–Leverett Equation). Finally, we solve thetwo-dimensional convection–diffusion equationu t + f (u)x + g(u) y = ε(u x x + u yy ),(6.19)

with ε = 0.01, the flux function of the formf (u) =u2
(u 2 + (1 − u)2 ),g(u) = f (u)(1 − 5(1 − u)2 ),and the initial data½u(x, y, 0) =1,0,x 2 + y 2 < 0.5,
otherwise.

Note that the above model includes gravitational effects in the y-direction.

FIG. 6.37. Problem (6.19)–(6.21)—solution at time T = 0.5 on a 200 ∗ 200 grid.

(6.20)(6.21)

280

KURGANOV AND TADMOR

FIG. 6.38. Problem (6.19)–(6.21)—solution at time T = 0.5 on a 100 ∗ 100 grid.

The solution, computed in the domain [−1.5, 1.5] × [−1.5, 1.5] by the SD2 scheme, isshown in Figs. 6.37 and 6.38. Our scheme also provides a highly satisfactory approximation
to this model with a nonlinear, degenerate diffusion.
ACKNOWLEDGMENT
The authors thank Dr. A. Medovikov (Institute of Numerical Mathematics, Russian Academy of Sciences,
Russia) for a number of useful discussions and for providing his ODE solver DUMKA3. The work of A. Kurganov
was supported in part by the NSF Group Infrastructure Grant. The work of E. Tadmor was supported in part by
NSF Grant DMS97-06827 and ONR Grant N00014-91-J-1076.

REFERENCES
1. A. M. Anile, V. Romano, and G. Russo, Extended hydrodynamical model of carrier transport in semiconductors, SIAM J. Appl. Math., in press.
2. P. Arminjon, D. Stanescu, and M.-C. Viallon, A two-dimensional finite volume extension of the Lax–Friedrichs
and Nessyahu–Tadmor schemes for compressible flow, in Proc. 6th Int. Symp. on CFD, Lake Tahoe, 1995,
M. Hafez and K. Oshima, editors, Vol. IV, pp. 7–14.
3. P. Arminjon and M.-C. Viallon, Généralisation du schéma de Nessyahu–Tadmor pour une équation hyperbolique à deux dimensions d’espace, C. R. Acad. Sci. Paris Sér. I 320, 85 (1995).
4. F. Bianco, G. Puppo, and G. Russo, High order central schemes for hyperbolic systems of conservation laws,
SIAM J. Sci. Comput. 21, 294 (1999).
5. P. Colella, private communication.
6. H. K. Dahle, Adaptive Characteristic Operator Splitting Techniques for Convection-Dominated Diffusion
Problems in One and Two Space Dimensions, Ph.D. thesis, University of Bergen, Norway, 1988.
7. B. Engquist and O. Runborg, Multi-phase computations in geometrical optics, J. Comput. Appl. Math. 74,
175 (1996).

HIGH-RESOLUTION CENTRAL SCHEMES

281

8. S. Evje, K. H. Karlsen, K.-A. Lie, and N. H. Risebro, Front tracking and operator splitting for nonlinear
degenerate convection–diffusion equations, Institut Mittag-Leffler Report, Stockholm, Sweden, 1997.
9. A. C. Fowler, Glaciers and ice sheets, in The Mathematics of Models for Climatology and Environment, edited
by J. I. Diaz, NATO ASI Series, Vol. 48 (Springer-Verlag, Berlin/New York, 1996), p. 302.
10. K. O. Friedrichs, Symmetric hyperbolic linear differential equations, Comm. Pure Appl. Math. 7, 345 (1954).
11. S. K. Godunov, A finite difference method for the numerical computation of discontinuous solutions of theequations of fluid dynamics, Mat. Sb. 47, 271 (1959).
12. J. Goodman, A. Kurganov, and P. Rosenau, Breakdown of Burgers-type equations with saturating dissipation
fluxes, Nonlinearity 12, 247 (1999).
13. A. Harten, High resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49, 357 (1983).
14. A. Harten, B. Engquist, S. Osher, and S. R. Chakravarthy, Uniformly high order accurate essentially nonoscillatory schemes, III, J. Comput. Phys. 71, 231 (1987).
15. A. Harten, P. D. Lax, and B. van Leer, On upstream differencing and Godunov-type schemes for hyperbolic
conservation laws, SIAM Rev. 25, 35 (1983).
16. A. Jameson, Internat. J. Comput. Fluid Dynamics 4, 171 (1995); 5, 1 (1995).
17. G.-S. Jiang, D. Levy, C.-T. Lin, S. Osher, and E. Tadmor, High-resolution non-oscillatory central schemeswith non-staggered grids for hyperbolic conservation laws, SIAM J. Numer. Anal. 35, 2147 (1998).
18. G.-S. Jiang and E. Tadmor, Non-oscillatory central schemes for multidimensional hyperbolic conservation
laws, SIAM J. Sci. Comput. 19, 1892 (1998).
19. S. Jin and Z. Xin, The relaxation schemes for hyperbolic systems of conservation laws in arbitrary spacedimensions, CPAM 48, 235 (1995).
20. K. H. Karlsen, K. Brusdal, H. K. Dahle, S. Evje, and K.-A. Lie, The corrected operator splitting approach
applied to a nonlinear advection–diffusion problem, Comput. Methods Appl. Mech. Eng. 167, 239 (1998).
21. K. H. Karlsen and K.-A. Lie, An unconditionally stable splitting for a class of nonlinear parabolic equations,
IMA J. Numer. Anal. 19, 609 (1999).
22. K. H. Karlsen and N. H. Risebro, An operator splitting method for nonlinear convection–diffusion equations,
Numer. Math. 77, 365 (1997).
23. K. H. Karlsen and N. H. Risebro, Corrected operator splitting for nonlinear parabolic equations, SIAM J.
Numer. Anal. 37, 980 (2000).
24. R. Kupferman, Simulation of viscoelastic fluids: Couette–Taylor flow, J. Comput. Phys., to appear.
25. R. Kupferman, A numerical study of the axisymmetric Couette–Taylor problem using a fast high-resolution
second-order central scheme, SIAM J. Sci. Comput. 20, 858 (1998).
26. R. Kupferman and E. Tadmor, A fast high-resolution second-order central scheme for incompressible flows,
Proc. Nat. Acad. Sci. 94, 4848 (1997).
27. A. Kurganov, Conservation Laws: Stability of Numerical Approximations and Nonlinear Regularization,
Ph.D. thesis, Tel-Aviv University, Israel, 1997.
28. A. Kurganov and P. Rosenau, Effects of a saturating dissipation in Burgers-type equations, Comm. Pure Appl.
Math. 50, 753 (1997).
29. P. D. Lax, Weak solutions of nonlinear hyperbolic equations and their numerical computation, Comm. Pure
Appl. Math. 7, 159 (1954).
30. B. van Leer, Towards the ultimate conservative difference scheme. III. Upstream-centered finite-difference
schemes for ideal compressible flow, J. Comput. Phys. 23, 263 (1977).
31. B. van Leer, Towards the ultimate conservative difference scheme. V. A second order sequel to Godunov’s
method, J. Comput. Phys. 32, 101 (1979).
32. D. Levy, G. Puppo, and G. Russo, Central WENO schemes for hyperbolic systems of conservation laws,
Math. Model. Numer. Anal. 33, 547 (1999).
33. D. Levy and E. Tadmor, Non-oscillatory central schemes for the incompressible 2-D Euler equations, Math.
Res. Lett. 4, 1 (1997).
34. X.-D. Liu and P. D. Lax, Positive schemes for solving multidimensional systems of hyperbolic conservation
laws, Comput. Fluid Dynamics J. 5, 1 (1996).

282

KURGANOV AND TADMOR

35. X.-D. Liu and S. Osher, Convex ENO high order multi-dimensional schemes without field by field decomposition or staggered grids, J. Comput. Phys. 142, 304 (1998).
36. X.-D. Liu and E. Tadmor, Third order nonoscillatory central scheme for hyperbolic conservation laws, Numer.
Math. 79, 397 (1998).
37. A. A. Medovikov, High order explicit methods for parabolic equations, BIT 38, 372 (1998).
38. H. Nessyahu and E. Tadmor, Non-oscillatory central differencing for hyperbolic conservation laws, J. Comput.
Phys. 87, 408 (1990).
39. H. Nessyahu and E. Tadmor, The convergence rate of approximate solutions for non-linear scalar conservation
laws, SIAM J. Numer. Anal. 29, 1505 (1992).
40. S. Osher, Convergence of generalized MUSCL schemes, SINUM 22, 947 (1984).
41. S. Osher and E. Tadmor, On the convergence of difference approximations to scalar conservation laws, Math.
Comput. 50, 19 (1988).
42. P. Roe, Approximate Riemann solvers, parameter vectors, and difference schemes, J. Comput. Phys. 43, 357(1981).
43. V. Romano and G. Russo, Numerical solution for hydrodynamical models of semiconductors, M3 AS, preprint.
44. C.-W. Shu, Total-variation-diminishing time discretizations, SISSC 6, 1073 (1988).
45. C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes,
J. Comput. Phys. 77, 439 (1988).
46. G. Sod, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws,
J. Comput. Phys. 22, 1 (1978).
47. E. Tadmor, Convenient total variation diminishing conditions for nonlinear difference schemes, SIAM J.
Numer. Anal. 25, 1002 (1988).