Sin-Chung Chang

Glenn Research Center, Cleveland, Ohio

A New Approach for Constructing 

Highly Stable High Order CESE Schemes 

NASA/TM—2010-216766

August 2010

AIAA–2010–543
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to help@

sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Telephone the NASA STI Help Desk at

 443–757–5802

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
Sin-Chung Chang

Glenn Research Center, Cleveland, Ohio

A New Approach for Constructing 

Highly Stable High Order CESE Schemes 

NASA/TM—2010-216766

August 2010

AIAA–2010–543

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

48th Aerospace Sciences Meeting

sponsored by the American Institute of Aeronautics and Astronautics

Orlando, Florida, January 4–7, 2010
Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

Available electronically at http://gltrs.grc.nasa.gov

This work was sponsored by the Fundamental Aeronautics Program 

at the NASA Glenn Research Center.

Level of Review: This material has been technically reviewed by technical management. 
A NEW APPROACH FOR CONSTRUCTING HIGHLY STABLE HIGH ORDER CESE SCHEMES

Sin-Chung Chang

NASA Glenn Research Center, Cleveland, OH 44135

E-mail: sin-chung.chang@nasa.gov

Abstract

A new approach is devised to construct high order CESE schemes which would avoid the common

shortcomings of traditional high order schemes including: (a) susceptibility to computational instabilities;

(b) computational inefficiency due to their local implicit nature (i.e., at each mesh points, need to solve a

system of linear/nonlinear equations involving all the mesh variables associated with this mesh point); (c)

use of large and elaborate stencils which complicates boundary treatments and also makes efficient parallel

computing much harder; (d) difficulties in applications involving complex geometries; and (e) use of problemspecific techniques which are needed to overcome stability problems but often cause undesirable side effects.

In fact it will be shown that, with the aid of a conceptual leap, one can build from a given 2nd-order CESE

scheme its 4th-, 6th-, 8th-,... order versions which have the same stencil and same stability conditions

of the 2nd-order scheme, and also retain all other advantages of the latter scheme. A sketch of

multidimensional extensions will also be provided.

1. Introduction

The space-time conservation element and solution element (CESE) method is a high-resolution and

genuinely multidimensional method for solving conservation laws [1–63]. Its nontraditional features include:

(i) a unified treatment of space and time; (ii) the introduction of conservation elements (CEs) and solution

elements (SEs) as the vehicles for enforcing space-time flux conservation; (iii) a time marching strategy

that has a space-time staggered stencil at its core and, as such, fluxes at an interface can be evaluated

without using any interpolation or extrapolation procedure (which, in turn, leads to the method’s ability

to capture shocks without using Riemann solvers); (iv) the requirement that each scheme be built from a

nondissipative core scheme and, as a result, the numerical dissipation can be controlled effectively; and (v)

the fact that mesh values of the physical dependent variables and their spatial derivatives are considered as

independent mesh variables to be solve for simultaneously. Note that CEs are non-overlapping space-time

subdomains introduced such that (i) the computational domain can be filled by these subdomains; and (ii)

flux conservation can be enforced over each of them and also over the union of any combination of them. On

the other hand, SEs are space-time subdomains introduced such that (i) the boundary of each CE can be

divided into component parts with each of them belonging to a unique SE; and (ii) within a SE, any physical

flux vector is approximated using simple smooth functions. In general, a CE does not coincide with a SE.

Without using flux-splitting, dimensional-splitting, mesh-alignment or other special techniques, since its

inception in 1991 [1], the unstructured-mesh compatible 2nd-order CESE method has been used to obtain

numerous accurate 1D, 2D and 3D steady and unsteady flow solutions with Mach numbers ranging from

0.0028 to 10 [15]. The physical phenomena modeled include traveling and interacting shocks, acoustic waves,

vortex shedding, viscous flows, detonation waves, cavitation, flows in fluid film bearings, heat conduction

with melting and/or freezing, MHD vortex, hydraulic jump, crystal growth, chromatographic problems and

solar wind [3–63]. In particular, its unexpected simple non-reflecting boundary conditions [12,50] and rather

unique capability to resolve both strong shocks and small disturbances (e.g., acoustic waves) simultaneously

[19,21,22] make the CESE method an effective tool for attacking computational aeroacoustics (CAA) problems. Note that the fact that the second-order CESE method is capable of solving CAA problems accurately

is an exception to the commonly-held belief that a second-order scheme is not adequate for modeling CAA

1
problems. Also note that, while numerical dissipation is needed for shock capturing, it may also result in

annihilation of small disturbances. Thus a solver that can handle both strong shocks and small disturbances

simultaneously must be able to overcome this difficulty.

In spite of its nontraditional features and robust capabilities, the core ideas of the CESE method are

simple. In fact, all of its key features are the results of a pursuit driven by these simple ideas. The first

and foremost is the belief that the method must be solid in physics. As such, in the CESE development,

conservation laws are enforced locally and globally in their natural space-time unity forms for 1D, 2D and 3D

cases. Moreover, because direct physical interaction generally occurs only among the immediate neighbors,

use of the simplest stencil also becomes a common CESE feature. Obviously, this requirement has the effect

of simplifying boundary-condition implementation.

The second idea emerges from the realization that stability and accuracy are two competing issues in

time-accurate computations, i.e., too much numerical dissipation will degrade accuracy while too little of

it will cause instability. In other words, to meet both accuracy and stability requirements, computation

must be performed away from the edge (“cliff”) of instability but not too far from it. This represents a

real dilemma in numerical method development. As an example, high order schemes generally have higher

accuracy and lower numerical dissipation. However, they are susceptible to computational instabilities. In

fact, in complicated real-world applications, not only they seldom live up to their nominal order of accuracy—

generally they possess only first-order accuracy when shocks are present, stability of these schemes often is

difficult to maintain without resorting to problem-specific treatments. To confront this issue head-on, in

the CESE development, generally it is required that a solver be built from a nondissipative (i.e., neutrally

stable) core scheme. By definition, computations involving a neutrally stable scheme are performed right on

the edge of instability and therefore the numerical results generated are nondissipative. As such numerical

dissipation can be controlled effectively if the deviation of a solver from its nondissipative core scheme can

be adjusted using some built-in parameters.

Moreover, because an accurate viscous flow simulation requires that the numerical dissipation be much

smaller than the physical dissipation which decreases as the Reynolds number increases, in principle, an

accurate and robust solver for high-Reynolds-number flows must be able to cut numerical dissipation as

the Reynolds number increases. Obviously this requirement can only be met by a solver built from a

nondissipative core solver.

Other CESE ideas are: (i) the flux at an interface be evaluated in a simple and consistent manner;

(ii) genuinely multidimensional schemes be built as simple, consistent and straightforward extensions of 1D

schemes so that multidimensional shocks can be captured without using any mesh alignment technique;

(iii) triangular and tetrahedral meshes be used in 2D and 3D cases, respectively, so that the method is

compatible with the simplest unstructured meshes and thus can be easily used to solve problems with

complex geometries; and (iv) logical structures and approximation techniques used be as simple as possible,

and special techniques that has only limited applicability and may cause undesirable side effects be avoided.

Fortunately for the CESE development, realization of the above lesser ideas (i)–(iv) follows naturally from

the first two core ideas.

Building on its initial successes, efforts to refine and improve the CESE method have continued in the

past few years. As an example, it was shown in [5] that the numerical dissipation of a dissipative extension

of a CESE core scheme may increase to an intolerable level as the value of the CF L number decreases from

its maximum stability bound. As such, in a case with a large CF L number disparity (e.g., a simulation

with a highly non-uniform spatial mesh and a spatially independent time step), the performance sensitivity

with respect to the CF L number can lead to a solution that is highly dissipative in a region where the

local CF L number ! 1. Even though a remedy was suggested in [5], a simple and robust solution to this

problem had not arrived until a family of new Courant number insensitive schemes [25,31,46–49,51,54,55,58]

was developed with fresh insights. Note that these new schemes have one important advantage, i.e., all

variables at each mesh point can be evaluated explicitly without resorting to solving a system

of linear/nonlinear algebraic equations involving these local mesh variables even in applications

where systems of multidimensional nonlinear PDEs are solved.

As another example, even though they are accurate enough to solve CAA problems, second-order CESE

2
solvers are not capable of resolving fine flow structures within a boundary layer without using relatively fine

meshes and/or meshes with large aspect ratios. To overcome this limitation, two neutrally stable 4th-order

schemes, referred to as the a(3) and a(4) schemes [59–61], were developed in 2006 and intended to serve as

the core schemes of other high order CESE schemes. Unfortunately, the dissipative extensions of these new

schemes turned out to have the same intractable stability problem that afflicts traditional high order schemes.

As such it became obvious to the author that overcoming this problem requires a fresh new prospective.

It will be shown in this paper that, with the aid of a conceptual leap, one can build from a given 2ndorder CESE scheme its 4th-, 6th-, 8th-,... order versions which have the same stability conditions of

the 2nd-order scheme and also retain all other advantages of the latter scheme. Thus the high

order CESE schemes so constructed would avoid the common shortcomings of traditional high order schemes

including: (a) susceptibility to computational instabilities; (b) computational inefficiency due to their local

implicit nature (i.e., at each mesh point, need to solve a system of linear/nonlinear equations involving all

the mesh variables associated with this mesh point); (c) use of large and elaborate stencils which complicates

boundary treatments and also makes efficient parallel computing much harder; (d) difficulties in applications

involving complex geometries; and (e) use of problem-specific techniques which are needed to overcome

stability problems but often cause undesirable side effects.

Because the new high order CESE schemes are built from 2nd-order CESE schemes, the remainder of

the paper includes an extensive review of 2nd-order CESE schemes (Sec. 2–4). It is followed by a thorough

discussion of several linear and nonlinear high order CESE schemes (Sec. 5–7). A sketch of multidimensional

extensions is provided in Sec. 8. The numerical results are presented in Sec. 9. Conclusions and discussions

are given in Sec. 10. As a final note, in Sec. 10, we also present a hybrid 3rd-order scheme built from the

Lax-Wendroff scheme and a 2nd-order CESE scheme.

2. Review of second-order 1D CESE schemes

For simplicity, we review existing 2nd-order CESE solvers for the simple partial differential equation

(PDE)

∂u

∂t + a

∂u

∂x = 0 (2.1)

where the advection speed a "= 0 is a constant. Let x1 = x, and x2 = t be the coordinates of a twodimensional Euclidean space E2. Then, because Eq. (2.1) can be expressed as ∇ · "

h = 0 with "h def = (au, u),

Gauss’ divergence theorem in the space-time E2 implies that Eq. (2.1) is the differential form of the integral

conservation law

!

S(V )

"

h · d"s = 0 (2.2)

As depicted in Fig. 1, here (i) S(V ) is the boundary of an arbitrary space-time region V in E2, and (ii)

d"s = dσ "n with dσ and "n, respectively, being the area and the unit outward normal vector of a surface element

on S(V ). Note that: (i) because "h · d"s is the space-time flux of "h leaving the region V through the surface

element d"s, Eq. (2.2) simply states that the total space-time flux of "h leaving V through S(V ) vanishes;

(ii) in E2, dσ is the length of a line segment on the simple closed curve S(V ); and (iii) all mathematical

operations can be carried out as though E2 were an ordinary two-dimensional Euclidean space.

To proceed, let Ω denote the set of all space-time staggered mesh points in E2 (dots in Fig. 2(a)),

where n = 0, ±1/2, ±1, ±3/2, ±2,..., and, for each n, j = n ± 1/2, n ± 3/2, n ± 5/2,.... Each (j, n) ∈ Ω

is associated with a solution element, i.e., SE(j, n). By definition, SE(j, n) is the interior of the space-time

region bounded by a dashed curve depicted in Fig. 2(b). It includes a horizontal line segment, a vertical line

segment, and their immediate neighborhood. By definition, the end points of the line segments referred to

above are excluded from SE(j, n) so that two SEs will not overlap.

Eq. (2.2) will be simulated numerically assuming that, for any (x, t) ∈ SE(j, n), u(x, t) and "h(x, t),

respectively, are approximated by

u∗(x, t; j, n) def

= un

j + (ux)

n

j (x − xj ) + (ut)

n

j (t − t

n) (2.3)

3
and "

h∗(x, t; j, n) def= "au∗(x, t; j, n), u∗(x, t; j, n)

# (x, t) ∈ SE(j, n) (2.4)

Note that: (i) un

j , (ux)nj , and (ut)nj are constants in SE(j, n), and (in a rough sense) they can be considered

as the numerical analogues of the values of u, ∂u/∂x, and ∂u/∂t at the mesh point (j, n), respectively, (ii)

(xj , tn) are the coordinates of the mesh point (j, n) with xj = j∆x and t

n = n∆t, and (iii) Eq. (2.4) is the

numerical analogue of the definition "h = (au, u).

Let u = u∗(x, t; j, n) satisfy Eq. (2.1) within SE(j, n). Then one has

(ut)

n

j = −a (ux)

n

j (2.5)

As a result, Eq. (2.3) reduces to

u∗(x, t; j, n) = un

j + (ux)

n

j

$

(x − xj ) − a (t − t

n)

% (2.6)

i.e., un

j and (ux)nj are the only independent mesh variables associated with (j, n). Note that, for the purpose

of uniquely defining the flux over a boundary of a CE, u∗(x, t; j, n) is defined here only for (x, t) ∈ SE(j, n)

for each (j, n) ∈ Ω. In some later developments, this function will also be used in applications other than

flux evaluation such as construction of finite-difference approximations. Hereafter, for such extended

applications, it is to be understood that u∗(x, t; j, n) is defined even for a point (x, t) ∈/ SE(j, n).

Let E2 be divided into non-overlapping rectangular regions (see Fig. 2(a)) referred to as basic conservation elements (BCEs). As depicted in Figs. 2(c)–2(e), (i) each (j, n) ∈ Ω is assigned with two BCEs, i.e.,

CE−(j, n) and CE+(j, n); (ii) each BCE has one and only one pair of diagonally opposite vertices which

belong to Ω; (iii) the space-time E2 can be filled by CE−(j, n) and CE+(j, n), (j, n) ∈ Ω; and (iv) CE(j, n),

which is the union of CE−(j, n) and CE+(j, n), is referred to as a compounded conservation element (CCE).

Given the above preliminaries, we are ready to describe the existing CESE solvers for Eq. (2.1).

2.1. The a scheme

Note that, among the line segments forming the boundary of CE−(j, n), AB and AD belong to SE(j, n),

while CB and CD belong to SE(j − 1/2, n − 1/2). Similarly, the boundary of CE+(j, n) belongs to either

SE(j, n) or SE(j + 1/2, n − 1/2). As a result, by imposing the two local flux conservation conditions at each

(j, n) ∈ Ω, i.e.,

!

S(CE+(j,n))

"

h∗ · d"s = 0 (j, n) ∈ Ω (2.7)

and

!

S(CE−(j,n))

"

h∗ · d"s = 0 (j, n) ∈ Ω (2.8)

and using Eqs. (2.4) and (2.6), one can obtain two equations for two unknowns un

j and (ux)nj . In fact, with

the aid of

ν def

= a

∆t

∆x and (ux ̄)

n

j

def

= ∆x

4 (ux)

n

j (j, n) ∈ Ω (2.9)

it can be shown that (see Comment (a) given at the end of this subsection) Eqs. (2.7) and (2.8) are equivalent

to

(1 − ν) [u + (1 + ν)ux ̄]

n

j = (1 − ν) [u − (1 + ν)ux ̄]

n−1/2

j+1/2 (j, n) ∈ Ω (2.7a)

and

(1 + ν) [u − (1 − ν)ux ̄]

n

j = (1 + ν) [u + (1 − ν)ux ̄]

n−1/2

j−1/2 (j, n) ∈ Ω (2.8a)

respectively. To simplify notation, in the above and hereafter we adopt a convention that can be explained

using an expression on the left side of Eq. (2.7a) as an example, i.e.,

[u + (1 + ν)ux ̄]

n

j = un

j + (1 + ν)(ux ̄)

n

j

4
By adding Eqs. (2.7a) and (2.8a) together, one has

un

j = 1

2

&

(1 − ν) [u − (1 + ν)ux ̄]

n−1/2

j+1/2 + (1 + ν) [u + (1 − ν)ux ̄]

n−1/2

j−1/2

'

(j, n) ∈ Ω (2.10)

Let 1 − ν2 "= 0, i.e., 1 − ν "= 0 and 1 + ν "= 0. Then Eqs. (2.7a) and (2.7b) can be divided by (1 − ν) and

(1 + ν), respectively. By subtracting the resulting equations from each other, one has

(ux ̄)

n

j = (uax ̄)

n

j (j, n) ∈ Ω (2.11)

with

(ua

x ̄)

n

j

def

= 1

2

&

[u − (1 + ν)ux ̄]

n−1/2

j+1/2 − [u + (1 − ν)ux ̄]

n−1/2

j−1/2

'

(j, n) ∈ Ω (2.12)

It has been shown that Eqs. (2.7) and (2.8) imply (i) Eq. (2.10) for all ν; and (ii) Eq. (2.11) if |ν| "= 1.

On the other hand, it can be shown that Eqs. (2.10) and (2.11) are equivalent to

[u + (1 + ν)ux ̄]

n

j = [u − (1 + ν)ux ̄]

n−1/2

j+1/2 (j, n) ∈ Ω (2.7b)

and

[u − (1 − ν)ux ̄]

n

j = [u + (1 − ν)ux ̄]

n−1/2

j−1/2 (j, n) ∈ Ω (2.8b)

for all ν. Because (i) Eqs. (2.7b) and (2.8b) imply Eqs. (2.7a) and (2.8a), respectively; and (ii) Eqs. (2.7a)

and (2.8a) are equivalent to Eqs. (2.7) and (2.8), respectively, one concludes that, for all ν, the last two

conservation conditions are enforced by the a scheme which, by definition, is formed by Eqs. (2.10) and

(2.11) for all ν.

Note that, because

∂u

∂x ̄ = ∆x4

∂u

∂x

if  ̄x

def

= x/(∆x/4), the normalized parameter (ux ̄)n

j can be interpreted as the value at (j, n) of the derivative

of u with respective to the normalized coordinate  ̄x. Also note that: (i) the superscript symbol “a” in the

parameter (ua

x ̄)n

j is introduced to remind the reader that Eq. (2.11) is valid for the a scheme, and (ii) the

a scheme is second-order accurate in space and time for both un

j and the non-normalized parameter (ux)nj

[1,5].

The review of the a scheme is concluded with the following remarks:

(a) Evaluation of Eqs. (2.7) and (2.8) can be facilitated by the following observations: because u∗(x, t; j, n) is

linear in x and t, the total flux of "h∗ leaving CE−(j, n) or CE+(j, n) through any of the four line segments

that form its boundary is equal to the scalar product of the vector "h∗ evaluated at the midpoint of the

line segment and the “surface” vector (i.e., the unit outward normal vector multiplied by the length) of

the line segment.

Note that, by definition, points B and D depicted in Fig. 2(c) do not belong to either SE(j, n) or

SE(j −1/2, n−1/2). This fact, however, does not pose a problem for flux evaluation over S (CE−(j, n))

because the values of "h∗ at isolated points do not contribute to the flux of "

h∗ over a finite line segment.

A similar comment applies to points D and F depicted in Fig. 2(d).

(b) Even though they have completely different origins and structures (one is a two-level scheme with two

independent variables per mesh point while another is a three-level scheme with a single independent

variable per mesh point), the a scheme and the classical leapfrog scheme are deeply related. In fact, it

can be shown that each of the two variables un

j and (ux)nj of the a scheme by itself satisfies the same

leapfrog scheme automatically [56]. Moreover, after sorting through the differences in the notations

used in typical CESE papers and traditional literatures, and also taking into account the decoupled

nature of a leapfrog solution, it is shown in [5, pp.301 and 313–315] that the two amplification factors

5
of the a scheme are identical to those of the neutrally stable leapfrog scheme. Thus, the a scheme is

also neutrally stable (i.e., nondissipative) if |ν| < 1 (see the additional discussions given in Sec. 2.2).

(c) Even though it is introduced to model a single PDE (i.e., Eq. (2.1)) with a single dependent variable

u, the a scheme is formed by two coupled discrete equations (i.e., Eqs. (2.10) and (2.11)) involving

two independent numerical variables un

j and (ux)nj . As such one would expect that the a scheme be

consistent with a pair of PDEs with one of them being Eq. (2.1). Indeed, it is shown in [1] that, regardless

how ∆t > 0 and ∆x > 0 are refined, Eqs. (2.10) and (2.11) are always consistent with Eq. (2.1) and

v − ∂u/∂x = 0 with the understanding that v being the analytical version of (ux)n

j .

(d) Because (i) the vector "h∗ at any surface element lying on any interface separating two neighboring BCE

is evaluated using the information from a single SE, and (ii) the unit outward normal vector on the

surface element pointing outward from one of these two neighboring BCEs is exactly the negative of that

pointing outward from another BCE, one concludes that the flux leaving one of these BCEs through the

interface is the negative of that leaving another BCE through the same interface. As a result of this flux

cancelation, the local flux conservation relations Eqs. (2.7) and (2.8) lead to a global flux conservation

relation, i.e., the total flux of "h∗ leaving the boundary of any space-time region that is the union of any

combination of BCEs will also vanish. In particular, because CE(j, n) is the union of CE−(j, n) and

CE+(j, n),

!

S(CE(j,n))

"

h∗ · d"s = 0 (j, n) ∈ Ω (2.13)

must follow from Eqs. (2.7) and(2.8). In fact, it can be shown that Eq. (2.13) is equivalent to Eq. (2.10).

(e) In addition to the nondissipative a scheme, as will be shown, there is a family of its dissipative extensions

in which only the less stringent conservation condition Eq. (2.13) is assumed [5]. Because Eq. (2.13) is

equivalent to Eq. (2.10), for each of these extensions, un

j is still evaluated using Eq. (2.10) while (ux ̄)nj

is evaluated using an equation different from Eq. (2.11).

2.2. The a-% scheme and the c scheme

Consider any (j, n) ∈ Ω. Then (j ± 1/2, n − 1/2) ∈ Ω. Let

u# n

j±1/2

def

= u∗ "

xj±1/2, tn; j ± 1/2, n − 1/2

#

= [u + (∆t/2)ut]

n−1/2

j±1/2 (2.14)

The validity of the last equality sign in Eq. (2.14) follows from Eqs. (2.5) and (2.6). Thus u# n

j±1/2 is the

first-order Taylor approximation at the point (xj±1/2, tn) (which ∈/ SE(j ±1/2, n−1/2)) evaluated using the

mesh variables at (j ±1/2, n−1/2). Moreover, with the aid of Eqs. (2.6) and (2.9), one can recast Eq. (2.14)

in the form

u# n

j±1/2 = (u − 2ν ux ̄)

n−1/2

j±1/2 (2.15)

Note that, by definition, (j ± 1/2, n) ∈/ Ω if (j, n) ∈ Ω. Thus u# n

j±1/2 is associated with a mesh point ∈/ Ω.

The reader is warned that similar situations may occur in the rest of this paper.

Because u# n

j±1/2 is a first-order Taylor’s approximation of u at (j ± 1/2, n),

(uc

x ̄)

n

j

def

= ∆x

4

(u# n

j+1/2 − u# nj−1/2

∆x

)

≡ u# n

j+1/2 − u# nj−1/2

4 (2.16)

is a central-difference approximation of ∂u/∂x at (j, n), normalized by the same factor ∆x/4 that appears

in Eq. (2.9). Note that: (i) the superscript “c” is used to remind the reader of the central-difference origin

of the term (uc

x ̄)n

j ; and (ii) by using Eqs. (2.15), (2.16) and (2.12), one has

(uc

x ̄)

n

j = 1

4

*

(u − 2νux ̄)

n−1/2

j+1/2 − (u − 2νux ̄)

n−1/2

j−1/2

+

(2.17)

6
and

(uc

x ̄ − uax ̄)

n

j = 1

4

*

(u + 2ux ̄)

n−1/2

j−1/2 − (u − 2ux ̄)

n−1/2

j+1/2

+

(2.18)

The a-% scheme is formed by Eq. (2.10) and

(ux ̄)

n

j = (uax ̄)

n

j + 2%(ucx ̄ − uax ̄)

n

j (2.19)

where % is an adjustable real number. Note that: (i) with the aid of Eqs. (2.12) and (2.18), Eq. (2.19) can

be rewritten as

(ux ̄)

n

j = 1

2

&

[(1 − %)u + (2% − 1 − ν)ux ̄]

n−1/2

j+1/2 − [(1 − %)u + (1 − ν − 2%)ux ̄]

n−1/2

j−1/2

'

(2.19a)

(ii) the a-% scheme reduces to the a scheme when % = 0; (iii) for the case % = 1/2, Eq. (2.19) reduces to

(ux ̄)

n

j = (ucx ̄)

n

j (2.20)

and (iv) because (uc

x ̄)n

j represents a central-difference approximation, hereafter, to simplify its frequent

references, the special a-% scheme with % = 1/2 will be referred to as the c scheme.

As a preliminary to later developments, we offer the following remarks on the a-% scheme:

(a) At each (j, n) ∈ Ω, Eqs. (2.10) and (2.11) imply Eqs. (2.7) and (2.8). Because Eq. (2.19) does not reduce

to Eq. (2.11) except in the special case % = 0, at each (j, n) ∈ Ω, generally the a-% scheme satisfies only

the single conservation condition Eq. (2.13) (which is equivalent to Eq. (2.10)) rather than the two

conservation conditions Eqs. (2.7) and (2.8). However, because (ua

x ̄)n

j generally is present on the right

side of Eq. (2.19), implementation of the a-% scheme generally will still be burdened with the cost of

solving two conservation conditions at each mesh point. The exception occurs only for the special case

% = 1/2 (i.e., the c scheme) in which Eq. (2.19) reduces to Eq. (2.20). As it turns out, implementation of

an Euler version of the c scheme does not require inverting any system of equations while that involving

an Euler version of any other a-% scheme (% "= 1/2) generally requires inverting, per mesh point and

per time step, a system of several linear equations. As an example, because the 3D Euler equations is

formed by a system of five PDEs with each of them being associated with three independent spatial

derivatives at each mesh point, a 3D Euler a-% (% "= 1/2) scheme requires inverting a system of 15

(= 5 × 3) equations. Partly because of this cost consideration, extensions of the c scheme have been

used extensively.

(b) For the a-% scheme, it is shown in [5] that the principal and spurious amplification factors per marching

step (i.e., per ∆t/2—recall that the solution is advanced by a time period ∆t/2 per marching step) are

λ+ and λ−, respectively, with

λ±(%, ν, θ) def

= % cos(θ/2) − iν sin(θ/2) ±

,

(1 − %)

$

(1 − %)cos2(θ/2) + (1 − ν2)sin2(θ/2)% (2.21)

Here (i) i def

= √−1, and (ii) θ, −∞ < θ < +∞, is the phase angle variation of a Fourier component over

the length ∆x (i.e., θ = k∆x where k denotes the phase angle variation per unit length). In addition, it

can be shown that the a-% scheme is stable ⇔ (i.e., “if and only if”) either

0 ≤ % ≤ 1 and |ν| < 1 (2.22a)

or

% = 1 and |ν| = 1 (2.22b)

(Note: Because Eq. (2.11) is a result of Eqs. (2.7) and (2.8) only if |ν| "= 1, the a and a-% schemes were

not defined for the case |ν| = 1 in [5]. As such the case Eq. (2.22a) was not presented in [5]). Also

one can show that the a-% scheme becomes more dissipative as the value of % increases from 0 to 1 [5].

7
Note that, unless specified otherwise, in the remainder of the paper the ranges of % and ν respectively,

are defined by Eq. (2.22a,b). Moreover, note that: (i) the mesh points used in the a-% scheme are

so arranged that those at any two time levels are staggered in space-time unless these two levels are

separated by a whole integer multiple of ∆t, and (ii) by applying Eqs. (2.10) and (2.19) successively

over two consecutive marching steps, one can derive a compounded form (Eq. (3.7) in [5]) of the a-%

scheme in which the mesh variables at any (j, n + 1) ∈ Ω are determined in terms of those at (j − 1, n),

(j, n) and (j + 1, n)—i.e., the space-time staggering nature of the original scheme is eliminated from

the compounded form. As such, (i) Hereafter two consecutive marching steps will be referred to as a

marching cycle, and (ii) in the following accuracy evaluation of the numerical amplification factors, we

choose to consider the amplification factors per marching cycle (i.e., per ∆t).

Note that, for several reasons which are related to the space-time staggered nature of the meshes used in

the CESE method (see pp. 313-315 in [5]), the same spatial-interval length and time-step size denoted

by ∆x and ∆t respectively in traditional literatures, are denoted by ∆x/2 and ∆t/2 respectively here.

Because the same definition θ def

= k∆x is used here and in traditional literatures, for the same Fourier

component (i.e., the same value of k), the corresponding value of θ in traditional literatures is only one

half of that in the current paper, i.e., the symbol θ/2 that appears on the right side of Eq. (2.21) would

be replaced by θ in traditional literatures.

(c) Let k be a constant. Then u = eik(x−at) represents a plane wave solution to Eq. (2.1). For this solution,

k is the phase angle variation per unit length. Thus, by using the relations θ = k∆x and ν = a∆t/∆x,

one concludes that

the analytical amplification factor per ∆t def

= eik[x−a(t+∆t)]

eik(x−at) = e−ika∆t = e−iνθ (2.23)

(d) The amplification factors of the a-% scheme per ∆t (i.e., per two consecutive marching steps) are given

by (λ±)2. According to Eq. (2.21), [λ±(0, ν, θ)]2, the amplification factors of the a scheme (which

corresponds to the case % = 0) per ∆t, have the following properties:

|[λ±(0, ν, θ)]2| = 1 if |ν| ≤ 1 (2.24)

lim ν→±1

[λ+(0, ν, θ)]2 = e∓iθ (2.25)

lim ν→±1

[λ−(0, ν, θ)]2 = e±iθ (2.26)

and

[λ±(0, 0, θ)]2 = 1 (2.27)

On the other hand, e−iνθ, the analytical amplification factor per ∆t, has the following properties:

|e−iνθ| = 1 (2.28)

lim ν→±1

e−iνθ = e∓iθ (2.29)

and

e−iνθ = 1 if ν = 0 (2.30)

For the a scheme, Eqs. (2.24)–(2.30) imply that: (i) the two amplification factor of the scheme, and

the analytical amplification factor all have the same constant absolute value (= 1) if |ν| ≤ 1, i.e.,

the scheme is nondissipative if |ν| ≤ 1; (ii) in the limit of |ν| → 1 (i.e., ν → 1 or ν → −1), the

principal amplification factor is identical to the analytical amplification factor and, thus, the former

has no dissipative or dispersive error in this limit; (iii) also in the limit of |ν| → 1, the phase angle

associated with the spurious amplification factor is exactly the negative of that associated with the

analytical amplification factor—thus the spurious amplification factor has a large dispersive error in

8
this limit except when |θ| ! 1 (i.e., when the wavelengths of the errors - 1); and (iv) when ν = 0, the

two amplification factors of the scheme, and the exact amplification factor are all equal to 1 and, thus,

the two amplification factors of the scheme have no dissipative or dispersive error if ν = 0. Because

the accuracy of a scheme is essentially hinged on the behaviors of the principal amplification factor [1],

according to the facts stated above, the a scheme tends to become very accurate when |ν| approaches

1 or 0. However, the dispersive errors associated with the nondissipative spurious amplification factor

(which could be introduced at t = 0 as a result of an inaccurate initial-value specification [1]) may

appear in a solution as persistent numerical wiggles when |ν| approaches 1.

(e) According to Eq. (2.21), [λ±(1/2, ν, θ)]2, the amplification factors of the c scheme (which corresponds

to the case % = 1/2) per ∆t, have the following properties:

lim ν→±1

[λ+(1/2, ν, θ)]2 = e∓iθ (2.31)

lim ν→±1

[λ−(1/2, ν, θ)]2 = − sin2(θ/2) (2.32)

and

[λ±(1/2, 0, θ)]2 = 1

2

*

1 ± cos(θ/2)-2 − cos2(θ/2)+

(2.33)

Because (i) eiπ = −1 and (ii)

1 −

*

cos(θ/2)-2 − cos2(θ/2)+2

≡ $1 − cos2(θ/2)%2 ≥ 0

i.e.,

1 ≥

.

.

.cos(θ/2)-2 − cos2(θ/2)

.

.

.

Eqs. (2.32) and (2.33) can be recast as

lim ν→±1

[λ−(1/2, ν, θ)]2 = sin2(θ/2)eiπ (2.32a)

and

[λ±(1/2, 0, θ)]2 = 1

2

.

.

.

1 ± cos(θ/2)-2 − cos2(θ/2)

.

.

. (2.33a)

respectively. By comparing Eqs. (2.29) and (2.30) with Eqs. (2.31), (2.32a) and (2.33a), one arrives at

the following conclusions for the c scheme: (i) in the limit of |ν| → 1, the principal amplification factor

is identical to the analytical amplification factor and, thus, the former has no dissipative or dispersive

errors in this limit; (ii) also in the limit of |ν| → 1, the spurious amplification factor generally has large

dissipative and dispersive errors; and (iii) when ν = 0, both the principal and spurious amplification

factors generally have large dissipative errors but no dispersive errors. According to the facts stated

above, like the a scheme, the c scheme also tends to become very accurate when |ν| approaches 1.

However, unlike the a scheme, the errors associated with the spurious amplification factor of the c

scheme do die out rapidly when |ν| approaches 1 if the value of |θ| is not too close to π (see Eq. (2.32a)),

i.e., if the wavelength of the discrete Fourier component is not close to 2∆x—the shortest possible

wavelength. Also, in sharp contrast to the a scheme, the c scheme becomes highly dissipative when ν

approaches 0.

From the above discussions, one concludes that:

(a) The advantages of the a scheme include: (i) it is nondissipative for all θ if |ν| ≤ 1; and (ii) when the value

of |ν| is close to 0 or 1, the scheme is very accurate. On the other hand, its disadvantages include: (i)

because it is nondissipative, its extensions for nonlinear equations generally are unstable; (ii) when the

value of |ν| is close to 1, the dispersive errors associated with the nondissipative spurious amplification

9
factors will appear in a solution as persistent numerical wiggles; and (iii) comparing with the c scheme,

it costs more to implement.

(b) The advantages of the c scheme include: (i) when the value of |ν| is close to 1, it is very accurate and

the errors associated with the spurious amplification factor generally die out rapidly; (ii) because of its

dissipative property, its extensions for nonlinear equations can be stable; and (iii) in terms of ease of

implementation and computer cost, it is much more superior than any other a-% scheme. On the other

hand, the c scheme has a serious disadvantage, i.e., it is very dissipative when ν approaches 0.

Based on the remarks given above, it becomes clear that a new solver of Eq. (2.1) possesses all the

advantages of the a and c schemes but none of their disadvantages if it is formed by Eq. (2.10) and a new

equation in which (i) (ux ̄)n

j is evaluated using a simple central-differencing procedure similar to that used

to obtain (uc

x ̄)n

j ; and (ii) the (ux ̄)nj so obtained is identical to (uc

x ̄)n

j in the limit of |ν| → 1 and identical to

(ua

x ̄)n

j in the limit of ν → 0. In other words, such an ideal solver (i) is comparable to the c scheme in ease

of implementation; (ii) becomes the c scheme in the limit of |ν| → 1; and (iii) becomes the a scheme in the

limit of ν → 0. In Sec. 3, it will be shown how such a solver can be constructed.

2.3. The w-α scheme—a special wiggle-suppressing scheme

If discontinuities are present in a numerical solution, any a-% scheme such as the c scheme is not equipped

to suppress numerical wiggles that generally appear near these discontinuities. To serve as a preliminary

for future development, here we shall briefly review an extension of the c scheme which was introduced as a

remedy for this deficiency [5].

To proceed, let

(ux ̄−)

n

j

def

= ∆x

4

(un

j − u# n

j−1/2

∆x/2

)

≡ un

j − u# n

j−1/2

2 (2.34)

and

(ux ̄+)

n

j

def

= ∆x

4

(u# n

j+1/2 − un

j

∆x/2

)

≡ u# n

j+1/2 − un

j

2 (2.35)

i.e., (ux ̄−)n

j and (ux ̄+)nj are normalized numerical analogues of ∂u/∂x at (j, n) evaluated from the left and

the right, respectively. By using Eqs. (2.16), (2.34), and (2.35), it can be shown that

(uc

x ̄)

n

j = 1

2 (ux ̄− + ux ̄+)

n

j (2.36)

i.e., (uc

x ̄)n

j is the simple average of (ux ̄−)nj and (ux ̄+)nj . As such, the c scheme can be extended by replacing

(uc

x ̄)n

j in Eq. (2.20) with an weighted average of (ux ̄−)nj and (ux ̄+)nj . In other words, the resulting extension

is formed by Eq. (2.10) and

(ux ̄)

n

j = (w−)

n

j (ux ̄−)

n

j + (w+)

n

j (ux ̄+)

n

j (2.37)

where (w−)n

j and (w+)nj , the weight factors associated with (ux ̄−)nj and (ux ̄+)nj respectively, must satisfy the

condition

(w−)

n

j + (w+)

n

j = 1 (2.38)

at all (j, n) ∈ Ω. In addition, it is assumed that the expression on the right side of Eq. (2.37) represents an

interpolation (rather than an extrapolation) of (ux ̄−)n

j and (ux ̄+)nj . This assumption ⇔

(w−)

n

j ≥ 0 and (w+)

n

j ≥ 0 (2.39)

For real variables x−, x+, and α ≥ 0, let W− and W+ be the functions defined by: (i) W−(x−, x+; α) =

W+(x−, x+; α) = 1/2 if x− = x+ = 0; and (ii)

W−(x−, x+; α) = |x+|

α

|x−|

α + |x+|α (|x−| + |x+| > 0) (2.40)

10
and

W+(x−, x+; α) = |x−|

α

|x−|

α + |x+|α (|x−| + |x+| > 0) (2.41)

if either x− "= 0 or x+ "= 0. Furthermore, let

(w±)

n

j = W±((ux ̄−)

n

j ,(ux ̄+)

n

j , α) (2.42)

Then (w−)n

j and (w+)nj so defined satisfy Eqs. (2.38) and (2.39) and have the property that

(w−)

n

j = (w+)

n

j = 1/2 if α = 0 or |(ux ̄−)

n

j | = |(ux ̄+)

n

j | (2.43)

Note that: (i) to avoid dividing by zero, in practice either a small positive number such as 10−20 is added to

each of the denominators in Eqs. (2.40) and (2.41), or one must make sure that the definition W−(x−, x+; α) =

W+(x−, x+; α) = 1/2 instead of Eqs. (2.40) and (2.41) is applied when x− = x+ = 0; and (ii) the special

cases of Eqs. (2.40) and (2.41) with α = 1 and α = 2 are used in the slope-limiter proposed by van Leer [64]

and van Albada et al. [65].

An extension of the c scheme is formed by Eqs. (2.10) and (2.37) with (w−)n

j and (w+)nj being defined

by Eq. (2.42). Because it involves an weighted average which is dependent on a parameter α, hereafter the

scheme is referred to as the w-α scheme. In a straightforward manner, the w-α scheme has been extended

to become the Euler w-α scheme. The details are given in [5].

Assuming α > 0 and |(ux ̄−)n

j | "= |(ux ̄+)nj |, then Eqs. (2.37), (2.40), and (2.41) imply that, of (ux ̄−)nj and

(ux ̄+)n

j , the one with smaller absolute value is assigned an weight factor > 1/2. This observation coupled

with Eqs. (2.36)–(2.39) leads to the conclusion that, of (ux ̄−)n

j and (ux ̄+)nj , (ux ̄)nj will have an algebraic

value closer to the one with smaller absolute value if (ux ̄)n

j is evaluated as an weighted average of (ux ̄−)nj

and (ux ̄+)n

j according to Eq. (2.37). As a result, (ux ̄)nj so evaluated has a smaller absolute value than that

evaluated using Eq. (2.20). In turn, numerical wiggles or overshoots can be annihilated by the additional

numerical dissipation introduced as a result of this local “flattening” of (ux ̄)n

j . It has been shown numerically

that the extension is stable if |ν| < 1 and α ≥ 0. Moreover, as a result of Eqs. (2.20), (2.36), (2.37) and

(2.43), (i) the extension reduces to the c scheme when α = 0; and (ii) even if α > 0, the extension behaves

very much like the c scheme in any smooth solution region (where the condition (ux ̄−)n

j = (ux ̄+)nj more or

less prevails) or at a solution extremum (where the condition (ux ̄−)n

j = −(ux ̄+)nj more or less prevails). As

such, the wiggle-suppressing power of the extension takes effect only if α > 0 and only in a solution region

where |(ux ̄−)n

j | and |(ux ̄+)nj | differ substantially.

This section is concluded with the following discussion of CESE non-reflecting boundary conditions.

2.4. Non-reflecting boundary conditions

An unique feature of the CESE method is its simple and robust non-reflecting boundary conditions

[12,50]. Let (j, n) be a mesh point on the right spatial boundary of the computational domain. Then (i)

(j − 1/2, n − 1/2) represents an interior mesh point at the (n − 1/2)th time level; and (ii)

un

j = un−1/2

j−1/2 and (ux ̄)

n

j = (ux ̄)

n−1/2

j−1/2 (2.44)

form a set of non-reflecting boundary conditions. On the other hand, for a mesh point (j, n) on the left

boundary, (i) (j + 1/2, n − 1/2) represents an interior mesh point at the (n − 1/2)th time level; and (ii)

un

j = un−1/2

j+1/2 and (ux ̄)

n

j = (ux ̄)

n−1/2

j+1/2 (2.45)

form a set of non-reflecting boundary conditions. The above non-reflecting conditions are applicable to all

the second-order schemes described in Sec. 2-4. A rigorous justification of these and other simple CESE

non-reflecting boundary conditions is given in [50].

3. The c-τ and c-τ ∗ schemes

11
The ideal solvers of Eq. (2.1) defined at the end of Sec. 2.2 will be developed here. As a preliminary, we

shall show that (ua

x ̄)n

j can be cast into a central-difference form when ν = 0.

According to Figs. 2(d), |EF| = |AD| = 0 when ∆t = 0. Moreover, because of the assumption a "= 0,

ν = 0 ⇔ ∆t = 0. Thus, in the case ν = 0, Eq. (2.7) reduce to the condition that the sum of the fluxes leaving

CE+(j, n) through its top face AF and bottom face ED is zero. Because (i) "n = (0, 1) at AF and "n = (0, −1)

at ED ("n denotes the outward unit normal vector at any boundary of CE+(j, n)), (ii) |AF| = |ED|, and

(iii) "h∗ = (au∗, u∗), by using the flux evaluation rule referred to in Remark (a) given at the end of Sec. 2.1,

one concludes that Eq. (2.7) reduces to

(u + ux ̄)

n

j = (u − ux ̄)

n−1/2

j+1/2 (ν = 0) (3.1)

if ν = 0. Similarly, Eq. (2.8) reduces to

(u − ux ̄)

n

j = (u + ux ̄)

n−1/2

j−1/2 (ν = 0) (3.2)

if ν = 0. Because Eqs. (2.10) and (2.11) (which form the a scheme) are equivalent to Eqs. (2.7) and (2.8) if

|ν| "= 1, they must be equivalent to Eqs. (3.1) and (3.2) when ν = 0. In fact, by subtracting Eq. (3.2) from

Eq. (3.1), one obtains Eq. (2.11) with (ua

x ̄)n

j being the reduced form of Eq. (2.12) when ν = 0, i.e.,

(ua

x ̄)

n

j = 1

2

*

(u − ux ̄)

n−1/2

j+1/2 − (u + ux ̄)

n−1/2

j−1/2

+

(ν = 0) (3.3)

On the other hand, by summing over Eqs. (3.1) and (3.2), one has the reduced form of Eq. (2.10) for the

case ν = 0.

With the aid of Eq. (3.3) and the facts that: (i) (u − ux ̄)

n−1/2

j+1/2 represents the approximation of u at

the midpoint of ED while (u + ux ̄)

n−1/2

j−1/2 represents that at the midpoint of CD (see Fig. 2(c,d)); and (ii)

the distance between the two midpoints is ∆x/2, one concludes that, for the special case ν = 0, (ua

x ̄)n

j is

indeed a central-difference approximation of ∂u/∂x at (j, n − 1/2) (which coincides with (j, n) when ν = 0),

normalized by the factor ∆x/4. QED.

According to the above discussions, construction of the ideal solvers is hinged on finding central-difference

approximations for (ux ̄)n

j such that each approximation (i) becomes (uc

x ̄)n

j in the limit of |ν| → 1, and (ii)

reduces to the expression on the right side of Eq. (3.3) when ν = 0. As a result, these new solvers will be

constructed as the subschemes of the c-τ scheme, a new class of CESE solvers for Eq. (2.1) to be described

immediately.

3.1. The c-τ scheme

To proceed, refer to Fig. 3. Here M+ and M− denote the midpoints of AF and AB, respectively.

Also P + and P − are two points on BF that satisfy the following conditions: (i) P + is to the right (left)

of M+ ⇔ P − is to the left (right) of M−; and (ii) |M+P +| = |M−P −|. In addition, the parameter τ is

defined by: (i) τ = 0 if P + coincides with M+; (ii) τ ∆x/4 = |M+P +| > 0 if P + is to the right of M+; and

(iii) τ ∆x/4 = −|M+P +| < 0 if P + is to the left of M+. Obviously, it follows from the above definitions

that (i) τ = 0 if P − coincides with M−; (ii) τ ∆x/4 = |M−P −| > 0 if P − is to the left of M−; and (iii)

τ ∆x/4 = −|M−P −| < 0 if P − is to the right of M−.

Next let

u(P +) def

= u∗ "

x(P +), tn; j + 1/2, n − 1/2

# (3.4)

and

u(P −) def

= u∗ "

x(P −), tn; j − 1/2, n − 1/2

# (3.5)

where x(P +) and x(P −) denote the x-coordinates of P + and P −, respectively. Then, according to Eqs. (2.5)

and (2.6), u(P +) is the first-order Taylor’s approximation of u at P + evaluated using the mesh variables

12
at (j + 1/2, n − 1/2) (i.e., point E) while u(P −) is that at P − evaluated using the mesh variables at

(j − 1/2, n − 1/2) (i.e., point C). Moreover, with the aid of Eqs. (2.6) and (2.9) along with the relations

x(P +) = xj + (1 + τ)

∆x

4 = xj+1/2 − (1 − τ)

∆x

4 , −∞ < τ < +∞ (3.6)

and

x(P +) = xj − (1 + τ)

∆x

4 = xj−1/2 + (1 − τ)

∆x

4 , −∞ < τ < +∞ (3.7)

Eqs. (3.4) and (3.5), respectively, can be simplified as

u(P +) = [u − (1 + 2ν − τ)ux ̄]

n−1/2

j+1/2 (3.8)

and

u(P −) = [u + (1 − 2ν − τ)ux ̄]

n−1/2

j−1/2 (3.9)

To proceed, note that: (i) the mesh point (j, n) (i.e., point A depicted in Fig. 3) is the midpoint of

P −P +, and (ii) |P −P +| = (1 + τ)∆x/2. Thus

(ˆux ̄)

n

j

def

= ∆x

4

/

u(P +) − u(P −)

(1 + τ)∆x/2

0

≡ u(P +) − u(P −)

2(1 + τ) (τ "= −1) (3.10)

represents a central-difference approximation of ∂u/∂x at the mesh point (j, n), normalized by the factor

∆x/4. Thus the new scheme defined by Eq. (2.10) and (ux ̄)n

j = (ˆux ̄)nj represents a solver for Eq. (2.1).

Moreover, by using Eqs. (3.8)–(3.10), the second component of the new scheme can be explicitly expressed

as

(ux ̄)

n

j = (ˆux ̄)

n

j ≡ 1

2(1 + τ)

&

[u − (1 + 2ν − τ)ux ̄]

n−1/2

j+1/2 − [u + (1 − 2ν − τ)ux ̄]

n−1/2

j−1/2

'

(τ "= −1) (3.11)

In view of the facts that (i) (ˆux ̄)n

j represents a central-difference approximation of (ux ̄)nj , and (ii) the approximation is associated with the parameter τ, hereafter the new scheme will be referred to as the c-τ

scheme.

Moreover, by using Eqs. (2.12), (2.18) and (3.11), one has

(ˆux ̄)

n

j = (uax ̄)

n

j +

2τ

1 + τ

(uc

x ̄ − uax ̄)

n

j − ν(1 − τ)

2(1 + τ)

*

(ux ̄)

n−1/2

j+1/2 − (ux ̄)

n−1/2

j−1/2

+

(τ "= −1) (3.12)

By comparing Eq. (3.12) with (2.19), one concludes that the c-τ scheme generally is different from the a-%

scheme. In fact, the c-τ scheme becomes the a-% scheme ⇔ either (i) τ = 1 or (ii) ν = 0. For the case τ = 1,

Eq. (3.12) implies that (ˆux ̄)n

j = (uc

x ̄)n

j . In other words, the c scheme is the special case of the c-τ scheme

with τ = 1, a fact that can also be deduced from the observation that the points P + and P − depicted in

Fig. 3, respectively, coincide with points F and B (i.e., the mesh points (j + 1/2, n) and (j −1/2, n)) if τ = 1.

On the other hand, when ν = 0, the c-τ scheme become the a-% scheme with % = 2τ/(1 + τ). In fact one can

further deduce that c-τ scheme reduces to the a scheme ⇔ ν = τ = 0.

Because the c-τ scheme is formed by two rather complicated equations (i.e., Eqs. (2.10) and (3.11))

involving two parameters ν and τ, it were not expected that its von Neumann stability conditions could be

cast into an explicit analytical form. However, by using the Jordan canonical form theorem [66], it is shown

rigorously in [46] that the c-τ scheme is von Neumann stable ⇔

ν2 ≤ 1, τ ≥ τo(ν2), and (ν2, τ) "= (1, 1) (3.13)

13
where

τo(s)

def

=







0 if s = 0

4 − s − 2

-2(2 − s − s2)

s

if 0 < s ≤ 3

11

s − 1 + √1 − 2s + 5s2

2s

if 3

11 ≤ s ≤ 1

(3.14)

Note that:

(a) The current stability conditions agree completely with those generated through numerical experiments

[31] (Note: the reader is warned that the condition τ ≥ τo(ν2) that appears in Eq. (3.13) here is

expressed as τ ≥ τo(|ν|) in Eq. (3.12) of [31], i.e., the function τo defined here is different from the

function τo introduced in [31]).

(b) It can be shown that [46]: (i) τo(s) is continuous at s = 0, i.e.,

lim

s→0+ τo(s) = τo(0) = 0

(ii) τo(s) is consistently defined at s = 3/11, i.e., both the second and third expressions on the right

side of Eq. (3.14) implies that τo(3/11) = 1/3; (iii) τo(s) is continuously differentiable at s = 3/11, i.e.,

lim

s→ 3

11

−

τ#

o(s) = lim s→ 3

11

+

τ#

o(s) = 121/90

where τ#

o(s) def

= dτo(s)/ds; (iv) τo(s) is strictly monotonically increasing in the interval 0 <s< 1; (v)

τo(1) = 1; and (vi)

s < τo(s) < √s, 0 <s< 1 (3.15)

(c) For any given fixed value of |ν| < 1, the c-τ scheme tends to become more dissipative as the value of τ

increases from its minimum stability bound τo(ν2).

With the above preliminaries, the ideal solvers of Eq. (2.1) will be constructed in Sec. 3.2.

3.2. The c-τ ∗ schemes

The value of τ used in the c-τ scheme generally can be chosen independent of ν. Here we will introduce

a subset of the c-τ scheme in which τ is a function of ν2 for each member of this subset. As a preliminary,

note that, by using the properties of τo(s) presented earlier, it can be shown that there exist infinitely many

choices of a strictly monotonically increasing smooth function h(s), 0 ≤ s < 1, which possesses the following

properties: (i)

h(0) = 0; and lims→1−

h(s) = 1 (3.16)

and (ii)

h(s) ≥ τo(s) if 0 < s < 1 (3.17)

By definition, the subset contains every special c-τ scheme with

τ = h(ν2) (|ν| < 1) (3.18)

where h(s) is one of the special functions defined above. By using Eqs. (3.16)–(3.18) and the fact that h(s)

is a strictly monotonically increasing function in the domain 0 ≤ s < 0, one concludes that, for each member

of the subset, (i)

τ = 0 if ν = 0 (3.19)

(ii)

lim

ν2→1− τ = 1 (3.20)

14
(iii) τ increases with |ν| in the domain 0 ≤ |ν| < 1; and (iv)

τ ≥ τo(ν2) (|ν| < 1) (3.21)

By using the results presented in the above items (i)–(iii), one can easily infer from Fig. 3 a simple relation

between the value of |ν| and the locations of P + and P −, i.e., as the value of |ν| increases from 0 to 1, P +

will move away from M+ and edge toward the mesh point (j + 1/2, n) while P − will move away from point

M− and edge toward the mesh point (j − 1/2, n).

Recall that (i) (ˆux ̄)n

j = (ua

x ̄)n

j if τ = ν = 0; and (ii) (ˆux ̄)nj = (uc

x ̄)n

j if τ = 1. As such, Eqs. (3.19) and

(3.20) imply that, for each member in the subset, (i) (ˆux ̄)n

j = (ua

x ̄)n

j if ν = 0; and (ii) (ˆux ̄)nj = (uc

x ̄)n

j in the

limit of |ν| → 1−. In other words, in the domain |ν| < 1, all members in the subset are the ideal solvers

defined at the end of Sec. 2.2. Moreover, by using Eq. (3.13) and (3.21), one can also show that these ideal

solvers are also stable in the same domain. Hereafter, each of these ideal solvers will be referred to as a c-τ ∗

scheme.

Corresponding to infinitely many choices of h, there are infinitely many different c-τ ∗ schemes. In

particular, with the aid of Eq. (3.15), one can see easily that the special choice

h(s) = √s (0 ≤ s < 1) (3.22)

represents a strictly monotonically increasing smooth function which also satisfies Eqs. (3.16) and (3.17).

Thus Eq. (3.18) implies that a c-τ ∗ scheme can be defined using the relation

τ = |ν| (|ν| < 1) (3.23)

Any c-τ ∗ scheme described above meets all the requirements of an ideal solver. However, because of

stability problem, a c-τ ∗ scheme and its multidimensional extensions may not be robust enough for some

complicated real-world applications. To overcome this difficulty, we introduce a special c-τ scheme which is

defined by Eqs. (2.10) and (3.11) with

τ = β |ν| (β ≥ 1; |ν| < 1) (3.24)

where β ≥ 1 is a constant parameter. For this scheme, Eq. (3.19) is also valid. Thus it becomes the

nondissipative a scheme when ν = 0 and, therefore, it will not become overly dissipative as ν → 0 (i.e., the

scheme’s numerical dissipation is CF L number insensitive). Moreover, because β ≥ 1, Eqs. (3.13) and (3.15)

along with Comment (c) given following Eq. (3.15) imply that the scheme is stable if |ν| < 1 and it becomes

more dissipative (i.e., more stable) as the value of β increases. However, also because β ≥ 1,

lim

|ν|→1−

τ = β ≥ 1 (3.25)

As such, except the special case β = 1, the scheme does not become the c scheme in the limit of |ν| → 1−,

i.e., it does not meet all the requirement of an ideal solver.

Note that, unless specified otherwise, in the remainder of this paper we consider only the

c-τ ∗ schemes or the special c-τ scheme defined by Eq. (3.25).

Next we will introduce weighted-averaging extensions of the c-τ schemes along with two classes of

advanced weighted-averaging techniques.

4. Other c-τ ∗ extensions and new weighted-averaging techniques

To proceed, let

(ˆux ̄−)

n

j

def

= ∆x

4

( un

j − u(P −)

(1 + τ)∆x/4

)

≡ un

j − u(P −)

(1 + τ) (4.1)

15
and

(ˆux ̄+)

n

j

def

= ∆x

4

( u(P +) − un

j

(1 + τ)∆x/4

)

≡ u(P +) − un

j

(1 + τ) (4.2)

Because |AP−| = |AP+| = (1 + τ)∆x/4 (see Fig. 3), it is easy to see that (ˆux ̄−)n

j and (ˆux ̄+)nj are two

normalized one-sided difference approximations of ∂u/∂x at the mesh point (j, n) with one being evaluated

from the left and another from the right. Also, it follows immediately from Eqs. (3.10), (4.1) and (4.2) that

(ˆux ̄)

n

j = 1

2

$

(ˆux ̄−)

n

j + (ˆux ̄+)

n

j

% (4.3)

Moreover, by (i) substituting Eqs. (2.10), (3.8) and (3.9) into Eqs.(4.1) and (4.2), and (ii) using Eqs. (3.11)

and (3.19), one arrives at the conclusion that

(ˆux ̄−)

n

j = (ˆux ̄+)

n

j = (ˆux ̄)

n

j (ν = 0) (4.4)

when ν = 0.

With the above preliminaries, several extensions of the c-τ ∗ scheme will be constructed in the following

subsections.

4.1. Scheme w-1

A comparison of Eqs. (4.1)–(4.3) with Eqs. (2.34)–(2.36) reveals that an obvious extension of the c-τ ∗

scheme can be obtained by replacing (ux ̄−)n

j and (ux ̄+)nj in Eqs. (2.37) and (2.42) with (ˆux ̄−)nj and (ˆux ̄+)nj ,

respectively. In other words, the new extension is formed by Eq. (2.10) and

(ux ̄)

n

j = (w−)

n

j (ˆux ̄−)

n

j + (w+)

n

j (ˆux ̄+)

n

j (4.5)

with

(w±)

n

j = W±((ˆux ̄−)

n

j ,(ˆux ̄+)

n

j , α) (4.6)

Because the scheme is the first extension of the c-τ ∗ scheme in which (ux ̄)n

j is expressed as an weighted

average of (ˆux ̄−)n

j and (ˆux ̄+)nj , for simplicity, hereafter it will be referred to as Scheme w-1. It has been

shown numerically that Scheme w-1 is stable if |ν| < 1 and α ≥ 0. However, as will be shown, Scheme w-1

is deficient in one key aspect.

Eqs. (4.4) and (4.6) along with the definition of W± imply that, for any given α ≥ 0, (w−)n

j → 1/2

and (w+)n

j → 1/2 as ν → 0. In other words, for Scheme w-1, the “weighted” average on the right side of

Eq. (4.5) will approach the simple average as ν → 0. According to an explanation given in the last paragraph

of Sec. 2.3, this implies that Scheme w-1 will lose its capability to suppress wiggles or overshoots when ν

becomes small. For this reason, even though the Euler version of Scheme w-1 performs much better than

that of the w-α scheme in its ability to resolve shocks and contact discontinuities crisply in a wide range

(from 1 to less than 0.001) of the global CF L number (i.e., the maximal value of local CF L numbers), it is

handicapped by the fact that wiggles or overshoots can appear near a discontinuity in a generated solution

when the local CF L number there becomes extremely small. In the following, it will be shown that this

weakness can be overcome by simple modifications of Eq. (4.6).

4.2. Scheme w-2

A new scheme, referred to as Scheme w-2 is formed by Eqs. (2.10) and (4.5) with (w±)n

j being given by

Eq. (2.42). In other words, although (ux ̄)n

j is still constructed as an weighted average of (ˆux ̄−)nj and (ˆux ̄+)nj ,

the associated weight factors (w±)n

j are evaluated using (ux ̄−)nj and (ux ̄+)nj (which are defined in Eqs. (2.34)

and (2.35), respectively). Because the last two parameters, respectively, are identical to the special cases of

(ˆux ̄−)n

j and (ˆux ̄+)nj with τ ≡ 1 (see Eqs. (2.15), (2.34), (2.35), (3.8), (3.9), (4.1), and (4.2)), their values do

not vary with ν. As such, (w±)n

j "= 1/2 and therefore the weighted average on the right side of Eq. (4.5) will

not turn into a simple average when ν = 0. In other words, Scheme w-2 is still capable of annihilating the

16
numerical wiggles near a discontinuity even if ν becomes small. It has been shown numerically that Scheme

w-2 again is stable if |ν| < 1 and α ≥ 0.

Note that a possible drawback of Scheme w-2 is that the relation |(ux ̄−)n

j | < |(ux ̄+)nj | (|(ux ̄−)nj | >

|(ux ̄+)n

j |) does not automatically follow from |(ˆux ̄−)nj | < |(ˆux ̄+)nj | (|(ˆux ̄−)nj | > |(ˆux ̄+)nj |) and vice versa.

As a result, at some local mesh points, it may happen that, of (ˆux ̄−)n

j and (ˆux ̄+)nj , the one with smaller

absolute value may not be associated with a weight factor > 1/2. According to a discussion given in the last

paragraph of Sec. 2.3, this implies that there is no guarantee that, at all localities, the weighted-averaging

induced numerical dissipation will be available to suppress wiggles or overshoots. Despite this possible failing,

fortunately it has been demonstrated numerically that, not only are they capable of suppressing wiggles or

overshoots robustly, Scheme w-2 and its Euler extensions are also highly accurate.

In the following, schemes that overcome the weakness of Scheme w-1 and also avoid the theoretically

possible failing associated with Scheme w-2 will be constructed using two new classes of weighted-averaging

techniques more advanced than that given in Eqs. (2.40) and (2.41).

4.3. New classes of weighted-averaging techniques

To pave the way, first we shall discuss a limitation of Eqs. (2.40) and (2.41) as a generator of weight

factors.

Let x± "= 0. Then, for a given α > 0, obviously W− → 1/2 and W+ → 1/2 as |x+/x−| → 1. As such,

when |x+/x−| is very close to 1, then both W− and W+ will be very close to 1/2 unless α - 1. As a result,

in case that (i) (ˆux ̄±)n

j "= 0, (ii) |(ˆux ̄+)nj /(ˆux ̄−)nj | is very close to 1 (e.g., when ν is very close to 0); and (iii)

Eqs. (4.6) is assumed, then the only way to prevent the weighted average that appears on the right side

of Eq. (4.5) from becoming almost a simple average (and thus causing the overshoot problem discussed in

Sec. 4.1) is to increase the value of α used. However, this approach may be impracticable because numerical

evaluation of a quantity such as xα for any real number x generally is hampered by round-off errors and

thus becomes highly inaccurate if the value of α becomes too large, say 100.

To overcome the limitation discussed above, two classes of weighted-averaging techniques much more

potent and flexible than that discussed in Sec. 2.3 were presented in Sec. 4.3 of [31]. These new techniques

can be used to constructed weighted averages of any set of N (> 1) parameters. By applying these new

techniques for the special case N = 2 [62], one can construct two more extensions of the c-τ ∗ scheme. These

extensions, referred to as Schemes w-3 and w-4, are presented in Sec. 4.4 and 4.5, respectively.

4.4. Scheme w-3

Let

σ = σo

|ν| (σo > 0 and ν "= 0) (4.7)

ξ(j, n; α) def

= |(ˆux ̄+)n

j |

α − |(ˆux ̄−)n

j |

α

|(ˆux ̄+)n

j |

α + |(ˆux ̄−)n

j |

α (α > 0 and |(ˆux ̄−)

n

j | + |(ˆux ̄+)

n

j | > 0) (4.8)

wn

j

def

= 1

2 [1 + σ ξ(j, n; α)] (4.9)

(w#

−)

n

j

def

=







wn

j if 0 ≤ wnj ≤ 1

1 if wn

j > 1

0 if wn

j < 0

(4.10)

and

(w#

+)

n

j

def = 1 − (w#

−)

n

j (4.11)

where σo > 0 and α > 0 are preset parameters in the order of 1. Note that: (i) to avoid dividing by zero, in

practice a small positive number such as 10−20 should be added to each of the denominators of the fractions

that appear in Eqs. (4.7) and (4.8); and (ii) Eqs. (4.10) and (4.11) imply that

(w#

−)

n

j + (w#+)

n

j = 1 (4.12)

17
and

(w#

−)

n

j ≥ 0 and (w#+)

n

j ≥ 0 (4.13)

i.e., (w#

−)n

j and (w#

+)n

j represent a pair of interpolating weight factors. With the aid of these weight factors,

Scheme w-3 is formed by Eq. (2.10) and

(ux ̄)

n

j = (w#−)

n

j (ˆux ̄−)

n

j + (w#+)

n

j (ˆux ̄+)

n

j (4.14)

To explore Scheme w-3, note that an immediate result of Eqs. (4.7)–(4.9) is

wn

j







> 1/2 if |(ˆux ̄+)n

j | > |(ˆux ̄−)nj |

= 1/2 if |(ˆux ̄+)n

j | = |(ˆux ̄−)nj |

< 1/2 if |(ˆux ̄+)n

j | < |(ˆux ̄−)nj |

(ν "= 0) (4.15)

In turn, by using Eq. (4.15), Eqs. (4.10) and (4.11) imply that







0 ≤ (w#

+)n

j < 1/2 < (w#

−)n

j ≤ 1 if |(ˆux ̄+)nj | > |(ˆux ̄−)nj |

(w#

+)n

j = (w#

−)n

j = 1/2 if |(ˆux ̄+)nj | = |(ˆux ̄−)nj |

0 ≤ (w#

−)n

j < 1/2 < (w#

+)n

j ≤ 1 if |(ˆux ̄+)nj | < |(ˆux ̄−)nj |

(ν "= 0) (4.16)

As such, of (ˆux ̄−)n

j and (ˆux ̄+)nj , the one with smaller absolute value is assigned an weight factor > 1/2 in

Eq. (4.14). As explained at the end of Sec. 2.3, it follows that Scheme w-3 also has the capability to suppress

wiggles or overshoots. In fact, as will be shown, it would retain this capability even when |ν| ! 1.

According to Eqs. (4.4) and (4.8), (i) ξ(j, n; α) implicitly is a function of ν, and (ii) ξ(j, n; α) → 0 as

ν → 0. Let ξ(j, n; α) varies linearly with |ν| in the immediate neighborhood of ν = 0, i.e.,

lim

ν→0+(ξ(j, n; α)/|ν|) = ξ+

o (j, n; α) and lim ν→0−(ξ(j, n; α)/|ν|) = ξ−o (j, n; α) (4.17)

where ξ+

o (j, n; α) and ξ−o (j, n; α) are parameters independent of ν. Then, by substituting Eq. (4.7) into

Eq. (4.9) and using Eq. (4.17), one has

lim

ν→0+

wn

j = 1

2

$

1 + σo ξ+

o (j, n; α)

% and lim ν→0−

wn

j = 1

2

$

1 + σo ξ−

o (j, n; α)

% (4.18)

By assumption, σo > 0. Thus, one concludes that

lim

|ν|→0+

wn

j "= 1

2 (4.19)

if

ξ+

o (j, n; α) "= 0 and ξ−o (j, n; α) "= 0 (4.20)

In turn, with the aid of Eq. (4.19), Eqs. (4.10) and (4.11) imply that

lim

|ν|→0+(w#

+)

n

j "= 1

2

and lim

|ν|→0+(w#

−)

n

j "= 1

2 (4.21)

if Eq. (4.20) is assumed. In other words, the expression on the right side of Eq. (4.14) generally would not

reduce to the simple average of (ˆux ̄−)n

j and (ˆux ̄+)nj , and as such, Scheme w-3 would retain its capability to

suppress wiggles and overshoots even when |ν| ! 1. This conclusion was supported by the numerical results

of Euler solvers [31] which were constructed using the weighted-average technique similar to that presented

here.

4.4. Scheme w-4

18
Scheme w-4 is originally presented in [25]. As will be shown in [62], the weight factors used in this

scheme can be constructed using the second class of weighted-averaging techniques reported in Sec. 4.3 of

[31]. Let

(η±)

n

j

def

= |(ˆux ̄±)n

j |

α

min{|(ˆux ̄−)n

j |

α, |(ˆux ̄+)n

j |

α} − 1 (α > 0 and |(ˆux ̄−)

n

j | + |(ˆux ̄+)

n

j | > 0) (4.22)

(  ̃w+)

n

j

def

= 1 + σ (η−)n

j

2 + σ[(η+)n

j + (η−)nj ] (ν "= 0) (4.23)

and

(  ̃w−)

n

j

def

= 1 + σ (η+)n

j

2 + σ[(η+)n

j + (η−)nj ] (ν "= 0) (4.24)

where α > 0 is a preset parameter in the order of 1 while σ is that defined in Eq. (4.7). Note that: (i) to

avoid dividing by zero, in practice a small positive number such as 10−20 should be added to the denominator

of the fraction that appears in Eq. (4.22); and (ii) Eq. (4.22) implies that

(η+)

n

j

5 > 0 if |(ˆux ̄+)n

j | > |(ˆux ̄−)nj |

= 0 if |(ˆux ̄+)n

j | ≤ |(ˆux ̄−)nj | (4.25)

and

(η−)

n

j

5 > 0 if |(ˆux ̄−)n

j | > |(ˆux ̄+)nj |

= 0 if |(ˆux ̄−)n

j | ≤ |(ˆux ̄+)nj | (4.26)

Eqs. (4.23)–(4.26) imply that: (i)

(  ̃w+)

n

j + (  ̃w−)

n

j = 1 (4.27)

and

(  ̃w+)

n

j > 0 and (  ̃w−)

n

j > 0 (4.28)

i.e., (  ̃w+)n

j and (  ̃w−)nj represent a pair of interpolating weight factors; and (ii)







0 < (  ̃w+)n

j < 1/2 < (  ̃w−)nj < 1 if |(ˆux ̄+)nj | > |(ˆux ̄−)nj |

(  ̃w+)n

j = (  ̃w−)nj = 1/2 if |(ˆux ̄+)nj | = |(ˆux ̄−)nj |

0 < (  ̃w−)n

j < 1/2 < (  ̃w+)nj < 1 if |(ˆux ̄+)nj | < |(ˆux ̄−)nj |

(ν "= 0) (4.29)

Scheme w-4 is formed by Eq. (2.10) and

(ux ̄)

n

j = (  ̃w−)

n

j (ˆux ̄−)

n

j + (  ̃w+)

n

j (ˆux ̄+)

n

j (4.30)

According to Eq. (4.29), of (ˆux ̄−)n

j and (ˆux ̄+)nj in Eq. (4.30), the one with the smaller absolute value is

also assigned an weight factor > 1/2. As such, Scheme w-4 also has the capability to suppress wiggles or

overshoots. In fact, if (η±)n

j varies linearly with |ν| in the immediate neighborhood of ν = 0, one would

conclude from the argument that lead to Eq. (4.21) that in general

lim

|ν|→0+(  ̃w+)

n

j "= 1

2

and lim

|ν|→0+(  ̃w−)

n

j "= 1

2 (4.31)

As such the expression on the right side of Eq. (4.30) generally would not reduce to the simple average of

(ˆux ̄−)n

j and (ˆux ̄+)nj , i.e., Scheme w-4 would also retain its capability to suppress wiggles and overshoots even

when |ν| ! 1. This conclusion was supported by the numerical results of Euler solvers [31] which were

constructed using the weighted-average technique similar to that presented here.

A rather thorough review of several 2nd-order 1D CESE schemes has been provided in Sec. 2–4. In the

following, it will be shown that each of these second-order scheme can be “naturally” extended to become its

19
4th-order version. In Sec. 5, first we will construct the a-4 and a-%-4 schemes, i.e., the 4th-order extensions

of the a and a-% schemes, respectively.

5. The a-4 and a-%-4 schemes

To construct the a-4 and a-%-4 schemes, again we consider Eq. (2.1). Also we will use (i) the same mesh

structure and the same CEs and SEs defined in Sec. 2 and (ii) the same definitions of "h∗(x, t; j, n), ν, and

(ux ̄)n

j defined in Eqs. (2.4) and (2.9). However, Eq. (2.3) is now replaced by

u∗(x, t; j, n) def

= un

j + (ux)

n

j (x − xj )+(ut)

n

j (t − t

n) + 1

2

(uxx)

n

j (x − xj )

2

+ (uxt)

n

j (x − xj )(t − t

n) + 1

2

(utt)

n

j (t − t

n)2 +

1

6

(uxxx)

n

j (x − xj )

3

+

1

2

(uxxt)

n

j (x − xj )

2(t − tn) + 1

2

(uxtt)

n

j (x − xj )(t − t

n)2 +

1

6

(uttt)

n

j (t − t

n)3

(5.1)

Here un

j , (ux)nj , (ut)nj , (uxx)nj , (uxt)nj , (utt)nj (uxxx)nj , (uxxt)nj , (uxtt)nj , and (uttt)nj are constants in SE(j, n).

They can be considered as the numerical analogues of the values of u, ∂u/∂x, ∂u/∂t, ∂2u/∂x2, ∂2u/∂x∂t,

∂2u/∂t2, ∂3u/∂x3, ∂3u/∂x2∂t, ∂3u/∂x∂t

2, and ∂3u/∂t3 at the mesh point (j, n), respectively. As such

u∗(x, t; j, n) represents a 3rd-order Taylor’s approximation of u.

For any (j, n) ∈ Ω, let u = u∗(x, t; j, n) be a solution to Eq. (2.1) in SE(j, n), i.e.,

∂u∗(x, t; j, n)

∂t + a

∂u∗(x, t; j, n)

∂x ≡ 0 (x, t) ∈ SE(j, n) (5.2)

Then one has

(ut)

n

j = −a(ux)

n

j , (uxt)

n

j = −a(uxx)

n

j , (utt)

n

j = a2(uxx)

n

j ,

(uxxt)

n

j = −a(uxxx)

n

j , (uxtt)

n

j = a2(uxxx)

n

j , (uttt)

n

j = −a3(uxxx)

n

j

(j, n) ∈ Ω (5.3)

Substituting Eq. (5.3) into Eq. (5.1), one obtains

u∗(x, t; j, n) = un

j + (ux)

n

j

$

(x − xj ) − a (t − t

n)

%

+

1

2

(uxx)

n

j

$

(x − xj ) − a (t − t

n)

%2

+

1

6

(uxxx)

n

j

$

(x − xj ) − a (t − t

n)

%3 (5.4)

Thus, (i) u∗(x, t; j, n) is a 3rd-order polynomial function of $

(x−xj )−a (t−t

n)

%

; and (ii) un

j , (ux)nj , (uxx)nj ,

and (uxxx)n

j are the only independent mesh variables associated with (j, n). Note that, for the same reason

presented immediately following Eq. (2.6), except for flux evaluation, Eq. (5.4) and Eq. (5.7) (see below) are

applicable even for (x, t) ∈/ SE(j, n).

Let

u∗

xx(x, t; j, n) def

= ∂2u∗(x, t; j, n)

∂x2 (5.5)

Then Eqs. (5.2) and (5.4), respectively, imply that

∂u∗

xx(x, t; j, n)

∂t + a

∂u∗

xx(x, t; j, n)

∂x ≡ 0 (x, t) ∈ SE(j, n) (5.6)

and

u∗

xx(x, t; j, n) = (uxx)

n

j + (uxxx)

n

j

$

(x − xj ) − a (t − t

n)

% (x, t) ∈ SE(j, n) (5.7)

i.e., (i) u = u∗

xx(x, t; j, n) is also a solution to Eq. (2.1) in SE(j, n); and (ii) u∗xx(x, t; j, n) is a first-order

polynomial function of $(x − xj ) − a (t − tn)

%

. Note that the fact that Eq. (5.2) implies Eq. (5.6) is but the

special case of the more general result that Eq. (2.1) implies

∂

∂t

6∂2u

∂x2

7

+ a

∂

∂x

6∂2u

∂x2

7

= 0 (5.8)

20
In fact, for a smooth function u(x, t), the conservation condition Eq. (2.2) implies the conservation condition

!

S(V )

"

hxx · d"s = 0 (5.9)

where

"

hxx

def

= ∂2"

h

∂x2 ≡

6

a

∂2u

∂x2 ,

∂2u

∂x2

7

(5.10)

As a result of the above observations on the “equivalence” relation between u and ∂2u/∂x2, the reader

should no longer be surprised by the fact that, by replacing the symbols un

j , (ux)nj , u∗(x, t; j, n), and

"

h∗(x, t; j, n) with (uxx)n

j , (uxxx)nj , u∗

xx(x, t; j, n), and "

h∗

xx(x, t; j, n), respectively, Eqs. (2.4) and (2.6) would

turn into "

h∗

xx(x, t; j, n) def = (au∗xx(x, t; j, n), u∗xx(x, t; j, n)) (x, t) ∈ SE(j, n) (5.11)

(i.e., the numerical analogue of Eq. (5.10)) and Eq. (5.7), respectively. This structural similarity strongly

suggests that the mesh variables (uxx)n

j and (uxxx)nj can be advanced in time by using one of the 2nd-order

CESE schemes described earlier with the understanding that the role played by un

j , (ux)nj , u∗(x, t; j, n), and

"

h∗(x, t; j, n) in a 2nd-order CESE scheme be replaced by (uxx)n

j , (uxxx)nj , u∗

xx(x, t; j, n), and "

h∗

xx(x, t; j, n),

respectively. Note that this simplified approach is justified by the observation that (uxx)n

j and (uxxx)nj ,

respectively, appear in the 2nd-order and 3rd-order terms on the right side of Eq. (5.4) and therefore the

4th-order accuracy requirement of u∗(x, t; j, n) demands only that (uxx)n

j and (uxxx)nj , respectively, be 2ndorder and 1st-order in accuracy, i.e., the same orders of accuracy achieved by un

j and (ux)nj in a 2nd-order

CESE scheme.

For the reason presented above, equations such as Eqs. (5.11) and (5.7) hereafter will be designated as

the substitution forms of Eqs. (2.4) and (2.6), respectively. With the aid of this designation, in the following

(uxx)n

j and (uxxx)nj will be advanced in time using the a and a-% schemes.

5.1. Time marching of (uxx)n

j and (uxxx)nj

By imposing the substitution forms of Eqs. (2.7) and (2.8), i.e.,

!

S(CE+(j,n))

"

h∗

xx · d"s = 0 (j, n) ∈ Ω (5.12)

and

!

S(CE−(j,n))

"

h∗

xx · d"s = 0 (j, n) ∈ Ω (5.13)

and using Eqs. (5.7) and (5.11), one can show that (i)

(ux ̄x ̄)

n

j = 1

2

&

(1 − ν) [ux ̄x ̄ − (1 + ν)ux ̄x ̄x ̄]

n−1/2

j+1/2 + (1 + ν) [ux ̄x ̄ + (1 − ν)ux ̄x ̄x ̄]

n−1/2

j−1/2

'

(j, n) ∈ Ω (5.14)

and, assuming |ν| "= 1, (ii)

(ux ̄x ̄x ̄)

n

j = (uax ̄x ̄x ̄)

n

j (j, n) ∈ Ω (5.15)

with

(ua

x ̄x ̄x ̄)

n

j

def

= 1

2

&

[ux ̄x ̄ − (1 + ν)ux ̄x ̄x ̄]

n−1/2

j+1/2 − [ux ̄x ̄ + (1 − ν)ux ̄x ̄x ̄]

n−1/2

j−1/2

'

(5.16)

where ν, (ux ̄)n

j , (ux ̄x ̄)nj , and (ux ̄x ̄x ̄)nj are defined by Eq. (2.9) and

(ux ̄x ̄)

n

j

def

=

8 ∆x

4

92

(uxx)

n

j and (ux ̄x ̄x ̄)

n

j

def

=

8 ∆x

4

93

(uxxx)

n

j (j, n) ∈ Ω (5.17)

21
It can be shown that Eqs. (5.14)–(5.15) are equivalent to the substitution form of Eqs. (2.10) and (2.11).

As such, the time marching scheme formed by Eqs. (5.14) and (5.15) literally is the same a scheme defined

earlier. Also, one can shown easily that Eq. (5.14) is equivalent to the substitution form of Eq. (2.13), i.e.,

!

S(CE(j,n))

"

h∗

xx · d"s = 0 (j, n) ∈ Ω (5.18)

Note that, as explained in Comment (d) given at the end of Sec. 2.1, Eq. (5.18) is a direct result of Eqs. (5.12)

and (5.13).

Next let

u∗

x ̄x ̄(x, t; j, n) def

=

8 ∆x

4

92

u∗

xx(x, t; j, n) and (ux ̄xt ̄ )(j, n) def

=

8 ∆x

4

92

(uxxt)

n

j (j, n) ∈ Ω (5.19)

Then the substitution forms of Eqs. (2.14)–(2.18) can be expressed as

(u#

x ̄x ̄)

n

j±1/2

def

= u∗

x ̄x ̄

"

xj±1/2, tn; j ± 1/2, n − 1/2

#

= [ux ̄x ̄ + (∆t/2)ux ̄xt  ̄ ]

n−1/2

j±1/2 (5.20)

(u#

x ̄x ̄)

n

j±1/2 = (ux ̄x ̄ − 2ν ux ̄x ̄x ̄)

n−1/2

j±1/2 (5.21)

(uc

x ̄x ̄x ̄)

n

j

def

= ∆x

4

((u#

x ̄x ̄)

n

j+1/2 − (u#

x ̄x ̄)

n

j−1/2

∆x

)

≡ (u#

x ̄x ̄)

n

j+1/2 − (u#

x ̄x ̄)

n

j−1/2

4 (5.22)

(uc

x ̄x ̄x ̄)

n

j = 1

4

*

(ux ̄x ̄ − 2νux ̄x ̄x ̄)

n−1/2

j+1/2 − (ux ̄x ̄ − 2νux ̄x ̄x ̄)

n−1/2

j−1/2

+

(5.23)

and

(uc

x ̄x ̄x ̄ − uax ̄x ̄x ̄)

n

j = 1

4

*

(ux ̄x ̄ + 2ux ̄x ̄x ̄)

n−1/2

j−1/2 − (ux ̄x ̄ − 2ux ̄x ̄x ̄)

n−1/2

j+1/2

+

(5.24)

respectively. Note that the validity of the last equality sign in Eq. (5.20) can be proved using Eqs. (5.7),

(5.17) and (5.19) along with the relation (uxxt)n

j = −a(uxxx)nj which is part of Eq. (5.3).

Moreover, let the symbol % be replaced by %1. Then the substitution form of Eq. (2.19) can be expressed

as

(ux ̄x ̄x ̄)

n

j = (uax ̄x ̄x ̄)

n

j + 2%1(ucx ̄x ̄x ̄ − uax ̄x ̄x ̄)

n

j (5.25)

Obviously Eq. (5.25) reduces to (i) (ux ̄x ̄x ̄)n

j = (ua

x ̄x ̄x ̄)n

j if %1 = 0; and (ii) (ux ̄x ̄x ̄)nj = (uc

x ̄x ̄x ̄)n

j if %1 = 1/2. With

the aid of Eqs. (5.16) and (5.24), Eq. (5.25) can be rewritten as

(ux ̄x ̄x ̄)

n

j = 1

2

&

[(1 − %1)ux ̄x ̄ + (2%1 − 1 − ν)ux ̄x ̄x ̄]

n−1/2

j+1/2 − [(1 − %1)ux ̄x ̄ + (1 − ν − 2%1)ux ̄x ̄x ̄]

n−1/2

j−1/2

'

(5.26)

Obviously Eq. (5.26) is the substitution form of Eq. (2.19a) if the symbol % in the latter equation is replaced

by %1.

Because (i) (ux ̄x ̄)n

j and (ux ̄x ̄x ̄)nj can be explicitly evaluated in terms of the mesh variables at the (n −

1/2)th time level using Eqs. (5.14) and (5.26), and (ii) the latter equations are equivalent to the substitution

forms of Eqs. (2.10) and (2.19a), it has been shown that indeed (ux ̄x ̄)n

j and (ux ̄x ̄x ̄)nj can be advanced in time

using the a-% scheme.

This completes the description of the first step in the construction of the a-%-4 scheme. In the following,

we describe the second step, i.e., the details of how un

j and (ux ̄)nj can be advanced in time. As a preliminary,

we conclude this subsection with the following remarks:

(a) As a result of its construction, the scheme formed by Eqs. (5.14) and (5.26) observes both the conservation relations Eqs. (5.12) and (5.13) if and only if %1 = 0.

22
(b) Because Eq. (5.14) is equivalent to Eq. (5.18), the scheme formed by Eq. (5.14) and (5.26) still observes

the less stringent conservation relation Eq. (5.18) if %1 "= 0.

(c) Because of how they are derived, by replacing the symbols %1, ux ̄x ̄, and ux ̄x ̄x ̄ with %, u, and ux ̄, respectively, Eqs. (5.14) and (5.26), respectively, will become Eq. (2.10) and (2.19a), the defining equations of

the a-% scheme.

5.2. Time marching of un

j and (ux)nj

Let

(s+)

n

j

def

=

/

u − (1 + ν)ux ̄ + 2(1 + ν + ν2)

3

ux ̄x ̄ − (1 + ν)(1 + ν2)

3

ux ̄x ̄x ̄

0n

j

(j, n) ∈ Ω (5.27)

(s−)

n

j

def

=

/

u + (1 − ν)ux ̄ + 2(1 − ν + ν2)

3

ux ̄x ̄ + (1 − ν)(1 + ν2)

3

ux ̄x ̄x ̄

0n

j

(j, n) ∈ Ω (5.28)

Then, assuming Eqs. (2.4), (2.9), (5.4), and (5.17), it can be shown that the conservation relations Eqs. (2.7)

and (2.8) are equivalent to

(1 − ν)

:/

u + (1 + ν)ux ̄ + 2(1 + ν + ν2)

3

ux ̄x ̄ + (1 + ν)(1 + ν2)

3

ux ̄x ̄x ̄

0n

j

− (s+)

n−1/2

j+1/2

;

= 0 (5.29)

and

(1 + ν)

:/

u − (1 − ν)ux ̄ + 2(1 − ν + ν2)

3

ux ̄x ̄ − (1 − ν)(1 + ν2)

3

ux ̄x ̄x ̄

0n

j

− (s−)

n−1/2

j−1/2

;

= 0 (5.30)

respectively.

By adding Eqs. (5.29) and (5.30) together, one has

un

j = 1

2

*

(1 − ν)(s+)

n−1/2

j+1/2 + (1 + ν)(s−)

n−1/2

j−1/2

+

− 2

3

(ux ̄x ̄)

n

j (j, n) ∈ Ω (5.31)

Let 1 − ν2 "= 0, i.e., 1 − ν "= 0 and 1 + ν "= 0. Then Eqs. (5.29) and (5.30) can be divided by (1 − ν) and

(1 + ν), respectively. By subtracting the resulting equations from each other, one has

(ux ̄)

n

j = (ua-4 x ̄ )

n

j (j, n) ∈ Ω (5.32)

where

(ua-4 x ̄ )

n

j

def

= 1

2

*

(s+)

n−1/2

j+1/2 − (s−)

n−1/2

j−1/2

+

− 2ν

3 (uxx)

n

j − 1

3

"

1 + ν2#

(ux ̄x ̄x ̄)

n

j (j, n) ∈ Ω (5.33)

Hereafter the superscript “4” is used to distinguish the parameters introduced in the current construction

of the 4th-order schemes from their counterparts that appear in the 2nd-order schemes.

At this juncture, note that:

(a) The zero-order and first-order terms in the third-order expansion on the right side of Eq. (5.4) are

identical to those in the first-order expansion on the right side of Eq. (2.6), and (ii) Eqs. (5.31) and

(5.32) were obtained by imposing the conservation relations Eqs. (2.7) and (2.8) and assuming Eqs. (2.4)

and (5.4) while Eqs. (2.10) and (2.11) were obtained by imposing the same conservation relations and

assuming Eqs. (2.4) and (2.6). As such, one would expect that Eqs. (5.31) and (5.32) differ from

Eqs. (2.10) and (2.11) only in the appearance of the extra 2nd-order and 3rd-order terms in the former

equations. With the aid of Eqs. (2.12), (5.27), (5.28), and (5.33), the reader can prove for himself that

this is indeed true.

23
(b) Assuming Eqs. (2.4), (2.9), (5.4), (5.17), (5.27), and (5.28), it has been shown that Eqs. (2.7) and (2.8)

imply (i) Eq. (5.31) for all ν; and (ii) Eq. (5.32) if |ν| "= 1. On the other hand, it can be shown that

Eqs. (5.31) and (5.32) are equivalent to

/

u + (1 + ν)ux ̄ + 2(1 + ν + ν2)

3

ux ̄x ̄ + (1 + ν)(1 + ν2)

3

ux ̄x ̄x ̄

0n

j

− (s+)

n−1/2

j+1/2 = 0 (5.29a)

and

/

u − (1 − ν)ux ̄ + 2(1 − ν + ν2)

3

ux ̄x ̄ − (1 − ν)(1 + ν2)

3

ux ̄x ̄x ̄

0n

j

− (s−)

n−1/2

j−1/2 = 0 (5.30a)

for all ν. Because (i) Eqs. (5.29a) and (5.30a) imply Eqs. (5.29) and (5.30), respectively; and (ii)

Eqs. (5.29) and (5.30) are equivalent to Eqs. (2.7) and (2.8), respectively, one concludes that, for all ν,

the last two conservation conditions are enforced by Eqs. (5.31) and (5.32).

(c) Under the current assumptions, it can be shown that Eq. (5.31) is equivalent to the conservation relation

Eq. (2.13).

For Eqs. (5.31) and (5.32) to be considered as a complete marching scheme for un

j and (ux ̄)nj , the

mesh variables (ux ̄x ̄)n

j and (ux ̄x ̄x ̄)nj which appear on right sides of Eqs. (5.31) and (5.33) must be expressed

explicitly in terms of the mesh variables at the (n − 1/2)th time level. In the current construction, this is

done by demanding that the former mesh variables be evaluated in terms of the latter mesh variables using

Eqs. (5.14) and (5.26). By using this assumption and Eqs. (5.27) and (5.28), Eqs. (5.31)–(5.33) can now be

cast into the forms by which un

j and (ux ̄)nj can be evaluated in terms of the mesh variables at the (n−1/2)th

time level, i.e.,

un

j = 1 − ν

2

/

u − (1 + ν)ux ̄ +

2ν(1 + ν)

3

ux ̄x ̄ + (1 + ν)(1 − ν2)

3

ux ̄x ̄x ̄

0n−1/2

j+1/2

+

1 + ν

2

/

u + (1 − ν)ux ̄ − 2ν(1 − ν)

3

ux ̄x ̄ − (1 − ν)(1 − ν2)

3

ux ̄x ̄x ̄

0n−1/2

j−1/2

(j, n) ∈ Ω (5.34)

and

(ux ̄)

n

j = $

ua-4 x ̄ (%1)

%n

j (j, n) ∈ Ω (5.35)

where

$

ua-4 x ̄ (%1)

%n

j

def

= 1

2

5

u − (1 + ν)ux ̄ +

1+3ν2 + %1(1 + ν2)

3

ux ̄x ̄ +

2

$

ν(1 − ν2) − %1(1 + ν2)

%

3

ux ̄x ̄x ̄

<n−1/2

j+1/2

− 1

2

5

u + (1 − ν)ux ̄ +

1+3ν2 + %1(1 + ν2)

3

ux ̄x ̄ +

2

$

ν(1 − ν2) + %1(1 + ν2)

%

3

ux ̄x ̄x ̄

<n−1/2

j−1/2

(j, n) ∈ Ω

(5.36)

At this juncture, note that:

(a) As a result of how Eqs. (5.34)–(5.36) are derived, one concludes that: (i) Eqs. (5.14) and (5.34) satisfy

Eq. (5.31) automatically; and (ii) Eqs. (5.14), (5.26), and (5.35) satisfy Eqs. (5.32) automatically.

Because the conservation relations Eq. (2.7) and (2.8), in turn, are enforced by Eqs. (5.31) and (5.32)

for all ν (see Comment (b) given following Eq. (5.33)), the two conservation relations are observed by

the scheme formed by Eqs. (5.14), (5.26), (5.34), and (5.35) for all ν and %1.

(b) It is pointed out earlier (see comments (a) and (b) given at the end of Sec. 5.1) that the scheme formed

by Eqs. (5.14) and (5.26) observes (i) the conservation relations Eqs. (5.12) and (5.13) if %1 = 0; and

(ii) only the conservation relation Eq. (5.18) if %1 "= 0.

24
(c) As a result of the comments made in the above items (a) and (b), one concludes that the special scheme

formed by Eqs. (5.14), (5.26), (5.34), and (5.35) with %1 = 0 observes the four conservation relations

Eqs. (2.7), (2.8), (5.12), and (5.13). In fact, for this special scheme, the time evolution of the four

independent mesh variables un

j , (ux ̄)nj , (ux ̄x ̄)nj , and (ux ̄x ̄x ̄)nj is completely determined by these four

conservation relations. As such, this special scheme is referred to as the a-4 scheme.

(d) On the other hand, for the case %1 "= 0, the scheme defined in the above item (a) satisfies only three

conservation relations, i.e., Eqs. (2.7), (2.8), and (5.18). In the following, this scheme will be extended

so that the extension satisfies only the two least stringent conservation relations Eq. (2.13) and (2.18).

To proceed, first we will construct the current version of (uc

x ̄)n

j which was defined in Eq. (2.16). As a

preliminary, Taylor’s expansion is invoked to obtain

u

"

xj ± [(1 + τ)∆x/4], tn#

=

&

1 ± [(1 + τ)∆x/4]∂/∂x + (1/2)[(1 + τ)∆x/4]2∂2/∂x2

± (1/6)[(1 + τ)∆x/4]3∂3/∂x3 + (1/24)[(1 + τ)∆x/4]4∂4/∂x4

'

u(xj , tn) + O

$

(∆x)

5%

(τ > −1) (5.37)

Here (i) τ > −1 is an adjustable parameter, and (ii) (xj + [(1 + τ)∆x/4], tn) and (xj − [(1 + τ)∆x/4], tn), are

the coordinates of points P + and P − depicted in Fig. 3, respectively. By subtracting the two equations in

Eq. (5.37) from each other and then rearranging the resulting equation, one has

∂u

∂x(xj , tn) = u (xj + [(1 + τ)∆x/4], tn) − u (xj − [(1 + τ)∆x/4], tn)(1 + τ)∆x/2 − [(1 + τ)∆x]

2

96

∂3u

∂x3 (xj , tn) + O $

(∆x)

4%

(5.38)

Next, by using Eqs. (3.6) and (3.7), one has

u(xj ± [(1 + τ)∆x/4], tn) = u(xj±1/2 ∓ [(1 − τ)∆x/4], tn−1/2 + ∆t/2) (5.39)

The expression on the right side of Eq. (5.39) can be expressed as a Taylor’s expansion about the point

(xj±1/2, tn−1/2). Substituting the expansion into Eq. (5.39), one has

u (xj ± [(1 + τ)∆x/4], tn) = u4(P ±) + 4th-order terms involving ∆t and ∆x (5.40)

where

u4(P ±) def

=

&

1 + =(∆t/2)∂/∂t ∓ [(1 − τ)∆x/4]∂/∂x

>

+ (1/2)=(∆t/2)∂/∂t ∓ [(1 − τ)∆x/4]∂/∂x

>2

+ (1/6)=(∆t/2)∂/∂t ∓ [(1 − τ)∆x/4]∂/∂x

>3

'

u(xj±1/2, tn−1/2)

(5.41)

Hereafter, we consider only a mesh refinement procedure in which the ratio ∆t/∆x and the parameter

τ are held constant as ∆t → 0 and ∆x → 0. Because the advection speed a is a constant, in such a mesh

refinement procedure, ν is also held constant as ∆t → 0 and ∆x → 0. Moreover, with this assumption, the

4th-order terms which appear on the right side of Eq. (5.40) can be expressed as O $(∆x)4%

, i.e., Eq. (5.40)

can be expressed as

u (xj ± [(1 + τ)∆x/4], tn) = u4(P ±) + O $

(∆x)

4% (5.42)

In turn, by substituting Eqs. (5.42) into Eq. (5.38), one has

∂u

∂x(xj , tn) = u4(P +) − u4(P −)(1 + τ)∆x/2 − [(1 + τ)∆x]

2

96

∂3u

∂x3 (xj , tn) + O $

(∆x)

3% (5.43)

25
For any (j, n) ∈ Ω, u4(P +) and u4(P −) will be approximated by

u4∗(P +) def

= u∗(xj + [(1 + τ)∆x/4], tn; j + 1/2, n − 1/2)

≡ $u − (1 − τ + 2ν)ux ̄ + (1/2)(1 − τ + 2ν)

2 ux ̄x ̄ − (1/6)(1 − τ + 2ν)3 ux ̄x ̄x ̄

%n−1/2

j+1/2

(5.44)

and

u4∗(P −) def

= u∗(xj − [(1 + τ)∆x/4], tn; j − 1/2, n − 1/2)

≡ $u + (1 − τ − 2ν)ux ̄ + (1/2)(1 − τ − 2ν)

2 ux ̄x ̄ + (1/6)(1 − τ − 2ν)3 ux ̄x ̄x ̄

%n−1/2

j−1/2

(5.45)

respectively. According to Eqs. (3.6), (3.7), and (5.1)–(5.4), u4∗(P +) is the 3rd-order Taylor’s approximation

(i.e., the error is O $(∆x)4%

) of u at P + evaluated using the mesh variables at (j+ 1/2, n−1/2) while u4∗(P −)

is that at P − evaluated using the mesh variables at (j − 1/2, n − 1/2). Note that the validity of the last

equality signs on the right sides of Eqs. (5.44) and (5.45) follows from Eqs. (2.9), (3.6), (3.7), (5.4) and

(5.17). Moreover, according to Eq. (5.43), the parameter

$

u4

x(τ)

%n

j

def

= u4∗(P +) − u4∗(P −)

(1 + τ)∆x/2 − [(1 + τ)∆x]

2

96 (uxxx)

n

j (5.46)

represents a numerical analogue of the value of ∂u/∂x at (j, n) with 3rd-order accuracy. This is consistent

with the 4th-order accuracy requirement of u∗(x, t; j, n) because (ux)n

j (the numerical analogue of ∂u/∂x)

appears in the first-order term on the right side of Eq. (5.4). By using Eq. (5.17) and

$

u4

x ̄(τ)

%n

j

def

= ∆x

4

$

u4

x(τ)

%n

j (5.47)

Eq. (5.46) can be cast into its normalized form, i.e.,

$

u4

x ̄(τ)

%n

j

def

= u4∗(P +) − u4∗(P −)

2(1 + τ) − (1 + τ)26 (ux ̄x ̄x ̄)

n

j (5.48)

Because points P + and P − coincide with the mesh points (j + 1/2, n) and (j − 1/2, n), respectively,

when τ = 1, the current version of (uc

x ̄)n

j is defined by

(uc-4 x ̄ )

n

j

def

= $u4

x ̄(1)%n

j = 1

4

*

u − 2νux ̄ + 2ν2ux ̄x ̄ − 4ν3

3

ux ̄x ̄x ̄

+n−1/2

j+1/2

− 1

4

*

u − 2νux ̄ + 2ν2ux ̄x ̄ − 4ν3

3

ux ̄x ̄x ̄

+n−1/2

j−1/2 − 2

3

(ux ̄x ̄x ̄)

n

j

(j, n) ∈ Ω (5.49)

The validity of the last equality sign in Eq. (5.49) can be proved by using Eqs. (5.44), (5.45), and (5.48) with

τ = 1. Substituting Eq. (5.26) into the last term on the right side of Eq. (5.49) and denoting the resulting

expression as $uc-4 x ̄ (%1)

%n

j , one concludes that

$

uc-4 x ̄ (%1)

%n

j = 1

4

5

u − 2ν ux ̄ +

*

2ν2 − 4(1 − %1)

3

+

ux ̄x ̄ + 4(1 + ν − ν3 − 2%1)

3

ux ̄x ̄x ̄

<n−1/2

j+1/2

− 1

4

5

u − 2ν ux ̄ +

*

2ν2 − 4(1 − %1)

3

+

ux ̄x ̄ − 4(1 − ν + ν3 − 2%1)

3

ux ̄x ̄x ̄

<n−1/2

j−1/2

(j, n) ∈ Ω (5.50)

26
By definition, the a-%-4 scheme is formed by Eqs. (5.14), (5.26), (5.34), and

(ux ̄)

n

j = $

ua-4 x ̄ (%1)

%n

j + 2%2

$

uc-4 x ̄ (%1) − ua-4 x ̄ (%1)

%n

j (j, n) ∈ Ω (5.51)

where %2 is another adjustable parameter. With the aid of Eqs. (5.36) and (5.50), Eq. (5.51) can be expressed

as

(ux ̄)

n

j = 1

6

&

3(1 − %2)u − 3(1 + ν − 2%2)ux ̄ + $

1+3ν2 − 6%2 + (1 + ν2)%1 + 2(1 − ν2)%1%2

%

ux ̄x ̄

+ 2$ν(1 − ν2)+2%2 − (1 + ν2)%1 − 2(1 − ν2)%1%2

%

ux ̄x ̄x ̄

'n−1/2

j+1/2

− 1

6

&

3(1 − %2)u + 3(1 − ν − 2%2)ux ̄ + $

1+3ν2 − 6%2 + (1 + ν2)%1 + 2(1 − ν2)%1%2

%

ux ̄x ̄

+ 2$ν(1 − ν2) − 2%2 + (1 + ν2)%1 + 2(1 − ν2)%1%2

%

ux ̄x ̄x ̄

'n−1/2

j−1/2 (j, n) ∈ Ω

(5.52)

We conclude this subsection with the following remarks:

(a) Eqs. (5.51) reduces to (i) (ux ̄)n

j = $

ua-4 x ̄ (%1)

%n

j if %2 = 0; and (ii) (ux ̄)n

j = $

uc-4 x ̄ (%1)

%n

j if %2 = 1/2. Because

the a-4 scheme is the special scheme formed by Eqs. (5.14), (5.26), (5.34) and (5.35) with %1 = 0, the

a-%-4 scheme reduces to the a-4 scheme when %1 = %2 = 0.

(b) Because of how they are derived, by replacing the symbol %2 with % and dropping all the 2nd-order and

3rd-order terms, Eqs. (5.34) and (5.52), respectively, will become Eqs. (2.10) and (2.19a), the defining

equations of the a-% scheme.

In the following, we will discuss the stability of the a-%-4 scheme. As a preliminary, the scheme will be

cast into a matrix form.

5.3. Matrix form of the a-%-4 scheme

The component equations of the a-%-4 scheme, i.e., Eqs. (5.14), (5.26), (5.34) and (5.52), can be cast

into the matrix equation

"q(j, n) = Q−(%1, %2, ν) "q(j + 1/2, n − 1/2) + Q+(%1, %2, ν) "q(j − 1/2, n − 1/2) (j, n) ∈ Ω (5.53)

where

"q(j, n) def

=





un

j

(ux ̄)n

j

(ux ̄x ̄)n

j

(ux ̄x ̄x ̄)n

j





(j, n) ∈ Ω (5.54)

and

Q±(%1, %2, ν) def

=





q±

11 q±12 q±13 q±14

q±

21 q±22 q±23 q±24

0 0 q±

33 q±34

0 0 q±

43 q±44





(5.55)

with

q±

11

def

= 1 ± ν

2 , q±

12

def

= ±

1 − ν2

2 , q±

13

def

= ∓ν(1 − ν2)

3 , q±

14

def

= ∓(1 − ν2)2

6 (5.56)

27
q±

21

def

= ∓

1 − %2

2 , q±

22

def

= 2%2 − 1 ± ν

2 , q±

23

def

= ∓

1+3ν2 − 6%2 + (1 + ν2)%1 + 2(1 − ν2)%1%2

6

q±

24

def

= ∓ν(1 − ν2) ∓ 2%2 ± (1 + ν2)%1 ± 2(1 − ν2)%1%2

3

(5.57)

q±

33

def

= 1 ± ν

2 , q±

34

def

= ±

1 − ν2

2 , q±

43

def

= ∓

1 − %1

2 , q±

44

def

= 2%1 − 1 ± ν

2 (5.58)

For simplicity, hereafter Q+(%1, %2, ν) and Q−(%1, %2, ν) may be abbreviated as Q+ and Q−, respectively.

Eq. (5.53) represents the form of the a-%-4 scheme by which mesh variables can be advanced over ∆t/2.

By applying Eq. (5.53) successively over two consecutive marching steps, one obtain a compounded form by

which mesh variables can be advanced over ∆t, i.e.,

"q(j, n + 1) = (Q+)

2"q(j − 1, n) + (Q+Q− + Q−Q+)"q(j, n)+(Q−)2"q(j + 1, n) (j, n) ∈ Ω (5.59)

Eq. (5.59) has the advantage that the mesh variables involved are associated with mesh points which are not

staggered in space-time.

To further explore the relations between the a-%-4 and a-% schemes, note that the component equations

of the a-% scheme, i.e., Eqs. (2.10) and (2.19a), can be cast into the matrix form

"q(j, n) = M−(%, ν) "q(j + 1/2, n − 1/2) + M+(%, ν) "q(j − 1/2, n − 1/2) (j, n) ∈ Ω (5.60)

where

"q(j, n) def

=





un

j

(ux ̄)n

j



 (j, n) ∈ Ω (5.61)

and

M±(%, ν) def

= 1

2





1 ± ν ±(1 − ν2)

∓(1 − %) 2% − 1 ± ν



 − ∞ < %< +∞ (5.62)

With the aid of Eqs. (5.56)–(5.58) and (5.62), it can be shown that





q±

11 q±12

q±

21 q±22



 = M±(%2, ν) (5.63)

and





q±

33 q±34

q±

43 q±44



 = M±(%1, ν) (5.64)

By comparing Eqs. (5.63) and (5.64) with Eq. (5.55), one now arrives at the conclusion that the 2 × 2

submatrices in the top-left and bottom right corners of the 4 × 4 matrix Q±(%1, %2, ν) are M±(%2, ν) and

M±(%1, ν), respectively.

Because of the relations Eqs. (5.63) and (5.64), as will be shown in the next subsection, the a-%-4 and

a-% schemes share essentially the same stability conditions.

5.4. Stability conditions of the a-%-4 scheme

Note that a scheme such as Eq. (5.53) would become computationally unstable if the round-off errors

introduced initially at any time level are amplified without bound as they propagate down the subsequent

time levels. To be specific, consider the effect of contamination of the initial values by round-off errors. Let

"q(j, n) denote the solution associated with the given exact initial values while "q(j, n) + δ"q(j, n) denote the

solution associated with the contaminated initial values. Then "q(j, n) must satisfy Eq. (5.53). Also, because

28
the coefficient matrices Q+ and Q− in the linear scheme Eq. (5.53) have constant elements independent of

mesh variables, the contaminated solution satisfies the equation

"q(j, n) + δ"q(j, n) = Q−

$

"q(j + 1/2, n − 1/2) + δ"q(j + 1/2, n − 1/2)%

+ Q+

$

"q(j − 1/2, n − 1/2) + δ"q(j − 1/2, n − 1/2)% (j, n) ∈ Ω (5.65)

By subtracting Eq. (5.53) from (5.65), one has

δ"q(j, n) = Q− δ"q(j + 1/2, n − 1/2) + Q+ δ"q(j − 1/2, n − 1/2) (j, n) ∈ Ω (5.66)

Thus the time evolution of δ"q(j, n), which represents the solution errors induced by the initial-value contamination, is also governed by the scheme Eq. (5.53). For this reason, stability of the linear scheme Eq. (5.53)

can be determined by the time-evolution behavior of "q(j, n) itself.

The stability of the a-%-4 scheme will be studied using the von Neumann method. In this method (see

Sec. 4 in [1] for a rigorous treatment based on the discrete Fourier analysis), the time-evolution behavior of

each Fourier component (characterized by a phase angle θ) of "q(j, n) is studied separately. Because these

components propagate in time independently for a linear scheme such as Eq. (5.53), the scheme is stable if

only if each component is bounded as it evolves in time.

Let

"q(j, n) = "q ∗(n, θ)eijθ (j, n) ∈ Ω; −∞ < θ < +∞; i def

= √−1 (5.67)

Here (i) θ is the phase angle variation of a Fourier component over the length ∆x, and (ii) "q ∗(n, θ) is a 2 × 1

column matrix. Substituting Eq. (5.67) into Eq. (5.59), one has

"q ∗(n + 1, θ) = [Q(%1, %2, ν, θ)]2 "q ∗(n, θ) n = 0, ±1/2, ±1, ±3/2,... (5.68)

where

Q(%1, %2, ν, θ) = eiθ/2Q−(%1, %2, ν) + e−iθ/2Q+(%1, %2, ν)

=





eiθ/2q−

11 + e−iθ/2q+11 eiθ/2q−12 + e−iθ/2q+12 eiθ/2q−13 + e−iθ/2q+13 eiθ/2q−14 + e−iθ/2q+14

eiθ/2q−

21 + e−iθ/2q+21 eiθ/2q−22 + e−iθ/2q+22 eiθ/2q−23 + e−iθ/2q+23 eiθ/2q−24 + e−iθ/2q+24

0 0 eiθ/2q−

33 + e−iθ/2q+33 eiθ/2q−34 + e−iθ/2q+34

0 0 eiθ/2q−

43 + e−iθ/2q+43 eiθ/2q−44 + e−iθ/2q+44





(5.69)

Because of Eq. (5.68), [Q(%1, %2, ν, θ)]2 is referred to as the amplification matrix of the a-%-4 scheme per

marching cycle (i.e., per ∆t) while Q(%1, %2, ν, θ) is referred to as the amplification factor per marching step

(i.e., per ∆t/2). Also by using Eq. (5.68), one has

"q ∗(n + m, θ) = [Q(%1, %2, ν, θ)]2m"q ∗(n, θ) m = 1, 2, 3 ... ; n = 0, ±1/2, ±1, ± = 3/2,... (5.70)

There is an important relation between the matrix Q(%1, %2, ν, θ) and the amplification matrix per marching

step of the a-% scheme, i.e.,

M(%, ν, θ) def

= eiθ/2M−(%, ν) + e−iθ/2M+(%, ν) ≡





cos(θ/2) − iν sin(θ/2) −i(1 − ν2) sin(θ/2)

i(1 − %) sin(θ/2) (2% − 1) cos(θ/2) − iν sin(θ/2)





− ∞ < % < +∞ (5.71)

29
In fact, it follow from Eqs. (5.63), (5.64), (5.69), and (5.71) that





eiθ/2q−

11 + e−iθ/2q+11 eiθ/2q−12 + e−iθ/2q+12

eiθ/2q−

21 + e−iθ/2q+21 eiθ/2q−22 + e−iθ/2q+22



 = M(%2, ν, θ) (5.72)

and





eiθ/2q−

33 + e−iθ/2q+33 eiθ/2q−34 + e−iθ/2q+34

eiθ/2q−

43 + e−iθ/2q+43 eiθ/2q−44 + e−iθ/2q+44



 = M(%1, ν, θ) (5.73)

i.e., the 2×2 submatrices in the top-left and bottom-right corners of the matrix Q(%1, %2, ν, θ) are M(%2, ν, θ)

and M(%1, ν, θ), respectively. As such, the amplification matrix of the a-% scheme is identical to each of the

above 2 × 2 submatrices of Q(%1, %2, ν, θ) except that the parameter % in the former matrix is replaced by %2

or %1 in the latter two submatrices. By using this property, it will be shown that, for each set of %1, %2, ν,

and θ, the eigenvalues of Q(%1, %2, ν, θ) can be expressed as λ±(%1, ν, θ) and λ±(%2, ν, θ) where the function

λ±(%, ν, θ) was defined in Eq. (2.21).

To prove the above proposition, note that, for each set of %1, %2, ν, and θ, the eigenvalues λ of

Q(%1, %2, ν, θ) are the four roots of the characteristic equation

det [Q(%1, %2, ν, θ) − λI4] = 0 (5.74)

Hereafter In denote the n × n identity matrix. Let 02 denote the 2 × 2 zero matrix. Then according to

Eqs. (5.69), (5.72), and (5.73), the 2 × 2 submatrices in the top-left, bottom-left, and bottom-right corners

of the 4 × 4 matrix [Q(%1, %2, ν, θ) − λI4] are [M(%2, ν, θ) − λI2], 02, and [M(%1, ν, θ) − λI2], respectively. As

such, by using elementary theory of determinant, we have

det [Q(%1, %2, ν, θ) − λI4] ≡ det [M(%1, ν, θ) − λI2] · det [M(%2, ν, θ) − λI2] (5.75)

It now follows from Eqs. (5.74) and (5.75) that, for each set of %1, %2, and ν, the eigenvalues of Q(%1, %2, ν, θ)

are the four roots of the characteristic equations

det [M(%1, ν, θ) − λI2] = 0 (5.76)

and

det [M(%2, ν, θ) − λI2] = 0 ((5.77)

Because the roots of Eq. (5.76) and ((5.77) are λ±(%1, ν, θ) and λ±(%2, ν, θ), respectively, the proof is completed. QED.

Because of the relation Eq. (5.70), for a given set of %1, %2, and ν, the a-%-4 scheme is said to be stable

relative to this set if and only if, for all θ, −∞ < θ < +∞, all the elements of the matrix [Q(%1, %2, ν, θ)]m

remain bounded as m → +∞. By using the Jordan canonical form theorem [66] and the proposition we have

just proved, it can be shown that, for each given set of %1, %2, ν, and θ, all the elements of [Q(%1, %2, ν, θ)]m

remain bounded as m → +∞ if and only if

max {|λ+(%1, ν, θ)| , |λ−(%1, ν, θ)| , |λ+(%2, ν, θ)| , |λ−(%2, ν, θ)|} ≤ 1 if Q(%1, %2, ν, θ) is nondefective (5.78)

and

max {|λ+(%1, ν, θ)| , |λ−(%1, ν, θ)| , |λ+(%2, ν, θ)| , |λ−(%2, ν, θ)|} < 1 if Q(%1, %2, ν, θ) is defective (5.79)

Note that an n × n matrix is said to be defective if and only if the number of its linearly independent

eigenvectors is less than n [66].

30
From the above discussions, it is seen that, for a given set of %1, %2, and ν, a necessary condition for the

stability of a-%-4 scheme is

max −∞<θ<+∞ {|λ+(%1, ν, θ)| , |λ−(%1, ν, θ)| , |λ+(%2, ν, θ)| , |λ−(%2, ν, θ)|} ≤ 1 (5.80)

By using Eq. (2.21) and a line of arguments which was used in establishing the stability condition Eq. (3.14)

of [5] (Note: In [5], the a-% scheme is not defined for the case |ν| = 1), one can show that Eq. (5.80) ⇔

0 ≤ %1 ≤ 1, 0 ≤ %2 ≤ 1 and |ν| ≤ 1 (5.81)

Thus Eq. (5.81) represents a set of necessary stability conditions for the a-%-4 scheme.

Because Q(%1, %2, ν, θ) does become defective for some special values of %1 %2, ν and θ, according to

Eq. (5.79), Eq. (5.80) and therefore Eq. (5.81) generally is not sufficient to insure the stability of the a-%-4

scheme. However, it is expected that the necessary and sufficient stability conditions (not yet established)

for the a-%-4 scheme would be formed by Eq. (5.81) and some exclusion conditions involving some boundaries

of the stability domain defined in Eq. (5.81).

From the above discussions, one concludes that the 2nd-order and 3rd-order terms which appear on the

right sides of Eqs. (5.34) and (5.52) (which are associated with the 2×2 submatrices in the top-right corners

of the matrices Q+ and Q−) play only a marginal role in the stability of the a-%-4 scheme even though they

play a decisive role in the 4th-order accuracy of the scheme

Next we will present a set of non-reflecting boundary conditions with 4th-order accuracy.

5.5. Non-reflecting boundary conditions

Let (j, n) be a mesh point on the right spatial boundary. It has been shown numerically the current

extension of Eq. (2.44) i.e.,

un

j = un−1/2

j−1/2 , (ux ̄)

n

j = (ux ̄)

n−1/2

j−1/2 , (ux ̄x ̄)

n

j = (ux ̄x ̄)

n−1/2

j−1/2 , and (ux ̄x ̄x ̄)

n

j = (ux ̄x ̄x ̄)

n−1/2

j−1/2 (5.82)

and a similar extension of Eq. (2.45) also form sets of non-reflecting boundary conditions. However, these

simple extensions are not 4th order in accuracy. In the following, we present non-reflecting boundary

conditions with 4th-order accuracy.

Let (j, n) be a mesh point on the right spatial boundary. Again we assume

(ux ̄x ̄x ̄)

n

j = (ux ̄x ̄x ̄)

n−1/2

j−1/2 (5.83)

However, (ux ̄x ̄)n

j will be determined in terms of the mesh variables at (j − 1/2, n − 1/2) by imposing the

conservation condition Eq. (5.13) which links the mesh variables at the two diagonally opposite mesh points

(j, n) and (j − 1/2, n − 1/2). Because Eq. (5.13) ⇔

(1 + ν) [ux ̄x ̄ − (1 − ν)ux ̄x ̄x ̄]

n

j = (1 + ν) [ux ̄x ̄ + (1 − ν)ux ̄x ̄x ̄]

n−1/2

j−1/2 (5.84)

we assume

[ux ̄x ̄ − (1 − ν)ux ̄x ̄x ̄]

n

j = [ux ̄x ̄ + (1 − ν)ux ̄x ̄x ̄]

n−1/2

j−1/2 (5.85)

(Eq. (5.85) is a condition stronger than Eq. (5.84)). By combining Eqs. (5.83) and (5.85), one has

(ux ̄x ̄)

n

j = [ux ̄x ̄ + 2(1 − ν)ux ̄x ̄x ̄]

n−1/2

j−1/2 (5.86)

Next (ux ̄)n

j will be evaluated using the Taylor’s expansion

(ux)

n

j =

E

ux +

∆x

2

uxx +

∆t

2

uxt +

1

2

8 ∆x

2

92

uxxx +

∆x

2

∆t

2

uxxt +

1

2

6 ∆t

2

72

uxttFn−1/2

j−1/2

(5.87)

31
With the aid of Eqs. (2.9), (5.3), and (5.17), Eq. (5.87) can be rewritten as

(ux ̄)

n

j = $

ux ̄ + 2(1 − ν)ux ̄x ̄ + 2(1 − ν)

2ux ̄x ̄x ̄

%n−1/2

j−1/2 (5.88)

The mesh variable un

j will be evaluated in terms of the mesh variables at (j − 1/2, n − 1/2) by imposing

the conservation condition Eq. (2.8). Because Eq. (2.8) ⇔ Eq. (5.30), we will assume Eq. (5.30a) (a stronger

condition than Eq. (5.30)). As such one has

un

j =

/

(1 − ν)ux ̄ − 2(1 − ν + ν2)

3

ux ̄x ̄ + (1 − ν)(1 + ν2)

3

ux ̄x ̄x ̄

0n

j

+ (s−)

n−1/2

j−1/2 (5.89)

where (s−)

n−1/2

j−1/2 is defined in Eq. (5.28). With the aid of Eqs. (5.83), (5.86), and (5.88), un

j can now be

evaluated in terms of the mesh variables at (j − 1/2, n − 1/2).

A similar treatment can also be applied at the left spatial boundary. As will be shown, these new

non-reflecting boundary conditions are 4th-order in accuracy.

6. The c-τ-4 and c-τ ∗-4 schemes

The c-τ-4 and c-τ ∗-4 scheme will be built from the c-τ scheme in a procedure similar to that was used

to built the a-%-4 scheme from a-% scheme. Thus the basic equations including Eqs. (2.4), (2.9), (5.1)–(5.7),

(5.11), and (5.17) will be assumed here. In addition, both the conservation conditions Eq. (2.13) and (5.18)

will also be assumed here.

6.1. The c-τ-4 scheme

Because Eq. (5.18) is equivalent to Eq. (5.14) assuming Eqs. (5.7), (5.11), and (5.17), in the current

development, the mesh variable (ux ̄x ̄)n

j will be evaluated using Eq. (5.14). Moreover, because (i) Eq. (5.31)

is equivalent to Eq. (2.13) assuming Eqs. (2.4), (2.9), (5.4), and (5.17); and (ii) Eq. (5.34) can be obtained

from Eqs. (5.31) and (5.14), the mesh variable un

j will be evaluated using Eq. (5.34). As such, the current

schemes share with the a-%-4 scheme two component equations, i.e., Eqs. (5.14) and (5.34).

In the a-%-4 scheme, (ux ̄x ̄x ̄)n

j is evaluated using Eq. (5.26) which is the substitution form of Eq. (2.19a)

(by which (ux ̄)n

j is evaluated in the a-% scheme) except that the symbol % in Eq. (2.19a) is replaced by %1

in Eq. (5.26). Similarly, in the c-τ-4 scheme, (ux ̄x ̄x ̄)n

j is evaluated using the substitution form of Eq. (3.11)

(by which (ux ̄)n

j is evaluated in the c-τ scheme) with the understanding that the symbol τ in Eq. (3.11) is

replaced by τ1 in the substitution form. Thus we have

(ux ̄x ̄x ̄)

n

j = 1

2(1 + τ1)

&

[ux ̄x ̄ − (1 + 2ν − τ1)ux ̄x ̄x ̄]

n−1/2

j+1/2 − [ux ̄x ̄ + (1 − 2ν − τ1)ux ̄x ̄x ̄]

n−1/2

j−1/2

'

(j, n) ∈ Ω; τ1 "= −1

(6.1)

In the c-τ-4 scheme, (ux ̄)n

j will be evaluated using an equation which is the result of the following

assumptions: (i) (ux ̄)n

j = $

u4

x ̄(τ2)

%n

j (see Eq. (5.48)); (ii) the mesh variable (ux ̄x ̄x ̄)n

j which appears on the

right side of Eq. (5.48) is evaluated using Eq. (6.1); and (iii) the parameter τ which appears on the right

sides of Eqs. (5.44) and (5.45) is replaced by τ2. Then, by using Eqs. (5.44), (5.46), (5.48), and (6.1), one

32
has

(ux ̄)

n

j = 1

2

5 u

1 + τ2

− 1 − τ2 + 2ν

1 + τ2

ux ̄ +

1

2

* (1 − τ2 + 2ν)2

1 + τ2

− (1 + τ2)2

3(1 + τ1)

+

ux ̄x ̄

+

1

6

* (1 + τ2)2(1 − τ1 + 2ν)

1 + τ1

− (1 − τ2 + 2ν)3

1 + τ2

+

ux ̄x ̄x ̄

<n−1/2

j+1/2

− 1

2

5 u

1 + τ2

+

1 − τ2 − 2ν

1 + τ2

ux ̄ +

1

2

* (1 − τ2 − 2ν)2

1 + τ2

− (1 + τ2)2

3(1 + τ1)

+

ux ̄x ̄

− 1

6

* (1 + τ2)2(1 − τ1 − 2ν)

1 + τ1

− (1 − τ2 − 2ν)3

1 + τ2

+

ux ̄x ̄x ̄

<n−1/2

j−1/2

(j, n) ∈ Ω; τ1 "= −1; τ2 "= −1

(6.2)

By definition, the c-τ-4 scheme is formed by Eqs. (5.14), (5.34), (6.1), and (6.2). It enforces the local

conservation relations Eqs. (2.13) and (5.18). According to an argument given in Comment (d) presented at

the end of Sec. 2.1, it will also enforce a global conservation relation for each of the fluxes of "h∗ and "

h∗

xx over

any space-time region that is the union of any combination of CE(j, n), (j, n) ∈ Ω. Moreover, the scheme is

completely explicit, i.e., the four independent mesh variables at each (j, n) can be evaluated directly in terms

of those at the (n − 1/2)th time level without resorting to solving a system of linear equations involving the

independent mesh variables at (j, n) such as Eqs. (2.7a) and (2.8a).

The c-τ and a schemes become the same scheme when ν = τ = 0 (see Eq. (3.12)). It happens that

the c-τ-4 and a-4 schemes also become the same scheme when ν = τ1 = τ2 = 0. In fact, for the case

ν = τ1 = τ2 = 0, Eqs. (6.1) and (6.2) reduce to

(ux ̄x ̄x ̄)

n

j = 1

2

*

(ux ̄x ̄ − ux ̄x ̄x ̄)

n−1/2

j+1/2 − (ux ̄x ̄ + ux ̄x ̄x ̄)

n−1/2

j−1/2

+

(6.3)

and

(ux ̄)

n

j = 1

2

E6

u − ux ̄ +

1

3

ux ̄x ̄

7n−1/2

j+1/2

−

6

u + ux ̄ +

1

3

ux ̄x ̄

7n−1/2

j−1/2

F

(6.4)

respectively. It can be shown that Eqs. (5.26) and (5.52) (the counterparts of Eqs. (6.1) and (6.2) in the

a-%-4 scheme) also reduce to Eqs. (6.3) and (6.4), respectively, when ν = %1 = %2 = 0. Because (i) the a-4

scheme is a special case of a-%-4 scheme with %1 = %2 = 0 and (ii) Eqs. (5.14) and (5.34) are shared by the

c-τ-4 and a-%-4 schemes, the proof is completed. QED.

Because the conservation relations Eqs. (2.7), (2.8), (5.12), and (5.13) are satisfied by the a-4 scheme for

all ν, it now follows that these conservation relations are satisfied by the c-τ-4 scheme when ν = τ1 = τ2 = 0.

The c-τ and c schemes (the special case of the a-% scheme with % = 1/2) become the same scheme when

τ = 1 (see Eq. 3.12). Let the c-4 scheme be the special case of the a-%-4 scheme with %1 = %2 = 1/2. Then

one can also show that the c-τ-4 and c-4 schemes become the same scheme when τ1 = τ2 = 1. In fact, for

the case τ1 = τ2 = 1, Eqs. (6.1) and (6.2) reduce to

(ux ̄x ̄x ̄)

n

j = 1

4

*

(ux ̄x ̄ − 2ν ux ̄x ̄x ̄)

n−1/2

j+1/2 − (ux ̄x ̄ − 2ν ux ̄x ̄x ̄)

n−1/2

j−1/2

+

(6.5)

and

(ux ̄)

n

j = 1

12

&$

3u − 6ν ux ̄ + 2(3ν2 − 1)ux ̄x ̄ + 4ν(1 − ν2)ux ̄x ̄x ̄

%n−1/2

j+1/2

− $3u − 6ν ux ̄ + 2(3ν2 − 1)ux ̄x ̄ + 4ν(1 − ν2)ux ̄x ̄x ̄

%n−1/2

j−1/2

'

(6.6)

33
respectively. It can be shown that Eqs. (5.26) and (5.52) also reduce to Eqs. (6.5) and (6.6), respectively,

when %1 = %2 = 1/2 (i.e., when the a-%-4 scheme becomes the c-4 scheme). QED.

6.2. Matrix form of the c-τ-4 scheme

The component equations of the c-τ-4 scheme, i.e., Eqs. (5.14), (5.34), (6.1) and (6.2), can be cast into

the matrix equation

"q(j, n) = QG−(τ1, τ2, ν) "q(j + 1/2, n − 1/2) + QG+(τ1, τ2, ν) "q(j − 1/2, n − 1/2) (j, n) ∈ Ω (6.7)

where "q(j, n) is defined in Eq. (5.54), and

QG±(τ1, τ2, ν) def

=





qˆ

±

11 qˆ

±

12 qˆ

±

13 qˆ

±

14

qˆ

±

21 qˆ

±

22 qˆ

±

23 qˆ

±

24

0 0 ˆq±

33 qˆ

±

34

0 0 ˆq±

43 qˆ

±

44





(6.8)

with

qˆ

±

11

def

= 1 ± ν

2 , qˆ

±

12

def

= ±(1 − ν2)

2 , qˆ

±

13

def

= ∓ν(1 − ν2)

3 , qˆ

±

14

def

= ∓(1 − ν2)2

6 (6.9)

qˆ

±

21

def

= ∓1

2(1 + τ2)

, qˆ

±

22

def

= −(1 − τ2 ∓ 2ν)

2(1 + τ2) , qˆ

±

23

def

= ± (1 + τ2)2

12(1 + τ1) ∓ (1 − τ2 ∓ 2ν)24(1 + τ2)

qˆ

±

24

def

= (1 + τ2)2(1 − τ1 ∓ 2ν)

12(1 + τ1) − (1 − τ2 ∓ 2ν)312(1 + τ2)

(6.10)

and

qˆ

±

33

def

= 1 ± ν

2 , qˆ

±

34

def

= ±(1 − ν2)

2 , qˆ

±

43

def

= ∓1

2(1 + τ1)

, qˆ

±

44

def

= −(1 − τ1 ∓ 2ν)

2(1 + τ1) (6.11)

To further explore the relations between the c-τ-4 and c-τ schemes, note that the component equations

of the c-τ scheme, i.e., Eqs. (2.10) and (3.11), can be cast into the matrix form

"q(j, n) = M

H−(τ, ν) "q(j + 1/2, n − 1/2) + MH+(τ, ν) "q(j − 1/2, n − 1/2) (j, n) ∈ Ω (6.12)

where "q(j, n) is defined in Eq. (5.61), and

M

H±(τ, ν) def

= 1

2





1 ± ν ±(1 − ν2)

∓1

1 + τ −1 − τ ∓ 2ν1 + τ





− ∞ < τ < +∞ (6.13)

With the aid of Eqs. (6.9)–(6.11) and (6.13), it can be shown that





qˆ

±

11 qˆ

±

12

qˆ

±

21 qˆ

±

22



 = M

H±(τ2, ν) (6.14)

and





qˆ

±

33 qˆ

±

34

qˆ

±

43 qˆ

±

44



 = M

H±(τ1, ν) (6.15)

34
By comparing Eqs. (6.14) and (6.15) with Eq. (6.8), one now arrives at the conclusion that the 2 × 2

submatrices in the top-left and bottom right corners of the 4 × 4 matrix QG±(τ1, τ2, ν) are M

H±(τ2, ν) and

M

H±(τ1, ν), respectively.

6.3. Stability conditions of the c-τ-4 scheme

By using similar arguments presented in Sec. 5.4, one concludes that the amplification matrices per

marching step of the c-τ-4 and c-τ schemes are

QG(τ1, τ2, ν, θ) = eiθ/2QG−(τ1, τ2, ν) + e−iθ/2QG+(τ1, τ2, ν)

=





eiθ/2qˆ

−

11 + e−iθ/2qˆ

+

11 eiθ/2qˆ

−

12 + e−iθ/2qˆ

+

12 eiθ/2qˆ

−

13 + e−iθ/2qˆ

+

13 eiθ/2qˆ

−

14 + e−iθ/2qˆ

+

14

eiθ/2qˆ

−

21 + e−iθ/2qˆ

+

21 eiθ/2qˆ

−

22 + e−iθ/2qˆ

+

22 eiθ/2qˆ

−

23 + e−iθ/2qˆ

+

23 eiθ/2qˆ

−

24 + e−iθ/2qˆ

+

24

0 0 eiθ/2qˆ

−

33 + e−iθ/2qˆ

+

33 eiθ/2qˆ

−

34 + e−iθ/2qˆ

+

34

0 0 eiθ/2qˆ

−

43 + e−iθ/2qˆ

+

43 eiθ/2qˆ

−

44 + e−iθ/2qˆ

+

44





(6.16)

and

M

H(τ, ν, θ) def

= eiθ/2M

H−(τ, ν) + e−iθ/2MH+(τ, ν) ≡

=





cos(θ/2) − iν sin(θ/2) −i(1 − ν2) sin(θ/2)

isin(θ/2)

1 + τ −

/

(1 − τ) cos(θ/2) + 2iν sin(θ/2)

1 + τ

0





(6.17)

respectively. Moreover, by using Eqs. (6.14)–(6.17), one concludes that





eiθ/2qˆ

−

11 + e−iθ/2qˆ

+

11 eiθ/2qˆ

−

12 + e−iθ/2qˆ

+

12

eiθ/2qˆ

−

21 + e−iθ/2qˆ

+

21 eiθ/2qˆ

−

22 + e−iθ/2qˆ

+

22



 = M

H(τ2, ν, θ) (6.18)

and





eiθ/2qˆ

−

33 + e−iθ/2qˆ

+

33 eiθ/2qˆ

−

34 + e−iθ/2qˆ

+

34

eiθ/2qˆ

−

43 + e−iθ/2qˆ

+

43 eiθ/2qˆ

−

44 + e−iθ/2qˆ

+

44



 = M

H(τ1, ν, θ) (6.19)

i.e., the 2×2 submatrices in the top-left and bottom-right corners of the matrix QG(τ1, τ2, ν, θ) are M

H(τ2, ν, θ)

and M

H(τ1, ν, θ), respectively. As such, the amplification matrix of the c-τ scheme is identical to each of the

above 2 × 2 submatrices of QG(τ1, τ2, ν, θ) except that the parameter τ in the former matrix is replaced by τ2

and τ1 in the latter two submatrices. By using this property and using a line of arguments similar to that

presented in Sec. 5.4, one can show that, for each set of τ1, τ2, ν, and θ, the eigenvalues of QG(τ1, τ2, ν, θ) are

the union of those of M

H(τ1, ν, θ) and M

H(τ2, ν, θ). As a result, by using a line of arguments presented in the

last part of Sec. 5.4, one can infer from the stability conditions Eqs. (3.13) of the c-τ scheme that

ν2 ≤ 1, τ1 ≥ τo(ν2), τ2 ≥ τo(ν2), (ν2, τ1) "= (1, 1), and (ν2, τ2) "= (1, 1) (6.20)

form a set of necessary stability conditions for the c-τ-4 scheme.

6.4. The c-τ ∗-4 schemes

Let h1(s) and h2(s) be any two functions that meet the same requirements imposed on the function h(s)

introduced in Sec. 3.2. Then, because of the close relations between the c-τ and c-τ-4 schemes established

above, here a c-τ ∗-4 scheme is defined to be a special c-τ-4 scheme with

τ1 = h1(ν2) and τ2 = h2(ν2) (|ν| < 1) (6.21)

35
According to Eq. (3.22), one can assume h1(s) = h2(s) = √s. As such a special c-τ ∗-4 scheme can be defined

by the relation

τ1 = τ2 = |ν| (|ν| < 1) (6.22)

For robust real-world applications, we may also introduce a special c-τ-4 scheme which is defined by an

extension of Eq. (3.24), i.e.,

τ1 = β1 |ν| and τ2 = β2 |ν| (β1 ≥ 1; β2 ≥ 1; |ν| < 1) (6.23)

where β1 ≥ 1 and β2 ≥ 1 are two adjustable parameters.

7. Fourth-order solvers for the inviscid Burgers equation

Consider the inviscid Burgers equation

∂u

∂t +

∂f(u)

∂x = 0 (7.1)

where

f(u) def

= u2

2 (7.2)

Because, in E2, Eq. (7.1) can be expressed as ∇ · "

h = 0 with

"

h def = (f, u) (7.3)

Gauss’ divergence theorem in the space-time E2 implies that Eq. (7.1) is the differential form of the integral

conservation law Eq. (2.2) if "h is defined by Eq. (7.3). Moreover, for a smooth function u(x, t), Eq. (2.2)

implies the conservation condition Eq. (5.9) if

"

hxx

def

= ∂2"

h

∂x2 ≡

6∂2f

∂x2 ,

∂2u

∂x2

7

(7.4)

7.1. The inviscid Burgers c-τ-4 scheme

In this section, the c-τ-4 scheme will be extended to become a solver for the inviscid Burgers equation. Here again Eq. (5.1) is assumed. However, for any (x, t) ∈ SE(j, n), the flux function f(u) will be

approximated by

f ∗(x, t, j, n) def

= f n

j + (fx)

n

j (x − xj )+(ft)

n

j (t − t

n) + 1

2

(fxx)

n

j (x − xj )

2

+ (fxt)

n

j (x − xj )(t − t

n) + 1

2

(ftt)

n

j (t − t

n)2 +

1

6

(fxxx)

n

j (x − xj )

3

+

1

2

(fxxt)

n

j (x − xj )

2(t − tn) + 1

2

(fxtt)

n

j (x − xj )(t − t

n)2 +

1

6

(fttt)

n

j (t − t

n)3

(x, t) ∈ SE(j, n) (7.5)

Here f n

j , (fx)nj , (ft)nj , (fxx)nj , (fxt)nj , (ftt)nj (fxxx)nj , (fxxt)nj , (fxtt)nj , and (fttt)nj are constants in SE(j, n).

They can be considered as the numerical analogues of the values of f, ∂f /∂x, ∂f /∂t, ∂2f /∂x2, ∂2f /∂x∂t,

∂2f /∂t

2, ∂3f /∂x3, ∂3f /∂x2∂t, ∂3f /∂x∂t2, and ∂3f /∂t3 at the mesh point (j, n), respectively. Thus

f ∗(x, t, j, n) represents a 3rd-order Taylor’s approximation of f. Moreover, as will be shown, by imposing proper assumptions, the coefficients on the right sides of Eqs. (5.1) and (7.5) can be explicity expressed

as the functions of the independent unknowns un

j , (ux)nj , (uxx)nj and (uxxx)nj .

For any (j, n) ∈ Ω, let u = u∗(x, t; j, n) and f(u) = f ∗(x, t, j, n) satisfy Eq. (7.1) uniformly in SE(j, n),

i.e.,

∂u∗(x, t; j, n)

∂t +

∂f ∗(x, t, j, n)

∂x ≡ 0 (j, n) ∈ Ω (7.6)

36
It can be shown that Eq. (7.6) ⇔, for all (j, n) ∈ Ω,

(ut)

n

j = −(fx)

n

j (7.7)

(uxt)

n

j = −(fxx)

n

j (7.8)

(utt)

n

j = −(fxt)

n

j (7.9)

(uxxt)

n

j = −(fxxx)

n

j (7.10)

(uxtt)

n

j = −(fxxt)

n

j (7.11)

and

(uttt)

n

j = −(fxtt)

n

j (7.12)

Moreover, by taking the numerical analogues of all the differential forms of Eq. (7.2) up to the third order,

one has, for all (j, n) ∈ Ω,

f n

j = 1

2

(u2)

n

j (7.13)

(fx)

n

j = (uux)

n

j (7.14)

(ft)

n

j = (uut)

n

j (7.15)

(fxx)

n

j = $

(uuxx + (ux)

2%n

j (7.16)

(fxt)

n

j = (uuxt + uxut)

n

j (7.17)

(ftt)

n

j = $

(uutt + (ut)

2%n

j (7.18)

(fxxx)

n

j = (uuxxx + 3uxuxx)

n

j (7.19)

(fxxt)

n

j = (uuxxt + utuxx + 2ux uxt)

n

j (7.20)

(fxtt)

n

j = (uuxtt + uxutt + 2utuxt)

n

j (7.21)

(fttt)

n

j = (uuttt + 3ututt)

n

j (7.22)

By substituting Eqs. (7.14), (7.16), and (7.19) into Eqs. (7.7), (7.8), and (7.10), respectively, one has

(ut)

n

j = −(uux)

n

j (7.23)

(uxt)

n

j = − $

uuxx + (ux)

2%n

j (7.24)

and

(uxxt)

n

j = −(uuxxx + 3uxuxx)

n

j (7.25)

respectively. In turn, by substituting Eqs. (7.23)–(7.25) into Eqs. (7.15), (7.17), and (7.20), one has

(ft)

n

j = −(u2ux)

n

j (7.26)

(fxt)

n

j = − $

2u(ux)

2 + u2uxx%n

j (7.27)

and

(fxxt)

n

j = − $

2(ux)

3 + 6uuxuxx + u2uxxx%n

j (7.28)

Moreover, by substituting Eqs. (7.27) and (7.28) into Eqs. (7.9) and (7.11), respectively, one has

(utt)

n

j = $

2u(ux)

2 + u2uxx%n

j (7.29)

37
and

(uxtt)

n

j = $

2(ux)

3 + 6uuxuxx + u2uxxx%n

j (7.30)

respectively.

Next, by substituting Eqs. (7.23), (7.24), (7.29), and (7.30) into Eqs. (7.18) and (7.21), one obtains

(ftt)

n

j = $

3u2(ux)

2 + u3uxx%n

j (7.31)

and

(fxtt)

n

j = $

6u(ux)

3 + 9u2uxuxx + u3uxxx%n

j (7.32)

In turn, by substituting Eq. (7.32) into Eq. (7.12), one has

(uttt)

n

j = − $

6u(ux)

3 + 9u2uxuxx + u3uxxx%n

j (7.33)

Finally, by substituting Eqs. (7.23), (7.29), and (7.33) into Eq. (7.22), one has

(fttt)

n

j = − $

12u2(ux)

3 + 12u3uxuxx + u4uxxx%n

j (7.34)

According to the set of sixteen equations formed by Eqs. (7.13), (7.14), (7.16), (7.19), and (7.23)–

(7.34), all the twenty coefficients on the right sides of Eqs. (5.1) and (7.5) other than un

j , (ux)nj , (uxx)nj , and

(uxxx)n

j are explicit functions of the latter coefficients. Thus there are only four independent mesh variables

associated with each (j, n) ∈ Ω. In fact, by substituting these sixteen equations into Eqs. (5.1) and (7.5),

one has

u∗(x, t; j, n) = un

j + (ux)

n

j

&

1 − (ux)

n

j (t − t

n) + $

(ux)

n

j

%2

(t − t

n)2

'$

x − xj − un

j (t − t

n)

%

+ (uxx)n

j

2

$

1 − 3(ux)

n

j (t − t

n)

%$x − xj − un

j (t − t

n)

%2

+ (uxxx)n

j

6

$

x − xj − un

j (t − t

n)

%3

(7.35)

and

f ∗(x, t, j, n) =(un

j )2

2 + (ux)

n

j

&

un

j + (ux)n

j

2

$

x − xj − 3un

j (t − t

n)

%

− $(ux)

n

j

%2

(t − t

n)

$

x − xj − 2un

j (t − t

n)

%'$

x − xj − un

j (t − t

n)

%

+ (uxx)n

j

2

&

un

j + (ux)

n

j

$

x − xj − 4un

j (t − t

n)

%'$

x − xj − un

j (t − t

n)

%2

+

un

j (uxxx)nj

6

$

x − xj − un

j (t − t

n)

%3

(7.36)

In turn, by substituting Eqs. (7.35) and (7.36) into Eq. (5.5) and

f ∗

xx(x, t; j, n) def

= ∂2f ∗(x, t, j, n)

∂x2 (7.37)

one has

u∗

xx(x, t; j, n) = (uxx)

n

j + (uxxx)

n

j (x − xj ) − (uuxxx + 3uxuxx)

n

j (t − t

n) (7.38)

and

f ∗

xx(x, t; j, n) = [uuxx + (ux)

2]n

j + (uuxxx + 3uxuxx)

n

j (x−xj )−$

2(ux)

3 + 6uuxuxx + u2uxxx%n

j (t−t

n) (7.39)

38
Note that, for the same reason presented immediately following Eq. (2.6), except for flux evaluation,

Eqs. (7.35) and (7.38) are applicable even for (x, t) ∈/ SE(j, n).

Next let the numerical analogues of Eqs. (7.3) and (7.4) be

"

h∗(x, t; j, n) def = (f ∗(x, t, j, n), u∗(x, t; j, n)) (x, t) ∈ SE(j, n) (7.40)

and "

h∗

xx(x, t; j, n) def = (f ∗xx(x, t; j, n), u∗xx(x, t; j, n)) (x, t) ∈ SE(j, n) (7.41)

respectively. Also, let

(ux ̄)

n

j

def

= ∆x

4 (ux)

n

j , (ux ̄x ̄)

n

j

def

=

8 ∆x

4

92

(uxx)

n

j , and (ux ̄x ̄x ̄)

n

j

def

=

8 ∆x

4

93

(uxxx)

n

j (j, n) ∈ Ω (7.42)

and

νn

j

def

= un

j ∆t

∆x , (νx ̄)

n

j

def

= (ux ̄)n

j ∆t

∆x , (νx ̄x ̄)

n

j

def

= (ux ̄x ̄)n

j ∆t

∆x , and (νx ̄x ̄x ̄)

n

j

def

= (ux ̄x ̄x ̄)n

j ∆t

∆x (j, n) ∈ Ω (7.43)

Obviously, (ux ̄)n

j , (ux ̄x ̄)nj , and (ux ̄x ̄x ̄)nj defined here are identical to those defined in Eqs. (2.9) and (5.17).

However, the “Courant number” νn

j and its “derivatives” (νx ̄)nj , (νx ̄x ̄)nj , and (νx ̄x ̄x ̄)nj are mesh-point dependent in the current development.

Like the c-τ-4 scheme, the current extension will be construct to enforce the least stringent conservation

conditions Eqs. (2.13) and (5.18). In fact, by using Eqs. (7.38), (7.39), (7.41), (7.42), and (7.43), it can be

shown that Eq. (5.18) ⇔

(ux ̄x ̄)

n

j = 1

2

&$

νx ̄(2νx ̄ − 1)ux ̄ + [1 − ν(1 − 6νx ̄)]ux ̄x ̄ − (1 − ν2)ux ̄x ̄x ̄

%n−1/2

j+1/2

+ $νx ̄(1 − 2νx ̄)ux ̄ + [1 + ν(1 − 6νx ̄)]ux ̄x ̄ + (1 − ν2)ux ̄x ̄x ̄

%n−1/2

j−1/2

'

(7.44)

Moreover, by using Eqs. (7.35), (7.36), (7.40), (7.42), and (7.43), it can be shown that Eq. (2.13) ⇔

2

*

u +

2ux ̄x ̄

3

+n

j

=

&"

1 − ν

2

#

u − $1 − ν2 + 2νx ̄ν2 − 4(νx ̄)

2ν2%

ux ̄ + $2(1 − ν3)

3 + 4νx ̄ν3%

ux ̄x ̄ − (1 − ν4)

3

ux ̄x ̄x ̄

'n−1/2

j+1/2

+

&"

1 + ν

2

#

u + $1 − ν2 + 2νx ̄ν2 − 4(νx ̄)

2ν2%

ux ̄ + $ 2(1 + ν3)

3 − 4νx ̄ν3%

ux ̄x ̄ + (1 − ν4)

3

ux ̄x ̄x ̄

'n−1/2

j−1/2

(7.45)

Expressing the mesh variable (ux ̄x ̄)n

j on the left side of Eq. (7.45) in terms of the mesh variables at the

(n − 1/2)th time level by using Eq. (7.44), one obtains

un

j =

&(2 − ν)

4 u −

$

3(1 − ν2) + 2(3ν2 − 1)(1 − 2νx ̄)νx ̄]

6

ux ̄ + ν(1 − ν2)(1 − 6νx ̄)

3

ux ̄x ̄ + (1 − ν2)2

6

ux ̄x ̄x ̄

'n−1/2

j+1/2

+

&(2 + ν)

4 u +

$

3(1 − ν2) + 2(3ν2 − 1)(1 − 2νx ̄)νx ̄]

6

ux ̄ − ν(1 − ν2)(1 − 6νx ̄)

3

ux ̄x ̄ − (1 − ν2)2

6

ux ̄x ̄x ̄

'n−1/2

j−1/2

(7.46)

39
It has been shown that, by imposing the conservation conditions Eqs. (2.13) and (5.18), one can obtain

Eqs. (7.44) and (7.46) by which both un

j and (ux ̄x ̄)nj can be explicitly evaluated in terms of the mesh

variables at the (n − 1/2)th time level. In the following, by using techniques developed earlier, similar

explicit expressions will be developed for (ux ̄)n

j and (ux ̄x ̄x ̄)nj .

To proceed, let P + and P − be the points depicted in Fig. 3, except that the parameter τ shown in the

figure be replaced by the mesh-point dependent parameter (τ1)n

j . Then the x-coordinates of P + and P − are

x(P +) = xj + $

1+(τ1)

n

j

% ∆x

4 = xj+1/2 − $

1 − (τ1)

n

j

% ∆x

4 , −∞ < (τ1)

n

j < +∞ (7.47)

and

x(P −) = xj − $

1+(τ1)

n

j

% ∆x

4 = xj−1/2 + $

1 − (τ1)

n

j

% ∆x

4 , −∞ < (τ1)

n

j < +∞ (7.48)

respectively. In the current development, we assume

(uxxx)

n

j = uxx(P +) − uxx(P −)

$

1+(τ1)n

j

%

∆x/2 (7.49)

where

uxx(P +) def

= u∗

xx "

x(P +), tn; j + 1/2, n − 1/2

# (7.50)

and

uxx(P −) def

= u∗

xx "

x(P −), tn; j − 1/2, n − 1/2

# (7.51)

Because |P −P +| = $

1+(τ1)n

j

%

∆x/2, the expression on the right side of Eq. (7.49) represents a centraldifference approximation of (uxxx)n

j . With the aid of Eqs. (7.38), (7.42), (7.43), (7.47), (7.48), (7.50), and

(7.51), it can be shown that Eq. (7.49) ⇔

(ux ̄x ̄x ̄)

n

j =

$

(1 − 6νx ̄)ux ̄x ̄ − (1 − τ1 + 2ν)ux ̄x ̄x ̄

%n−1/2

j+1/2 − $

(1 − 6νx ̄)ux ̄x ̄ + (1 − τ1 − 2ν)ux ̄x ̄x ̄

%n−1/2

j−1/2

2

$

1+(τ1)n

j

% (7.52)

To simplify notation, in Eq. (7.52) and hereafter in this section we adopt a convention that can be explained

using two expression on the right side of Eq. (7.52) as examples, i.e.,

$

(1 − 6νx ̄)ux ̄x ̄ − (1 − τ1 + 2ν)ux ̄x ̄x ̄

%n−1/2

j+1/2 = $

(1 − 6νx ̄)ux ̄x ̄

%n−1/2

j+1/2 − $

1 − (τ1)

n

j + 2νn−1/2

j+1/2

%

(ux ̄x ̄x ̄)

n−1/2

j+1/2

and

$

(1 − 6νx ̄)ux ̄x ̄ + (1 − τ1 − 2ν)ux ̄x ̄x ̄

%n−1/2

j−1/2 = $

(1 − 6νx ̄)ux ̄x ̄

%n−1/2

j−1/2 + $

1 − (τ1)

n

j − 2νn−1/2

j−1/2

%

(ux ̄x ̄x ̄)

n−1/2

j−1/2

In other words, the parameter τ1 which appears on the right side of an equation such as Eq. (7.52)

and the mesh variable which appears on the left side are associated with the same mesh point

(j, n).

Next, let P + and P − be the points depicted in Fig. 3, except that the parameter τ shown in the figure

be replaced by the mesh-point dependent parameter (τ2)n

j . Then the x-coordinates of P + and P − are

x(P +) = xj + $

1+(τ2)

n

j

% ∆x

4 = xj+1/2 − $

1 − (τ2)

n

j

% ∆x

4 , −∞ < (τ2)

n

j < +∞ (7.53)

and

x(P −) = xj − $

1+(τ2)

n

j

% ∆x

4 = xj−1/2 + $

1 − (τ2)

n

j

% ∆x

4 , −∞ < (τ2)

n

j < +∞ (7.54)

40
respectively. In the current development, (ux)n

j is approximated with the current version of the expression

on the right side of Eq. (5.46), i.e.,

(ux)

n

j = u4∗(P +) − u4∗(P −)

$

1+(τ2)n

j

%

∆x/2 −

=$1+(τ2)n

j

%

∆x

>2

96 (uxxx)

n

j (7.55)

where

u4∗(P ±) def

= u∗ "

x(P ±), tn; j ± 1/2, n − 1/2

# (7.56)

According to Eq. (7.35), u4∗(P +) is the 3rd-order Taylor’s approximation (i.e., the error is O $

(∆x)4%

) of u

at P + evaluated using the mesh variables at (j + 1/2, n − 1/2) while u4∗(P −) is that at P − evaluated using

the mesh variables at (j − 1/2, n − 1/2). As such, Eq. (5.43) implies that Eq. (7.55) is an approximation of

3rd-order accuracy for the value of ∂u/∂x at the mesh point (j, n).

Next, with the aid Eqs. (7.35), (7.42), (7.43), (7.53), and (7.54), it can be shown that

u4∗(P ±) =&

u ∓ (1 − τ2 ± 2ν)[1 − 2νx ̄ + 4(νx ̄)

2]ux ̄ +

1

2

(1 − τ2 ± 2ν)

2(1 − 6νx ̄)ux ̄x ̄

∓

1

6

(1 − τ2 ± 2ν)

3ux ̄x ̄x ̄

'n−1/2

j±1/2

(7.57)

In turn, by using Eqs. (7.42),(7.52) and (7.57), it can be shown that Eq. (7.55) ⇔

(ux ̄)

n

j =

5 u

2(1 + τ2) − (1 − τ2 + 2ν)

$

1 − 2νx ̄ + 4(νx ̄)2

%

ux ̄

2(1 + τ2) + (1 − 6νx ̄)

* (1 − τ2 + 2ν)2

4(1 + τ2) − (1 + τ2)212(1 + τ1)

+

ux ̄x ̄

+

* (1 + τ2)2(1 − τ1 + 2ν)

12(1 + τ1) − (1 − τ2 + 2ν)312(1 + τ2)

+

ux ̄x ̄x ̄

<n−1/2

j+1/2

−

5 u

2(1 + τ2) + (1 − τ2 − 2ν)

$

1 − 2νx ̄ + 4(νx ̄)2

%

ux ̄

2(1 + τ2) + (1 − 6νx ̄)

* (1 − τ2 − 2ν)2

4(1 + τ2) − (1 + τ2)212(1 + τ1)

+

ux ̄x ̄

−

* (1 + τ2)2(1 − τ1 − 2ν)

12(1 + τ1) − (1 − τ2 − 2ν)312(1 + τ2)

+

ux ̄x ̄x ̄

<n−1/2

j−1/2

(7.58)

Note that the simplifying convention established earlier for the parameter τ1 in Eq. (7.52)

applies here for both τ1 and τ2, i.e., everywhere on the right sides of Eqs. (7.57) and (7.58),

the symbols τ1 and τ2 are the shorthands for (τ1)n

j and (τ2)nj , respectively.

By definition, the Burgers c-τ-4 scheme is formed by Eqs. (7.44), (7.46), (7.52), and (7.58). It enforces

the local conservation relations Eqs. (2.13) and (5.18). According to an argument given in Comment (d)

presented at the end of Sec. 2.1, it will also enforce a global conservation relation for each of the fluxes

of "h∗ and "

h∗

xx over any space-time region that is the union of any combination of CE(j, n), (j, n) ∈ Ω.

Moreover, even though it is a solver for a nonlinear equation Eq. (7.1), by explicitly specifying (τ1)n

j and

(τ2)n

j as functions of the mesh variables at (j − 1/2, n − 1/2) and (j + 1/2, n − 1/2), the scheme will become

completely explicit, i.e., the four independent mesh variables at each (j, n) can be evaluated in terms of

those at the (n − 1/2)th time level without resorting to solving a system of nonlinear equations involving

the independent mesh variables at (j, n). Also note that, because νn

j = (νx ̄)nj = 0 when ∆t = 0, Eqs. (7.52)

and (7.58) reduce to Eqs. (6.3) and (6.4) respectively when ∆t = τ1 = τ2 = 0 In fact, it can be shown that

the Burgers c-τ-4 scheme satisfies Eqs. (2.7), (2.8), (5.12), and (5.13) (i.e., the conservation relations over

each CE±(j, n)) when ∆t = τ1 = τ2 = 0.

41
To turn the current scheme into an explicit scheme, as an example, we can assume that

(τ1)

n

j = β1|νnj | and (τ2)

n

j = β2|νnj | (β1 ≥ 1; β2 ≥ 1) (7.59)

where β1 ≥ 1 and β2 ≥ 1 are two adjustable constant parameters. According to Eqs. (7.43) and (7.46),

Eq. (7.59) effectively defines (τ1)n

j and (τ2)nj in terms of the mesh variables at (j + 1/2, n − 1/2) and

(j − 1/2, n − 1/2). As another example, one can assume

(τ1)

n

j = β1 max&

|νn−1/2

j−1/2 |, |νn−1/2j+1/2 |

'

and (τ2)

n

j = β2 max&

|νn−1/2

j−1/2 |, |νn−1/2j+1/2 |

'

(β1 ≥ 1; β2 ≥ 1) (7.60)

where β1 ≥ 1 and β2 ≥ 1 are again two adjustable constant parameters. Other possible choices of (τ1)n

j and

(τ2)n

j are specified in Eqs. (7.108)–(7.110).

In the following, it will be shown that the stability conditions of the Burgers c-τ-4 scheme are essentially

similar to those of the c-τ-4 scheme.

7.2. Stability conditions of the inviscid Burgers c-τ-4 scheme

As explained in Sec. 5.4, for linear solvers such as the a-%-4 or c-τ-4 scheme, propagation of roundoff errors is governed by the same linear system of equations governing the associated independent mesh

variables. As will be shown, the above simple relation is not valid for a nonlinear solver such as the Burgers

c-τ-4 scheme.

To proceed, note that Eqs. (7.1) and (7.2) imply that the dependent variable u has the dimension of

∆t/∆x. Thus νn

j , (νx ̄)nj , (νx ̄x ̄)nj , and (νx ̄x ̄x ̄)nj defined in Eq. (7.43) are all dimensionless. As a result, by

multiplying both sides of each of Eqs. (7.44), (7.46), (7.52), and (7.58) with the factor ∆t/∆x, one can turn

the latter equations into their dimensionless versions. The dimensionless versions can be obtained from

the original versions by replacing the symbols u, ux ̄, ux ̄x ̄ and ux ̄x ̄x ̄ with ν, νx ̄, νx ̄x ̄, and νx ̄x ̄x ̄, respectively.

Explicitly, the dimensionless versions of Eqs. (7.46), (7.58), (7.44), and (7.52) are

νn

j =

&(2 − ν)

4 ν −

$

3(1 − ν2) + 2(3ν2 − 1)(1 − 2νx ̄)νx ̄]

6 νx ̄ + ν(1 − ν2)(1 − 6νx ̄)3 νx ̄x ̄ + (1 − ν2)26 νx ̄x ̄x ̄

'n−1/2

j+1/2

+

&(2 + ν)

4 ν +

$

3(1 − ν2) + 2(3ν2 − 1)(1 − 2νx ̄)νx ̄]

6 νx ̄ − ν(1 − ν2)(1 − 6νx ̄)3 νx ̄x ̄ − (1 − ν2)26 νx ̄x ̄x ̄

'n−1/2

j−1/2

(7.61)

(νx ̄)

n

j =

5 ν

2(1 + τ2) − (1 − τ2 + 2ν)

$

1 − 2νx ̄ + 4(νx ̄)2

%

νx ̄

2(1 + τ2) + (1 − 6νx ̄)

* (1 − τ2 + 2ν)2

4(1 + τ2) − (1 + τ2)212(1 + τ1)

+

νx ̄x ̄

+

* (1 + τ2)2(1 − τ1 + 2ν)

12(1 + τ1) − (1 − τ2 + 2ν)312(1 + τ2)

+

νx ̄x ̄x ̄

<n−1/2

j+1/2

−

5 ν

2(1 + τ2) + (1 − τ2 − 2ν)

$

1 − 2νx ̄ + 4(νx ̄)2

%

νx ̄

2(1 + τ2) + (1 − 6νx ̄)

* (1 − τ2 − 2ν)2

4(1 + τ2) − (1 + τ2)212(1 + τ1)

+

νx ̄x ̄

−

* (1 + τ2)2(1 − τ1 − 2ν)

12(1 + τ1) − (1 − τ2 − 2ν)312(1 + τ2)

+

νx ̄x ̄x ̄

<n−1/2

j−1/2

(7.62)

(νx ̄x ̄)

n

j = 1

2

&$

(2νx ̄ − 1)(νx ̄)

2 + [1 − ν(1 − 6νx ̄)]νx ̄x ̄ − (1 − ν2)νx ̄x ̄x ̄

%n−1/2

j+1/2

+ $(1 − 2νx ̄)(νx ̄)

2 + [1 + ν(1 − 6νx ̄)]νx ̄x ̄ + (1 − ν2)νx ̄x ̄x ̄

%n−1/2

j−1/2

'

(7.63)

42
and

(νx ̄x ̄x ̄)

n

j =

$

(1 − 6νx ̄)νx ̄x ̄ − (1 − τ1 + 2ν)νx ̄x ̄x ̄

%n−1/2

j+1/2 − $

(1 − 6νx ̄)νx ̄x ̄ + (1 − τ1 − 2ν)νx ̄x ̄x ̄

%n−1/2

j−1/2

2

$

1+(τ1)n

j

% (7.64)

respectively.

In the current study, we assume that

|(ux ̄x ̄x ̄)

n

j | ! |(ux ̄x ̄)

n

j | ! |(ux ̄)

n

j | !| unj | (j, n) ∈ Ω (7.65)

|νn

j | ≤ 1 (j, n) ∈ Ω (7.66)

and

|(τ1)

n

j | = O(|νnj |) and |(τ2)

n

j | = O(|νnj |), (j, n) ∈ Ω (7.67)

Note that, by using Eqs. (7.35) and (7.42), one has

u∗(xj ±

∆x

4 , tn; j, n) = un

j ± (ux ̄)

n

j +

1

2

(ux ̄x ̄)

n

j ±

1

6

(ux ̄x ̄x ̄)

n

j (j, n) ∈ Ω (7.68)

Thus, Eq. (7.65) simply states that the magnitude of a term on the right side of Eq. (7.68) is much larger

than that of the next higher-order term to its right. In view of Eqs. (7.42), obviously, this assumption can

be made valid by choosing a small enough value of ∆x. Moreover, by using Eqs. (7.43), it can be shown that

Eqs. (7.65) and (7.66) ⇔

|(νx ̄x ̄x ̄)

n

j | !| (νx ̄x ̄)

n

j | !| (νx ̄)

n

j | !| νnj | ≤ 1 (j, n) ∈ Ω (7.69)

With the aid of Eqs. (7.42) and (7.43), it is seen that the assumption |(νx ̄)n

j | !| νnj | in Eq. (7.69) ⇔

∆x(ux)

n

j ! unj (j, n) ∈ Ω (7.70)

Within a smooth solution region, (ux)n

j is a good approximation of the value of ∂u/∂x at (j, n). Thus

Eq. (7.70) implies that

|un−1/2

j+1/2 − un−1/2j−1/2 | !| un−1/2j±1/2 | (j, n) ∈ Ω (7.71)

By using Eq. (7.43), one can cast Eq. (7.71) into its dimensionless form, i.e.,

|νn−1/2

j+1/2 − νn−1/2j−1/2 | ! |νn−1/2j±1/2 | (j, n) ∈ Ω (7.72)

In turn, Eq. (7.72) implies that

|νn−1/2

j±1/2 − ν ̃n

j | ! |ν ̃nj | (j, n) ∈ Ω (7.73)

with

ν ̃n

j

def

= 1

2

8

νn−1/2

j+1/2 + νn−1/2j−1/2

9

(j, n) ∈ Ω (7.74)

Eq. (7.73) implies that

νn−1/2

j±1/2

.

=  ̃νn

j (j, n) ∈ Ω (7.75)

Hereafter the symbol “ .=” denote the fact that the magnitude of the difference of the quantities on the both

sides of the symbol is much less than the magnitude of either quantity itself. In other words, either is a close

approximation of another.

43
The stability study will be carried out using Eqs. (7.61)–(7.64). In these dimensionless equations, the

independent mesh variables at any mesh point (j, n) ∈ Ω are νn

j , (νx ̄)nj , (νx ̄x ̄)nj , and (νx ̄x ̄x ̄)nj . Hereafter the

round-off errors associated with these variables are denoted by δνn

j , δ(νx ̄)nj , δ(νx ̄x ̄)nj , and δ(νx ̄x ̄x ̄)nj , respectively. Because the magnitude of the round-off error of a quantity with a smaller magnitude

generally is not necessary smaller than that of a quantity with a larger magnitude, in contrast

to Eq. (7.69), we assume that

|δ(τ1)

n

j | ≈ |δ(τ2)

n

j | ≈ |δνnj | ≈ |δ(νx ̄)

n

j | ≈ |δ(νx ̄x ̄)

n

j | ≈ |δ(νx ̄x ̄x ̄)

n

j | (j, n) ∈ Ω (7.76)

Hereafter (i) the symbol “≈” implies that the quantities on both sides of this symbol have the same order

of magnitude, and (ii) δ(τ1)n

j and δ(τ2)nj denote the round-off errors of (τ1)nj and (τ2)nj , respectively. Also

because round-off errors could vary erratically from one mesh point to another, in contrast to

Eq. (7.72), in general, we have

|δνn−1/2

j+1/2 − δνn−1/2j−1/2 | ≈ |δνn−1/2j±1/2 | ≈ |δ(τ1)

n

j | ≈ |δ(τ2)

n

j | (j, n) ∈ Ω (7.77)

Let x%, / = 1, 2, 3, 4, be independent mesh variables at a mesh point. Let δx%, / = 1, 2, 3, 4 be their

round-off errors, respectively. Then because each δx% is a very small perturbation of x%, the round-off errors

of a differentiable function f(x1, x2, x3, x4) can be approximated by

δf(x1, x2, x3, x4) def

= I

4

%=1

∂f

∂x%

δx% (7.78)

As an example, let

x1

def

= νn−1/2

j+1/2 , x2

def = (νx ̄)

n−1/2

j+1/2 , x3

def = (νx ̄x ̄)

n−1/2

j+1/2 , and x4

def = (νx ̄x ̄x ̄)

n−1/2

j+1/2 (j, n) ∈ Ω (7.79)

Then the expression on the right side of Eq. (7.61) associated with the mesh point (j + 1/2, n − 1/2) can be

expressed as g(x1, x2, x3, x4) where

g(x1, x2, x3, x4) def

= 1

2

"

1 − x1

2

#

x1 − 1

2

$

1 − (x1)

2%

x2 − 1

3

$

3(x1)

2 − 1

%

(1 − 2x2)(x2)

2

+

x1

3

$

1 − (x1)

2%

(1 − 6x2)x3 +

1

6

$

1 − (x1)

2%2

x4

(7.80)

Eq. (7.80) implies that

∂g

∂x1

= 1 − x1

2 + x1x2 − 2x1(1 − 2x2)(x2)

2 +

$

1 − 3(x1)2%

3 (1 − 6x2)x3 − 2x1x43

$

1 − (x1)

2% (7.81)

∂g

∂x2

= −

$

1 − (x1)2%

2 − 23

$

3(x1)

2 − 1

%

x2(1 − 3x2) − 2x1x3

$

1 − (x1)

2% (7.82)

∂g

∂x3

= x1

$

1 − (x1)2%

3 − 2x1x2

$

1 − (x1)

2% (7.83)

and

∂g

∂x4

=

$

1 − (x1)2

%2

6 (7.84)

Note that Eq. (7.69) implies that

|x4| ! |x3| !| x2| !| x1| ≤ 1 (7.85)

44
As such, by excluding the case

|1 − (x1)

2| ! 1 (7.86)

one concludes that the expression on the right side of each of Eqs. (7.81)–(7.83) is dominated by the first

term which involves only the variable x1. Moreover, because Eq. (7.76) implies that

|δx4| ≈ |δx3| ≈ |δx2| ≈ |δx1| (7.87)

With the aid of Eqs. (7.78) and (7.79), one concludes that

δg(x1, x2, x3, x4) .=

81 − x1

2

9

δx1 −

*1 − (x1)2

2

+

δx2 +

x1

$

1 − (x1)2%

3

δx3 +

$

1 − (x1)2%2

6

δx4

=

/81 − ν

2

9

δν −

81 − ν2

2

9

δνx ̄ + ν(1 − ν2)

3

δνx ̄x ̄ + (1 − ν2)2

6

δνx ̄x ̄x ̄

0n−1/2

j+1/2

(7.88)

Next, by combining Eqs. (7.75) and (7.88), one arrives at the conclusion

δg(x1, x2, x3, x4) .

=

81 − ν ̃

2

9n

j

δνn−1/2

j+1/2 −

81 − ν ̃2

2

9n

j

δ(νx ̄)

n−1/2

j+1/2

+

*ν ̃(1 − ν ̃2)

3

+n

j

δ(νx ̄x ̄)

n−1/2

j+1/2 +

* (1 − ν ̃2)2

6

+n

j

δ(νx ̄x ̄x ̄)

n−1/2

j+1/2

(7.89)

Note that, according to Eqs. (7.75) and (7.79), Eq. (7.86) ⇔

"

ν ̃n

j

#2 .

= 1 (7.90)

As such the validity of Eq. (7.89) has been established except the case Eq. (7.90).

By excluding the case Eq. (7.90) and using similar arguments, one can also derive the counterpart of

Eq. (7.89) for the expression on the right side of Eq. (7.61) associated with the mesh point (j −1/2, n−1/2).

After combining Eq. (7.89) and its counterpart, one concludes that, excluding the case Eq. (7.90),

δνn

j

.

=

81 − ν ̃

2

9n

j

δνn−1/2

j+1/2 +

81 +  ̃ν

2

9n

j

δνn−1/2

j−1/2 −

81 − ν ̃2

2

9n

j

*

δ(νx ̄)

n−1/2

j+1/2 − δ(νx ̄)

n−1/2

j−1/2

+

+

*ν ̃(1 − ν ̃2)

3

+n

j

*

δ(νx ̄x ̄)

n−1/2

j+1/2 − δ(νx ̄x ̄)

n−1/2

j−1/2

+

+

* (1 − ν ̃2)2

6

+n

j

*

δ(νx ̄x ̄x ̄)

n−1/2

j+1/2 − δ(νx ̄x ̄x ̄)

n−1/2

j−1/2

+

(7.91)

The proof of Eq. (7.91) is completed if one can show that it is valid even for the case Eq. (7.90). Note

that Eq. (7.90) ⇔ either

ν ̃n

j

.

= 1 (7.92)

or

ν ̃n

j

.

= −1 (7.93)

but not both. On the other hand, it can be shown easily that (i)

δνn

j

.

= δνn−1/2

j−1/2 (7.94)

for the case Eq. (7.92), and (ii)

δνn

j

.

= δνn−1/2

j+1/2 (7.95)

45
for the case Eq. (7.93). By comparing Eq. (7.91) with Eqs. (7.94) and (7.95), it is clear that Eq. (7.91) is

also valid for the case Eq. (7.90). QED.

Except for a possible complication to be discussed later, in a similar manner one can derive from

Eqs. (7.62)–(7.64) the following results:

δ(νx ̄)

n

j

.

=

* 1

2(1 + τ2)

+n

j

8

δνn−1/2

j+1/2 − δνn−1/2j−1/2

9

−

*1 − τ2 + 2 ̃ν

2(1 + τ2)

+n

j

δ(νx ̄)

n−1/2

j+1/2

−

* 1 − τ2 − 2 ̃ν

2(1 + τ2)

+n

j

δ(νx ̄)

n−1/2

j−1/2 +

* (1 − τ2 + 2 ̃ν)2

4(1 + τ2) − (1 + τ2)212(1 + τ1)

+n

j

δ(νx ̄x ̄)

n−1/2

j+1/2

−

* (1 − τ2 − 2 ̃ν)2

4(1 + τ2) − (1 + τ2)212(1 + τ1)

+n

j

δ(νx ̄x ̄)

n−1/2

j−1/2

+

* (1 + τ2)2(1 − τ1 + 2 ̃ν)

12(1 + τ1) − (1 − τ2 + 2 ̃ν)312(1 + τ2)

+n

j

δ(νx ̄x ̄x ̄)

n−1/2

j+1/2

+

* (1 + τ2)2(1 − τ1 − 2 ̃ν)

12(1 + τ1) − (1 − τ2 − 2 ̃ν)312(1 + τ2)

+n

j

δ(νx ̄x ̄x ̄)

n−1/2

j−1/2

(7.96)

δ(νx ̄x ̄)

n

j

.

= 1

2

5

(1 − ν ̃)

n

j δ(νx ̄x ̄)

n−1/2

j+1/2 + (1 +  ̃ν)

n

j δ(νx ̄x ̄)

n−1/2

j−1/2 − (1 − ν ̃2)

n

j

*

δ(νx ̄x ̄x ̄)

n−1/2

j+1/2 − δ(νx ̄x ̄x ̄)

n−1/2

j−1/2

+<

(7.97)

and

δ(νx ̄x ̄x ̄)

n

j

.

=

* 1

2(1 + τ1)

+n

j

*

δ(νx ̄x ̄)

n−1/2

j+1/2 − δ(νx ̄x ̄)

n−1/2

j−1/2

+

−

*1 − τ1 + 2 ̃ν

2(1 + τ1)

+n

j

δ(νx ̄x ̄x ̄)

n−1/2

j+1/2 −

* 1 − τ1 − 2 ̃ν

2(1 + τ1)

+n

j

δ(νx ̄x ̄x ̄)

n−1/2

j−1/2

(7.98)

In deriving Eqs. (7.96) and (7.98) from Eqs. (7.62) and (7.64), the reader should bear in mind the

simplifying convention established earlier for the parameters τ1 and τ2. As an example, we have

/ ν

2(1 + τ2)

0n−1/2

j±1/2

def

= νn−1/2

j±1/2

2

$

1+(τ2)n

j

% (7.99)

Thus, by using Eqs. (7.75) and (7.78), Eq. (7.99) implies that

δ

/ ν

2(1 + τ2)

0n−1/2

j±1/2

.

=

/ 1

2(1 + τ2)

0n

j

δνn−1/2

j±1/2 −

/ ν ̃

2(1 + τ2)2

0n

j

δ(τ2)

n

j (7.100)

In turn, Eq. (7.100) implies that

δ

/ ν

2(1 + τ2)

0n−1/2

j+1/2

− δ

/ ν

2(1 + τ2)

0n−1/2

j−1/2

.

=

/ 1

2(1 + τ2)

0n

j

8

δνn−1/2

j+1/2 − δνn−1/2j−1/2

9

(7.101)

The expression on the right side of Eq. (7.101) is the first term on the right side of Eq. (7.96).

Eqs. (7.91) and (7.96)–(7.98) can be cast into the matrix form

δ"q(j, n) .= $

QG−(τ1, τ2, ν ̃)

%n

j δ"q(j + 1/2, n − 1/2) + $

QG+(τ1, τ2, ν ̃)

%n

j δ"q(j − 1/2, n − 1/2) (j, n) ∈ Ω (7.102)

46
where

δ"q(j, n) def=





δνn

j

δ(νx ̄)n

j

δ(νx ̄x ̄)n

j

δ(νx ̄x ̄x ̄)n

j





(j, n) ∈ Ω (7.103)

and

$

QG±(τ1, τ2, ν ̃)

%n

j

def

=





qˆ

±

11 qˆ

±

12 qˆ

±

13 qˆ

±

14

qˆ

±

21 qˆ

±

22 qˆ

±

23 qˆ

±

24

0 0 ˆq±

33 qˆ

±

34

0 0 ˆq±

43 qˆ

±

44





n

j

(j, n) ∈ Ω (7.104)

with

qˆ

±

11

def

= 1 ± ν ̃

2 , qˆ

±

12

def

= ±(1 − ν ̃2)

2 , qˆ

±

13

def

= ∓ν ̃(1 − ν ̃2)

3 , qˆ

±

14

def

= ∓(1 − ν ̃2)2

6 (7.105)

qˆ

±

21

def

= ∓1

2(1 + τ2)

, qˆ

±

22

def

= −(1 − τ2 ∓ 2 ̃ν)

2(1 + τ2) , qˆ

±

23

def

= ± (1 + τ2)2

12(1 + τ1) ∓ (1 − τ2 ∓ 2 ̃ν)24(1 + τ2)

qˆ

±

24

def

= (1 + τ2)2(1 − τ1 ∓ 2 ̃ν)

12(1 + τ1) − (1 − τ2 ∓ 2 ̃ν)312(1 + τ2)

(7.106)

and

qˆ

±

33

def

= 1 ± ν ̃

2 , qˆ

±

34

def

= ±(1 − ν ̃2)

2 , qˆ

±

43

def

= ∓1

2(1 + τ1)

, qˆ

±

44

def

= −(1 − τ1 ∓ 2 ̃ν)

2(1 + τ1) (7.107)

A comparison of Eqs. (7.102)–(7.107) with Eqs. (6.7), (5.54), and (6.8)–(6.11) reveals a striking similarity

in the structures of these two sets of equations. On the other hand, it was explained in Sec. 5.4 that: (i)

computational instability would occur if round-off errors are amplified without bound as they propagate down

the subsequent time levels; and (ii) for linear schemes like the a-%-4 and c-τ-4 schemes, the time evolution of

the mesh variables and that of the associated round-off errors are governed by the same system of equations.

As a result of the facts stated above, one now arrives at the conclusion that the stability conditions of

the nonlinear inviscid Burgers c-τ-4 scheme can be obtained from those of the linear c-τ scheme

by replacing the constant coefficients ν, τ1, and τ2 with the mesh-point dependent coefficients

ν ̃n

j , (τ1)nj , and (τ2)nj , respectively.

7.3. The inviscid Burgers c-τ ∗-4 schemes

Let h1(s) and h2(s) be any two functions that meet the same requirements imposed on the function

h(s) introduced in Sec. 3.2. Then, because of the close relations among the c-τ, c-τ-4 and current schemes

established above, here an “ideal” inviscid Burgers c-τ ∗-4 scheme is defined by an extension of Eq. (6.21),

i.e.,

(τ1)

n

j = h1

"

( ̃νn

j )

2# and (τ2)n

j = h2

"

( ̃νn

j )

2# (|ν ̃n

j | < 1) (7.108)

According to Eq. (3.22), one can assume h1(s) = h2(s) = √s. As such a special inviscid Burgers c-τ ∗-4

scheme can be defined by the relation

(τ1)

n

j = (τ2)

n

j = |ν ̃nj | (|ν ̃nj | < 1) (7.109)

47
For robust real-world applications, we may also introduce a special inviscid Burgers c-τ-4 scheme which

is defined by an extension of Eq. (6.23), i.e.,

(τ1)

n

j = β1 |ν ̃nj | and (τ2)

n

j = β2 |ν ̃nj | (β1 ≥ 1; β2 ≥ 1; |ν ̃nj | < 1) (7.110)

where β1 ≥ 1 and β2 ≥ 1 are two adjustable parameters. In the rest of this section, we consider only the

special schemes which satisfy at least one of Eqs. (7.59), (7.60), and (7.108)–(7.110). As a result of Eq. (7.43),

for these schemes,

(τ1)

n

j → 0 and (τ2)

n

j → 0 as ∆t → 0 (7.111)

7.4. Shock capturing schemes

As explained in Sec. 2.3, the current Burgers solver needed to be modified so that it is capable of

resolving shocks. Specifically, we need to turn the expression on the right side of Eq. (7.52) and/or that of

Eq. (7.58) into proper weighted averages discussed in Sec. 2.3 and Sec. 4.

First we consider Eq. (7.52). Let

ux ̄x ̄(P +) def

=

8 ∆x

4

92

uxx(P +) and ux ̄x ̄(P −) def=

8 ∆x

4

92

uxx(P −) (7.112)

where uxx(P +) and uxx(P −) are defined by Eqs. (7.47), (7.48), (7.50), and (7.51). Then it can be shown

that

ux ̄x ̄(P +) = $

(1 − 6νx ̄)ux ̄x ̄ − (1 − τ1 + 2ν)ux ̄x ̄x ̄

%n−1/2

j+1/2 (7.113)

and

ux ̄x ̄(P −) = $

(1 − 6νx ̄)ux ̄x ̄ + (1 − τ1 − 2ν)ux ̄x ̄x ̄

%n−1/2

j−1/2 (7.114)

Thus Eq. (7.52) implies that

(ux ̄x ̄x ̄)

n

j = ux ̄x ̄(P +) − ux ̄x ̄(P −)

2

$

1+(τ1)n

j

% = 1

2 (ux ̄x ̄x ̄− + ux ̄x ̄x ̄+)

n

j (7.115)

where

(ux ̄x ̄x ̄−)

n

j

def

= (ux ̄x ̄)n

j − ux ̄x ̄(P −)

1+(τ1)n

j

and (ux ̄x ̄x ̄+)

n

j

def

= ux ̄x ̄(P +) − (ux ̄x ̄)n

j

1+(τ1)n

j

(7.116)

with the understanding that (ux ̄x ̄)n

j is that defined in Eq. (7.44). Eqs. (7.115) and (7.116) can be considered

as the substitution forms of Eqs. (3.10) and (4.1)–(4.3).

According to Eq. (7.115), (ux ̄x ̄x ̄)n

j evaluated using Eq. (7.52) is the simple average of (ux ̄x ̄x ̄−)nj and

(ux ̄x ̄x ̄+)n

j . Moreover, by using Eqs. (7.111), (7.113)–(7.116) and

νn

j → 0, (νx ̄)

n

j → 0, ,(νx ̄x ̄)

n

j → 0, and (νx ̄x ̄x ̄)

n

j → 0 as ∆t → 0 (7.117)

(which follows from Eq. (7.43)), one can show that

lim

∆t→0

(ux ̄x ̄x ̄−)

n

j = lim

∆t→0

(ux ̄x ̄x ̄+)

n

j = lim

∆t→0

(ux ̄x ̄x ̄)

n

j = 1

2

*

(ux ̄x ̄ − ux ̄x ̄x ̄)

n−1/2

j+1/2 − (ux ̄x ̄ + ux ̄x ̄x ̄)

n−1/2

j−1/2

+

(7.118)

Because (ux ̄x ̄x ̄−)n

j and (ux ̄x ̄x ̄+)nj satisfy Eq. (7.118), the simple average in Eq. (7.115) can be turned into

an weighted average using any of the techniques described in Sec. 4. As an example, consider the technique

presented in Sec. 4.4. The current counterparts of Eqs. (4.22)–(4.24) are

$

(η1)±

%n

j

def

= |(ux ̄x ̄x ̄±)n

j |

α1

min{|(ux ̄x ̄x ̄−)n

j |

α1 , |(ux ̄x ̄x ̄+)n

j |

α1 } − 1 (α1 > 0 and |(ux ̄x ̄x ̄−)

n

j | + |(ux ̄x ̄x ̄+)

n

j | > 0) (7.119)

48
$

(  ̃w1)+

%n

j

def

=

5 1 + σ1 (η1)−

2 + σ1[(η1)+ + (η1)−]

<n

j

( ̃νn

j "= 0) (7.120)

and

$

(  ̃w1)−

%n

j

def

=

5 1 + σ1 (η1)+

2 + σ1[(η1)+ + (η1)−]

<n

j

( ̃νn

j "= 0) (7.121)

respectively. Here (i) α1 > 0 is a preset parameter in the order of 1, and (ii)

(σ1)

n

j

def

= (σ1)o

|ν ̃n

j | ((σ1)o > 0 and  ̃νn

j "= 0) (7.122)

with (σ1)o being a preset parameter in the order of 1. With the above definitions, the current counterpart

of Eq. (4.30) can be expressed as

(ux ̄x ̄x ̄)

n

j = $

(  ̃w1)−

%n

j (ux ̄x ̄x ̄−)

n

j + $

(  ̃w1)+

%n

j (ux ̄x ̄x ̄+)

n

j (7.123)

Next consider Eq. (7.58). Note that, with the aid of Eq. (7.52) and (7.57), it can be shown that Eq. (7.58)

⇔

(ux ̄)

n

j = u4∗(P +) − u4∗(P −)

2

$

1+(τ2)n

j

% −

$

1+(τ2)n

j

%2

6 (ux ̄x ̄x ̄)

n

j (7.124)

In the following, with the aid of Eq. (7.124), (ux ̄±)n

j will be defined so that

(ux ̄)

n

j = 1

2

(ux ̄− + ux ̄+)

n

j (7.125)

and

lim

∆t→0

(ux ̄−)

n

j = lim

∆t→0

(ux ̄+)

n

j = lim

∆t→0

(ux ̄)

n

j = 1

2

E6

u − ux ̄ +

1

3

ux ̄x ̄

7n−1/2

j+1/2

−

6

u + ux ̄ +

1

3

ux ̄x ̄

7n−1/2

j−1/2

F

(7.126)

Note that the validity of the last equality sign follows directly from Eqs. (7.58), (7.111), and (7.117).

To facilitate the development, note that here (i) P + and P − are the points depicted in Fig. 3, except

that the parameter τ shown in the figure be replaced by the mesh-point dependent parameter (τ2)n

j ; and

(ii) u4∗(P +) and u4∗(P −) defined in Eq. (7.56) are 4th-order Taylor’s approximations of the values of u at

points P + and P −, respectively. Alternatively, these two values can also be approximated by

u ̃4∗(P ±) def= u∗ "

x(P ±), tn; j, n# (7.127)

respectively. With the aid of Eqs. (7.35), (7.42), (7.53), and (7.54), Eq. (7.127) implies that

u ̃4∗(P ±) = un

j ± $

1+(τ2)

n

j

%

(ux ̄)

n

j +

$

1+(τ2)n

j

%2

2 (ux ̄x ̄)

n

j ±

$

1+(τ2)n

j

%3

6 (ux ̄x ̄x ̄)

n

j (7.128)

In turn, Eq. (7.128) implies that

u ̃4∗(P +) − un

j

1+(τ2)n

j

− 1+(τ2)n

j

2 (ux ̄x ̄)

n

j = un

j − u ̃4∗(P −)

1+(τ2)n

j

+ 1+(τ2)n

j

2 (ux ̄x ̄)

n

j (7.129)

Moreover, by using Eqs. (7.44), (7.46), (7.52), (7.57), (7.58), (7.111), (7.117) and (7.128), one concludes that

lim

∆t→0

u ̃4∗(P ±) = lim ∆t→0u4∗(P ±) = 6u ∓ ux ̄ +

1

2

ux ̄x ̄ ∓

1

6

ux ̄x ̄x ̄

7n−1/2

j±1/2

(7.130)

49
An immediate result of Eqs. (7.129) and (7.130) is

lim

∆t→0

(u#

x ̄+)

n

j = lim

∆t→0

(u#

x ̄−)

n

j (7.131)

where

(u#

x ̄+)

n

j

def

= u4∗(P +) − un

j

1+(τ2)n

j

− 1+(τ2)n

j

2 (ux ̄x ̄)

n

j (7.132)

and

(u#

x ̄−)

n

j

def

= un

j − u4∗(P −)

1+(τ2)n

j

+ 1+(τ2)n

j

2 (ux ̄x ̄)

n

j (7.133)

As a result of Eqs. (7.132) and (7.133), one has

(u#

x ̄+)

n

j + (u#x ̄−)

n

j = u4∗(P +) − u4∗(P −)

1+(τ2)n

j

(7.134)

Combining Eqs. (7.124), (7.131) and (7.134), one concludes that the requirements Eqs. (7.125) and (7.126)

are met by the definition

(ux ̄±)

n

j

def = (u#

x ̄±)

n

j −

$

1 + (τ2)n

j

%2

6 (ux ̄x ̄x ̄)

n

j (7.135)

Because (ux ̄+)n

j and (ux ̄−)nj meet the condition that they have the same value in the limit of ∆t → 0,

the simple average on the right side of Eq. (7.125) can also be converted into an weighted average using

any of the techniques described in Sec. 4. As an example, consider the technique presented in Sec. 4.4. The

current counterparts of Eqs. (7.119)–(7.121) are

$

(η2)±

%n

j

def

= |(ux ̄±)n

j |

α2

min{|(ux ̄−)n

j |

α2 , |(ux ̄+)n

j |

α2 } − 1 (α2 > 0 and |(ux ̄−)

n

j | + |(ux ̄+)

n

j | > 0) (7.136)

$

(  ̃w2)+

%n

j

def

=

5 1 + σ2 (η2)−

2 + σ2[(η2)+ + (η2)−]

<n

j

( ̃νn

j "= 0) (7.137)

and

$

(  ̃w2)−

%n

j

def

=

5 1 + σ2 (η2)+

2 + σ2[(η2)+ + (η2)−]

<n

j

( ̃νn

j "= 0) (7.138)

respectively. Here (i) α2 > 0 is a preset parameter in the order of 1, and (ii)

(σ2)

n

j

def

= (σ2)o

|ν ̃n

j | ((σ2)o > 0 and  ̃νn

j "= 0) (7.139)

with (σ2)o being a preset parameter in the order of 1. With the above definitions, the current counterpart

of Eq. (7.123) can be expressed as

(ux ̄)

n

j = $

(  ̃w2)−

%n

j (ux ̄−)

n

j + $

(  ̃w2)+

%n

j (ux ̄+)

n

j (7.140)

Note that a result of Eqs. (7.125) and (7.135) is

(ux ̄)

n

j = 1

2

"

u#

x ̄+ + u#x ̄−

#n

j −

$

1+(τ2)n

j

%2

6 (ux ̄x ̄x ̄)

n

j (7.141)

Thus an alternative shock-capturing form of Eq. (7.125) can be obtained by replacing the simple-average

term on the right side of Eq. (7.141) with an weighted average of (u#

x ̄+)n

j and (u#

x ̄−)n

j .

50
8. A sketch of multidimensional extensions

Consider the 2D advection equation

∂u

∂t + ax

∂u

∂x

+ ay

∂u

∂y = 0 (8.1)

where ax, and ay are constants. Let x1 = x, x2 = y, and x3 = t be the coordinates of a three-dimensional

Euclidean space E3. By using Gauss’ divergence theorem in the space-time E3, it can be shown that Eq. (8.1)

is the differential form of the integral conservation law

!

S(V )

"

h · d"s = 0 (8.2)

Here (i) S(V ) is the boundary of an arbitrary space-time region V in E3, (ii)

"

h def = (axu, ayu, u) (8.3)

is a current density vector in E3, and (iii) d"s = dσ "n with dσ and "n, respectively, being the area and the

outward unit normal of a surface element on S(V ).

In addition, by using an argument similar to that was used to obtain the conservation condition Eq. (5.9),

it can be shown that, for a smooth solution u(x, y, t) of Eq. (8.1), we have

!

S(V )

"

hxx · d"s = 0,

!

S(V )

"

hxy · d"s = 0, and !

S(V )

"

hyy · d"s = 0 (8.4)

where

"

hxx

def

= ∂2"

h

∂x2 ≡

6

ax

∂2u

∂x2 , ay

∂2u

∂x2 ,

∂2u

∂x2

7

, "hxy

def

= ∂2"

h

∂x∂y ≡

6

ax

∂2u

∂x∂y

, ay

∂2u

∂x∂y

, ∂2u

∂x∂y

7

and "hyy

def

= ∂2"

h

∂y2 ≡

6

ax

∂2u

∂y2 , ay

∂2u

∂y2 ,

∂2u

∂y2

7

(8.5)

In a 2D CESE scheme [13], a spatial domain can be divided into triangles of arbitrary shape with the

understanding that any two neighboring triangles share a common side. Each triangle and the associated

spatial solution point is assigned with an integer cell number j. The x- and y- coordinates of a spatial solution

point j are denoted by xj and yj , respectively. In addition, a space-time solution point (j, n) is defined to

be the point with x = xj , y = yj , and t = t

n def

= n∆t. As in the 1D case, each (j, n) is associated with a

solution element SE(j, n). However, corresponding to the three sides of a triangle, each (j, n) is associated

with three basic space-time conservation elements CEr(j, n), r = 1, 2, 3. Again these BCEs are constructed

so that the boundary of each BCE can be divided into component parts with each of them belonging to a

unique SE. The details are described in [13].

To construct a 4th-order scheme, for any (x, y, t) ∈ SE(j, n), the dependent variable u is approximated

by a numerical analogue of the third-order Taylor’s expansion, i.e.,

u∗(x, y, t; j, n) def= un

j + (ux)

n

j (x − xj )+(uy)

n

j (y − yj)+(ut)

n

j (t − t

n) + 1

2

(uxx)

n

j (x − xj )

2

+

1

2

(uyy)

n

j (y − yj )

2 +

1

2

(utt)

n

j (t − t

n)2 + (uxy)n

j (x − xj )(y − yj ) + (uxt)

n

j (x − xj )(t − t

n)

+ (uyt)

n

j (y − yj )(t − t

n) + 1

6

(uxxx)

n

j (x − xj )

3 +

1

6

(uyyy)

n

j (y − yj)

3 +

1

6

(uttt)

n

j (t − t

n)3

+

1

2

(uxxy)

n

j (x − xj )

2(y − yj) + 1

2

(uxyy)

n

j (x − xj )(y − yj)

2 +

1

2

(uxxt)

n

j (x − xj )

2(t − tn)

+

1

2

(uyyt)

n

j (y − yj )

2(t − tn) + 1

2

(uxtt)

n

j (x − xj )(t − t

n)2 +

1

2

(uytt)

n

j (y − yj)(t − t

n)2

+ (uxyt)

n

j (x − xj )(y − yj)(t − t

n)

(8.6)

51
Let u = u∗(x, y, t; j, n) satisfy Eq. (8.1) uniformly in SE(j, n), i.e.,

∂u∗(x, y, t; j, n)

∂t + ax

∂u∗(x, y, t; j, n)

∂x

+ ay

∂u∗(x, y, t; j, n)

∂y ≡ 0 (x, y, t) ∈ SE(j, n) (8.7)

Then it can be shown that

(ut)

n

j = −(axux + ayuy)

n

j (8.8)

(uxt)

n

j = −(axuxx + ayuxy)

n

j (8.9)

(uyt)

n

j = −(axuxy + ayuyy)

n

j (8.10)

(utt)

n

j = −(axuxt + ayuyt) = $

(ax)

2uxx + 2axayuxy + (ay)2uyy%n

j (8.11)

(uxxt)

n

j = −(axuxxx + ayuxxy)

n

j (8.12)

(uxyt)

n

j = −(axuxxy + ayuxyy)

n

j (8.13)

(uyyt)

n

j = −(axuxyy + ayuyyy)

n

j (8.14)

(uxtt)

n

j = −(axuxxt + ayuxyt)

n

j = $

(ax)

2uxxx + 2axayuxxy + (ay)2uxyy%n

j (8.15)

(uytt)

n

j = −(axuxyt + ayuyyt)

n

j = $

(ax)

2uxxy + 2axayuxyy + (ay)2uyyy%n

j (8.16)

and

(uttt)

n

j = −(axuxtt + ayuytt)

n

j = − $

(ax)

3uxxx + 3(ax)2ayuxxy + 3ax(ay)2uxyy + (ay)3uyyy%n

j (8.17)

Substituting Eqs. (8.8)–(8.17) into Eq. (8.6), one has

u∗(x, y, t; j, n) = un

j + (ux)

n

j

$

(x − xj ) − ax(t − t

n)

%

+ (uy)

n

j

$

(y − yj ) − ay(t − t

n)

%

+

1

2

(uxx)

n

j

$

(x − xj ) − ax(t − t

n)

%2

+

1

2

(uyy)

n

j

$

(y − yj ) − ay(t − t

n)

%2

+ (uxy)

n

j

$

(x − xj ) − ax(t − t

n)

%$(y − yj ) − ay(t − t

n)

%

+

1

6

(uxxx)

n

j

$

(x − xj ) − ax(t − t

n)

%3

1

6

(uyyy)

n

j

$

(y − yj ) − ay(t − t

n)

%3

+

1

2

(uxxy)

n

j

$

(x − xj ) − ax(t − t

n)

%2$

(y − yj) − ay(t − t

n)

%

+

1

2

(uxyy)

n

j

$

(x − xj ) − ax(t − t

n)

%$(y − yj ) − ay(t − t

n)

%2

(8.18)

Thus u∗(x, y, t; j, n) is a 3rd-order polynomial function of $

(x − xj ) − ax(t − tn)

%

and $(y − yj) − ay(t − tn)

%

specified by the coefficients un

j , (ux)nj , (uy)nj , (uxx)nj , (uyy)nj , (uxy)nj , (uxxx)nj , (uyyy)nj , (uxxy)nj , and (uxyy)nj .

For a reason to be given later, in the current extension, only eight of these ten coefficients are considered

as the independent mesh variables at (j, n) while two of them (i.e., (uxxy)n

j and (uxyy)nj ) are considered as

functions of another four independent mesh variables at (j, n). In other words, there are twelve independent

mesh variables at each solution point (j, n).

To construct the a-4 scheme for Eq. (8.1), we will impose the following three conservation conditions:

!

S(CEr(j,n))

"

h∗ · d"s = 0, r = 1, 2, 3 (8.19)

where "

h∗(x, y, t; j, n) def

= $axu∗(x, y, t; j, n), ayu∗(x, y, t; j, n), u∗(x, y, t; j, n)

% (8.20)

Eqs. (8.19) and (8.20) are the numerical analogues of Eqs. (8.2) and (8.3), respectively.

52
Note that Eq. (8.18) implies that

∂2u∗(x, y, t; j, n)

∂x2 = (uxx)

n

j + (uxxx)

n

j

$

(x − xj ) − ax(t − t

n)

%

+ (uxxy)

n

j

$

(y − yj ) − ay(t − t

n)

% (8.21)

∂2u∗(x, y, t; j, n)

∂x∂y = (uxy)

n

j + (uxxy)

n

j

$

(x − xj ) − ax(t − t

n)

%

+ (uxyy)

n

j

$

(y − yj ) − ay(t − t

n)

% (8.22)

and

∂2u∗(x, y, t; j, n)

∂y2 = (uyy)

n

j + (uxyy)

n

j

$

(x − xj ) − ax(t − t

n)

%

+ (uyyy)

n

j

$

(y − yj) − ay(t − t

n)

% (8.23)

Thus, at each (j, n), imposition of conservation conditions Eq. (8.19) and the numerical analogue of Eq. (8.4),

i.e.,

!

S(CEr(j,n))

"

h∗

xx · d"s = 0, r = 1, 2, 3 (8.24)

!

S(CEr(j,n))

"

h∗

xy · d"s = 0, r = 1, 2, 3 (8.25)

and

!

S(CEr(j,n))

"

h∗

yy · d"s = 0, r = 1, 2, 3 (8.26)

would result in twelve conditions for the ten coefficients that appear on the right side of Eq. (8.18) if one

assumes that (i)

"

h∗

xx(x, y, t; j, n) def

= $axu∗

xx(x, y, t; j, n), ayu∗xx(x, y, t; j, n), u∗xx(x, y, t; j, n)

% (8.27)

"

h∗

xy(x, y, t; j, n) def

= $axu∗

xy(x, y, t; j, n), ayu∗xy(x, y, t; j, n), u∗xy(x, y, t; j, n)

% (8.28)

and "

h∗

yy(x, y, t; j, n) def

= $axu∗

yy(x, y, t; j, n), ayu∗yy(x, y, t; j, n), u∗yy(x, y, t; j, n)

% (8.29)

and (ii)

u∗

xx(x, y, t; j, n) = ∂2u∗(x, y, t; j, n)

∂x2 (8.30)

u∗

xy(x, y, t; j, n) = ∂2u∗(x, y, t; j, n)

∂x∂y (8.31)

and

u∗

yy(x, y, t; j, n) = ∂2u∗(x, y, t; j, n)

∂y2 (8.32)

To overcome the above problem of over determinacy, we replace Eqs. (8.30)–(8.32) with

u∗

xx(x, y, t; j, n) = (uxx)

n

j + (uxxx)

n

j

$

(x − xj ) − ax(t − t

n)

%

+ (uxxy)

n

j

$

(y − yj ) − ay(t − t

n)

% (8.33)

u∗

xy(x, y, t; j, n) = (uxy)

n

j + (uxyx)

n

j

$

(x − xj ) − ax(t − t

n)

%

+ (uxyy)

n

j

$

(y − yj ) − ay(t − t

n)

% (8.34)

and

u∗

yy(x, y, t; j, n) = (uyy)

n

j + (uyyx)

n

j

$

(x − xj ) − ax(t − t

n)

%

+ (uyyy)

n

j

$

(y − yj ) − ay(t − t

n)

% (8.35)

53
respectively. Note that: (i) the expression on the right sides of Eq. (8.33) is identical to that of Eq. (8.21)

except that the symbol (uxxy)n

j in Eq. (8.21) is replaced by (uxxy)nj in Eq. (8.33); (ii) the expression on

the right side of Eq. (8.34) is identical to that of Eq. (8.22) except that the symbols (uxxy)n

j and (uxyy)nj

in Eq. (8.22) are replaced by (uxyx)n

j and (uxyy)nj , respectively; and (iii) the expression on the right side

of Eq. (8.35) is identical to that of Eq. (8.23) except that the symbol (uxyy)n

j in Eq. (8.23) is replaced by

(uyyx)n

j . In spite of the facts that (i) (uxxy)nj , (uxxy)nj , and (uxyx)nj are the “numerical analogues” of the

value of the same mixed spatial derivative ∂3u/∂x2∂y at (j, n); and (ii) (uxyy)n

j , (uxyy)nj , and (uyyx)nj are the

“numerical analogues” of the value of the same mixed spatial derivative ∂3u/∂x∂y2 at (j, n), here (i) (uxxy)n

j ,

(uxxy)n

j , and (uxyx)nj are treated as different numerical parameters and thus can assume different

values; and (ii) (uxyy)n

j , (uxyy)nj , and (uyyx)nj are also treated as different numerical parameters

and thus can assume different values. In fact, (uxxy)n

j , (uxyx)nj , (uxyy)nj , and (uyyx)nj are considered

as independent mesh variables at (j, n) while (i) (uxxy)n

j is considered as a function of (uxxy)nj

and (uxyx)n

j ; and (ii) (uxyy)nj is considered as a function of (uxyy)nj and (uyyx)nj . In fact, one may

assume that

(uxxy)

n

j = 1

2

(uxxy + uxyx)

n

j and (uxyy)

n

j = 1

2

(uxyy + uyyx)

n

j (8.36)

As such, at each (j, n), there are twelve independent mesh variables, i.e., un

j , (ux)nj , (uy)nj , (uxx)nj ,

(uxxx)n

j , (uxxy)nj , (uxy)nj , (uxyx)nj , (uxyy)nj , (uyy)nj , (uyyx)nj , and (uyyy)nj . They can be determined in terms

of the known mesh variables by imposing the twelve conservations conditions given in Eqs. (8.19) and (8.24)–

(8.26). In fact, with the aid of Eqs. (8.27) and (8.33), the three mesh variables (uxx)n

j , (uxxx)nj , and (uxxy)nj

can be determined in terms of the known mesh variables by imposing the three conditions given in Eq. (8.24).

Similarly, (uxy)n

j , (uxyx)nj , (uxyy)nj can be solved by using Eq. (8.25) while (uyy)nj , (uyyx)nj , (uyyy)nj can be

solved by using Eq. (8.26). After the nine 2nd-order and 3rd-order independent mesh variables at (j, n)

are determined, the last three independent mesh variables, i.e., un

j , (ux)nj , and (uy)nj can be solved using

Eq. (8.19). Thus the a-4 scheme for Eq. (8.1) can be constructed easily. Moreover, with the aid of Eqs. (8.33)–

(8.36), the dissipative extensions of the 2D a-4 scheme can also be constructed easily using well-established

CESE techniques.

We conclude this section with the following remarks:

(a) The identities

∂2

∂x2

6∂u

∂y

7

= ∂2

∂x∂y

6∂u

∂x

7

= ∂2

∂y∂x

6∂u

∂x

7

(8.37)

and

∂2

∂y2

6∂u

∂x

7

= ∂2

∂x∂y

6∂u

∂y

7

= ∂2

∂y∂x

6∂u

∂y

7

(8.38)

are derived from a limiting process in which ∆x, ∆y → 0. Because ∆x and ∆y are finite in a numerical

simulation, logically the numerical analogues of the third-order mixed spatial derivatives are not required

to satisfy the numerical analogues of Eqs. (8.37) and (8.38).

(b) Because (uxxy)n

j and (uxyy)nj appear as 3rd-order terms on the right side of Eq. (8.18), the 4th-order

accuracy requirement of u∗(x, t; j, n) demands only that (uxxy)n

j and (uxyy)nj be 1st-order in accuracy,

i.e., the same orders of accuracy achieved by (ux)n

j and (uy)nj in a 2nd-order 2D CESE scheme. This

rather low requirement for accuracy further justifies the above introduction of (uxxy)n

j , (uxyx)nj , (uxyy)nj ,

and (uyyx)n

j .

(c) For the 3D CESE method, tetrahedrons are the most natural building blocks for structured or unstructured spatial meshes. As it turns out, the high order CESE framework described here can be naturally

extended in a 3D space using spatial meshes built from tetrahedrons.

9. Numerical results

54
To assess the stability and accuracy of the a-%-4 and c-τ-4 schemes, consider the model problem defined

by the PDE Eq. (2.1) and the exact solution

u = ue(x, t) def = sin [2π(x − t)] (9.1)

We have

a = λ = T = 1 (9.2)

where λ = wavelength and T = period. Let (i)

uxe(x, t) def

= ∂ue(x, t)

∂x

, uxxe(x, t) def

= ∂2ue(x, t)

∂x2 , and uxxxe(x, t) def

= ∂3ue(x, t)

∂x3 (9.3)

and (ii) the spatial domain of unit length be divided into K uniform intervals. Thus

∆x = 1/K, ∆t = ν∆x and t = n∆t (9.4)

where n = number of marching cycles (each cycle is formed by two consecutive marching steps), and t =

total marching time.

In Tables 1-14, the numerical errors of several computations using the a-%-4 and c-τ-4 schemes are

presented in terms of the L2 error norms

E(K, n, ν) def

=

JKKL

1

K

I

K

j=1

[un

j − ue(xj , tn)]2 (9.5)

Ex(K, n, ν) def

=

JKKL

1

K

I

K

j=1

[(ux)n

j − uxe(xj , tn)]2 (9.6)

Exx(K, n, ν) def=

JKKL

1

K

I

K

j=1

[(uxx)n

j − uxxe(xj , tn)]2 (9.7)

and

Exxx(K, n, ν) def

=

JKKL

1

K

I

K

j=1

[(uxxx)n

j − uxxxe(xj , tn)]2 (9.8)

Note that (ux)n

j , (uxx)nj , and (uxxx)nj which appear in Eqs. (9.5)–(9.8) are not the normalized mesh variables

defined in Eq. (7.42).

The numerical errors of the a-%-4 and c-τ-4 schemes with periodic boundary conditions are shown in

Tables 1–7. In all simulations, the values of n are adjusted so that t = 5.4 for all the results presented.

The results with ν = 1 generated using the a-%-4 scheme are given in Table 1 (%1 = %2 = 0) and Table 2

(%1 = %2 = 0.5). For these cases where the value of ν is right at the stability boundary, the errors for un

j and

(uxx)n

j are machine zero. On the other hand, the values of Ex and Exxx are reduced by factors of about 16

and 4, respectively, as both K and n double their values. Thus, for both the cases (a) ν = 1 and %1 = %2 = 0;

and (b) ν = 1 and %1 = %2 = 0.5, the a-%-4 scheme is 4th-order and 2nd-order in accuracy for (ux)n

j and

(uxxx)n

j , respectively.

The results with ν = %1 = %2 = 0.5 generated using the a-%-4 scheme are shown in Table 3. It is seen

that, for this case, the a-%-4 scheme is 4th-order in accuracy for both un

j and (ux)nj while it is 2nd-order in

accuracy for both (uxx)n

j and (uxxx)nj .

The results with ν = 0.05 and %1 = %2 = 0.5 generated using the a-%-4 scheme are shown in Table 4.

It is seen that, for this case with a very small value of ν, the a-%-4 scheme has an accuracy higher than

55
4th-order for both un

j and (ux)nj when the meshes used are relatively coarse. For the coarse mesh cases, it

also has an accuracy higher than 2nd-order for both (uxx)n

j and (uxxx)nj . According to the discussions given

in Sec. 2 and 5, the special scheme with %1 = %2 = 0.5 (referred to as the c scheme) has a huge advantage

in computational efficiency over other a-%-4 schemes. As such the c scheme is emphasized in our numerical

study of the a-%-4 scheme.

The numerical results of the three cases with (a) ν = τ1 = τ2 = 1, (b) ν = τ1 = τ2 = 0.5, and (c)

ν = τ1 = τ2 = 0.05, are presented in Tables 5–7, respectively. According to a discussion given in Sec. 6,

all these cases are the so-called c-τ ∗-4 schemes which form a special class of the c-τ-4 schemes and are

computational efficient and also Courant number insensitive. For case (a), the errors for un

j and (uxx)nj are

again machine zero while it is 4th-order and 2nd-order in accuracy for (ux)n

j and (uxxx)nj , respectively. For

case (b), again the scheme is 4th-order in accuracy for un

j . Also it is 3rd-order, 2nd-order and 1st-order in

accuracy for (ux)n

j , (uxx)nj , and (uxxx)nj when the meshes used become finer and finer. In fact, for this case

we have (i) E = 2.45 × 10−8, Ex = 2.14 × 10−6, Exx = 7.31 × 10−3, and Exxx = 0.299 when K = 200 and

n = 2160; and (ii) E = 1.53 × 10−9, Ex = 2.76 × 10−7, Exx = 1.83 × 10−3, and Exxx = 0.161 when K = 400

and n = 4320. For case (c), as the meshes used become finer and finer, the scheme converges to 4th-order,

3rd-order, 2nd-order, and 1st-order in accuracy for un

j , (ux)nj , (uxx)nj , and (uxxx)nj , respectively.

According to discussions presented in Sec. 2, 3, 5, and 6, in the case where |ν| ! 1, the c-τ ∗-4 schemes

should be much more accurate than the 4th-order c scheme. This conclusion is indeed consistent with the

numerical results presented in Tables 4 and 7.

The results presented in Tables 8–14 are for the cases which are identical to those associated with the

results presented in Tables 1–7 except that the periodic boundary conditions used in the latter cases are

replaced in the former cases by the boundary conditions in which (i) the exact solution Eq. (9.1) is specified

at the left spatial boundary x = 0, and (ii) the 4th-order non-reflecting boundary conditions Eqs. (5.83),

(5.86), (5.88), and (5.89) are imposed on the right spatial boundary x = 1. As shown by the results presented

in Tables 8–14, this changes of boundary conditions has no impact on the orders of accuracy for the a-%-4 and

c-τ ∗-4 schemes. In fact, this is true even for the numerical results associated with the right spatial boundary

x = 1.

The last numerical results to be presented is a shock solution for the problem defined by Eq. (7.1) and

the initial condition

u =

: 1 if x < 0

0 if x > 0

(t = 0) (9.9)

The exact solution of this problem is

ue(x, t) =







1 if x < t/2

0 if x > t/2

(9.10)

The numerical solution is generated using a scheme described in Sec. 7. Specifically, we assume

Eqs. (7.44), (7.46), (7.110), (7.119)–(7.123) and the alternative shock-capturing form of Eq. (7.125) described

following Eqs. (7.141) (the same procedure Eqs. (7.136)–(7.140) is used to obtain the weighted average of

(u#

x ̄+)n

j and (u#

x ̄−)n

j ). We also assume that

β1 = β2 = α1 = α2 = 2 and (σ1)o = (σ2)o = 1 (9.11)

The spatial computational domain is defined by −0.505 ≤ x ≤ 0.505 and divided into 101 uniform

interval, i.e ∆x = 0.01. We also assume that ∆t = 0.009 and n = 50, i.e., t = n∆t = 0.45. According to

Eq. (9.10),

ue(x, 0.45) = : 1 if x < 0.225

0 if x > 0.225

(9.12)

56
The numerical solution is shown against the exact solution in Fig. 4. Note that even though the maximal

CF L number = 0.9 for the exact solution, numerically, the maximal CF L number reached during the

simulation is 0.967.

10. Conclusions and discussions

A thorough and rigorous discussion of a new approach for constructing highly stable high-order CESE

schemes has been presented. Each of the linear and nonlinear 4th-order schemes described in Sec. 5–7 has

the same stencil and same stability conditions of its 2nd-order version and still retains all other advantages of

the latter scheme. In particular, the c-τ-4 scheme and its extensions for solving the inviscid Burgers equation

Eq. (7.1) are completely explicit, i.e., implementation of these schemes does not require inverting systems

of linear/nonlinear algebraic equations locally or globally. As such the new approach can overcome several

common shortcomings of traditional high order schemes that were listed in the abstract.

Because the new 4th-order schemes are built from their 2nd-order versions, in addition to providing a

review of the CESE core ideas and advantages in Sec. 1, a rather extensive review of several key 2nd-order

CESE schemes were provided in Sec. 2–4. The special CESE techniques (e.g., the advanced weighted-average

techniques described in Sec. 4) and the stability conditions presented in this review form a foundation for

the 4th-order schemes discussed in Sec. 5–7.

To explain clearly how a 4th-order scheme can have the same stencil and same stability conditions of

the 2nd-order scheme from which it is built, the conceptual leap that underlines the current development

is clearly explained early in Sec. 5. This key idea is then thoroughly validated in the rest of the paper.

The validation includes rigorous discussions of the stability conditions of the new 4th-order schemes. In

particular, in Sec. 7.2, we present a stability study of a nonlinear scheme based on a rather exhaustive

analysis of the propagation of round-off errors. As a result of this study, the reader can understand why the

underlying idea is applicable even in a nonlinear case.

To complement the new 4th-order schemes, a set of simple 4th-order non-reflecting boundary conditions

for linear schemes is presented in Sec. 5.5. Its potency is validated in Sec. 9. The extensions of these

boundary conditions for nonlinear schemes are under development.

By introducing new independent mesh variables (e.g., those appear in Eqs. (8.33)–(8.35)) and the

concept of dependent mesh variables (e.g., those defined in Eq. (8.36)), it is shown in Sec. 8 that the

framework described in Sec. 5–7 can be extended to construct multidimensional high order schemes.

Accuracy and stability of the new 4th-order CESE schemes are rigorously assessed using the numerical

results presented in Sec. 9. It is shown that the theoretical predictions made in Sec. 5–7 are consistent with

the numerical results.

Even though they are not presented in the paper, from the details provided here, in a straightforward

manner, the reader can build from a given 2nd-order CESE scheme its 6th-, 8th-,... order versions which

have the same stencil and same stability conditions of the 2nd-order scheme.

Finally note that the current approach can also be used to construct high-order schemes from classical

2nd-order schemes. In fact, by using a “regular” space-time mesh in which (i) the set of mesh points used

include all (j, n) with j, n = 0, ±1, ±2,..., and (ii) xj = j∆x and t

n = n∆t, a 3rd-order hybrid scheme has

been built from the Lax-Wendroff and the second-order a scheme. To construct this scheme, Eq. (5.4) is

replaced by

u∗(x, t; j, n) = un

j + (ux)

n

j

$

(x − xj ) − a (t − t

n)

%

+

1

2

(uxx)

n

j

$

(x − xj ) − a (t − t

n)

%2 (10.1)

Also, to account for the new mesh structure, we have

ν def= a

∆t

∆x

, (ux ̄)

n

j

def

= ∆x

2 (ux)

n

j , and (ux ̄x ̄)

n

j

def

=

8 ∆x

2

92

(uxx)

n

j , j, n = 0 ± 1, ±2,... (10.2)

Then the new hybrid scheme is formed by

(ux ̄x ̄)

n

j = ν(ν + 1)

2 (ux ̄x ̄)

n−1

j−1 + (1 − ν2)(ux ̄x ̄)

n−1

j + ν(ν − 1)

2 (ux ̄x ̄)

n−1

j+1 (10.3)

57
(ux ̄)

n

j = 1

2

$

(s+)

n−1

j+1 − (s−)

n−1

j−1

%

− 2ν

3 (ux ̄x ̄)

n

j (10.4)

and

un

j = 1

2

$

(1 − ν)(s+)

n−1

j+1 + (1 + ν)(s−)

n−1

j−1

%

− 2

3

(ux ̄x ̄)

n

j (10.5)

where

(s+)

n

j

def

=

/

u − (1 + ν)ux ̄ + 2(1 + ν + ν2)

3

ux ̄x ̄

0n

j

j, n = 0, ±1, ±2,... (10.6)

and

(s−)

n

j

def

=

/

u + (1 − ν)ux ̄ + 2(1 − ν + ν2)

3

ux ̄x ̄

0n

j

j, n = 0, ±1, ±2,... (10.7)

As expected, the hybrid scheme is third order in accuracy and is stable if ν2 < 1.

References

1. S.C. Chang and W.M. To, A New Numerical Framework for Solving Conservation Laws–The Method

of Space-Time Conservation Element and Solution Element, NASA TM 104495, August 1991.

2. S.C. Chang, On An Origin of Numerical Diffusion: Violation of Invariance under Space-Time Inversion,

in Proceedings, 23rd Conference on Modeling and simulation, April 30-May 1, 1992, Pittsburgh, PA,

USA, edited by W.G. Volt and M.H. Mickle, Part 5, p. 2727. Also published as NASA TM 105776.

3. S.C. Chang and W.M. To, A Brief Description of a New Numerical Framework for Solving Conservation Laws—The Method of Space-Time Conservation Element and Solution Element, in Proceedings

of the Thirteenth International Conference on Numerical Methods in Fluid Dynamics, Rome, Italy,

1992, edited by M. Napolitano and F. Sabetta, Lecture Notes in Physics 414, (Springer-Verlag, New

York/Berlin, 1992), p. 396.

4. S.C. Chang, New Developments in the Method of Space-Time Conservation Element and Solution

Element—Applications to the Euler and Navier-Stokes Equations, NASA TM 106226.

5. S.C. 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, 295 (1995).

6. S.T. Yu and S.C. Chang, Treatments of Stiff Source Terms in Conservation Laws by the Method of

Space-Time Conservation Element and Solution Element, AIAA Paper 97-0435 (1997).

7. S.C. Chang, S.T. Yu, A. Himansu, X.Y. Wang, C.Y. Chow, and C.Y. Loh, The method of space-time

conservation element and solution element—A new paradigm for numerical solution of conservation

laws, in Computational Fluid Dynamics Review 1998 edited by M.M. Hafez and K. Oshima (World

Scientific, Singapore), Vol. 1, p. 206.

8. T. Molls and F. Molls, Space-Time Conservation Method Applied to Saint Venant Equations, J. of

Hydraulic Engr., 124(5), 501 (1998).

9. S.T. Yu, S.C. Chang, P.C.E. Jorgenson, S.J. Park and M.C. Lai, “Treating Stiff Source Terms in Conservation Laws by the Space-Time Conservation Element and Solution Element Method,” in Proceedings

of the 16th International Conference on Numerical Method in Fluid Dynamics, Arcachon, France, 6-10

July, 1998, edited by C.H. Bruneau, (Springer-Verlag Berlin Heidelberg 1998), p. 433.

10. S.C. Chang, X.Y. Wang, and C.Y. Chow, The Space-Time Conservation Element and Solution Element

Method—A New High-Resolution and Genuinely Multidimensional Paradigm for Solving Conservation

Laws. I. The Two-Dimensional Time Marching Schemes, NASA TM—1998-208843, December, 1998.

11. S.C. Chang, X.Y. Wang, and C.Y. Chow, The Space-Time Conservation Element and Solution Element

Method—A New High-Resolution and Genuinely Multidimensional Paradigm for Solving Conservation

Laws. II. Numerical Simulation of Shock Waves and Contact Discontinuities, NASA TM—1998-208844,

December, 1998.

12. S.C. Chang, X.Y. Wang, and C.Y. 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, 89 (1999).

58
13. X.Y. Wang, and S.C. Chang, A 2D non-splitting unstructured triangular mesh Euler solver based on the

space-time conservation element and solution element method, Computational Fluid Dynamics Journal,

8(2), 309 (1999).

14. C. Zoppou and S. Roberts, Space-Time Conservation Method Applied to Saint Venant Equations: A

Discussion, J. of Hydraulic Engr., 125(8), 891 (1999).

15. N.S. Liu and K.H. Chen, Flux: An Alternative Flow Solver for the National Combustion Code, AIAA

Paper 99-1079.

16. X.Y. Wang and S.C. Chang, A 3D structured/unstructured Euler solver based on the space-time conservation element and solution element method, in A Collection of Technical Papers, 14th AIAA CFD

Conference, June 28–July 1, 1999, Norfolk, Virginia, AIAA Paper 99-3278.

17. S.C. Chang, X.Y. Wang and W.M. To, Application of the space-time conservation element and solution

element method to one-dimensional convection-diffusion problems, J. Comput. Phys., 165, 189 (2000).

18. J. Qin, S.T. Yu, Z.C. Zhang, and M.C. Lai, Direct Calculations of Cavitating Flows by the Space-Time

CESE Method, J. Fuels & Lubricants, SAE Transc., 108(4), 1720 (2000).

19. C.Y. Loh, L.S. Hultgren and S.C. Chang, Wave computation in compressible flow using the space-time

conservation element and solution element method, AIAA J., 39(5), 794 (2001).

20. Z.C. Zhang, S.T. Yu, and S.C. 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, 168 (2002).

21. K.B.M.Q. Zaman, M.D. Dahl, T.J. Bencic, and C.Y. Loh, Investigation of A ‘Transonic Resonance’

with Convergent-Divergent Nozzles, J. Fluid Mech., 463, 313 (2002).

22. C.Y. Loh and K.B.M.Q. Zaman, Numerical Investigation of ‘Transonic Resonance’ with A ConvergentDivergent Nozzle, AIAA J., 40(12), 2393 (2002).

23. S. Motz, A. Mitrovic, and E.-D. Gilles, Comparison of Numerical Methods for the Simulation of Dispersed Phase Systems, Chemical Engineering Science, 57, 4329 (2002).

24. S. Cioc and T.G. Keith, Application of the CESE Method to One-Dimensional Flow in Fluid Film

Bearings, STLE Tribology Transactions, 45, 167 (2002).

25. S.C. Chang, Courant Number Insensitive CESE Schemes, AIAA Paper 2002-3890 (2002).

26. S. Cioc and T.G. Keith, Application of the CESE Method to Two-Dimensional Flow in Fluid Film

Bearings, International J. of Numer. Methods for Heat & Fluid Flow, 13(2), 216 (2003).

27. S. Cioc, F. Dimofte, T.G. Keith, and D.P. Fleming, Computation of Pressurized Gas Bearings Using

the CESE Method, STLE Tribology Transactions, 46(1), 128 (2003).

28. S. Cioc, F. Dimofte, and T.G. Keith, Application of the CESE Method to Wave Journal Bearings, STLE

Tribology Transactions, 46(2), 179 (2003).

29. A. Ayasoufi and T.G. Keith, Application of the Conservation Element and Solution Element Method

in Numerical Modeling of Heat Conduction with Melting and/or Freezing, International J. of Numer.

Methods for Heat & Fluid Flow, 13(4), 448 (2003).

30. I.S. Chang, Unsteady Underexpanded Jet Flows, AIAA Paper 2003-3885.

31. S.C. Chang and X.Y. Wang, Multidimensional Courant Number Insensitive CESE Euler Solvers for

Applications Involving Highly Nonuniform Meshes, AIAA Paper 2003-5280.

32. A. Ayasoufi and T.G. Keith, Application of the Conservation Element and Solution Element Method

in Numerical Modeling of Axisymmetric Heat Conduction with Melting and/or Freezing, JSME International J. Series B, 47(1), 115 (2004).

33. A. Ayasoufi and T.G. Keith, Application of the Conservation Element and Solution Element Method in

Numerical Modeling of Three-dimensional Heat Conduction with Melting and/or Freezing, Transactions

of the ASME, J. of Heat Transfer, 126(6), 937 (2004).

34. Y.I. Lim, S.C. Chang, and S.B. Jorgensen, A Novel Partial Differential Algebraic Equation (PDAE)

Solver: Iterative Space-Time Conservation Element/Solution Element (CESE) Method, Computers and

Chemical Engineering, 28, 1309 (2004)

35. Y.I. Lim and S.B. Jorgensen, A Fast and Accurate Numerical Method for Solving Simulated Moving

Bed (SMB) Chromatographic Separation Problems, Chemical Engineering Science, 59, 1931 (2004).

59
36. Y.I. Lim, An Optimization Strategy for Nonlinear Simulated Moving Bed Chromatography: Multi-level

Optimization Procedure (MLOP), Korea J. Chem. Eng., 21(4), 836 (2004).

37. C.K. Kim, S.T. John Yu, and Z.C. Zhang, Cavity Flow in Scramjet Engine by the Space-Time Conservation Element and Solution Element Method, AIAA J., 42(5), 912 (2004).

38. M. Zhang, S.T. John Yu, S.C. Lin, S.C. Chang, and I. Blankson, Solving Magnetohydrodynamic Equations Without Special Treatment for Divergence-Free Magnetic Field, AIAA J., 42(12), 2605 (2004).

39. K.S. Im, M.C. Lai, S.T. John Yu, and Robert R. Matheson, Jr., Simulation of Spray Transfer Process

in Electrostatic Rotary Bell Sprayer, ASME J. of Fluid Engineering, 126(3), 449 (2004).

40. S. Jerez, J.V. Romero, and M.D. Rosello, A Semi-Implicit Space-Time CESE Method to Improve Mass

Conservation through Tapered Ducts in Internal Combustion Engines, Math. and Computer Modeling,

40, 941 (2004).

41. B.S. Venkatachari, G.C. Cheng, and S.C. Chang, Development of A Transient Viscous Flow Solver

Based on Conservation Element-Solution Element Framework, AIAA Paper 2004-3413.

42. S.C. Chang, Y. Wu, V. Yang, and X.Y. Wang, Local Time Stepping Procedures for the Space-Time

Conservation Element and Solution Element Method, International J. Comput. Fluid Dynamics, 19(5),

359 (2005).

43. Y.I. Lim, S.B. Jorgensen, and I.H. Kim, Computer-Aided Model Analysis for Ionic Strength-Dependent

Effective Charge of Protein in Ion-Exchange Chromatography, Bio-Chem. Eng. J., 25(2), 125 (2005).

44. T.I. Tseng and R.J. Yang, Simulation of the Mach Reflection in Supersonic Flows by the CESE Method,

Shock Waves, 14(4), 307 (2005).

45. B. Wang, H. He, and S.-T.J. Yu, Direct Calculation of Wave Implosion for Detonation Initiation, AIAA

J., 43(10), 2157 (2005).

46. S.C. Chang, Explicit von Neumann Stability Conditions for the c-τ Scheme—A Basic Scheme in the

Development of the CESE Courant Number Insensitive Schemes, NASA TM 2005-213627, April 2005.

47. J.C. Yen and D.A. Wagner, Computational Aeroacoustics Using a Simplified Courant Number Insensitive CESE Method, AIAA Paper 2005-2820.

48. I.S. Chang, C.L. Chang, and S.C. Chang, Unsteady Navier-Stokes Rocket Nozzle Flows, AIAA Paper

2005-4353.

49. S.C. Chang, Courant Number and Mach Number Insensitive CESE Euler Solvers, AIAA Paper 2005-

4355

50. S.C. Chang, A. Himansu, C.Y. Loh, X.Y. Wang, and S.T. Yu, Robust and Simple Non-Reflecting Boundary Conditions for the Euler Equations–A New Approach Based on the Space-Time CESE Method, in

Proceedings, NSF-CBMS Regional Research Conference on Mathematical Methods in Nonlinear Wave

Propagation, North Carolina A&T State University, Greensboro, North Carolina, May 15-19, 2002,

edited by D.P. Clemence and G. Tang, p. 155-190, Vol. 379 in Contemporary Mathematics, American

Mathematical Society (2005). Also published as NASA/TM-2003-212495/Rev1.

51. B.S. Venkatachari, G.C. Cheng, and S.C. Chang, Courant Number Insensitive Transient Viscous Flow

Solver Based on CESE Framework, AIAA Paper 2005-00931.

52. T.I. Tseng and R.J. Yang, Numerical Simulation of Vorticity Production in Shock Diffraction, AIAA

J., 44(5), 1040 (2006).

53. M. Zhang, S.-T. Yu, S.C.H. Lin, S.C. Chang, and I. Blankson, Solving the MHD Equations By the

Space-Time Conservation Element and Solution Element Method, J. Comput. Phys., 214, 599 (2006).

54. I.S. Chang, C.L. Chang, and S.C. Chang, 3D Unsteady Navier-Stokes Rocket Nozzle Flows, AIAA Paper

2006-4775.

55. C.L. Chang, Time-accurate, Unstructured-Mesh Navier-Stokes Computations with the Space-time CESE

Method, AIAA Paper 2006-4780.

56. S.C. Chang, On Space-Time Inversion Invariance and Its Relation to Non-Dissipativeness of a CESE

Core Scheme, AIAA Paper 2006-4779. Also published as NASA TM 2006-214385.

57. X. Feng, Y. Zhou, and S.T. Wu, A Novel Numerical Implementation for Solar Wind Modeling by the

Modified Conservation Element/Solution Element Method, Astrophys. J., 655(1), 1110 (2007).

60
58. C.L. Chang, Three-Dimensional Navier-Stokes Calculations Using the Modified Space-Time CESE

Method, AIAA Paper 2007-5818.

59. S.C. Chang, The a(3) Scheme—A Fourth-Order Neutrally Stable CESE Solver, AIAA Paper 2007-4321.

60. S.C. Chang, The a(4) Scheme—A High-Order Neutrally Stable CESE Solver, AIAA Paper 2007-5820.

61. S.C. Chang, The a(3) Scheme—A Fourth-Order Space-Time Flux-Conserving and Neutrally Atable

CESE Solver (A revised and expanded version of AIAA Paper 2007-4321), NASA TM 2008-215138,

April 2008.

62. S.C. Chang, Courant Number Insensitive Schemes and New Classes of Wiggle-Suppressing WeightedAveraging Techniques, in preparation.

63. Other CESE references are posted on: http://www.grc.nasa.gov/www/microbus.

64. B. van Leer, Toward the Ultimate Conservative Difference Scheme. IV. A New Approach to Numerical

Convection, J. Comput. Phys., 23, 276 (1977).

65. G.D. van Albada, B. van Leer, and W.W. Robert, A Comparative Study of Computational Methods in

Cosmic Gas Dynamics, Astronom. and Astrophys., 108, 76 (1982).

66. R.A. Horn and C.R. Johnson, Matrix Analysis, Cambridge University Press, 1985.

61
Table 1. Numerical results of the a -ε -4 scheme with periodic boundary conditions

K=25, n=135 K=50, n=270 K=100, n=540

E 1.07E-15 1.38E-15 1.43E-15

E X 2.91E-05 1.81E-06 1.13E-07

E xx 4.21E-14 5.49E-14 5.66E-14

E xxx 1.09 0.272 6.79E-02

Table 2. Numerical results of the a -ε -4 scheme with periodic boundary conditions

K=25, n=135 K=50, n=270 K=100, n=540

E 1.07E-15 1.38E-15 1.43E-15

E X 2.48E-05 1.54E-06 9.62E-08

E xx 4.21E-14 5.49E-14 5.66E-14

E xxx 0.929 0.231 5.77E-02

Table 3. Numerical results of the a -ε -4 scheme with periodic boundary conditions

K=25, n=270 K=50, n=540 K=100, n=1080

E 1.59E-03 9.95E-05 6.21E-06

E X 1.00E-02 6.27E-04 3.91E-05

E xx 1.86 0.467 0.117

E xxx 11.7 2.94 0.735

ν = 1.0, ε 1

 = ε 2 = 0, t = 5.4

ν = 1.0, ε 1 = ε 2 = 0.5, t = 5.4

ν = 0.5, ε 1

 = ε 2 = 0.5, t = 5.4

Table 4. Numerical results of the a -ε -4 scheme with periodic boundary conditions

K=25, n=2700 K=50, n=5400 K=100, n=10800

E 1.14E-02 3.30E-04 1.35E-05

E X 7.30E-02 2.14E-03 8.69E-05

E xx 4.83 0.842 0.171

E xxx 30.6 5.37 1.09

Table 5. Numerical results of the c -τ-4 scheme with periodic boundary conditions

K=25, n=135 K=50, n=270 K=100, n=540

E 1.07E-15 1.38E-15 1.43E-15

E X 2.48E-05 1.54E-06 9.62E-08

E xx 4.21E-14 5.49E-14 5.66E-14

E xxx 0.929 0.231 5.77E-02

ν = 0.05, ε 1 = ε 2 = 0.5, t = 5.4

ν = τ 1

 = τ 2 = 1.0, t = 5.4
Table 6. Numerical results of the c -τ-4 scheme with periodic boundary conditions

K=25, n=270 K=50, n=540 K=100, n=1080

E 1.05E-04 6.37E-06 3.94E-07

E X 5.58E-04 1.08E-04 1.59E-05

E xx 4.76E-01 1.18E-01 2.93E-02

E xxx 0.366 0.64 0.505

Table 7. Numerical results of the c-τ-4 scheme with periodic boundary conditions

K=25, n=2700 K=50, n=5400 K=100, n=10800

E 5.10E-04 3.17E-05 1.98E-06

E X 3.67E-03 2.61E-04 2.02E-05

E xx 1.06 0.265 6.61E-02

E xxx 7.16 1.92 0.546

Table 8. Numerical results of the a -ε -4 scheme with non-reflecting boundary conditions

K=25, n=135 K=50, n=270 K=100, n=540

E 2.18E-15 2.34E-15 2.48E-15

E X 4.21E-05 2.67E-06 1.67E-07

E xx 8.57E-14 9.21E-14 9.79E-14

E xxx 1.58 0.4 0.1

ν = τ 1

 = τ 2 = 0.5, t = 5.4

ν = τ 1 = τ 2 = 0.05, t = 5.4

ν = 1.0, ε 1

 = ε 2 = 0, t = 5.4

Table 9. Numerical results of the a -ε -4 scheme with non-reflecting boundary conditions

K=25, n=135 K=50, n=270 K=100, n=540

E 2.18E-15 2.34E-15 2.48E-15

E X 2.46E-05 1.54E-06 9.60E-08

E xx 8.57E-14 9.21E-14 9.79E-14

E xxx 0.923 0.23 5.76E-02

Table 10. Numerical results of the a -ε -4 scheme with non-reflecting boundary conditions

K=25, n=270 K=50, n=540 K=100, n=1080

E 2.34E-05 1.30E-06 7.76E-08

E X 5.30E-04 4.98E-05 4.58E-06

E xx 0.227 5.63E-02 1.40E-02

E xxx 1.44 0.507 0.18

ν = 1.0, ε 1 = ε 2 = 0.5, t = 5.4

ν = ε 1

 = ε 2 = 0.5, t = 5.4
Table 11. Numerical results of the a -ε -4 scheme with non-reflecting boundary conditions

K=25, n=2700 K=50, n=5400 K=100, n=10800

E 2.99E-04 8.25E-06 2.76E-07

E X 3.21E-03 1.06E-04 7.88E-06

E xx 0.522 9.20E-02 1.98E-02

E xxx 5.28 1.59 0.586

Table 12. Numerical results of the c -τ-4 scheme with non-reflecting boundary conditions

K=25, n=135 K=50, n=270 K=100, n=540

E 2.18E-15 2.34E-15 2.48E-15

E X 2.46E-05 1.56E-06 9.60E-08

E xx 8.57E-14 9.21E-14 9.80E-14

E xxx 0.923 0.23 5.76E-02

Table 13. Numerical results of the c-τ-4 scheme with non-reflecting boundary conditions

K=25, n=270 K=50, n=540 K=100, n=1080

E 2.64E-06 1.49E-07 8.84E-09

E X 2.35E-04 2.66E-05 3.09E-06

E xx 6.53E-02 1.57E-02 3.84E-03

E xxx 2.94 1.45 0.709

ν = τ 1 = τ 2 = 0.5, t = 5.4

ν = 0.05, ε 1

 = ε 2 = 0.5, t = 5.4

ν = τ 1 = τ 2 = 1.0, t = 5.4

Table 14. Numerical results of the c-τ-4 scheme with non-reflecting boundary conditions

K=25, n=2700 K=50, n=5400 K=100, n=10800

E 9.11E-06 5.46E-07 3.27E-08

E X 8.25E-04 6.75E-05 5.69E-06

E xx 0.125 3.12E-02 7.83E-03

E xxx 2.86 1.07 0.393

ν = τ 1 = τ 2 = 0.05, t = 5.4
t

∆x/2 ∆x/2

∆t/2

n = 0

n = 1/2

n = 1

n = 3/2

n = 2

j = 0

j = 1/2

j = 1

j = 3/2

∆t/2

x

j = –3/2

j = –1

j = –1/2

(j,n)

A

A

B

B

C

C

D

A

D

D

E

E

F

F

(j-1/2,n-1/2) (j+1/2,n-1/2)

(j,n)

(j–1/2,n–1/2) (j+1/2,n–1/2)

(j,n) (j,n)

2(e) CE(j,n)

2(c) CE–(j,n). 2(d) CE+(j,n).

Figure 2.—The SEs and CEs.

2(b) SE(j,n).

∆t

∆t

∆x ∆x

2(a)—The space-time mesh.

Figure 1.—A surface element on the boundary

 S(V) of an arbitrary space-time volume V.

x

ds

→

→

→

dr

t

V

S(V)

r = (x,t)
Figure 3.—Definition of P—

, M—, M+, and P+. 

Figure 4.—A Burgers solution. 

-0.2

0

0.2

0.4

0.6

0.8

1

1.2

-0.6 -0.4 -0.2 0 0.2 0.4 0.6

u

x

Numerical

Analytical
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-08-2010 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

A New Approach for Constructing Highly Stable High Order CESE Schemes 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Chang, Sin-Chung 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 599489.02.07.03.04.03.01 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-17391 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2010-216766 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Category: 64 

Available electronically at http://gltrs.grc.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

A new approach is devised to construct high order CESE schemes which would avoid the common shortcomings of traditional high order 

schemes including: (a) susceptibility to computational instabilities; (b) computational inefficiency due to their local implicit nature (i.e., at 

each mesh points, need to solve a system of linear/nonlinear equations involving all the mesh variables associated with this mesh point); (c) 

use of large and elaborate stencils which complicates boundary treatments and also makes efficient parallel computing much harder; (d) 

difficulties in applications involving complex geometries; and (e) use of problem-specific techniques which are needed to overcome stability 

problems but often cause undesirable side effects. In fact it will be shown that, with the aid of a conceptual leap, one can build from a given 

2nd-order CESE scheme its 4th-, 6th-, 8th-,... order versions which have the same stencil and same stability conditions of the 2nd-order 

scheme, and also retain all other advantages of the latter scheme. A sketch of multidimensional extensions will also be provided. 

15. SUBJECT TERMS

CESE method; High-order schemes 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

72 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18