arXiv:1212.5315v2 [math.NA] 28 Nov 2016

A HYBRID FD-FV METHOD FOR FIRST-ORDER HYPERBOLIC
CONSERVATION LAWS ON CARTESIAN GRIDS: THE SMOOTH PROBLEM
CASE
XIANYI ZENG

Abstract. We present a class of hybrid FD-FV (finite difference and finite volume) methods for solving general hyperbolic conservation laws written in first-order form. The presentation focuses on one- and two-dimensional Cartesian grids; however, the generalization to higher dimensions is straightforward. These methods use both cell-averaged values and nodal values as dependent variables to discretize the governing partial differential
equation (PDE) in space, and they are combined with method of lines for integration in
time. This framework is absent of any Riemann solvers while it achieves numerical conservation naturally. This paper focuses on the accuracy and linear stability of the proposed
FD-FV methods, thus we suppose in addition that the solutions are sufficiently smooth.
In particular, we prove that the spatial-order of the FD-FV method is typically one-order
higher than that of the discrete differential operator, which is involved in the construction of
the method. In addition, the methods are linearly stable subjected to a Courant-FriedrichLewy condition when appropriate time-integrators are used. The numerical performance of
the methods is assessed by a number of benchmark problems in one and two dimensions.
These examples include the linear advection equation, nonlinear Euler equations, the solid
dynamics problem for linear elastic orthotropic materials, and the Buckley-Leverett equation.

1. Introduction
Hyperbolic conservation laws [6] are used to describe many problems in the areas of
fluid dynamics such as compressible flows, aeroacoustics, magnetohydrodynamics. More
recently, they are also used to formulate the problems of solid dynamics in both infinitesimal strain regime [18] and finite strain regime [1].
Conventional numerical methods for conservation laws are typically constructed in a
way such that the same discretization concept is used to construct the discrete formula of
each dependent variable. In the case of many conventional finite difference (FD) methods,
such as [16], this discrete formula represents a discretization of governing equation in its
strong form; whereas in the case of classical finite volume (FV) methods, such as those
in [27] and [19], the discrete formula represents a numerical realization of the integral
form of the governing equation to achieve natural conservation.
One disadvantage of these methods, however, is that they typically require large (and
overlapping) stencils to achieve higher-order accuracy. For example, to achieve secondorder accuracy in space, an explicit one-dimensional finite volume method requires information of five cells to update one particular cell average (i.e., two additional cells to each
side); and a similar situation holds for the finite difference methods. A more compact
Date: November 29, 2016.
2010 Mathematics Subject Classification. 65M12 and 35L65.
Key words and phrases. First-order hyperbolic conservation laws; High-order accuracy; FD-FV method;
Cartesian grids; Linear stability.
1

2

X. ZENG

formulation can be achieved by the discontinuous Galerkin (DG) formulation (see, for example [5]), however, at the cost of having duplicate degrees of freedom on the interface
between elements.
In recent years, there is continuous interest in developing numerical methods that are
more compact than the conventional FD or FV methods, but still require fewer unknowns
than the DG-type methods, for example, by using different discretization formula for different unknowns. To the best knowledge of the author, numerical methods of this type date
back to the 1970s: van Leer [26] constructed a third-order accurate method (scheme V of
the reference) that used both cell averages and nodal values at cell faces for 1D conservation laws. This work, however, is not further explored. A non-exhausted list of similar
works that employ variables with different meanings includes the multi-moment methods
([28] and [13]), the space-time control element and solution element methods ([3], [4], and
[29]), the staggered mesh method [22] and dual-mesh method [12], the PN P M methods [7],
and more recently the active flux scheme [8].
The present work also falls into this category in that we choose the discrete variables
to approximate both cell averages and nodal values that collocate at cell interfaces. In
particular, we focus on discretization in space and use the method of lines for integration
in time. In our computational framework, on the one hand, the semi-discretization formula
for the cell averages is constructed using the integral form of the governing equation, just
like the FV methods, but with the numerical fluxes replaced by physical fluxes evaluated
at nodal values. On the other hand, the semi-discretization formula for the nodal values is
constructed by local linearization and applying discrete differential operators that involve
both nodal values and cell averages. In the one-dimensional case, the latter construction
has direct connection to the Hermite interpolation polynomials. Due to these reasons, we
call the present method the hybrid FD-FV methods (or simply FD-FV methods). The FDFV methods are more compact than the conventional finite difference and finite volume
methods. For example in one space dimension, we have in average two variables for each
mesh cell, which enable us to construct higher-order discrete differential operators than the
finite difference methods given the same stencil. Furthermore, we will show that using a
discrete differential operator that is constructed to pth -order operator in space, the resulting
FD-FV method is formally (p + 1)th -order accurate. Thus we are able to construct a fourthorder accurate FD-FV method using a stencil of three consecutive cells (see the method
3,up−biased
given by [D x
] in Section 2.3), whereas in the case of conventional FD or FV
methods the spatial-order is only first-order.
In this paper, we only consider explicit time-integrators for simplicity. In addition, because we focus on the analysis of the spatial order of accuracy of the proposed FD-FV
methods and their linear stability properties, we also assume that the solutions are sufficiently smooth. Enhancement of the nonlinear stability of the methods and their performance to solve discontinuous solutions will be addressed in future work. To this end, the
remainder of the paper is organized as follows. We describe the general formulation of
the FD-FV methods for one-dimensional first-order hyperbolic system in Section 2. The
proofs of the existence of the coefficients and the theorem concerning the spatial order of
accuracy of the FD-FV methods are provided in Appendix A and Appendix B, respectively.
Section 3 extends the FD-FV methodology to two-dimensional Cartesian grids, with evidence of the conjectured order of spatial accuracy provided in Appendix C. Numerical
examples are offered in Section 4, and finally Section 5 concludes this paper.

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

3

2. The hybrid FD-FV approach for one-dimensional problems
We describe in this section the fundamental framework of the hybrid FD-FV method
for the one-dimensional hyperbolic conservation law in first-order form:
(2.1)

wt + f (w
w) x = 0,

which is defined on the rectangular space-time domain (x, t) ∈ Ω × [0, T ], and supplemented
by appropriate initial and boundary conditions. We assume that w(x, t) takes values in Rd ,
and f : Rd → Rd is a continuously differentiable function of w. Here and throughout the
remainder of the paper, we use normal font variables such as t and x to denote scalars,
whereas we use bold font to designate vectors or matrices.
2.1. General FD-FV formulation for scalar problems. We first divide the computaS
tional domain Ω into cells with non-overlapping interiors: Ω = j C j , where C j = [x j−1/2 , x j+1/2 ].
We denote by h j = x j+1/2 − x j−1/2 and x j = (x j−1/2 + x j+1/2 )/2 the length and the center of
the cell Ω j , respectively. In this paper, we assume uniform grids (h j ≡ h), but the methodology applies equally to nonuniform 1D grids.
For simplicity, we first consider scalar problems, in which case w and f of (2.1) are
replaced by scalar functions w and f :
(2.2)

wt + f (w) x = 0.

The target of the hybrid FD-FV methods is to compute the numerical approximations to
the nodal values at cell faces:
(2.3)

wnj+1/2 ≈ w(x j+1/2 , tn ),

and those to cell-averaged values for each cell:
Z
1
w(x, tn )dx.
wnj ≈
(2.4)
h Cj
Because we use the method of lines to decouple the discretization schemes in space and in
time, we only consider the semi-discretized variables in this section, namely:
Z
1
w j+1/2 (t) ≈ w(x j+1/2 , t), w j (t) ≈
w(x, t)dx.
h Cj
n

Here and throughout the remainder of the paper, we use † j , †nj+1/2 and † j , † j+1/2 to designate the approximate cell averages and nodal values to the fully discretized system and
semi-discretized system, respectively. In the latter case, we often omit expressing the dependence in t for convenience.
Remark 2.1. For the purpose of analysis, we suppose that the cell averages are initialized
exactly by the initial condition. But this is not necessary in practical computations. In the
numerical examples provided in Section 4, we use appropriate quadrature rules to initialize
the cell averages w0j = w j (0).
On the one hand, in order to derive the semi-discretized formula for cell averages, we
integrate equation (2.2) over each cell C j to obtain:

 Z


i
 1 h 
d  1

w(x, t)dx + f w(x j+1/2 , t) − f w(x j−1/2 , t) = 0,
dt h C j
h
which then leads to the following discretization scheme in space for w j :


i
dw j 1 h 
(2.5)
+ f w j+1/2 − f w j−1/2 = 0,
dt
h

4

X. ZENG

where w j±1/2 are the nodal variables.
On the other hand, we linearize the governing equation at a node x j+1/2 to derive the
def

semi-discretized formula for nodal variables. Supposing that J(w) == f 0 (w) is the Jacobian
of the flux function, we seek a numerical scheme for constructing the ODE for w j+1/2 in
the form:
dw j+1/2
+ J(w j+1/2 )[D x w] j+1/2 = 0.
(2.6)
dt
Here we use [D x ] to denote a discrete differential operator that uses both adjacent cellaveraged values and nodal values, to distinguish it from the conventional finite difference
operators D x , which only involves the nodal variables. A general (linear) form of this
operator is:
q
q
1 X
1X
[D x w] j+1/2 =
(2.7)
αl w j+l +
βl w j+1/2+l ,
h l=−q+1
h l=−q
where q is a positive integer that indicates the stencil of the method, and αl and βl are
constants to be determined later.
Remark 2.2. One may recover a conventional finite difference method, by setting all constants αl to zero. In this case, the equations for the nodal variables decouple from the
cell averages. Thus in theory, the conventional finite difference methods are a subset of
the present hybrid FD-FV framework. However, this is not encouraged, for the reason
discussed in the accuracy analysis below.
2.2. Spatial accuracy of FD-FV operators. There are mainly two considerations in design the constants αl and βl : the first one is that we want to achieve a certain order of
accuracy, and the second one is that we need at least linear stability when these operators
are applied to solve an simple advection problem. To this end, we first focus on the accuracy analysis of the resulting FD-FV method by using (2.5) and (2.6) with a set of constants
in this section, and move on to the linear stability analysis in Section 2.3.
We apply the Taylor series expansion to the nodal variables and cell-averaged variables
to study the accuracy of the FD-FV operators [D x ] and the resulting FD-FV method. Considering performing the expansion about a node x j+1/2 , we have for a nearby nodal variable
w(x j+l+1/2 ) (here we omit the dependence on t for simplicity):
∞
X
(lh)m m
(2.8)
w(x j+l+1/2 ) = w(x j+1/2 ) +
∂ w(x j+1/2 ),
m! x
m=1
and for the average of the cell C j+l :


Z
Z
∞
X

(x − x j+1/2 )m m
1 x j+l+1/2
1 x j+l+1/2 
w(x j+1/2 ) +
w(x)dx =
∂ x w(x j+1/2 ) dx
h x j+l−1/2
h x j+l−1/2
m!
m=1
(2.9)
∞
m+1
m+1
X
[l
− (l − 1) ]hm m
= w(x j+1/2 ) +
∂ x w(x j+1/2 ),
(m + 1)!
m=1
where we use ∂mx to denote the mth derivative in x. To this end, we may define the order of
the operator [D x ] by Definition 2.3, and relate it to the coefficients by Theorem 2.5.
Definition 2.3. The FD-FV operator [D x ] defined by (2.7) is pth -order accurate for some
integer p ≥ 1 if:
(2.10)

p+1

[D x w] j+1/2 = ∂ x w j+1/2 + c p ∂ x w j+1/2 h p + O(h p+1 )

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

5

for any sufficiently smooth function w(x), where each subscript j+1/2 denotes the exact data
at x j+1/2 . In this formula, c p , 0 is a constant that is independent of the function w and the
cell size h.
Example 2.4. The operator [D x w] j+1/2 = (w j+1/2 − w j )/h is a first-order operator with
c1 = −1/3; and the operator [D x w] j+1/2 = 2(−w j+3/2 + 3w j+1 − 2w j+1/2 )/h is a second-order
operator with c2 = −1/12.
Theorem 2.5. The FD-FV operator (2.7) is pth -order accurate for some integer p ≥ 1 if
and if the constants αl , −q + 1 ≤ l ≤ q and βl , −q ≤ l ≤ q satisfy:

(2.11)

q
q
X
X
lm+1 − (l − 1)m+1
lm
αl +
βl = δm1 ,
(m + 1)!
m!
l=−q+1
l=−q
q
q
X
X
l p+1
l p+2 − (l − 1) p+2
αl +
βl = c p .
(p + 2)!
(p + 1)!
l=−q
l=−q+1

m = 0, 1, · · · , p,

Here δm1 is the Kronecker delta that takes the value 1 when m = 1 and 0 otherwise.
The proof of Theorem 2.5 is straightforward: one simply needs to plug (2.8–2.9) into
(2.7) and use Definition 2.3. What is more interesting is the existence of the coefficients
that satisfy (2.11). It turns out that at least in the 1D case, these coefficients have close
relation to the construction of Hermite interpolation polynomials, see Appendix A. Thus
in most situations, we can construct the coefficients αl and βl that leads to an arbitrary
order operator [D x ]. Note that the preceding results extend naturally to nonuniform grids,
seeing their close relation to the Hermite interpolation theory.
Our next target is to establish the connection of the order of [D x ] and that of the resulting
semi-discretized FD-FV scheme, between which the latter is to be defined. Considering
applying the FD-FV method (2.5–2.6) to solve the Cauchy problem for the scalar advection
equation:
wt + cw x = 0,

(2.12)

(x, t) ∈ R × [0, T ],

for some constant c , 0 and the simple wave initial condition w(x, 0) = eikx , where k is the
wavenumber. The exact solutions for the cell averages and nodal variables at t = T are:
Z


1 jh ik(x−cT )
1
(2.13)
w∗j =
e
dx = e−ickT ei jθ 1 − e−iθ , w∗j+1/2 = e−ickT ei jθ .
h ( j−1)h
iθ
Here i denotes the imaginary unit, θ = kh is the numerical wavenumber, and we suppose
x j−1/2 = ( j − 1)h as well as x j = ( j − 1/2)h for the ease of expression; we also use the
superscript ∗ to denote the exact solutions.
The semi-discretized solutions to (2.12) by using (2.5–2.6) leads to a coupled ODE
system for w j (t) and w j+1/2 (t):
dw j c
+ (w j+1/2 − w j−1/2 ) = 0,
dt
h
+ c[Aw] j+1/2 + c[Nw] j+1/2 = 0,

(2.14)

dw j+1/2
dt
where we split the operator [D x ] into two components by [D x ] = [A] + [N]:
(2.15)

def

[Aw] j+1/2 ==

q
1 X
αl w j+l ,
h l=−q+1

def

[Nw] j+1/2 ==

q
1X
βl w j+l+1/2 .
h l=−q

6

X. ZENG

It is not difficult to see that if the initial data is exact, the solutions to (2.14) are given by
the form:
1
(2.16)
w j = A(t)ei jθ (1 − e−iθ ), w j+1/2 = N(t)ei jθ , ∀ j.
iθ
Here A(·) and N(·) are two scalar functions that solve the ODE system:
(2.17)

A0 (t) + ickN(t) = 0,
N (t) + cka(θ)A(t) + ckb(θ)N(t) = 0,
0

with initial condition A(0) = N(0) = 1. The two functions a(θ) and b(θ) of (2.17) are:
q
q
−iθ X
X
def 1
def 1 − e
ilθ
e
α
,
and
b(θ)
==
(2.18)
a(θ) ==
eilθ βl .
l
θ l=−q
iθ2 l=−q+1
Comparing (2.13) and (2.16), we have the following definition of the order of accuracy of
the semi-discretized FD-FV scheme.
Definition 2.6. The semi-discretized FD-FV scheme defined by [D x ] is at least pth -order
accurate, if the following is true. If the method is applied to solve the Cauchy problem
specified by (2.12) with initial condition w(x, 0) = eikx , in which case the approximated
solutions are in the form of (2.16), then the following identities hold:




(2.19)
A(t) = e−ickt 1 + ca (t)θ p + O(θ p+1 ) , N(t) = e−ickt 1 + cn (t)θ p + O(θ p+1 ) ,
where ca and cn are functions of time and they are independent of θ.
Remark 2.7. Although we express the two scalars A(·) and N(·) as functions of time, they
actually also depend on the spatial parameter θ = hk. However, in our proof of spatialorder of accuracy in next theorem we suppose that the mesh size is fixed; thus we omit the
dependence of these functions on θ for simplicity.
The main result concerning the spatial-order of the semi-discretized FD-FV method and
the designed-order of the operator [D x ] is given by the next theorem.
def Pq
Theorem 2.8. Given a pth -order accurate [D x ] such that b0 == l=−q βl , 0, it leads to a
semi-discretized FD-FV method that is at least (p + 1)th -order accurate provided also that
cb0 > 0, where c is the constant advection speed of eq. (2.12).
The proof of this theorem is based on solving (2.17) exactly and performing Taylor
series expansions, which is provided as Appendix B. The important implication of Theorem 2.8 is that by allowing a mixed form of cell-averages and nodal values to design the
differential operator [D x ], one achieves one-order higher accuracy compared to conventional finite difference methods. To see why this is true, we note that conventional finite
difference operators fit into the hybrid FD-FV framework by setting all αl to be zero (Remark 2.2), in which case b0 = 0. By slight modification of the proof of Appendix B to
b0 = 0, it is easy to see that the maximum possible order of accuracy is pth , i.e., the designed order of the differential operator. It is also noted that the upwind-like condition
cb0 > 0 is very important to prevent a second mode from blowing up as the mesh size goes
to zero, see eq. (B.5) and the discussion that follows.
To conclude this section, we listed in Table 1 several operators [D x ] as well as their
designed order and corresponding b0 . In the table, we use “F-biased” and “B-biased” to
denote forward biased or backward biased stencils, respectively. Due to the upwind-like
condition cb0 > 0, we also use the terminology “upwind” (designated by the superscript up )
or “upwind-biased” (designated by the superscripts up−b or up−biased ) to denote an operator.

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

7

Table 1. FD-FV operators [D x ] with b0 , 0.
[D
Order
b0
 x w] j+1/2

1st Forward: 2 w j+1 − w j+1/2 /h
−2


Backward: 2 w j+1/2 − w j /h
2


2nd Forward: 2 −w j+3/2 + 3w j+1 − 2w j+1/2 /h
−6


Backward: 2 2w j+1/2 − 3w j + w j−1/2 /h
6


3rd Forward: w j+2 − 8w j+3/2 + 17w j+1 − 10w j+1/2 /(2h)
−9


F-biased: −2w j+3/2 + 7w j+1 − 4w j+1/2 − w j /(2h)
−3


B-biased: w j+1 + 4w j+1/2 − 7w j + 2w j−1/2 /(2h)
3


Backward: 10w j+1/2 − 17w j + 8w j−1/2 − w j−1 /(2h)
9


4th Forward: −2w j+5/2 + 7w j+2 − 16w j+3/2 + 23w j+1 − 12w j+1/2 /(2h) −15


F-biased: w j+2 − 12w j+3/2 + 31w j+1 − 18w j+1/2 − 2w j /(6h)
−5


B-biased: 2w j+1 + 18w j+1/2 − 31w j + 12w j−1/2 − w j−1 /(6h)
5


Backward: 12w j+1/2 − 23w j + 16w j−1/2 − 7w j−1 + 2w j−3/2 /(2h)
15
2,up

For example, [D x ] denotes the second-order backward operator if c > 0, and the secondorder forward operator if c < 0.
2.3. Linear stability analysis. We study the linear stability of both the semi-discretized
and fully-discretized FD-FV methods in this section. In the semi-discretized case, we
apply the von Neumann analysis to the same advection problem as in previous section. In
the view of (2.13) and (2.16), we can plot the dissipation error and the dispersion error by
comparing A(t) and N(t) with e−ickt , respectively. For simplicity, we suppose c > 0 and
only consider operators of Table 1 with positive b0 (Theorem 2.8).
For a simple wave problem, the numerical solutions to the semi-discretized system is
obtained by solving (2.17), which gives rise to two eigenvalues λ1 and λ2 ; these values correspond to the two superposed traveling waves in the numerical solutions A(t) and N(t), as
demonstrated by (B.3). Thus contrary to the von Neumann analysis to conventional finite
difference systems, we study the dispersion property of the FD-FV methods by plotting
def

θ̃1,2 against θ = kh, where θ̃1,2 == k̃1,2 h and k̃1,2 are the wave numbers of the corresponding
traveling wave, i.e., =(kλ1,2 ), the imaginary part of kλ1,2 . Similarly, we study the dissipation property by plotting the numerical dissipation rate ε̃1,2 against ε ≡ 0, where ε̃1,2 are
the real parts of θλ1,2 , <(θλ1,2 ), and it needs to be non-negative for linear stability. To this
end, we can rewrite (B.3) in the form:
A(t) = a1 e−(ε̃1 +icθ̃1 )t/h + a2 e−(ε̃2 +icθ̃2 )t/h and N(t) = n1 e−(ε̃1 +icθ̃1 )t/h + n2 e−(ε̃2 +icθ̃2 )t/h ,
where a1,2 and n1,2 are constants that are independent of t. These quantities are important
for the accuracy analysis, as shown in Appendix B, but they do not affect the linear stability
of the semi-discretized FD-FV schemes.
In Figure 2.1 we plot the dispersion relations in the range θ ∈ [0, 2π], in which the right
end of the interval corresponds to resolving the wave length by one cell. Clearly, we see
that the wave speed associated with λ1 is closer to the true wave speed as the order of the
operator increases; whereas the wave associated with λ2 travels in the opposite direction as
the physical wave. Note, however, the second wave has a magnitude diminishing very fast
as h → 0, see the accuracy analysis in Appendix B. Similarly, the dissipation properties

8

X. ZENG

5

θ̃1

4

0

Exact
[D x1,up]
[D x2,up]
[D x3,up−b]
[D x3,up]
[D x4,up−b]
[D x4,up]

−1
−2
−3

θ̃2

6

3

−4

2

−5

1
0
0

2

4

6

−6
0

2

4

θ

θ

(a) θ̃1 vs. θ.

(b) θ̃2 vs. θ.

θ − 2π
[D x1,up]
[D x2,up]
[D x3,up−b]
[D x3,up]
[D x4,up−b]
[D x4,up]
6

Figure 2.1. Dispersion property of various FD-FV operators [D x ]:
(2.1a) θ=(λ1 ) vs. θ, (2.1b) θ=(λ2 ) vs. θ. Here “up ” and “up−biased ”
denote the upwind and the upwind-biased stencils, respectively.

ε̃ 1

10

5

15

ε ≡1,up
0
[D x ]
2,up
[D x ]
[D x3,up−b]
[D x3,up]
[D x4,up−b]
[D x4,up]

0

−5
0

ε ≡1,up
0
[D x ]
2,up
[D x ]
[D x3,up−b]
[D x3,up]
[D x4,up−b]
[D x4,up]

10

ε̃ 2

15

5

0

2

4

6

−5
0

2

4

θ

θ

(a) ε̃1 vs. θ.

(b) ε̃2 vs. θ.

6

Figure 2.2. Dissipation property of various FD-FV operators [D x ]:
(2.2a) θ<(λ1 ) vs. θ, (2.2b) θ<(λ2 ) vs. θ.
of the FD-FV operators are plotted in the same range in Figure 2.2. From the figures, the
4,up
dissipation rate of [D x ] is positive around θ ≈ π, or two cells per wave length. Thus the
method is linearly unstable in this regime on the semi-discretized level. This observation is
consistent with a class of leapfrog methods [15], which claims that the number of upwind
points in the stencil cannot exceed the number of downwind points by more than two.
Furthermore, the dispersion and dissipation relations associated with λ1,2 show certain
symmetry properties about θ = π; this is further discussed in Remark B.1 of Appendix B.
Remark 2.9. An alternative way to investigate the linear stability of the FD-FV methods is
to write the semi-discretized system in matrix form, and plot the eigenvalue distribution of
the resulting linear ODE system. This approach not only applies to periodic boundary conditions but also works well with other types of boundary conditions, such as the examples
given in Section 2.5.

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

9

Table 2. Pair of FD-FV operator and Runge-Kutta time-integrator, and
the corresponding maximum allowed Courant numbers.
Operator Time-integrator Formal order of the FD-FV method νmax
1,up
[D x ]
RK2
2
1.0
2,up
[D x ]
RK3
3
0.409
3,up−biased
[D x
]
RK4
4
0.808
3,up
[D x ]
RK4
4
0.309
4,up−biased
[D x
]
RK5
5
0.494
Finally, we consider the fully-discretized FD-FV methods by pairing each [D x ] operator with an appropriate time-integrator. For simplicity, only explicit Runge-Kutta type
methods are considered in this paper. In the view of Theorem 2.8, for each [D x ] that is
a pth -order accuracy operator we choose a (p + 1)th -order accurate scheme to integrate the
resulting ODE system in time. In particular, we consider the 2nd -order and 3rd -order total variational diminishing Runge-Kutta schemes by S. Gottlieb et al. [9], “the” 4th -order
method as described by E. Hairer et al. [11], and the 5th -order method by H. A. Luther [20].
These methods are denoted by RK2, RK3, RK4, and RK5, respectively; and their corresponding FD-FV operators [D x ] as well as the maximum allowed Courant number (denoted by νmax ) are summarized in Table 2.
2.4. Extension to systems of equations. In this section, we extend the methods presented before to a system of hyperbolic conservation laws. Considering the governing
equation (2.1), we define the nodal variable w j+1/2 and the cell average w j similarly as
(2.3–2.4). Then in analogy to (2.5), the semi-discretized formula for w j is given by:


i
wj 1 h 
dw
+ f w j+1/2 − f w j−1/2 = 0.
dt
h
To construct the spatial discretization for nodal values, we use the hyperbolic assumption
def
and compute the characteristic decomposition of the flux Jacobian J = R Λ L , where J ==
w. We denote the diagonal entries of the diagonal matrix Λ by λ1 ≤ λ2 ≤ · · · ≤ λd , and
∂ f /∂w
let r k and l k be the right and left eigenvalues associated with λk , 1 ≤ k ≤ d, respectively,
such that R = [rr 1 , r 2 , · · · , r d ] and R −1 = L = [ll1 , l 2 , · · · , l d ]T . Using these notations,
equation (2.1) is rewritten locally as decoupled scalar problems as:

(2.20)

L w t + Λ L w x = 0.
To this end, the semi-discretized formula for w j+1/2 is given by:
Lw j+1/2 )
d(L
d(ωk, j+1/2 )
Lw )] j+1/2 = 0 or
+ Λ [D x (L
+ λk [D x ω] j+1/2 = 0, k = 1, · · · , q,
dt
dt
def

w j+1/2 ), and ωk, j+1/2 == l Tk w j+1/2 . Thus
where L and Λ are computed using flux Jacobian J (w
1,up
by using a particular FD-FV operator, such as [D x ], to solve the hyperbolic system (2.1),
we mean to apply the scalar operator to each characteristics.
Example 2.10. Consider the hyperbolic equation:
" # " #
"
# "
#"
#
u
v
u+v
1 0 u+v
+
= 0 =⇒
+
= 0,
v t u x
u − v t 0 −1 u − v x

10

X. ZENG

1,up

the spatial discretization for (u j+1/2 , v j+1/2 ) using [D x ] is given by:

i

d(u j+1/2 + v j+1/2 ) 2 h



(u
+
v
)
−
(u
+
+
v
)
=0

j+1/2
j+1/2
j
j

dt
h
.

h
i

d(u
−
v
)
2

j+1/2
j+1/2


− (u j+1 − v j+1 ) − (u j+1/2 − v j+1/2 ) = 0

dt
h
2.5. Boundary conditions. To conclude this section, we describe the numerical treatments of the Dirichlet boundary condition for incoming characteristics, and briefly study
the linear stability property by applying the methods to solve simple model problems. For
simplicity, we consider the advection problem:
wt + cw x = 0,

(x, t) ∈ [0, x0 ] × [0, T ],

where

c > 0,

x0  0,

with the Dirichlet boundary condition:
w(0, t) = g(t)
for some known data g(t), t ∈ [0, T ]. Then we follow the next two-step procedure to enforce
p
this boundary condition in the spatial discretization of the FD-FV method given by [D x ],
where p denotes the order of this operator. Supposing x1/2 = 0, we first set w1/2 = g. Next,
p
for all small indices j such that calculating [D x w] j+1/2 requires data collocated to the left
q
of x = 0, we use an operator [D x w] of order q and with a biased stencil at this node, such
that the corresponding stencil stays within [0, x0 ]. In order not to degrade the overall order
of accuracy of the method, q should be at least p−1 [10].
3,up

Example 2.11. Consider the FD-FV method given by [D x ] (i.e., the 3rd -order backward
operator of Table 1), then applying this operator at x1/2 and x3/2 requires a stencil that
penetrates into outside the computational domain. To enforce the boundary condition at
2,up
3,up−biased
x = 0, we set w1/2 = g and apply [D x ] or [D x
] at x3/2 to discretize the equation
for the nodal variable in space.
Although the hyperbolic problems are in general not well-posed with Neumann boundary conditions, in practice we may use the numerical boundary condition of Neumann type
for nonlinear hyperbolic systems. For example, we consider the 1D Euler equations:

  
 ρ   ρu 

   2
(2.21)
 ρu  +  ρu + p  = 0,
E t (E + p)u x
which is defined on (x, t) ∈ [0, x0 ] × [0, T ]. Here ρ, u, and p are density, velocity, and
pressure, respectively; and E is the total energy per unit volume that is determined by the
ideal gas law: E = p/(γ − 1) + ρu2 /2 with γ ≡ 1.4. Then the wall-boundary condition at
x = 0:
u(0, t) ≡ 0
leads to the Neumann boundary condition for the total energy at this location: E x (0, t) ≡ 0.
To see this, we first have u(0, t) = 0 and ut (0, t) = 0, then from the momentum and energy
equations we obtain:
px
1
0 = (ρu)t + (ρu2 + p) x = p x =⇒ E x =
+ ρ x u2 + ρuu x = 0 at x = 0.
γ−1 2
Hence we describe here the numerical Neumann boundary condition, such as w x (0, t) = g(t)
at the left boundary x1/2 = 0. This procedure is also composed of two steps: First we
q, f orward
compute w1/2 such that [D x
]1/2 = g, then for all small indices j such that calculating
p
[D x w] j+1/2 requires data collocated to the left of x1/2 we use the same strategy as in the

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

[D x1,up]
[D x2,up]
[D x3,up−b]
[D x3,up]
[D x4,up−b]

Imaginary part

4
2
0
−2

4

Imaginary part

6

11

[D x1,up]
[D x2,up]
[D x3,up−b]
[D x3,up]
[D x4,up−b]

2

0

−2

−4
−6
−10

−8

−6

−4

Real part

−2

0

(a) Periodic boundary condition.

2

−4
−7

−6

−5

−4

−3

−2

Real part

−1

0

1

2

(b) Dirichlet boundary condition.

Figure 2.3. The spectra of various FD-FV schemes to solve wt + w x = 0
on x ∈ [0, 1] using N = 50 cells: (2.3a) periodic boundary condition, and
(2.3b) Dirichlet boundary condition at x = 0.
Dirichlet case. Similarly as before, q should be at least p−1 in order not to degrade the
formal order of accuracy.
3,up

Example 2.12. Consider again the FD-FV method given by [D x ] and homogeneous
3, f orward
boundary condition w x = 0, then we first compute w1/2 by setting [D x
w]1/2 = 0 to
have:
1
w1/2 = (wn2 − 8wn3/2 + 17w1 ).
10
2,up
3,up−biased
Next we apply [D x ] or [D x
] at x3/2 as in Example 2.11.
Finally, we study the spectral property of the resulting ODE systems obtained by applying the FD-FV methods to solve a model problem with various boundary conditions. In
particular, we solve the advection problem wt + w x = 0 on the domain x ∈ [0, 1] by meshing
the latter into N = 50 uniform cells. In the first scenario, we consider the periodic boundary
condition w(0, t) = w(1, t), and denote by W the array of unknowns:
W = [w1 , w2 , · · · , wN , w1/2 , w3/2 , · · · , wN−1/2 ]T .
Then the spatial discretization leads to a ODE system for W :
1
1
(2.22)
W t = DW , h = .
h
N
Here D is a 2N × 2N matrix that is determined by the operator [D x ].
1,up
matrix D corresponding to [D x ] is given by:



 0 1 · · · 0 0 
 0 0
 0 0 · · · 0 0 
 1 0
"
#



0 I − L1



D=
, L 1 =  ... ... . . . ... ...  and L 2 =  ... ...
L2 −2II
2L



 0 0 · · · 0 1 
 0 0



1 0 ··· 0 0
0 0

For example, the
···
···
..
.

0
0
..
.

1
0
..
.






 ,

· · · 0 0 
··· 1 0

with 0 and I being the N × N zero matrix and identity matrix, respectively. The semidiscretization of the FD-FV method is linearly stable if and only if all the eigenvalues of
the matrix D of (2.22) stay within the left complex domain (i.e., the real parts are nonpositive). To this end, we plot the eigenvalue distributions of the linearly stable operators
in Section 2.3, as shown in Figure 2.3a. Clearly, all the eigenvalues are in the left half
complex domain, indicating that the methods are linearly stable.

12

X. ZENG

In the second scenario, we consider the Dirichlet boundary condition at x1/2 = 0, namely
w(0, t) = g(t); and we denote by W the array of 2N − 1 unknowns:
W = [w1 , w2 , · · · , wN , w3/2 , · · · , wN−1/2 ]T .
Then applying the various FD-FV methods to solve this problem yields the ODE system:
1
W t = DW + g,
h
where g designates the contribution from boundary conditions. The eigenvalue distributions computed using several FD-FV operators are presented in Figure 2.3b, which again
3,up−biased
indicates linear stability for all operators. Note that we use [D x
] at x3/2 and use
3,up
3,up
[D x ] at the last node xN+1/2 in the cases of FD-FV methods that are defined by [D x ]
4,up−biased
and [D x
].
3. Extension to two space dimensions
The extension of the preceding methodology to two space dimensions is complicated by
the fact that for a two-dimensional cell, the face-averaged flux cannot be evaluated exactly
by using nodal values as in one-dimensional case. In particular, consider the general firstorder conservation law written in Cartesian coordinates:
w) x + g (w
w)y = 0,
w t + f (w

(3.1)

where w (x, y, t) takes value in Rd , and f , g : Rd → Rd are both continuously differentiable.
Let the rectangular computational domain Ω be divided into cells with non-overlapping
interiors: Ω = ∪ j1 , j2 C j1 , j2 , where C j1 , j2 = [x j1 −1/2 , x j1 +1/2 ]×[y j2 −1/2 , y j2 +1/2 ]; and similar to
the 1D case, we consider uniform grids for simplicity and let x j1 +1/2 − x j1 −1/2 ≡ h1 , y j2 +1/2 −
y j2 −1/2 ≡ h2 , and define x j1 = (x j1 −1/2 + x j1 +1/2 )/2, y j2 = (y j2 −1/2 + y j2 +1/2 )/2. Clearly, the
cell-averaged value is defined such that:
Z
1
w nj1 , j2 ≈
(3.2)
w (x, y, tn )dxdy.
h1 h2 C j , j
1 2

However, integrating (3.1) over C j1 , j2 yields:


Z

d  1
w (x, y, t)dxdy

dt h1 h2 C j , j
1 2
 Z

Z y j +1/2 
y j2 +1/2 

 
2
1  1
1
+ 
f w (x j1 +1/2 , y, t) dy −
f w (x j1 −1/2 , y, t) dy
h1 h2 y j −1/2
h2 y j −1/2
 Z 2

Z x2 j +1/2 
x j1 +1/2 

 

1
1
1  1
g w (x, y j2 +1/2 , t) dx −
g w (x, y j2 −1/2 , t) dx = 0,
+ 
h2 h1 x j −1/2
h1 x j −1/2
1

1

which contains the face-averaged fluxes that are not available directly. In order to obtain
a high-order scheme, these fluxes need to be evaluated to a certain order of accuracy. One
way to do this is to define the face-averaged flux itself as a variable, however it requires
additional governing equations to describe the dynamics of this variable, which increase the
complexity of the methodology significantly. Instead, we choose a more direct extension
of the 1D methodology: we select several nodes on a particular cell face and approximate
the face-averaged flux by a weighted sum of the fluxes evaluated at these nodes. That
is, we first choose a quadrature rule for integration on the interval [−1, 1] with quadrature
points −1 ≤ q1 , < · · · < qm ≤ 1 and corresponding weights ω1 , · · · , ωm > 0 such that ω1 + · · · +

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

13

ωm = 2. Then the nodal variables on the edge [x j1 −1/2 , x j1 +1/2 ] × {y j2 +1/2 } are collocated at
(x j1 +qk h1 /2 , y j2 +1/2 ), k = 1, · · · , m; and the face-averaged flux is:
Z

1
h1

x j1 +1/2
x j1 −1/2

m

g(w
w(x, y j2 +1/2 , t)dx ≈

1X
w(x j1 +qk h1 /2 , y j2 +1/2 , t)).
ωk g(w
2 k=1

Similarly, the nodal variables on the edge {x j1 +1/2 } × [y j2 −1/2 , y j2 +1/2 ] are collocated at the
points (x j1 +1/2 , y j2 +qk h2 /2 ), k = 1, · · · , m; and the face-averaged flux is:
1
h2

Z

y j2 +1/2

y j2 −1/2

m

w(x j1 +1/2 , y, t)dy ≈
f (w

1X
w(x j1 +1/2 , y j2 +1/2+qk h2 /2 , t)).
ωk f (w
2 k=1

By doing so, we lose the exactness of the cell-boundary fluxes as in the 1D case, hence
the order of this approximation must be chosen carefully without degrading the designed
order of the FD-FV scheme. In particular, we use at least pth -order quadrature rule to
design a FD-FV scheme with the same order.
1,up
2,up
For simplicity, we only consider two 1D FD-FV operators, namely [D x ] and [D x ],
and extend them to 2D problems. The extensions for the other operators, however, follow the same procedure. The common component of these FD-FV methods is the semidiscretized formula for the cell-averages, namely:
(3.3)

 P


Pm
m
w j1 , j2
dw
k=1 ωk f j1 +1/2, j2 ;k − f j1 −1/2, j2 ;k
k=1 ωk g j1 , j2 +1/2;k − g j1 , j2 −1/2;k
+
+
= 0,
dt
2h1
2h2
def

def

where we define x j1 +1/2, j2 ;k == (x j1 +1/2 , y j2 + qk h2 /2), x j1 , j2 +1/2;k == (x j1 + qk h1 /2, y j2 +1/2 ),
and for any such subscript combination x ∗ we denote the corresponding nodal value and
w∗ ) or g (w
w∗ ), respectively. Note that there is no ambiguity
nodal flux by w ∗ = w (xx∗ , t) and f (w
on which index should the last subscript k corresponds to, because it always acts on the
index that is not a half integer. We shall also use the notations x j1 ;k = x j1 + qk h1 /2 and
y j2 ;k = y j2 + qk h2 /2 in subsequent sections.
With (3.3) established, we focus on the semi-discretization formula for nodal values in
the following sections. For simplicity, we restrict ourself first to the case of simple scalar
advection equation:
(3.4)

wt + c1 w x + c2 wy = 0

for some constants c1 > 0 and c2 > 0. Extension to general conservation systems is similar
to that in Section 2.4, and the details will be provided at the end of this section.

1,up

3.1. A second-order scheme based on [D x ]. To construct a second-order FD-FV scheme,
we need at least second-order accurate quadrature rule. In particular, we select m = 1,
q1 = 0, and ω1 = 2. In this case, we have exactly one node on each cell face; hence instead
of using the notations x j1 +1/2, j2 ;k and x j1 , j2 +1/2;k introduced before, we will write x j1 +1/2, j2
and x j1 , j2 +1/2 for simplicity. The same rule holds for the other quantities associated with
the nodal values.

14

X. ZENG

w j1 +1/2, j2 +1

w j1 , j2 +1

w j1 , j2 +1/2
w j1 −1/2, j2

w j1 , j2

w j1 +1/2, j2

w j1 +1/2, j2

w j1 −1, j2 +1/2w j1 , j2 +1/2w j1 +1, j2 +1/2

w j1 , j2 w j1 +1, j2

w j1 , j2 −1/2

w j1 , j2

w j1 +1/2, j2 −1

(a) Stencil for w j1 , j2 .

(b) Stencil for w j1 +1/2, j2 .

(c) Stencil for w j1 , j2 +1/2 .

Figure 3.1. Stencil of the semi-discretization formula (3.5) for the (3.1a)
cell-averaged value, and (3.1b–3.1c) nodal values. ◦ – location of cell
averages; • – location of nodal values.
1,up

The semi-discretized formula for a second-order FD-FV scheme using [D x
by the following:

(3.5)

] is given

dw j1 , j2 c1 w j1 +1/2, j2 − c1 w j1 −1/2, j2 c2 w j1 , j2 +1/2 − c1 w j1 , j2 −1/2
+
+
= 0,
dt
h1
h2
dw j1 +1/2, j2 c1 (2w j1 +1/2, j2 − 2w j1 , j2 ) c2 (w j1 +1/2, j2 − w j1 +1/2, j2 −1 )
+
+
= 0,
dt
h1
h2
dw j1 , j2 +1/2 c1 (w j1 , j2 +1/2 − w j1 −1, j2 +1/2 ) c2 (2w j1 , j2 +1/2 − 2w j1 , j2 )
+
+
= 0.
dt
h1
h2
1,up

In these equations, the one-dimensional FD-FV operator [D x ] is used to discretize the
1,up
spatial derivative ∂ x w j1 +1/2, j2 in (3.5)2 ; and a similar operator [Dy ] is used to discretize
∂y w j1 , j2 +1/2 in (3.5)3 . The other two spatial derivatives in the last two equations of (3.5)
are discretized by a conventional finite-difference method, i.e., the first-order upwind order
1,up
1,up
that is denoted by D x or Dy . These notations differ from the FD-FV operators in the
absence of the bracket. These lead to the simplified semi-discretized formula for nodal
values:
dw j1 +1/2, j2
1,up
1,up
+ c1 [D x w] j1 +1/2, j2 + c2 Dy w j1 +1/2, j2 = 0,
dt
(3.6)
dw j1 , j2 +1/2
1,up
1,up
+ c1 D x w j1 , j2 +1/2 + c2 [Dy w] j1 , j2 +1/2 = 0,
dt
1,up

where the exact formula for D x
1,up

Dy

w j1 +1/2, j2 =

1,up

and Dy

w j1 +1/2, j2 − w j1 +1/2, j2 −1
,
h1

are:
1,up

Dx

w j1 , j2 +1/2 =

w j1 , j2 +1/2 − w j1 −1, j2 +1/2
.
h2

To this end, the stencil of the semi-discretization formula for w j1 , j2 , w j1 +1/2, j2 , and
w j1 , j2 +1/2 given in (3.5) are illustrated in Figure 3.1. Note that we also include the stencils
for situations c1 > 0 and/or c2 > 0. Similar to the one-dimensional case, although we use
discrete differential operators of first-order spatial accuracy, the method (3.5) leads to a
second-order accurate method. To show this, we consider the Cauchy problem of solving
(3.4) on R2 × [0, T ] using a uniform mesh with cell size h1 × h2 . In particular, let the initial
condition be a simple planar wave of the form:
(3.7)

w(x, y, 0) = ei(k1 x+k2 y) ,

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

15

Then we formalize the statement that (3.5) is second-order accurate in space by the following statement:
Conjecture 3.1. Fixing any k1 and k2 , then applying the method (3.5) to solve the preceding
Cauchy problem with exact initial data (3.7) leads to semi-discretized solutions w j1 , j2 (T ),
w j1 +1/2, j2 (T ) and w j1 , j2 +1/2 (T ), which differ from the exact values by O(h2 ).
Here we assume that O(h1 ) = O(h2 ) = O(h) always holds, which is particularly true if
we keep the aspect ratio of each cell the same when refining the mesh. This conjecture
is based on a similar argument to the proof of Theorem 2.8, and more details concerning
why this conjecture could be true is provided in Appendix C. Finally, we combine the
semi-discretization (3.5) and the second-order total variational diminishing Runge-Kutta
scheme in time to obtain a globally second-order accurate method. We omit the linear
stability analysis of the resulting fully-discrete scheme, and mention that the stability limit
given by the Courant-Friedrichs-Lewy condition is:
νmax
, νmax = 1.0.
(3.8)
∆t <
c1 /h1 + c2 /h2
2,up

3.2. A third-order scheme based on [D x ]. Next we consider constructing a third2,up
order FD-FV method by utilizing the operator [D x ]. In this case, the quadrature rule
to evaluate the edge-averaged flux needs to be at least third-order accurate.
√ To this end,
we choose the Gauss-Legendre quadrature points with m = 2, q1,2 = ±1/ 3, and ω1,2 = 1.
Then the stencil for the semi-discretization of a cell-averaged value w j1 , j2 is arranged as
Figure 3.2a, and the formula is:
dw j1 , j2 c1 (w j1 +1/2, j2 ;1 + w j1 +1/2, j2 ;2 ) − c1 (w j1 −1/2, j2 ;1 + w j1 −1/2, j2 ;2 )
+
dt
2h1
(3.9)
c2 (w j1 , j2 +1/2;1 + w j1 , j2 +1/2;2 ) − c2 (w j1 , j2 −1/2;1 + w j1 , j2 −1/2;2 )
+
= 0,
2h2
The semi-discretization formula for the nodal values, however, needs more elaboration.
2,up
For example, considering apply [D x ] to the nodal value at x j1 +1/2, j2 ;2 (Figure 3.2a), we
need the line-averaged data for the dashed line in the figure. In particular, we need an
approximation of the quantity:
Z x j +1/2
1
1
w xj1 , j2 ;2 (t) ≈
w(x, y j2 ;2 , t)dx
h1 x j −1/2
1

to at least third-order accuracy in space. For this purpose, we omit the variable t and define
the variable w xj1 (y) as:
Z x j +1/2
1
def 1
x
w j1 (y) ==
w(x, y)dx,
h1 x j −1/2
1

and our target is to find the approximations w xj1 , j2 ;k ≈ w xj1 (y j2 ;k ), k = 1, 2. Due to our choice
of the quadrature rule to locate the nodal values, w xj1 (y j1 ±1/2 ) can be evaluated to the thirdorder accuracy by:
1
1
w xj1 (y j2 ±1/2 ) = w j1 , j2 ±1/2;1 + w j1 , j2 ±1/2;2 + O(h31 ).
2
2
In addition, seeing that:
Z
Z y j +1/2
1
1
1
w(x, y)dxdy =
w x (y)dy,
h1 h2 C j , j
h2 y j −1/2 j1
1 2

2

16

X. ZENG

w xj1 , j2 ;2
w j1 +1/2, j2 ;2
w j1 +1/2, j2 ;2
w j1 , j2

w j1 , j2 +1/2;1
w j1 , j2 w j1 +1, j2

w j1 +1/2, j2 ;1

w j1 , j2

w j1 , j2 −1/2;1w j1 , j2 −1/2;2
(a) Stencil for w j1 , j2 .

w j1 , j2 +1

(b) Stencil for w j1 +1/2, j2 ;2 .

(c) Stencil for w j1 , j2 +1/2 ; 1.

Figure 3.2. Stencil of the semi-discretization for: (3.2a), the cellaveraged value; and (3.2b–3.2c), the nodal values. ◦ – location of cell
averages; • – location of nodal values;  – reconstructed line-averages.
in which the left hand side is represented by w j1 , j2 , we obtain a third-order accurate formula
for w xj1 , j2 ;k for k = 1, 2, namely:
 2

 3qk qk 1  w j1 , j2 +1/2;1 + w j1 , j2 +1/2;2
x
w j1 , j2 ;k = 
+ − 
+
4
2 4
2
 2



 3qk qk 1  w j1 , j2 −1/2;1 + w j1 , j2 −1/2;2  3 3q2k 
− − 
+  −

 w j , j .
4
2 4
2
2
2  1 2
2,up

This enables us to construct the discrete differential formula [D x w] j1 +1/2, j2 ;k by:

2 
2,up
[D x w] j1 +1/2, j2 ;k =
2w j1 +1/2, j2 ;2 − 3w xj1 , j2 ;k + w j1 −1/2, j2 ;k ,
h1
2,up

and similarly the formula for [Dy w] j1 , j2 +1/2;k , k = 1, 2. Finally, for the other two derivatives ∂y w j1 +1/2, j2 ;k and ∂ x w j1 , j2 +1/2;k , we use the similar strategy as in Section 3.1 and use
the second-order upwind finite-difference operators along the corresponding grid lines. We
2,up
2,up
denote these discrete operators by Dy w j1 +1/2, j2 ;k and D x w j1 , j2 +1/2;k , respectively, and
write the semi-discretized equations for all the nodal values as:
dw j1 +1/2, j2 ;k
2,up
2,up
+ c1 [D x w] j1 +1/2, j2 ;k + c2 Dy w j1 +1/2, j2 ;k = 0,
dt
(3.10)
dw j1 , j2 +1/2;k
2,up
2,up
+ c1 D x w j1 , j2 +1/2;k + c2 [Dy w] j1 , j2 +1/2;k = 0.
dt
The equations (3.9) and (3.10) constitute the whole set of semi-discretization formula for
our third-order FD-FV scheme to solve (3.4); and the stencil for each cell-averaged value
and nodal value (also taking into account of the possibility c1 < 0 and/or c2 < 0) are illustrated in Figure 3.2, in which the square designates the intermediate line-averaged quantities.
In principal, the method presented in this section could be extended to construct FDFV scheme of arbitrary order of accuracy; however, how to proof that the formal order of
accuracy of the resulting FD-FV scheme is one-order higher than the order of the discrete
differential operators [D x,y ] and D x,y that are used to construct the method is not clear.
Nevertheless, we make the conjecture that this is true, and verify it by numerical examples
in next section.

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

17

Finally, we combine the spatial discretization scheme presented with the third-order
total variational diminishing Runge-Kutta scheme in time to obtain a globally third-order
accurate method, and use the Courant-Friedrichs-Lewy condition:
νmax
(3.11)
∆t <
, νmax = 0.306
c1 /h1 + c2 /h2
in the numerical examples.
3.3. Extension to system of equations. The methods in previous sections can be written
in a more general form for scalar problem:
wt + f (w) x + g(w)y = 0
as follows:
P
P
dw j1 , j2
k ωk ( f j1 +1/2, j2 ;k − f j1 −1/2, j2 ;k )
k ωk (g j1 , j2 +1/2;k − g j1 , j2 −1/2;k )
+
+
= 0,
dt
h1
h2
dw j1 +1/2, j2 ;k
(3.12)
+ f 0 (w j1 +1/2, j2 ;k )[D x w] j1 +1/2, j2 ;k + g0 (w j1 +1/2, j2 ;k )Dy w j1 +1/2, j2 ;k = 0,
dt
dw j1 , j2 +1/2;k
+ f 0 (w j1 , j2 +1/2;k )D x w j1 , j2 +1/2;k + g0 (w j1 , j2 +1/2;k )[Dy w] j1 , j2 +1/2;k = 0,
dt
where [D x ] and [Dy ] denote one-dimensional FD-FV operators, and D x and Dy denote
conventional finite-difference operators.
The method (3.12) could be extended to the system (3.1) by using the characteristic
decomposition of the Jacobian matrices:
∂f
= J f = RfΛfLf,
w
∂w

∂gg
= J g = RgΛg Lg.
w
∂w

And we have:
P
P
w j1 , j2
dw
k ωk ( f j1 +1/2, j2 ;k − f j1 −1/2, j2 ;k )
k ωk ( f j1 , j2 +1/2;k − g j1 , j2 −1/2;k )
+
+
= 0,
dt
h1
h2
w j1 +1/2, j2 ;k
dw
L f w)] j1 +1/2, j2 ;k + R g Λ g Dy (L
Lg w) j1 +1/2, j2 ;k = 0,
+ R f Λ f [D x (L
dt
w j1 , j2 +1/2;k
dw
L f w) j1 , j2 +1/2;k + R g Λ g [Dy (L
Lg w)] j1 , j2 +1/2;k = 0,
+ R f Λ f D x (L
dt
where for simplicity we omit the subscripts on all the matrices associated with the Jacobian
matrices.
4. Numerical assessment
We apply the FD-FV methods developed in preceding sections to solve various problems in one and two space dimensions. We focus on problems with solutions that are sufficiently smooth in the paper, but show one simple example with non-convex flux function
and discontinuous initial condition to demonstrate the potential of the FD-FV methods to
capture correct solution structure for general fluxes. The latter issue, however, is not fully
investigated and is left for future work.
4.1. One-dimensional problems. In this section we focus on the smooth problems, such
as the advection equation and Euler equations that admit smooth solutions. We also solve
the Buckley-Leverett equation, an example of non-convex flux, with discontinuous initial
data. Note that although our methods are constructed based on Taylor series expansions,
they perform well in solving this problem and capture the nonlinear waves accurately.

18

X. ZENG

0

10

−2

10

2

−4

10

3

−6

10

4

5

−8

10

−2

10

D x1,up
D x2,up
D x3,up−biased
D x3,up
D x4,up−biased
−1

10

Mesh size (log scale)

L1 error in nodal value (log scale)

L1 error in cell−averaged value (log scale)

0

10

−2

10

2

−4

10

3

−6

10

4

5

−8

10

−2

10

D x1,up
D x2,up
D x3,up−biased
D x3,up
D x4,up−biased
−1

10

Mesh size (log scale)

Figure 4.1. L1 -errors computed by applying the FD-FV methods in Table 2 to solve (4.1): cell averages (left), and nodal values (right).
4.1.1. One-dimensional advection equations. First we consider an advection problem with
periodic boundary conditions:


ut + 2u x = 0,
(x, t) ∈ [−1, 1] × [0, 1],




1

(4.1)
u(x, 0) = 1 + sin(πx), x ∈ [−1, 1],





 u(−1, t) = u(1,2t)
t ∈ [0, 1].
The L1 -norms of the errors for both cell averages and nodal values are computed using the
exact solution at T = 1, and these errors are plotted in logarithmic scales in Figure 4.1. All
the five methods in Table 2 are tested, and the Courant numbers are selected to be 90%
the values of νmax in the last column. For each of the method, we compute the numerical
solutions using four uniform meshes with numbers of cells ranging from 20 to 160.
Next, we consider an advection problem with Dirichlet boundary condition:


ut + u x = 0,
(x, t) ∈ [−0.5, 0.5] × [0, 0.5],






 u(x, 0) = 1 + 1 x3 sin(2πx)1{x≤0} ,
x ∈ [−0.5, 0.5],
(4.2)


2


1
1
1



 u(−0.5, t) = 1 + (t + )3 sin(2π(t + )), t ∈ [0, 0.5],
2
2
2
where 1{x≤0} takes the value 1 if x ≤ 0 or 0 if x > 0. Using the same numerical methods and
the uniform grids with the same number of cells as before, the L1 -norms of the numerical
errors are computed and plotted in logarithmic scales in Figure 4.2. Note that in this
problem, no boundary condition is specified at the right end because the latter only has
3,up−biased
outgoing characteristics. From a numerical point of view, however, applying [D x
]
4,up−biased
and [D x
] at the last node requires data to the right of the boundary, which is not
3,up
available. For this reason, we locally use a purely upwind operator, namely [D x ], at the
last node in these two situations. Because this change of operator is local in space, it does
not affect the overall order of accuracy, especially for the fifth-order scheme constructed
4,up−biased
by [D x
].

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

−2

10

10

−3

10

2
−4

10

−5

10

3
−6

10

−7

4

−8

5

10

10

−9

10 −2
10

D x1,up
D x2,up
D x3,up−biased
D x3,up
D x4,up−biased
−1

10

Mesh size (log scale)

L1 error in nodal values (log scale)

L1 error in cell−averaged value (log scale)

−2

19

−3

10

2
−4

10

−5

10

3
−6

10

−7

4

−8

5

10

10

−9

10

−2

10

D x1,up
D x2,up
D x3,up−biased
D x3,up
D x4,up−biased
−1

10

Mesh size (log scale)

Figure 4.2. L1 -errors computed by applying the FD-FV methods in Table 2 to solve (4.2): cell averages (left), and nodal values (right).

Figures 4.1–4.2 show that the observed orders of accuracy of the FD-FV methods are
consistent with that predicted by Theorem 2.8, at least for linear problems. In particular,
we achieve the desired order of accuracy for both initial value problem and initial boundary
value problem.
Remark 4.1. We do not use any special treatment to enforce the Dirichlet boundary condition at intermediate Runge-Kutta stages to solve (4.2). As pointed out by Mark H. Carpenter et al. [2], this strategy may not be appropriate for high-order Runge-Kutta schemes.
However, there seems to be no fix of this issue in literature for general explicit Runge-Kutta
methods. Thus we do not apply the fix for the fourth-order Runge-Kutta in their work in
our computations; nevertheless our numerical solutions show no degrading of the spatial
orders of accuracy.
4.1.2. One-dimensional Euler equations. Our next example concerns the nonlinear Euler
equations in 1D. In particular, we consider solving the problem (2.21) on (x, t) ∈ [−1, 1] ×
[0, 0.3] with the following smooth initial conditions and periodic boundary conditions:
(4.3)

1
1
1



ρ(x, 0) = 1 + sin(πx), u(x, 0) = 2 + sin(πx), p(x, 0) = 1 + sin(πx), x ∈ [−1, 1];



2
2
2





 ρ(−1, t) = ρ(1, t), u(−1, t) = u(1, t), p(−1, t) = p(1, t),
t ∈ [0, 0.3].
The solutions to this problem at T = 0.3 is sufficiently smooth for our purpose to measure the order of accuracy by applying the FD-FV methods to solve a nonlinear problem.
4,up−biased
In particular, we compute the reference solution using [D x
] and the fifth-order
Runge-Kutta method on a very fine grid with 2,560 uniform cells.
We test the same set of numerical methods to solve this problem, using four successively
refined meshes with the number of uniform cells ranging from 40 to 320. The computed
L1 -norms in density are plotted in logarithmic scales in Figure 4.3; and the convergence

X. ZENG

−2

−2

10

10

−3

10

−4

10

2

−5

10

−6

10

3

−7

10

−8

10

4

−9

10

5

−10

10

−11

10

−2

D x1,up
D x2,up
D x3,up−biased
D x3,up
D x4,up−biased

10

Mesh size (log scale)

−1

10

L1 error in nodal density (log scale)

L1 error in cell−averaged density (log scale)

20

−3

10

−4

10

2

−5

10

−6

10

3

−7

10

−8

10

4

−9

10

5

−10

10

−11

10

−2

D x1,up
D x2,up
D x3,up−biased
D x3,up
D x4,up−biased

10

−1

10

Mesh size (log scale)

Figure 4.3. L1 -errors in density computed by applying the FD-FV methods in Table 2 to solve (2.21) and (4.3): cell-averaged density (left), and
nodal density (right).
curves clearly show that the desired spatial orders are achieved for the nonlinear Euler
equations.
For completeness, we also present the L1 -norms of the errors in velocity and pressure
in Table 3, and they also demonstrate the target orders of convergence.
Remark 4.2. In order to compute the numerical errors in the prime variables u and p, especially for the cell-averaged ones, we need to compute these quantities as functions of
the conservative variables. In particular, denoting the cell-averaged values for the conserw) and p(w
w) are
vative variables by w = [ρ ρu E]T , the “cell-averaged” prime variables u(w
computed using:
!
ρu
1 (ρu)2
w) = , p(w
w) = (γ − 1) E −
u(w
.
ρ
2 ρ
4.1.3. One-dimensional Buckley-Leverett equation. The last 1D example is a scalar con4w2
servation law (2.2) defined by the Buckley-Leverett flux function f (w) = 4w2 +(1−w)
2 [17].
We solve this problem on the domain (x, t) ∈ [−0.5, 2.5] × [0, 0.5], with the initial condition:


0.1 −0.5 < x < 0.0,



0.5 0.0 < x < 1.0,
(4.4)
w(x, 0) = 


 0.2 1.0 < x < 2.0.
The solution to this problem has more complex behavior: There is an inflection point inside
each of the upstream wave and downstream wave that is given by the solution 0 < w0 <
2 f (w )
0
= 0. This leads to the solution structure containing both shock and rarefaction in
1, d dx
2
each of the two waves.
The numerical solutions computed using the preceding numerical methods on a uniform
mesh using 150 cells are plotted in Figures 4.4. In these figures, only the nodal values are
used to plot the solutions, and the cell-averaged values are almost on top of the corresponding solution curves.

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

21

w)
w)
u
u(w
p
p(w
Cells L1 errors Rates L1 errors Rates L1 errors Rates L1 errors Rates
40 6.61e-3
4.72e-3
7.54e-3
5.98e-3
80 1.80e-3 1.87 1.33e-3 1.83 2.01e-3 1.91 1.61e-3 1.90
1,up
[D x ]
160 4.70e-4 1.94 3.50e-4 1.92 5.17e-4 1.96 4.14e-4 1.96
320 1.20e-4 1.97 8.98e-5 1.96 1.31e-4 1.98 1.05e-4 1.98
40 4.05e-4
3.54e-4
4.62e-4
4.10e-4
80 5.18e-5 2.97 4.54e-5 2.96 5.92e-5 2.96 5.23e-5 2.97
2,up
[D x ]
160 5.88e-6 3.14 5.11e-6 3.15 6.68e-6 3.15 5.83e-6 3.17
320 8.21e-7 2.84 7.22e-7 2.82 9.33e-7 2.84 8.26e-7 2.82
40 7.66e-5
6.61e-5
7.85e-5
6.49e-5
80 5.20e-6 3.88 4.46e-6 3.89 5.29e-6 3.89 4.26e-6 3.93
3,up−biased
[D x
]
160 3.36e-7 3.95 2.87e-7 3.96 3.39e-7 3.96 2.72e-7 3.97
320 2.13e-8 3.98 1.82e-8 3.98 2.13e-8 3.99 1.71e-8 3.99
40 4.39e-5
3.63e-5
4.32e-5
3.27e-5
80 2.82e-6 3.96 2.43e-6 3.90 2.82e-6 3.94 2.21e-6 3.88
3,up
[D x ]
160 1.79e-7 3.98 1.54e-7 3.98 1.78e-7 3.98 1.41e-7 3.97
320 1.12e-8 3.99 9.70e-8 3.99 1.12e-8 3.99 8.89e-8 3.99
40 8.12e-6
6.63e-6
8.40e-6
6.51e-6
80
2.94e-7
4.79
2.69e-7
4.62
3.09e-7
4.76
2.64e-7 4.62
4,up−biased
[D x
]
160 9.78e-9 4.91 9.09e-9 4.89 1.02e-8 4.92 9.13e-9 4.86
320 3.15e-10 4.96 2.93e-10 4.95 3.27e-10 4.97 2.94e-10 4.96
Table 3. L1 errors in velocity and pressure computed by applying the
FD-FV methods in Table 2 to solve (2.21) and (4.3).
Method

0.5

[D x1,up]
[D x2,up]
[D x3,up−b]

0.5

w

0.4

w

0.4
0.3

0.3

0.2

0.2

0.1
−0.5

[D x3,up]
[D x4,up−b]

0.1
0

0.5

1

x

1.5

2

2.5

−0.5

0

0.5

1

x

1.5

2

2.5

Figure 4.4. Solutions to the Buckley-Leverett problem with initial con1,up
2,up
3,up−biased
dition (4.4): (left) – [D x ], [D x ], and [D x
]; and (right) –
3,up
4,up−biased
[D x ] and [D x
].
Remark 4.3. Although the methods are constructed using the Taylor series expansions,
they happen to lead to correct solution structures for this particular problem. However, we
emphasize that without further modification to enhance the nonlinear stability, the methods
are in general not suitable to compute discontinuous solutions. For example, applying the
methods in Table 2 to solve the classical Sod’s shock tube problem [24] for Euler equations
leads to spurious oscillations known as Gibbs phenomenon. Nevertheless, we present the
simple problem with Buckley-Leverett flux function in this paper to demonstrate its potential to capture the correct solution structure for non-convex fluxes. A full study of the

22

X. ZENG

nonlinear stability and convergence to the entropy solution of the FD-FV methods, however, is beyond the scope of this paper, and its analysis and enhancement will be addressed
in future work.
4.2. Two-dimensional problems.
4.2.1. Two-dimensional Euler equations. The first 2D problem is the isentropic vortex
advection problem [23]. In particular, we solve the 2D Euler equations:

 
  
 ρ   ρu   ρv 
 ρu   ρu2 + p   ρuv 
 = 0
 + 
  + 
(4.5)
 ρv   ρuv   ρv2 + p 
E t (E + p)u x (E + p)v y
on the domain (x, y, t) ∈ [−5, 5]2 × [0, 10]. Here u and v are velocity components, and the
total energy is determined by E = p/(γ − 1) + ρ(u2 + v2 )/2 with γ ≡ 1.4. We specify periodic
boundary condition on all edges, and set the initial condition to be an isentropic vortex
given by:

1

2

! γ−1

(γ
−
1)


2
2

ρ(x, y, 0) = 1 −
exp 1 − x − y
,



8γπ2


!



y
1



u(x, y, 0) = 1 −
exp (1 − x2 − y2 ) ,



2π
2
!
(4.6)


1

x


2
2

exp
(1
−
x
−
y
)
,
v(x,
y,
0)
=
1
+



2π
2



γ


2


! γ−1

(γ
−
1)

2
2


,
exp 1 − x − y
 p(x, y, 0) = 1 −
8γπ2
where the constant is chosen to be  = 5. The exact solution at one period T = 10 is almost
identical to the initial data, which is used to compute the numerical errors.
We test the two numerical methods in Section 3 using four uniform grids with number
of cells ranging from 202 to 1602 , and summarize the L1 -norms of the numerical errors
in Table 4. From the results we see that the two conjectures on the order of the methods
are confirmed numerically, at least when the mesh is sufficiently refined. Note that the
discussion of Remark 4.2 holds here, which explains the computation of u(w), v(w), and
p(w) in Table 4 from the computed cell-averaged conservative variables.
4.2.2. Two-dimensional linear elasticity problem with orthotropic materials. In the last
example, we compute the wave propagation in an anisotropic medium that is modeled
by linear elastic orthotropic material. This is a linear problem that describes the solid
dynamics of linear elastic material with infinitesimal strain assumption and the plane strain
assumption. The governing equation is given by:
 


 
 σ11   C11 v1   C12 v2 
 σ   C21 v1   C22 v2 
 

 22  
(4.7)
 σ12  −  C33 v2  −  C33 v1  = 0,
 v   ρ−1 σ   ρ−1 σ 
11 
0
 1   −1
  0 12 
v2 t ρ0 σ12 x ρ−1
0 σ22 y
where ρ0 is a constant with the usual meaning of the density of the material, v1 and v2
denotes the velocity in x- and y-directions, respectively; and σi j and Ci j are components
of the stress tensor and the stiffness matrix (in Voigt notation), respectively.

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

23

Method of Section 3.1
Method of Section 3.2
Grid sizes L1 errors Rates L1 errors Rates L1 errors Rates L1 errors Rates
ρ
ρ
ρ
ρ
20 × 20 1.94e+0
1.87e+0
2.50e-1
2.38e-1
40 × 40 8.48e-1 1.19 8.12e-1 1.20 4.52e-2 2.47 4.34e-2 2.46
80 × 80 2.56e-1 1.73 2.41e-1 1.75 6.33e-3 2.83 6.10e-3 2.83
160 × 160 6.98e-2 1.88 6.26e-2 1.95 8.24e-4 2.94 7.83e-4 2.96
w)
w)
u
u(w
u
u(w
20 × 20 6.12e+0
6.36e+0
7.29e-1
6.93e-1
40 × 40 2.74e+0 1.16 2.79e+0 1.19 1.30e-1 2.49 1.29e-1 2.43
80 × 80 8.25e-1 1.73 8.19e-1 1.77 1.81e-2 2.84 1.87e-2 2.78
160 × 160 2.18e-1 1.92 2.14e-1 1.94 2.31e-3 2.97 2.50e-3 2.90
w)
w)
v
v(w
v
v(w
20 × 20 5.46e+0
5.76e+0
6.06e-1
5.70e-1
40 × 40 2.06e+0 1.40 2.08e+0 1.47 1.11e-1 2.45 1.10e-1 2.37
80 × 80 5.92e-1 1.80 5.83e-1 1.84 1.56e-2 2.83 1.60e-2 2.78
160 × 160 1.55e-1 1.93 1.50e-1 1.97 2.01e-3 2.96 2.16e-3 2.89
w)
w)
p
p(w
p
p(w
20 × 20 2.46e+0
2.41e+0
3.65e-1
3.53e-1
40 × 40 1.09e+0 1.17 1.03e+0 1.23 5.73e-2 2.67 5.71e-2 2.63
80 × 80 3.26e-1 1.74 2.91e-1 1.82 7.88e-3 2.86 8.08e-3 2.82
160 × 160 8.38e-2 1.89 7.56e-2 1.95 1.02e-3 2.94 1.21e-3 2.74
Table 4. L1 errors of solving the isentropic vortex problem by the two
methods of Section 3.

In particular, to derive (4.7), we start with the following classical equation for linear
elasticity with plane strain:
(4.8)

ρ0 u tt = ∇ · σ ,

where u = [u1 , u2 ]T is the displacement vector, and σ is the symmetric Cauchy stress tensor
that is given by:
"
#
σ σ
σ = 11 12 , σ12 = σ21 .
σ21 σ22
The stress tensor σ is connected to the strain tensor ∇uu by the constitutive equation, which
is given in Voigt notation by:

 


 σ11   C11 C12 0   u1,1 

 
 

(4.9)
 σ22  =  C21 C22 0   u2,2  ,
σ12
0 0 C33 u1,2 + u2,1
where we use ,1 and ,2 to denote ,x and ,y , respectively.
Then the hyperbolic equations (4.7) are obtained from (4.8) and (4.9) by choosing
σ11 , σ22 , σ12 , and v1 ≡ du1 /dt, v2 ≡ du2 /dt as the variables. Note that (4.7) is a linear
system, hence solving (4.7) is equivalent to solving the conservation law obtained from the
classical momentum equation.

24

X. ZENG

(a) Finite volume method.

(b) 2nd -order FD-FV method.

Figure 4.5. The velocity magnitude ||vv|| =

q

(c) 3rd -order FD-FV method.

v21 + v22 computed on a uni-

form 2012 mesh.
In this test, we choose ρ0 = 3.0 and the following stiffness matrix:


 
 C11 C12 0   2.0 0.99 0 
 C21 C22 0  =  0.99 0.5 0  ,


 
0 0 C33
0 0 2.0
which is only slightly positively definite and creates a strong anisotropic effect. We solve
the problem on the rectangular domain (x, y, t) ∈ [−1, 1]2 × [0, 0.7] with all initial conditions
set to zero.
To create the wave in the medium, we perturb the latter using a point pulse in the
x-momentum equation by adding a Ricker wavelet s(t) to the right hand side of the second
last equation of (4.7):
2 2 2
s(t) = (1 − 2π2 f 2 tˆ2 )e−π f tˆ ,
where f = 4.0 is the peak frequency and tˆ = t − (2π2 f 2 )−1/2 so that the pulse has zero
strength at t = 0.0. At T = 0.7, the fastest wave has not reached the domain boundary,
which allows the homogeneous Dirichlet numerical boundary condition for all variables in
our computation. This can be verified by the reference solution computed using a secondorder finite volume method [27] on a uniform 2012 mesh. In particular, we use the Rusanov
flux [21], the van Albada limiter [25], and the second-order total variational diminishing
Runge-Kutta method to compute the reference solution, which is plotted in Figure 4.5a.
The numerical solutions computed by the two FD-FV methods are plotted in Figures 4.5b–
4.5c, in which the cell-averaged data are used to make the figures. Clearly, the results
confirm that the FD-FV methods compute solutions with better resolution comparing to
that by the finite volume method on the same mesh. This observation is in consistent with
our prediction in that: (1) on the same mesh the FD-FV methods use more variables than
the finite volume method, and (2) the FD-FV method in Section 3.2 is formally one-order
higher than the second-order finite volume method.
5. Conclusions
We formulate a computational framework for the discretization of general first-order
hyperbolic conservation law in space using both cell-averages and nodal values. In the
one-dimensional case, we prove that the formal order of spatial accuracy of the proposed
FD-FV method is in general one-order higher than the order of the differential operator
involved in the construction of the method. This is unlike the conventional finite difference

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

25

or finite volume methods, in which case the order of the scheme is typically the same as
that of the discrete differential operator. We extend the methodology to Cartesian grids
in two space dimensions, and provide strong evidence that this property of spatial order
accuracy carries to 2D smooth problems. All the 1D and 2D methods are equipped with
explicit time-integrators with appropriate order, and linear stability analysis shows that
these methods are linearly stable when a Courant-Friedrich-Lewy condition is satisfied.
The numerical performance of the proposed methods is assessed by solving various
smooth problems in one and two space dimensions, and the theoretical orders of accuracy
are confirmed by the numerical tests. These problems include both linear problems, such
as scalar advection equation and solid dynamic equation for linear elastic orthotropic materials, and nonlinear problems, such as the 1D and 2D Euler equations. We also show that
the method is able to capture the correct solution structure for a problem with non-convex
Buckley-Leverett flux function, at least in the specific problem that is considered. A full investigation of enhancing the nonlinear stability and the convergence to the entropy solution
of the proposed FD-FV methods is work in progress and will be left for future publication;
nevertheless this Buckley-Leverett example shows the potential of the FD-FV methods for
computing the correct solutions for hyperbolic problems with general flux functions.

References
[1] Miquel Aguirre, Antonio J. Gil, Javier Bonet, and Aurelio Arranz Carre no. A vertex centred finite volume
jamesonschmidtturkel (jst) algorithm for a mixed conservation formulation in solid dynamics. J. Comput.
Phys., 259:672–699, February 2014.
[2] Mark H. Carpenter, David Gottlieb, Saul Abarbanel, and Wai-Sun Don. The theoretical accuracy of rungekutta time discretizations for initial boundary value problems: a study of the boundary error. SIAM J. Sci.
Comput., 16(6):1241–1252, November 1995.
[3] Sin-Chung Chang. The method of space-time conservation element and solution element – a new approach
for solving the navier-stokes and euler equations. J. Comput. Phys., 119(2):295–324, July 1995.
[4] Sin-Chung Chang, Xiao-Yen Wang, and Chuen-Yen Chow. The space-time conservation element and solution element method: A new high-resolution and genuinely multidimensional paradigm for solving conservation laws. J. Comput. Phys., 156(1):89–136, November 1999.
[5] Bernado Cockburn and Chi-Wang Shu. The runge-kutta discontinuous galerkin method of conservation laws
v: multidimensional systems. J. Comput. Phys., 141(2):199–224, April 1998.
[6] Constantine M. Dafermos. Hyperbolic Conservation Laws in Continuum Physics, volume 325 of
Grundlehren der mathematischen Wissenschaften. Springer, 3 edition, 2010.
[7] Michael Dumbser, Dinshaw S. Balsara, Eleuterio F. Toro, and Claus-Dieter Munz. A unified framework for
the construction of one-step finite volume and discontinuous galerkin schemes on unstructured meshes. J.
Comput. Phys., 227(18):8209–8253, September 2008.
[8] Timothy Eymann and Philip Roe. Active flux schemes. In 49th AIAA Aerospace Sciences Meeting
including the New Horizons Forum and Aerospace Exposition, January 2011. Orlando, Florida.
[9] Sigal Gottlieb and Chi-Wang Shu. Total variational diminishing runge-kutta schemes. Math. Comput.,
67(221):73–85, January 1998.
[10] Bertil Gustafsson. High Order Difference Methods for Time Dependent PDE, volume 38 of Springer Series
in Computational Mathematics. Springer, 1st edition, 2008.
[11] Ernst Hairer, Syvert P. Nørsett, and Gerhard Wanner. Solving Ordinary Differential Equations I: Nonstiff
Problems, volume 8 of Springer Series in Computational Mathematics. Springer, 2 edition, 1993.
[12] H. T. Huynh. A piecewise-parabolic dual-mesh method for the euler equations. In 12th AIAA Computational
Fluid Dynamics Conference, June 1995. San Diego, California.
[13] Yohsuke Imai, Takayuki Aoki, and Kenji Takizawa. Conservative form of interpolated differential operator
scheme for compressible and incompressible fluid dynamics. J. Comput. Phys., 227(4):2263–2285, February
2008.
[14] Eugene Isaacson and Herbert Bishop Keller. Analysis of Numerical Methods. Dover Publication, reprint
edition, June 1994.
[15] A. Iserles. Generalized leapfrog methods. IMA J. Numer. Anal., 6(4):381–392, October 1986.

26

X. ZENG

[16] Antony Jameson. Positive schemes and shock modelling for compressible flows. Int. J. Numer. Meth. Fl.,
20(8–9):743–776, April–May 1995.
[17] S. A. Karabasov and V. M. Goloviznin. Compact accurately boundary-adjusting high-resolution technique
for fluid dynamics. J. Comput. Phys., 228(19):7426–7451, October 2009.
[18] Grady I. Lemoine, M. Yvonne Ou, and Randall J. LeVeque. High-resolution finite volume modeling of wave
propagation in orthotropic poroelastic media. SIAM J. Sci. Comput., 35(1):B176–B206, 2013.
[19] Randall LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge Texts in Applied Mathematics. Cambridge University Press, 2002.
[20] H. A. Luther. Further explicit fifth-order runge-kutta formulas. SIAM Rev., 8(3):374–380, July 1966.
[21] V. V. Rusanov. The calculation of interaction of non-steady shock waves with obstacles. USSR Comput.
Math. Math. Phys., 1(2), 1962. 304–320.
[22] Richard Sanders and Alan Weiser. High resolution staggered mesh approach for nonlinear hyperbolic systems of conservation laws. J. Comput. Phys., 101(2):314–329, August 1992.
[23] Chi-Wang Shu. High-order finite difference and finite volume weno schemes and discontinuous galerkin
methods for cfd. Int. J. Comput. Fluid D., 17(2):107–118, 2003.
[24] Gary A. Sod. A survey of several finite difference methods for systems of nonlinear hyperbolic conservation
laws. J. Comput. Phys., 27(1):1–31, April 1978. Review.
[25] G. D. van Albada, B. van Leer, and Jr. W. W. Roberts. A comparative study of computational methods in
cosmic gas dynamics. Astron. Astrophys., 108(1):76–84, April 1982.
[26] Bram van Leer. Towards the ultimate conservative difference scheme iv. a new approach to numerical convection. J. Comput. Phys., 23(3):276–299, March 1977.
[27] Bram van Leer. Towards the ultimate conservative difference scheme v. a second-order sequel to godunov’s
method. J. Comput. Phys., 32(1):101–136, July 1979.
[28] Feng Xiao, Akio Ikebata, and Takashi Hasegawa. Numerical simulations of free-interface fluids by a multiintegrated momentmethod. Comput. Struct., 83(6–7):409–423, February 2005.
[29] Zeng-Chan Zhang, S. T. John Yu, and Sin-Chung Chang. A space-time conservation element and solution
element method for solving the two- and three-dimensional unsteady euler equations using quadrilateral and
hexahedral meshes. J. Comput. Phys., 175(1):168–199, January 2002.

Appendix A. The connection of the FD-FV operator and the Hermite interpolation
polynomials
We sketch here a proof of the existence of the constants αl and βl that satisfy the condition (2.11) providing p ≤ 4q, i.e., the stencil of the operator [D x ] is large enough. Furthermore, if p < 4q we can simply add more equations to (2.11), so that we obtain a system of
4q + 1 equations for 4q + 1 coefficients of the form:
q
q
X
X
lm
lm+1 − (l − 1)m+1
αl +
βl = γm h−m , m = 0, · · · 4q,
(A.1)
(m
+
1)!
m!
l=−q
l=−q+1
where γm are arbitrary constants. Instead of showing the matrix on the left hand side of
(A.1) is nonsingular, we link it to Hermite interpolation polynomials and show that the
solutions to this system exist for any right hand sides.
To this end, we consider the 2q + 1 cell faces in the stencil: x j+1/2−q < x j+1/2−(q−1) <
· · · < x j+1/2 < x j+1/2+1 < · · · < x j+1/2+q , and construct the 4q + 2 “basis” polynomials of
degree 4q + 1 using the Hermite interpolation theory [14]. We denote these polynomials
by Φ−q , Φ−q+1 , · · · , Φq−1 , Φq and Ψ−q , Ψ−q+1 , · · · , Ψq−1 , Ψq , and they are determined by the
conditions:
Φl (x j+1/2+k ) = δlk ,

Φ0l (x j+1/2+k ) = 0;

Ψl (x j+1/2+k ) = 0,

Ψ0l (x j+1/2+k ) = δlk

for all −q ≤ k, l ≤ q, where δlk is the Kronecker delta. Next, we define 4q + 1 polynomials
of degree 4q by:
q
1X 0
Φ (x), −q + 1 ≤ l ≤ q;
ψl (x) = Ψ0l (x), −q ≤ l ≤ q.
φl (x) =
h k=l k

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

Then it can be easily verified that:
(A.2)
Z
1
φl (x)dx = δlk , −q + 1 ≤ l, k ≤ q,
h C j+k
and
(A.3)
Z
1
ψl (x)dx = 0,
h C j+k

φl (x j+1/2+k ) = 0,

−q ≤ l ≤ q, −q + 1 ≤ k ≤ q,

27

−q + 1 ≤ l ≤ q, −q ≤ k ≤ q;

ψl (x j+1/2+k ) = δlk ,

−q ≤ l, k ≤ q.

Then we will verify that:
αl =

4q
X

γm ∂mx φl (x),

−q + 1 ≤ l ≤ q,

βl =

and

m=0

4q
X

γm ∂mx ψl (x),

−q ≤ l ≤ q

m=0

satisfy the system (A.1). Indeed, consider any polynomial p(x) of degree no larger than
4q, it may be uniquely written as:
q
q
X
X
p(x) =
cl φl (x) +
dl ψl (x).
l=−q+1

l=−q

Using (A.2–A.3), we obtain:
Z
1
p j+l =
p(x)dx = cl , −q + 1 ≤ l ≤ q;
h C j+l
Thus:
q
X
l=−q+1

αl p j+l +

q
X

βl p j+1/2+l =

l=−q

=

q
X

4q
X

γm ∂mx φl (x) +

q
X

dl

4q
X

−q ≤ l ≤ q.

γm ∂mx ψl (x)

l=−q m=0 
q
q
4q
 X
 X
X

m


γm ∂ x 
cl φl (x) +
dl ψl (x) =
γm ∂mx p(x).

l=−q+1
4q
X
m=0

cl

p j+1/2+l = p(x j+1/2+l ) = dl ,

m=0


l=−q+1

l=−q

m=0

Because this equality holds for any polynomial p(x) of degree equal to or less than 4q, by
applying the Taylor series expansions to the terms on the left hand side, we see that (A.1)
holds.
Remark A.1. The preceding procedure may be modified slightly for the case of nonuniform
grids, which then establish the general existence result for the FD-FV operator [D x ] of any
order of accuracy.
Appendix B. Proof of Theorem 2.8
Rewriting the system (2.17) in matrix form:
" #
" #
"
#
d A
A
0 i
def
+ ikC
= 0, C = C(θ) ==
,
N
a(θ) b(θ)
dt N
and using the initial condition A(0) = N(0) = 1, we may solve the ODE exactly to obtain:
"
#
" #
A(t)
1
(B.1)
= e−ckC(θ)t
.
N(t)
1
Let the two eigenvalues of C(θ) be:
!
q
1
2
(B.2) λ1 (θ) = b(θ) + b(θ) + 4ia(θ)
2

and

1
λ2 (θ) = b(θ) −
2

!
q
2
b(θ) + 4ia(θ) ,

28

X. ZENG

respectively, we compute from (B.1):
λ2 − i −ckλ1 t λ1 − i −ckλ2 t
e
−
e
,
λ2 − λ1
λ2 − λ1
λ1 (λ2 − i) −ckλ1 t λ2 (λ1 − i) −ckλ2 t
e
−
e
.
N(t) =
i(λ2 − λ1 )
i(λ2 − λ1 )
A(t) =

(B.3)

Next we compute the Laurent series of λ1,2 about θ = 0 (as we will see soon, if b0 =
β = 0, then the eigenvalues have no pole at 0 and we have the Taylor series expansion
l=−q l
as in conventional finite difference analysis). For this purpose, we define am and bm by:
Pq

q
X
1
(lm+1 − (l − 1)m+1 )αl ,
(m + 1)! l=−q+1

am =

bm =

q
1 X m
l βl ,
m! l=−q

which simplifies the expression for a(·) and b(·) to:
∞

a(θ) =

a0 X m+1
+
i am+1 θm
θ m=0

∞

b(θ) =

and

b0 X m+1
+
i bm+1 θm .
θ m=0

By the assumption that [D x ] is pth -order accurate and Theorem 2.5, we have:
a(θ) + b(θ) = i + r(θ)θ ,
p

where

def

r(θ) ==

∞
X

im+p+1 (am+p+1 + bm+p+1 )θm .

m=0

Then we are ready to compute the square root in (B.2). Because 2i − b = −b0 /θ + O(1) , 0
for small θ , 0, we have:
!2
2irθ p
4r2 θ2p
2
2
p
2
p
,
b + 4ia = b + 4i(i − b + rθ ) = (2i − b) + 4irθ = 2i − b +
+
2i − b
(2i − b)2
in which there is the estimate:
2irθ p 2i p+2 c p θ p + O(θ p+1 ) 2i p+2 c p θ p+1 + O(θ p+2 )
=
=
= O(θ p+1 ).
2i − b
2i − b0 /θ + O(1)
−b0 + O(θ)
It follows that:
p
2i p+2 c p p+1
2irθ p
b2 + 4ia = 2i − b +
+ O(θ2p+3 ) = 2i − b −
θ + O(θ p+2 ).
2i − b
b0
Note that we make the particular choice between the two possible roots so that using (B.2),
λ1 → i as θ → 0. Then the two eigenvalues are given by:
λ1 = i −

i p+2 c p p+1
θ + O(θ p+2 )
b0

and

λ2 =

b0
+ (b1 − 1)i + O(θ).
θ

Plugging these expressions into the first terms in the right hand sides of (B.3) leads to:


i p+2 cc p kt p+1
λ2 − i −ckλ1 t
−ickt 
p+2 
e
=e
θ + O(θ ) ,
1 +
λ2 − λ1
b0


(B.4)


i p+1 c p (ickt − 1) p+1
λ1 (λ2 − i) −ckλ1 t
e
= e−ickt 1 +
θ + O(θ p+2 ) ;
i(λ2 − λ1 )
b0
whereas plugging them into the second terms shows that:
(B.5)

e−ckλ2 t = e−

cbo t
h

e−ick(b1 −1)t (1 + O(θ)),

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

29

which decays to zero exponentially as t → ∞ or h → 0 provided that cb0 > 0 (we suppose
t, h > 0), otherwise this term grows without bound as the mesh is refined. Hence we see
that if cb0 > 0, as stated in the theorem, the following equations hold:


i p+2 cc p kt p+1
−ickt 
p+2 
A(t) = e
θ + O(θ ) ,
1 +
b0


i p+1 c p (ickt − 1) p+1
p+2 
−ickt 
θ + O(θ ) ;
N(t) = e
1 +
b0
and by Definition 2.6, the semi-discretized FD-FV method is at least (p + 1)th -order accurate, which completes the proof.
Remark B.1. By defining z1,2 = θλ1,2 (θ), it is not difficult to see that both z1 and z2 are
symmetric functions of θ about the value θ = π. Thus with our preceding choice of λ1,2 ,
λ1 has significant effect as the mesh is refined, i.e., θ → 0; whereas the other eigenvalue
λ2 has more significant effect when the mesh is approaching the under-resolved limit, i.e.,
θ → 2π or one cell per wave length.
Appendix C. Evidences for Conjecture 3.1
Similar to the analysis in one space dimension, we use the superscript ∗ to denote the
exact solutions. Suppose x j1 +1/2 = ( j1 + 1/2)h1 and y j2 +1/2 = ( j2 + 1/2)h2 , we have the
exact solutions at t = T :



1
w∗j1 , j2 = 2
e−i(c1 k1 +c2 k2 )T ei( j1 θ1 + j2 θ2 ) eiθ1 /2 − e−iθ1 /2 eiθ2 /2 − e−iθ2 /2 ,
i θ1 θ2
w∗j1 +1/2, j2 = e−i(c1 k1 +c2 k2 )T ei(( j1 +1/2)θ1 + j2 θ2 ) ,

w∗j1 , j2 +1/2 = e−i(c1 k1 +c2 k2 )T ei( j1 θ1 +( j2 +1/2)θ2 ) ,

where θ1 = k1 h1 and θ2 = k2 h2 .
It is not difficult to see that providing exact initial data at t = 0, the solutions to the
semi-discretized problem are



1
w j1 , j2 (t) = 2
A(t)ei( j1 θ1 + j2 θ2 ) eiθ1 /2 − e−iθ1 /2 eiθ2 /2 − e−iθ2 /2 ,
i θ1 θ2
w j1 +1/2, j2 (t) = N1 (t)ei(( j1 +1/2)θ1 + j2 θ2 ) ,

w j1 , j2 +1/2 (t) = N2 (t)ei( j1 θ1 +( j2 +1/2)θ2 ) ,

where A(t), N1 (t), and N2 (t) solve the ODE system:
(C.1)
  
 
−ic2 k2 τ−1
0
−ic1 k1 τ−1
 A  
  A 
1
2
d   
  
 N1  =  2c1 k1 τ1 τ2 θ1−1 e−iθ1 /2 −2c1 k1 θ1−1 − ic2 k2 τ2 e−iθ2 /2
0
  N1  ,
dt  N  
−1
−iθ1 /2   N 
−1 −iθ2 /2
e
0
−2c
k
θ
−
ic
k
τ
e
2c
k
τ
τ
θ
2
2
2 2 2
1 1 1
2 2 1 2 2
with τ j = (eiθ j /2 − e−iθ j /2 )/(iθ j ) = 1 + O(θ2j ), j = 1, 2. One way to compute the three eigenvalues associated with (C.1) is directly finding the roots of the cubic characteristic polynomial, just as what we have done in Appendix B. However, this method does not extend to
more general cases when the ODE system is larger; thus we do not perform the calculation
here but simply make the following conjecture and will verify it by numerical examples
afterwards:
2c j k j
+ O(1), j = 1, 2.
(C.2)
λ0 = −i(c1 k1 + c2 k2 ) + O(θ2 ), and λ j = −
θj
In particular, we consider three sets of parameters:

30

X. ZENG

7

0

Set (a)
Set (b)
Set (c)

(ℑλ0 + c 1 k 1 + c 2 k 2 )/θ22

6

−0.5

ℜλ0 /θ22

5
4
3

−1

Set (a)
Set (b)
Set (c)

−1.5

2
−2

1
0
0

0.05

(a)

0.1

θ2

=λ0 +c1 k1 +c2 k2
θ22

0.15

−2.5
0

0.2

0.05

vs. θ2 .

(b)

0.1

θ2

<λ0
θ22

0.15

0.2

vs. θ2 .

Figure C.1. The (C.1a) imaginary and (C.1b) real parts of the first eigenvalue λ0 as θ2 → 0.
1

θ1

−1.9

θ1

−1

ℜλj θj /(c j k j )

ℑλj θj /(c j k j )

0

−2

−1.8

θ1
θ2
θ2

Set (a)
Set (b)
Set (c)

−3

Set (a)
Set (b)
Set (c)

θ2

θ1

−2

θ2
θ1

−2.1
θ2

−4

θ2
θ1

−5
0

0.05

(a)

=λ j
c jk j

0.1

θ2

0.15

0.2

−2.2
0

vs. θ2 , j = 1, 2.

0.05

(b)

<λ j θ j
c jk j

0.1

θ2

0.15

0.2

vs. θ2 , j = 1, 2.

Figure C.2. The (C.2a) imaginary and (C.2b) real parts of the other
eigenvalue λ1,2 as θ2 → 0.
(a) c1 k1 = 1.0, c2 k2 = 1.0, θ1 /θ2 = 5.0;
(b) c1 k1 = 5.0, c2 k2 = 1.0, θ1 /θ2 = 1.0;
(c) c1 k1 = 1.0, c2 k2 = 1.2, θ1 /θ2 = 1.0.
For each set of the parameters, we plot the following curves to demonstrate (C.2): (1)
(=λ0 + c1 k1 + c2 k2 )/θ22 vs. θ2 , (2) (<λ0 )/θ22 vs. θ2 , (3) (<λ j )θ j /(c j k j ) vs. θ2 , j = 1, 2, and
(4) (=λ j )/(c j k j ) vs. θ2 , j = 1, 2. For the conjecture to be valid, all curves should converges
to a finite value as θ2 → 0; and in the case of (3), this constant value should be −2 for
j = 1, 2. The asymptotic behaviors of the imaginary and real parts of λ0 are plotted in
Figure C.1; and those of the other two eigenvalues λ1,2 are plotted in Figure C.2. Note
that it appears from Figure C.1b that <λ0 could be o(θ2 ) instead of O(θ2 ); nevertheless we
stick to the latter because it is sufficient for our purpose.

FD-FV FOR FIRST-ORDER HYPERBOLIC CONSERVATION LAWS

31

Supposing that (C.2) is true, it remains to show that the solutions of (C.1) are secondorder in θ1,2 . This can be done, by explicitly constructing the formula for A(t) and N1,2 (t)
similar to (B.3). We will only do this for A(t), and the procedure for the other quantities follows similarly. Denoting the matrix in (C.1) by C, and we compute its eigenvalue
decomposition as:
C = RDL,
where D is the diagonal matrix with diagonal entries λ0 , λ1 , λ2 ; and the other two matrices
are given by:


1
1
1


2b1 τ2
2b1 τ2
2b1 τ2


R =  θ1 λ0 +2a1 +iθ1 b2 θ1 λ1 +2a1 +iθ1 b2 θ1 λ2 +2a1 +iθ1 b2  ,


2b2 τ1
2b2 τ1
2b2 τ1
θ2 λ0 +2a2 +iθ2 b1 θ2 λ1 +2a2 +iθ2 b1 θ2 λ2 +2a2 +iθ2 b1
and
 1

−ia2 /τ1
−ia1 /τ2


 d0 d0 (λ0 +2a1 /θ1 +ib2 ) d0 (λ0 +2a2 /θ2 +ib1 ) 
−i
a
/τ
−i
a
/τ
1
L =  d d (λ +2a1 /θ2 +ib ) d (λ +2a2 /θ1 +ib )  ,
 1 1 1 1 1 2 1 1 2 2 1 
−ia1 /τ2
−ia2 /τ1
 1

d2 d2 (λ2 +2a1 /θ1 +ib2 ) d2 (λ2 +2a2 /θ2 +ib1 )
where a j = c j k j , b j = a j τ j e−iθ j /2 = a j (1 + O(θ j )) for j = 1, 2, and:
2ia2 b2 θ2
2ia1 b1 θ1
dj = 1−
−
, j = 0, 1, 2.
(λ j θ1 + 2a1 + ib2 θ1 )2 (λ j θ2 + 2a2 + ib1 θ2 )2
Using (C.2), we have:
d0 = 1 −

2ia21 θ1 (1 + O(θ1 ))
(2a1

Thus if we write A(t) =

+ O(θ))2

−

2ia22 θ2 (1 + O(θ2 ))
+ O(θ))2

(2a2
λ
t
λ
t
0
1
A0 e + A1 e + A2 eλ2 t ,

1
1
= 1 − iθ1 − iθ2 + O(θ2 ).
2
2

there is:


!
ia1 θ1 (1 + O(θ22 )) ia2 θ2 (1 + O(θ12 )) 
1
ia1 θ1 /τ2
1 
ia2 θ2 /τ1

A0 =
1−
=
−
−
1 −
d0
λ0 θ1 + 2a1 + ib2 θ1 ! λ0 θ2 + 2a2 + ib1 θ2
d0 
2a1 + O(θ)
2a2 + O(θ) 
1
1
1
=
1 − iθ1 − iθ2 + O(θ2 ) = 1 + O(θ2 ),
d0
2
2

which confirms that A0 eλ0 t = e−i(c1 k1 +c2 k2 )t + O(θ2 ). For the remaining two terms A j eλ j t ,
j = 1, 2, we note that A j grows (at most) at polynomial rate of 1/θ; whereas by (C.2), eλ j t
decays exponentially in 1/θ j providing that c j > 0, which is assumed. Hence we conclude
that A(t) = e−i(c1 k1 +c2 k2 )t + O(θ2 ), and similar result holds for N1 (t) and N2 (t).
Department of Mathematical Sciences, University of Texas at El Paso, El Paso TX 79902, United States.,
Tel.: +1-915-747-6759
E-mail address, X. Zeng: xzeng@utep.edu

