H.T. Huynh

Glenn Research Center, Cleveland, Ohio

On Formulations of Discontinuous Galerkin and

Related Methods for Conservation Laws

NASA/TM—2014-218135

June 2014
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 to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 443–757–5803

• Phone the NASA STI Information Desk at

 443–757–5802

• Write to:

 STI Information Desk

 NASA Center for AeroSpace Information

 7115 Standard Drive

 Hanover, MD 21076–1320
H.T. Huynh

Glenn Research Center, Cleveland, Ohio

On Formulations of Discontinuous Galerkin and

Related Methods for Conservation Laws

NASA/TM—2014-218135

June 2014

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

International Conference on Spectral and High-Order Methods (ICOSAHOM)

sponsored by the University of Utah and Arizona State University

Salt Lake City, Utah, June 23–27, 2014
Acknowledgments

The author was supported by the Aeronautical Sciences Project of NASA.

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://www.sti.nasa.gov

Level of Review: This material has been technically reviewed by technical management. 
NASA/TM—2014-218135 1 

On Formulations of Discontinuous Galerkin and 

Related Methods for Conservation Laws 

H.T. Huynh 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

huynh@grc.nasa.gov

Abstract 

A formulation for the discontinuous Galerkin (DG) method that leads to solutions using the 

differential form of the equation (as opposed to the standard integral form) is presented. The formulation 

includes (a) a derivative calculation that involves only data within each cell with no data interaction 

among cells, and (b) for each cell, corrections to this derivative that deal with the jumps in fluxes at the 

cell boundaries and allow data across cells to interact. The derivative with no interaction is obtained by a 

projection, but for nodal-type methods, evaluating this derivative by interpolation at the nodal points is 

more economical. The corrections are derived using the approximate (Dirac) delta functions. The 

formulation results in a family of schemes: different approximate delta functions give rise to different 

methods. It is shown that the current formulation is essentially equivalent to the flux reconstruction (FR) 

formulation. Due to the use of approximate delta functions, an energy stability proof simpler than that of 

Vincent, Castonguay, and Jameson (2011) for a family of schemes is derived. Accuracy and stability of 

resulting schemes are discussed via Fourier analyses. Similar to FR, the current formulation provides a 

unifying framework for high-order methods by recovering the DG, spectral difference (SD), and spectral 

volume (SV) schemes. It also yields stable, accurate, and economical methods. 

1. Introduction 

For conservation laws, e.g., those in the field of Computational Fluid Dynamics (CFD), low-order 

methods are generally robust and reliable; as a result, they are routinely employed in practical 

calculations. For the same computing cost, high-order methods can provide considerably more accurate 

solutions, but they are more complicated and less robust. The need to improve and develop new highorder methods with favorable properties has attracted the interest of many researchers. 

The discontinuous Galerkin (DG) method is currently among the most widely used high-order 

numerical methods for solving conservation laws on unstructured meshes. It was introduced for the 

neutron transport equation by Reed and Hill (1973), analyzed by LaSaint and Raviart (1974) and 

developed and made popular for fluid dynamics equations by Cockburn, Shu, Bassi, Rebay, and others 

(Cockburn, Karniadakis, and Shu 2000, Bassi and Rebay 1997a,b, 2000, Cockburn and Shu 2005, 2009, 

Shu 2012, and the references therein). Efficient DG schemes storing data at nodal points known as nodal 

DG methods were elaborated by Hesthaven and Warburton (2008). 

Alternative approaches to high-order accuracy employing the differential form (as opposed to DG 

which employs the integral form) have been proposed. Kopriva and Kolias (1996) pioneered this 

approach with the staggered-grid spectral method on quadrilateral meshes. The approach was extended to 

triangular meshes by Liu, Vinokur, and Wang (2006) and named spectral difference (SD), and solutions 

for a wide range of problems were obtained by Wang et al. (2007) and Liang et al. (2009a and 2009b). 

Another class of schemes called spectral volume (SV) based on the idea of subdividing each cell into 

subcells or control volumes in a structured manner was proposed by Wang, Zhang, and Liu (2004).
NASA/TM—2014-218135 2 

Recently, an approach to high-order accuracy with the advantage of simplicity and economy called 

flux reconstruction (FR) was introduced by the author (Huynh 2007, 2009). The approach amounts to 

evaluating the derivative of a discontinuous piecewise polynomial function by employing its 

straightforward derivative estimate together with a correction that accounts for the jumps at the interfaces. 

The FR framework unifies several existing schemes: with appropriate choices of correction terms, it 

recovers DG, SD, as well as SV, and the FR versions are generally simpler and more economical than the 

original versions. In addition, the approach results in numerous new methods that are stable and super 

accurate, i.e., more accurate than expected (also known as super convergent). Extensions of FR to 

unstructured triangular meshes was provided by Wang and Gao (2009), to the 2D Navier-Stokes 

equations on meshes of mixed elements by Gao and Wang (2009) and Gao et al. (2013), to the 3D Euler 

and Navier-Stokes equations on mixed meshes by Haga et al. (2010, 2011) and Wang et al. (2011), and to 

dynamic meshes by Yu et al. (2012). Whereas FR approach is also known as CPR (Correction Procedure 

via Reconstruction), due to the usage in references closely related to this work, the name FR is employed. 

A mathematical foundation for the FR approach was recently provided by Jameson (2010), who 

proved that a particular SD scheme (recovered via FR) is energy-stable for 1D linear advection. Vincent, 

Castonguay, and Jameson (2011a) subsequently extended this result, and proved that a one-parameter 

family of FR methods is energy-stable for linear advection. This family, referred to as Energy Stable Flux 

Reconstruction (ESFR), was extended to linear advection on 2D triangular grids by Castonguay, Vincent 

and Jameson (2012), to linear advection-diffusion in 1D by Castonguay et al. (2013), and to linear 

advection-diffusion on 2D triangular grids by Williams et al. (2011) and Williams et al. (2013). Von 

Neumann analyses for these schemes were investigated by Vincent, Castonguay and Jameson (2011b) and 

nonlinear stability was discussed by Jameson, Vincent and Castonguay (2012). The ESFR methods were 

shown to be equivalent a filtered DG scheme by Allaneau and Jameson (2012). 

A review of recent developments for the FR methods was presented by Huynh, Wang, and Vincent 

(2014). An open-source Python based framework for solving advection-diffusion type problems on 

streaming architectures using the FR approach can be found at www.pyfr.org/. 

In this paper, a formulation for the discontinuous Galerkin (DG) methods that leads to solutions using 

the differential form of the equation is presented. The formulation includes (a) a derivative calculation 

that involves only data within each cell with no data interaction among cells and (b) for each cell, 

corrections to this derivative that deal with the jumps in fluxes at the cell boundaries and allow data 

across cells to interact. The derivative with no interaction is obtained by a projection, but for nodal-type 

methods, evaluating of this derivative by interpolation at the nodal points is more economical. The 

corrections are derived using the approximate (Dirac) delta functions. The formulation results in a family 

of schemes: different approximate delta functions give rise to different methods. It is shown that the 

current formulation is essentially equivalent to the flux reconstruction (FR) formulation. Due to the use of 

approximate delta functions, an energy stability proof simpler than that of Vincent, Castonguay, and 

Jameson (2011) for a family of schemes is derived. Accuracy and stability of resulting schemes are 

discussed via Fourier analyses. 

The formulations to be discussed are basic in nature. Therefore, in the hope of attracting the interest of 

researchers who are not familiar with (some of) these approaches, this paper is written in an essentially 

self-contained manner. It is organized as follows. The DG method is reviewed in Section 2. Section 3 

introduces new strong forms for DG via the use of the approximate Dirac delta functions. The FR 

schemes, obtained by integrating the strong forms of Section 3, are discussed in Section 4. An energy 

stability proof for a family of schemes is presented in Section 5. Section 6 contains representative FR 

schemes and their properties obtained by Fourier analysis. Conclusions and discussion are presented in 

Section 7. A brief discussion of Fourier stability and accuracy analyses can be found in the Appendix. 

 
NASA/TM—2014-218135 3 

2. Discontinuous Galerkin Formulations 

2.1. Preliminaries. Consider the conservation law 

ݑ௧ ൅ ݂௫ ൌ 0 (2.1)

with initial condition ݑሺݔ, 0ሻ ൌ ݑ୧୬୧୲ሺݔሻ. The solution ݑ is assumed to be periodic or of compact support 

so that boundary conditions are trivial. The flux ݂ depends on ݑ. With ܽሺݑሻ ൌ ݂݀/݀ݑ, the above can be 

cast in nonconservation form: 

ݑ௧ ൅ ܽݑ௫ ൌ 0. (2.2)

Let the domain of calculation Ω be divided into possibly nonuniform cells or elements 	ܧ௝,	 ݆ ൌ 1, 2, ...

Denote the center of ܧ ௝by ݔ ௝and its width by ݄௝. Instead of dealing with the global elements, it is more 

convenient to deal with the local (or reference) element ܫ ൌ ሾെ1, 1ሿ. With ߦ varying on ܫ and ݔ on 	ܧ௝, the 

linear function mapping ܫ onto ܧ ௝and its inverse are 

ݔሺߦሻ ൌ ݔ ൅ ௝ߦ݄௝/2				 and				ߦሺݔሻ ൌ 2ሺݔ െ ݔ௝ሻ/݄௝ . (2.3)

A function ݎ௝ሺݔሻ on ܧ ௝results in a function on ܫ denoted by ݎ௝ሺߦሻ for simplicity of notation, where 

ݎ௝ሺߦሻ ൌ ݎ௝൫ݔሺߦሻ൯.

The derivatives in the global and local descriptions are related by the chain rule 

ሻݔሺ௝ݎ݀

݀ݔ ൌ ቆ 2

݄௝

ቇ

ሻߦሺ௝ݎ݀

(2.4) . ߦ݀

For simplicity, when the index is not essential, ܧ ௝is denoted by ܧ. Let ݒ and ݓ be two functions on ܧ ;

their inner product is given by 

ሺݒ, ݓሻா ൌ නݒሺݔሻݓሺݔሻ ݀ݔ

ா

.

The ܮଶ norm of ݒ on ܧ is (a different norm will be discussed later) 

‖ݒ‖ா ൌ ቆන൫ݒሺݔሻ൯

ݔ݀ ଶ

ா

ቇ

ଵ/ଶ

.

The inner product in the global description relates to that in the local description by, 

ሺݒ, ݓሻா ൌ නݒሺݔሻݓሺݔሻ ݀ݔ

ா

ൌ ቆ݄௝

2 ቇ නݒሺߦሻݓሺߦሻ ݀ߦ

ூ

ൌ ቆ݄௝

2 ቇ ሺݒ, ݓሻூ. (2.5)

For any nonnegative integer ݉, let ࡼ ௠be the space of polynomials of degree ݉ or less. The 

convention ࡼ ௠ൌ ሼ0ሽ for ݉൏0 is also used. 

Let the Legendre polynomial ܮ ௠be defined as the unique polynomial of degree ݉ that is orthogonal 

to ࡼ௠ିଵ and ܮ௠ሺ1ሻ ൌ 1. Then, it is well-known (Hildebrandt 1987, Hesthaven and Warburton 2008) that 
NASA/TM—2014-218135 4 

‖ܮ‖௠ଶ ൌ ሺܮ௠, ܮ௠ሻ ൌ

2

2݉ ൅ 1 . (2.6)

2.2. Discontinuous Galerkin formulations. Since there are several subtle differences among the 

various standard and new formulations, for ease of comparison, first, the weak and strong forms of the 

DG method are re-derived. 

Let ݇ be a nonnegative integer. At time ݐ, let the solution ݑሺݔ, ݐሻ be approximated on each cell ܧ ௝by a 

polynomial of degree ݇ in ݔ denoted by ݑ௝ሺݔ, ݐሻ. The collection of all ሼݑ௝ሺ. , ݐሻሽ as ݆ varies forms a 

function denoted by ݑ௛, which is generally discontinuous across cell interfaces. In the reference 

description for ܧൌܧ௝, with ݐ understood, ݑ ௝can be expressed in modal form as 

ሻߦሺ௜ܮ ௜ ,௝ݑ෍ൌ ሻߦሺ௝ݑ

௞

௜ୀ଴

. (2.7)

Assume that the data ݑ௝,௜ ሺݐ௡ሻ ൌ ݑ௝, ௜ are known for all cells. We wish to calculate ௗ

ௗ௧

 .௡ݐൌݐ at ௜ ,௝ݑ

Note first that since ݑ ௛is generally discontinuous across cell interfaces, ݂ሺݑ௛ሻ defined by ݂ሺݑ௛ሻሺݔሻ ൌ

݂ሺݑ௛ሺݔሻሻ also shares this property. For the data among cells to interact and, simultaneously, the method 

to be conservative, we need to define a unique flux for each interface that is common for the two adjacent 

cells. Also note that, on each cell, ݂ሺݑ௛ሻ can be a non-polynomial function and needs to be approximated 

as will be discussed later. 

Focusing on the cell ܧൌܧ௝, let ߶ be a test function, i.e., an element of ࡼ) ௞independent of 	ݐ(. The 

conservation law (2.1) leads to the following requirement 

߲

߲ݐ ሺݑ௛,߶ ሻா ൅ ቀ൫݂ሺݑ௛ሻ൯

௫, ߶ቁ

ா

ൌ 0. (2.8)

The above relation on ܧ, however, involves no interaction with the data of the neighboring cells. To 

account for data interaction, first, by integration by parts: 

߲

߲ݐ ሺݑ௛,߶ ሻா ൅ ሺ݂ሺݑ௛ሻ߶ሻడா െ ሺ݂ሺݑ௛ሻ, ߶௫ሻா ൌ 0.

For the boundary term ሺ݂ሺݑ௛ሻ߶ሻడா, at each interface ݆ ൅ 1/2, from the left and right ݑ values, namely, 

ݑି ൌ ݑ௝ାଵ/ଶ

ି ൌ ݑ௝ሺݔ௝ାଵ/ଶሻ and ݑା ൌ ݑ௝ାଵ/ଶ

ା ൌ ݑ௝ାଵሺݔ௝ାଵ/ଶሻ, (2.9)

a flux value 

݂௝ାଵ/ଶ

ூ ൌ ݂ூ

(2.10) ାሻݑ ି,ݑሺ

common for the two adjacent cells ܧ ௝and ܧ௝ାଵ is employed (the superscript ܫ represents ‘interface’). This 

common flux, also called numerical flux, is typically an upwind-biased flux for advection-type problems. 

In the case of an advection with constant speed ܽ, 

 
NASA/TM—2014-218135 5 

݂ூ ൌ 1

2 ܽሺݑି ൅ ݑାሻ െ 1

2 ߚ|ܽ|ሺݑା െ ݑିሻ (2.11)

where 0൑ߚ൑1) with ߚൌ0 recovering a centered scheme and ߚൌ1 recovering an upwind scheme). 

For a diffusion-type problem, the common flux is typically a centered quantity. As for ߶, the value from 

ܧ) no upwinding) is used. For the DG method, the solution is required to satisfy 

߲

߲ݐ ሺݑ௛,߶ ሻா ൅ ሺ݂ூ

߶ሻడா െ ሺ݂ሺݑ௛ሻ, ߶௫ሻா ൌ 0. (2.12)

The above is often referred to as the weak form. Denoting the two boundaries of ܧ by the indices ܮ and ܴ, 

߲

߲ݐ ሺݑ௛,߶ ሻா ൅ ݂ோ

ூ ߶ோ െ ݂௅ூ ߶௅ െ ሺ݂ሺݑ௛ሻ, ߶௫ሻா ൌ 0. (2.13)

Let the two boundary fluxes using the data inside the cell ܧൌܧ ௝be 

݂௅

ା ൌ ݂ሺݑ௝ሺݔ௝ିଵ/ଶሻሻ and ݂ோ

ି ൌ ݂ሺݑ௝ሺݔ௝ାଵ/ଶሻሻ. (2.14)

Again at the two boundaries of ܧ, set 

ሾ݂ሿ௅ ൌ ݂௅

ூ െ ݂௅

ା and ሾ݂ሿோ ൌ ݂ோ

ூ െ ݂ோି. (2.15)

That is, ሾ݂ሿ represents the jump or difference between the common interface flux and the flux value just 

inside the cell (and ሾ݂ሿ௅ is not ݂௅

ା െ ݂௅

ି). By applying integration by parts again to (2.12), 

߲

߲ݐ ሺݑ௛,߶ ሻா ൅ ቀ൫݂ሺݑ௛ሻ൯

௫, ߶ቁ

ா

൅ ሺሾ݂ሿ߶ሻడா ൌ 0. (2.16)

Comparing the above with (2.8), the difference is the boundary term ሺሾ݂ሿ߶ሻడா, which accounts for data 

interaction. Equivalently, 

߲

߲ݐ ሺݑ௛,߶ ሻா 	൅ ቀ൫݂ሺݑ௛ሻ൯

௫, ߶ቁ

ா

൅ ሾ݂ሿோ߶ோ െ ሾ݂ሿ௅߶௅ ൌ 0. (2.17)

Note that the weak and strong forms (2.13) and (2.17) are equivalent if the inner products are carried 

out exactly. In practice, however, the inner products are approximated by quadrature rules, and the two 

forms generally produce different results. 

3. New Strong Forms for the DG Method 

The next goal is to rewrite the strong form (2.17) in a manner involving no test function. The 

formulation below deals with ݂௫ directly and thus differs from the flux reconstruction (FR) formulation, 

which first reconstructs the flux ݂. As will be shown, the two formulations yield essentially the same 

schemes with FR typically formulated with the nodal form. The approximate Dirac delta functions below 

also simplify the energy-stability proof in Section 5. 

3.1. Approximate Dirac delta functions. These functions are motivated by the boundary terms in 

(2.17). From here on, unless otherwise stated, inner products take place in the reference frame, i.e., on 

ܫ ൌ ሾെ1,1ሿ. For a fixed ߙ on ܫ, let the approximate Dirac delta function at ߙ to degree ݇ be denoted by 

ߜఈ, ௞ and defined as the linear functional on ࡼ ௞that satisfies, for all ߶ in ࡼ௞, 
NASA/TM—2014-218135 6 

ߜఈ, ௞ሺ߶ሻ ൌ ߶ሺߙሻ. (3.1)

By a theorem in functional analysis, there exists a member ߛఈ, ௞ in ࡼ ௞such that for all ߶ in ࡼ௞, 

൫ߛఈ, ௞, ߶൯ ൌ ߜఈ, ௞ሺ߶ሻ ൌ ߶ሺߙሻ. (3.2)

In fact, ߛఈ, ௞ can be calculated. First, since it is of degree ݇, it can be expressed as 

ߛఈ, ௞ሺߦሻ ൌ෍ܾ௜ ܮ௜

௞

௜ୀ଴

ሺߦሻ. (3.3)

For any ݉൑݇, to obtain ܾ௠, apply the dot product with ܮ ௠to the above and, using the fact that ܮ ௠is 

orthogonal to all ܮ௜, ്݅݉, 

൫ߛఈ, ௞, ܮ௠൯ൌ ܾ௠ሺܮ௠, ܮ௠ሻ.

By (3.2), ൫ߛఈ, ௞, ܮ௠൯ൌܮ௠ሺߙሻ. This fact together with (2.6) imply 

ܮ௠ሺߙሻ ൌ 2

2݉ ൅ 1 ܾ௠.

That is, ܾ௠ ൌ ଶ௠ାଵ

ଶ

ܮ௠ሺߙሻ. Consequently, 

ߛఈ, ௞ሺߦሻ ൌ ෍2݅ ൅ 1

௜ܮ ሻߙሺ௜ܮ 2

௞

௜ୀ଴

ሺߦሻ. (3.4)

See e.g., (Morse and Feshbach 1953, p. 719) where an expression for the Dirac delta function in terms of 

general eigenvectors can be found. 

For simplicity, ߜఈ, ௞ is considered identical to ߛఈ, 

ߜఈ ൌ ߜఈ, ௞ ൌ ߛఈ, ௞ ൌ ߛఈ.

Note that for ߙ ൌ േ1, |ܮ௜ሺߙሻ| ൌ 1; therefore, ܮ௜ሺߙሻଶ ൌ 1, and 

ߜଵ, ௞ሺ1ሻ ൌ ߜିଵ, ௞ሺെ1ሻ ൌ ෍2݅ ൅ 1

2

௞

௜ୀ଴

.

As a result, 

ߜଵ, ௞ሺ1ሻ ൌ ߜିଵ, ௞ሺെ1ሻ ൌ ሺ݇൅1ሻଶ

2 . (3.5)

 
NASA/TM—2014-218135 7 

The above is the maximum value of ߜേଵ on ܫ and is related to the fact that in conjunction with an explicit 

time-stepping method, the CFL condition for the DG method is inversely proportional to ሺ݇൅1ሻଶ. In 

addition, the ݅-th mode of ߜേଵ satisfies, by (2.6), 

ฯ

2݅ ൅ 1

2 ܮ௜ሺേ1ሻ ܮ௜ฯ

ଶ

ൌ 2݅ ൅ 1

2 .

Therefore, after some algebra, 

ฮߜേଵฮ ൌ ฮߜേଵ, ௞ฮ ൌ ݇൅1

√2 . (3.6)

The graphs of ߜ଴,௞ , ݇ ൌ 0, 2, 4, 6, and 8 are shown in Figure 3.1(a). Note that since ܮଶ௞ାଵሺ0ሻ ൌ 0, 

ߜ଴, ଶ௞ାଵ ൌ ߜ଴, ଶ௞. For comparison, the graphs of ߜଵ, ௞, ݇ ൌ 4, ... , 8 are shown in Fig. 3.1(b). Note the 

difference in scales. 

3.2. Strong form S1. We can now derive the new strong forms. Again unless otherwise stated, we use 

the local description, and the indices ܮ and ܴ represent െ1 and 1, respectively. By (2.4) and (2.5), the 

strong form (2.17) implies 

ቆ

݄௝

2 ቇ ሺݑ௛,߶ ሻ௧ ൅ ቀ൫݂ሺݑ௛ሻ൯

క

, ߶ቁ ൅ ሾ݂ሿோ߶ோ െ ሾ݂ሿ௅߶௅ ൌ 0. (3.7)

Using (3.4), set 

1 ൅ 2݅ ෍ ൌ௅ ߜ

2 ሺെ1ሻ௜ ܮ௜

௞

௠ୀ଴

ሺߦሻ and ߜோ ൌ ෍ 2݅ ൅ 1

௜ܮ 2

௞

௠ୀ଴

ሺߦሻ. (3.8a,b)

(a) (b) 

Figure 3.1.—(a) Graphs of ߜ଴,௞ , ݇ ൌ 0, 2, 4, 6, and 8. Since ܮଶ௞ାଵሺ0ሻ ൌ 0, ߜ଴, ଶ௞ାଵ ൌ ߜ଴, ଶ௞. (b) Graphs of ߜଵ, ௞ , 

݇ ൌ 4, 5, ... , 8. 

 
NASA/TM—2014-218135 8 

Then ߜ ௅and ߜோ are of degree ݇ and, for all ߶ in ࡼ௞, 

ሺߜ௅,߶ ሻ ൌ ߶௅ and ሺߜோ, ߶ሻ ൌ ߶ோ. (3.9)

Thus, (3.7) and the above imply 

ቆ

݄௝

2 ቇ ሺݑ௛,߶ ሻ௧ ൅ ቀ൫݂ሺݑ௛ሻ൯

క

, ߶ቁ ൅ ሾ݂ሿோሺߜோ, ߶ሻ െ ሾ݂ሿ௅ሺߜ௅,߶ ሻ ൌ 0. (3.10)

For any integer ݉൒0, denote by ࣪ ௠the projection onto ࡼ௠, i.e., with any integrable function ݒ, 

࣪௠ሺݒሻ ൌ ෍ሺݒ, ܮ௜ሻ

௜ܮ ଶ௜‖ܮ‖

௠

௜ୀ଴

ൌ ෍2݅ ൅ 1

௜ܮ ሻ௜ܮ ,ݒሺ 2

௠

௜ୀ଴

. (3.11)

Due to the fact that ߶ is in ࡼ௞,) 3.10 (is equivalent to 

ቆ

݄௝

2 ቇ ሺݑ௛,߶ ሻ௧ ൅ ൫࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ, ߶൯ ൅ ሾ݂ሿோሺߜோ, ߶ሻ െ ሾ݂ሿ௅ሺߜ௅,߶ ሻ ൌ 0. (3.12)

The projection ࣪ ௞above serves the purpose of cancelling out the test function ߶. Indeed, since ݑ௛, 

࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ, ߜோ, and ߜ ௅are all in ࡼ ௞and, since ߶ spans ࡼ௞, the above implies the following strong 

form referred to as S1, 

S1 ቆ

݄௝

2 ቇ ሺݑ௛ሻ௧ ൅ ࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ ൅ ሾ݂ሿோߜோ െ ሾ݂ሿ௅ߜ ௅ൌ 0. (3.13)

Two remarks concerning the form S1 are in order. 

1. Generally, after integrating by parts twice, we get back the original equation. Here, due to the 

different flux values employed for the boundary terms in each of the two integrations by parts, the strong 

form S1 looks similar to but not identical with the original conservation laws in the local description. The 

differences are: (a) the derivative with no interaction is provided by ࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ and not just ሺ݂ሺݑ௛ሻሻక

and (b) the correction for the derivative, namely, ሾ݂ሿோߜோ െ 	ሾ݂ሿ௅ߜ௅, provides data interaction among cells 

and, at the same time, deals with the flux jumps at the interfaces. 

2. For numerical algorithms, the interaction of the data among cells involves the functions ߜ ௅and ߜோ, 

which are derived exactly and are given by (3.8). Concerning ࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ, which represents the 

derivative with no interaction, ݑ ௛is in ࡼ௞, but ݂ሺݑ௛ሻ needs not be a polynomial in general. For nodaltype schemes, ࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ can be approximated conveniently and effectively by using the chain rule at 

the nodal points ߦ௜, ݅ ൌ 0, ... , ݇: ൫݂ሺݑ௛ሻ൯

క

ሺߦ௜ሻ ൌ ܽ൫ݑ௛ሺߦ௜ሻ൯	ሺݑ௛ሻకሺߦ௜ሻ. 

3. The strong form (3.13) leads to a family of schemes by altering the criteria of approximation for the 

Dirac delta function. What is crucial is that the resulting scheme must be stable and accurate. 

3.3. Strong form S2. We next present an alternative to the strong form S1. Concerning again the 

derivative with no interaction, we can first approximate ݂ሺݑ௛ሻ by its projection ࣪௞ሺ݂ሺݑ௛ሻሻ and then 

calculate the derivative: the result is referred to as S2, 
NASA/TM—2014-218135 9 

S2 ቆ

݄௝

2 ቇ ሺݑ௛ሻ௧ ൅ ሾ࣪௞ሺ݂ሺݑ௛ሻሿక ൅ ሾ݂ሿோߜோ െ ሾ݂ሿ௅ߜ ௅ൌ 0. (3.14)

For this strong form, the two boundary fluxes interior to cell ܧൌܧ ௝are defined as, instead of (2.14), 

݂௅

ା ൌ ሾ࣪௞ሺ݂ሺݑ௝ሻሻሿሺݔ௝ିଵ/ଶሻ and ݂ோ

ି ൌ ሾ࣪௞ሺ݂ሺݑ௝ሻሻሿሺݔ௝ାଵ/ଶሻ (3.15)

and the jumps ሾ݂ሿ at the two interfaces are defined using the above interior values. 

The following remarks are in order. 

1. For this strong form, similar to the case of S1, concerning the derivative with no interaction 

ሾ࣪௞ሺ݂ሺݑ௛ሻሿక, with nodal-type schemes, we can first approximate ݂ሺݑ௛ሻ by its values at the nodal points 

and then evaluate the derivative of the resulting polynomial. Such an approximation is said to be via 

‘Lagrange polynomial’ whereas that in the second remark after (3.13) is via ‘chain rule’. These two ways 

of approximating the derivative with no interaction were mentioned by the author as alternative 

approaches (Huynh 2007). Later, Wang and Gao (2009) found that for nonlinear problems such as the 

Euler equations, the solution via the chain rule was more accurate. 

2. With the current strong forms S1 and S2, the conservation law for nodal-type methods is solved via 

a combination of interpolation and projection: the derivative with no interaction is approximated via 

interpolation, and the two boundary corrections are carried out exactly to degree ݇ via projection since the 

approximate Dirac delta functions ߜേଵ are exact up to degree ݇. The evaluation of ߜേଵሺߦሻ at the nodes 

incurs no additional (interpolation) errors. Also note that the approximation to the derivative with no 

interaction ൫݂ሺݑ௛ሻ൯

క

 can be carried out via projection as well but, typically, projection for such a function 

is costlier than interpolation; moreover, projection is typically approximated by a quadrature which, 

similar to interpolation, incurs errors. 

3. For the case of advection with constant speed, the strong forms S1 and S2 become identical and is 

called strong form S. The reason is that since ݂ሺݑ௝ሻ ൌ ܽݑ ௝is of degree ݇, both approaches for the 

derivative with no interaction result in ܽሺݑ௝ሻక, which is of degree ݇െ1. In the general case of nonlinear 

conservation laws, however, the two forms are different: the polynomial ሾ࣪௞ሺ݂ሺݑ௛ሻሿక for S2 is of degree 

݇െ1 as opposed to ࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ for S1, which is of degree ݇. The lower degree approximation of S2 

appears to be consistent with the numerical findings of Wang and Gao (2009) mentioned above. 

4. The key difference between the current strong forms and a standard DG strong form such as (2.17) 

is that by employing the approximate Dirac delta functions, the current ones involve no test functions. As 

will be shown, they also unify existing methods by recovering the spectral difference schemes via an 

appropriate modification to the approximate Dirac delta function; in addition, these strong forms result in 

a family of energy-stable methods that have the property of super accuracy (i.e., more accurate than 

expected via Fourier analysis as discussed by Huynh (2007) and Vincent et al. (2011)). 

3.4. Conservation property. Next, we show that the strong form S1 is conservative, i.e., 

߲

ݔ݀ ሻݔሺ௛ݑන ݐ߲

ா

ൌ ݂௅

୍ െ ݂ோ୍

. (3.16)

Note first that for any integrable function ݓ on ܫ, the projection preserves the cell average quantity, 
NASA/TM—2014-218135 10 

න ࣪௞ሺݓሻ

ଵ

ିଵ

ݓ න ൌ ߦ݀

ଵ

ିଵ

(3.17) .ߦ݀

Therefore, with ݒక ൌ ݓ, 

න ࣪௞൫ݒక൯

ଵ

ିଵ

݀ߦ ൌ න ݒక

ଵ

ିଵ

݀ߦ ൌ ݒሺ1ሻ െ ݒሺെ1ሻ. (3.18)

Applying the above to ࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ, we obtain 

න ࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ

ଵ

ିଵ

݀ߦ ൌ න ሺ݂ሺݑ௛ሻሻక

ଵ

ିଵ

݀ߦ ൌ ݂ோ

ି െ ݂௅ା (3.19)

where ݂௅

ା and ݂ோ

ି are the boundary fluxes interior to ܧ as defined in (2.14). Next, for the cell ܧ௝, in the 

local description, (3.12) with ߶ൌ1 implies 

ቆ

݄௝

2 ቇ ߲

௝ݑ න ݐ߲

ଵ

ିଵ

݀ߦ ൅ න ࣪௞ሾሺ݂൫ݑ௝൯ሻకሿ

ଵ

ିଵ

݀ߦ ൅ ሾ݂ሿோ െ ሾ݂ሿ௅ ൌ 0. (3.20)

By (3.23) and (2.15), Eq. (3.20) follows. This completes the proof. 

The conservative property of the scheme S2 is similar with (3.15) replacing (2.14) and, therefore, 

omitted. 

A modification to assure conservation property when using the derivative with no interaction via the 

chain rule on 2D unstructured meshes was provided by Gao and Wang (2013). 

4. DG Scheme Via Flux Reconstruction (FR) Approach 

We now present a formulation which deals with ݂ rather than ݂௫ as in the previous section. It will be 

shown that for the DG method, the slope ݂௫ for the conservation law is given by ሺܨ௝ሻ௫ where, on each cell 

ܧ௝, ܨ ௝is a polynomial of degree ݇൅1 defined by (a) ܨ ௝takes on the common interface flux values at the 

two cell boundaries and (b) ࣪௞ିଵ൫ܨ௝൯ൌ࣪௞ିଵ൫݂ሺݑ௛ሻ൯. The key idea here for constructing ܨ ௝is to 

integrate the strong form S1. 

4.1. Approximating the flux ݂ሺݑ௛ሻ with no interaction. In the form S1, the derivative with no 

interaction is approximated by ࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ. This approximation is equivalent to a certain approximation 

for the flux ݂ itself. In what manner is ݂ approximated? 

To answer the above question, first, by the fundamental theorem of Calculus, 

൫݂ሺݑ௛ሻ൯ሺߟሻ ൌ ݂௅

ା ൅ න ൫݂ሺݑ௛ሻ൯

క

ఎ

ିଵ

(4.1) .ߦ݀

Next, with subscript IPD short for ‘integrating the projection of the derivative’, set 

୍݂୔ୈሺߟሻ ൌ ݂௅

ା ൅ න ࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ

ఎ

ିଵ

(4.2) .ߦ݀
NASA/TM—2014-218135 11 

Then ୍݂୔ୈ is of degree ݇൅1 and, the second equation below is due to (3.19), 

୍݂୔ୈሺെ1ሻ ൌ ݂௅

ା and ୍݂୔ୈሺ1ሻ ൌ ݂ோ

ି. (4.3)

That is, ୍݂୔ୈ interpolates ݂ሺݑ௛ሻ at the two cell boundaries. Two conditions for ୍݂୔ୈ are known, ݇

conditions remain to describe how it approximates ݂ሺݑ௛ሻ. 

To obtain these remaining conditions, let ߶ be in ࡼ௞. Integration by parts yields 

൫୍݂୔ୈ, 	߶క൯ ൌ ݂ோ

ି߶ோ െ ݂௅

ା߶௅ െ ൫ሺ୍݂୔ୈሻక, ߶൯. (4.4)

By (4.2), ሺ୍݂୔ୈሻక ൌ ࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ. Thus, the above implies 

൫୍݂୔ୈ, 	߶క൯ ൌ ݂ோ

ି߶ோ െ ݂௅

ା߶௅ െ ൫࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ, ߶൯. (4.5)

Next, in (4.5), since ߶ is in ࡼ௞, we can replace ࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ by ሺ݂ሺݑ௛ሻሻక, and since ߶క is in ࡼ௞ିଵ, we 

can replace ୍݂୔ୈ by ࣪௞ିଵሺ୍݂୔ୈሻ, 

൫࣪௞ିଵሺ୍݂୔ୈሻ, 	߶క൯ ൌ ݂ோ

ି߶ோ െ ݂௅

ା߶௅ െ ሺሺ݂ሺݑ௛ሻሻక, ߶ሻ. (4.6)

By applying integration by parts to the inner product of the right hand side, 

൫࣪௞ିଵሺ୍݂୔ୈሻ, ߶క൯ ൌ ሺሺ݂ሺݑ௛ሻሻ, ߶కሻ. (4.7)

As ߶ spans ࡼ௞,߶ క spans ࡼ௞ିଵ; therefore, the above implies 

࣪௞ିଵሺ୍݂୔ୈሻ ൌ ࣪௞ିଵሺ݂ሺݑ௛ሻሻ. (4.8)

Thus, in the strong form S1, ࣪௞ሾ൫݂ሺݑ௛ሻ൯

క

ሿ approximates ሺ݂ሺݑ௛ሻሻక in the sense of projection by its 

definition. Equivalently, its integral ୍݂୔ୈ defined by (4.2), which is a polynomial of degree ݇൅1, 

approximates ݂ሺݑ௛ሻ in the sense of a combination of interpolation at the two boundaries as in (4.3), and 

projection as above. 

4.2. Correction functions. We now deal with the boundary terms of the form S1 or (3.13). With ߜோ by 

(3.8b), let the correction function for the right boundary be denoted by ݃ோ, ஽ீ or simply ݃ோ and defined by 

݃ோሺߦሻ ൌ න ߜோ

క

ିଵ

(4.9) .ߟ݀ሻߟሺ

(In Section 6, ݃ோ denotes the correction function for a more general scheme, not just DG.) Then, ݃ோ is of 

degree ݇൅1 and 

݃ோ′ൌߜோ. (4.10)

Moreover, 

݃ோሺെ1ሻ ൌ 0 and ݃ோሺ1ሻ ൌ 1, (4.11a,b)
NASA/TM—2014-218135 12 

where the second equation results from (4.9) with ߦൌ1 together with (3.2) where ߶ൌ1. Note that since 

ߜோ approximates the Dirac delta function at ߦൌ1, ݃ோ approximates the Heaviside or unit step (step up) 

function at the same location. 

Next, it is shown that 

࣪௞ିଵሺ݃ோሻ ൌ 0. (4.12)

Indeed, by integration by parts, for any polynomial ݒ of degree ݇, 

ሺ݃ோ, ݒ′ሻ ൌ ݒሺ1ሻ െ න ݃ோ

ᇱ ሺߦሻݒሺߦሻ݀ߦ

ଵ

ିଵ

ൌ ݒሺ1ሻ െ න ߜோሺߦሻݒሺߦሻ݀ߦ

ଵ

ିଵ

ൌ 0. (4.13)

As ݒ spans ࡼ௞, ݒ′ spans ࡼ௞ିଵ; therefore, (4.12) holds. 

Let the correction function for the left boundary be defined by (note the change in the integration 

limits) 

݃௅ሺߦሻ ൌ න ߜ௅

ଵ

క

(4.14) .ߟ݀ሻߟሺ

(Again a more precise notation is ݃௅, ஽ீ; however, the subscript DG is understood in this section.) Then, 

݃௅ is of degree ݇൅1 and 

݃௅

ᇱ ൌ െߜ௅.) 4.15(

Moreover, 

݃௅ሺെ1ሻ ൌ 1 and ݃௅ሺ1ሻ ൌ 0, (4.16a,b)

where the first equation results from (4.14) with ߦ ൌ െ1 together with (3.2) where ߶ൌ1. Note that since 

ߜ ௅approximates the Dirac delta function at ߦ ൌ െ1, ݃௅ approximates the step down function which takes 

on the value 1 at ߦ ൌ െ1 and the value 0 for െ1 ൏ ߦ ൑ 1. 

In addition, 

࣪௞ିଵሺ݃௅ሻ ൌ 0. (4.17)

The proof is similar to that of (4.12) and is omitted. 

4.3. Flux reconstruction. The reconstruction of the flux on ܧൌܧ ௝can now be carried out. At the two 

interfaces, recall that ݂௅

ା ൌ ݂ሺݑ௝ሺെ1ሻሻ, ݂ோ

ି ൌ ݂ሺݑ௝ሺ1ሻሻ, ሾ݂ሿ௅ ൌ ݂௅ூ െ ݂௅

ା, and ሾ݂ሿோ ൌ ݂ோ

ூ െ ݂ோି. 

From Subsection 4.1 above, ୍݂୔ୈ is of degree ݇൅1, involves no interaction among cells, takes on the 

values ݂௅

ା and ݂ோ

ି at the left and right interfaces, respectively, and ࣪௞ିଵሺ୍݂୔ୈሻ ൌ ࣪௞ିଵሺ݂ሺݑ௛ሻሻ. Next, the 

function 

୍݂୔ୈ ൅ ሾ݂ሿ௅ ݃௅
NASA/TM—2014-218135 13 

recovers the common flux value at the left interface, namely ݂௝ିଵ/ଶ

ூ , due to the condition ݃௅ሺെ1ሻ ൌ 1. It’s 

value at the right interface is unchanged from that of ୍݂୔ୈ, namely the value ݂ோ

ି; this is due to the 

condition ݃௅ሺ1ሻ ൌ 0. The function 

ܨ ௝ൌ ୍݂୔ୈ ൅ ሾ݂ሿ௅ ݃௅ ൅ ሾ݂ሿோ ݃ோ (4.18)

takes on the common interface flux values at the two boundaries, 

ܨ௝ሺെ1ሻ ൌ ݂௅

ூ

 and ܨ௝ሺ1ሻ ൌ ݂ோ

ூ

. (4.19)

In addition, on ܧ௝, by (4.8), (4.12), and (4.17), 

࣪௞ିଵ൫ܨ௝൯ ൌ ࣪௞ିଵ൫݂ሺݑ௝ሻ൯. (4.20)

Finally, by (4.2), (4.10), and (4.15), Eq. (4.18) implies 

ܨ௝′ൌ࣪௞ሾሺ݂ሺݑ௛ሻሻకሿ െ ሾ݂ሿ௅ߜ ൅ ௅ሾ݂ሿோߜோ. (4.21)

Note that the above right hand side is identical to that of the strong form S1 of (3.13). This completes the 

derivation of DG schemes via FR approach. 

A few remarks are in order. 

1. Due to (4.19), the collection of all ሼܨ௝ሽ as ݆ varies forms a function which is continuous across cell 

interfaces, and ܨ௝′ yields the DG approximation for ݂௫ for the conservation law. 

2. The flux reconstruction form (4.18) above leads to a family of schemes by using different criteria to 

define the correction function ݃௅ and ݃ோ. What is crucial is that the criteria must result in stable and 

accurate methods. 

3. In the weak form (2.13) for the DG scheme, ݂ሺݑ௛ሻ can be replaced by ࣪௞ିଵሺ݂ሺݑ௛ሻሻ since ߶௫ is of 

degree ݇െ1, 

߲

߲ݐ ሺݑ௛,߶ ሻா ൅ ݂ோ

ூ ߶ோ െ ݂௅ூ ߶௅ െ ሺ࣪௞ିଵሺ݂ሺݑ௛ሻሻ, ߶௫ሻா ൌ 0. (4.22)

Thus, due to (4.19) and (4.20), the definition for ܨ ௝by (4.18) is, loosely put, consistent with the weak 

form above. 

5. Energy Stable Flux Reconstruction Schemes (ESFR) 

Energy stability for a family of schemes can now be shown. Consider the case of advection with 

constant speed ܽ൒0 (the wind blows from left to right). We also assume that the initial data is periodic 

or of compact support so that boundary conditions are trivial. 

5.1. Energy stability of the DG scheme. First, it is well-known that the DG scheme with the upwind 

biased flux (2.11) is energy stable. The simple proof, which can be found in (Jameson 2011), plays a 

critical role below; therefore, it is included. 

The strong form (2.17) on ܧൌܧ ௝with ߶ൌݑ ௝implies 
NASA/TM—2014-218135 14 

൬

߲

൰௝ݑ ,௝ݑ ݐ߲

ா

	൅ ቄቀܽ	൫ݑ௝൯

௫

ቁ௝ݑ ,

ா

൅ ሺ݂௝ାଵ/ଶ

ூ െ ݂௝ାଵ/ଶି ሻݑ௝ାଵ/ଶି െ ሺ݂௝ିଵ/ଶூ െ ݂௝ିଵ/ଶ

ା ሻݑ௝ିଵ/ଶା ቅ ൌ 0. (5.1)

Note that ൫ܽሺݑ௝ሻక, ݑ௝൯ ൌ ௔

ଶ ቂ൫ݑ௝ାଵ/ଶ

ି ൯

ଶ െ ൫ݑ௝ିଵ/ଶ

ା ൯

ଶ

ቃ. With ݂௝ାଵ/ଶ

ூ ൌ ܽݑ௝ାଵ/ଶூ for each interface, the 

expression in the curly bracket above can be written as 

ܽሺݑ௝ାଵ/ଶ

ூ െ ଵ

ଶ

ݑ௝ାଵ/ଶ

ି ሻݑ௝ାଵ/ଶି െ ܽሺݑ௝ିଵ/ଶூ െ ଵ

ଶ

ݑ௝ିଵ/ଶ

ା ሻݑ௝ିଵ/ଶା . (5.2)

At each interface ݆ െ 1/2, collecting the contributions from the element ܧ௝ିଵ and then ܧ௝, there is a total 

contribution of 

ܽሺݑ௝ିଵ/ଶ

ூ െ ଵ

ଶ

ݑ௝ିଵ/ଶ

ି ሻݑ௝ିଵ/ଶି െ ܽሺݑ௝ିଵ/ଶூ െ ଵ

ଶ

ݑ௝ିଵ/ଶ

ା ሻݑ௝ିଵ/ଶା . (5.3)

Omitting the subscript ݆ െ 1/2, we obtain ܽሾሺݑூ െ ଵ

ଶ

ݑିሻݑି െ ሺݑூ െ ଵ

ଶ

ݑାሻݑାሿ. After some algebra, and 

using the upwind-biased flux (2.11), the result is 

ଵ

ଶ ߚ|ܽ|ሺݑା െ ݑିሻଶ. (5.4)

Recall that 0൑ߚ൑1. Thus, on the whole physical domain, with our periodic or of compact support 

assumption for the data, (5.1) and the above imply 

1

2

߲

௛ݑ	 න ݐ߲

ݔ݀  ଶ

ஐ

	ൌ െ෍1

2 ߚ|ܽ|൫ݑ௝ିଵ/ଶ

ା െ ݑ௝ିଵ/ଶ

ି ൯

ଶ

௝

൑ 0. (5.5)

That is, energy of the numerical solution can only decrease in time. This completes the proof that the DG 

scheme in semi-discrete form using the upwind-biased common flux is energy stable for the linear 

advection equation. 

What we need from the above proof is the fact that the quantity in the square bracket of (5.1) implies 

(5.4) at each interface after the square bracket quantities from the two adjacent cells are collected. 

5.2. Energy stability of a family of FR scheme. For the current case of advection, as discussed in the 

third remark after (3.15), the two strong forms S1 and S2 become identical and are referred to as strong 

form S. On ܧ௝, set 

݂௝ ൌ ݂൫ݑ௝൯ ൌ ܽݑ௝.) 5.6(

To obtain a family of schemes, we replace ߜ in the strong form S of (3.13) by ߜ൅߳, where ߳ provides the 

variation: 

ቆ

݄௝

2 ቇ ሺݑ௝ሻ௧ ൅	ሺ݂௝ሻక െ ሾ݂ሿ௅ሺߜ ൅ ௅߳௅ሻ ൅ ሾ݂ሿோሺߜோ ൅ ߳ோሻ ൌ 0. (5.7)

In other words, we approximate the Dirac delta function not by ߜ but by ߜ൅߳, which results in a family 

of schemes. What is crucial is that ߳ must be chosen in a manner that the resulting scheme is stable and 

accurate. In this section, we assume ߳ involves only the highest mode: 
NASA/TM—2014-218135 15 

߳௅ ൌ ሺെ1ሻ௞߳௞ܮ ௞and ߳ோ ൌ ߳௞ܮ௞.) 5.8a,b)

The coefficient ሺെ1ሻ௞ above is for being consistent with (3.8a). We wish to show that if 

߳௞ ൐ െ

2݇ ൅ 1

2 , (5.9)

then the resulting scheme is energy-stable under a certain energy norm to be defined later. 

Compared with the coefficient of the highest mode of ߜோ which, by (3.8b), equals ଶ௞ାଵ

ଶ , the above 

condition is equivalent to the fact that the highest mode of ߜோ ൅ ߳ோ has a strictly positive coefficient, i.e., 

1 ൅ 2݅෍ ൅ ௞ܮ ௞ߙ ൌ ோ߳ ൅ ோߜ

௜ܮ 2

௞ିଵ

௜ୀ଴

(5.10) . ଵି ௞,ோߜ ൅ ௞ܮ ௞ߙ ൌ

where 

1 ൅ 2݇ ൌ௞ ߙ

2 ൅ ߳௞ ൐ 0. (5.11)

For the left boundary, 

(5.12) . ଵି ௞,௅ߜ ൅ ௞ܮ ௞ߙ௞ሻ1ሺെ ൌ௅ ߳௅ ൅ ߜ

The following proof is similar to but simpler than that of Vincent et al. (2011) due to the use of the 

approximate Dirac delta functions. 

Using technique employed in (5.1), multiplying (5.7) by ݑ ௝and integrating on ܫ, we obtain 

݄௝

4 		

݀

௝ݑ න ݐ݀

ଶ

ଵ

ିଵ

݀ߦ ൌ െ න ݑ௝൫݂௝൯ߦ݀ క

ଵ

ିଵ

൅ ሾ݂ሿ௅ݑ௝, ௅ െ ሾ݂ሿோݑ௝, ோ ൅

																					 ሾ݂ሿ௅ න ݑ௝߳ ௅݀ߦ

ଵ

ିଵ

െ ሾ݂ሿோ න ݑ௝߳ோ ݀ߦ

ଵ

ିଵ

.

(5.13)

Set 

ܴ ൌ	ሾ݂ሿ௅ න ݑ௝߳ ௅݀ߦ

ଵ

ିଵ

െ ሾ݂ሿோ න ݑ௝߳ோ ݀ߦ

ଵ

ିଵ

. (5.14)

Eq. (5.13) implies 

݄௝

4 		

݀

௝ݑ න ݐ݀

ଶ

ଵ

ିଵ ݀ߦ 	െ 	ܴ ൌ െ ቊන ݑ௝൫݂௝൯

ߦ݀ క

ଵ

ିଵ

െ ሾ݂ሿ௅ݑ௝, ൅ ௅ ሾ݂ሿோݑ௝, ோቋ. (5.15)

The square bracket quantity above equals the square bracket quantity in (5.1), which will result in (5.4) at 

the interfaces. 
NASA/TM—2014-218135 16 

We now derive conditions on ߳௞ so that the above left hand side results in the time rate of change of a 

certain energy norm. Our goal is to relate ܴ of (5.14), which does not involve ௗ

ௗ௧, with ௗௗ௧ ൫	ݑ௝,௞ 

ଶ ሺݐሻ൯ where 

௞ ,௝ݑ

ଶ is, up to a constant factor, the energy of the ݇-th mode of ݑ௝. 

To this end, first, on ܫ, 

ฮݑ௝ฮ

ଶ ൌ න ݑ௝

ଶ

ଵ

ିଵ

௜ ,௝ݑ ෍ ൌ ߦ݀

ଶ ሺܮ௜, ܮ௜ሻ

௞

௜ୀ଴

ൌ ෍ 2

௜ ,௝ݑ 1 ൅ 2݅

ଶ

௞

௜ୀ଴

. (5.16)

Next, since ߳௅ and ߳ோ involve only the highest mode as expressed in (5.8), and since the ݇-th mode of 

 implies) 5.14 (.Eq ,௞ܮ௞ ,௝ݑ is௝ ݑ

ܴ ൌ 		 ߳௞ ݑ௝,‖ ௞ ܮ‖௞ଶ ሼ ሺെ1ሻ௞ሾ݂ሿ௅ െ ሾ݂ሿோ ሽ. (5.17)

Following (Jameson 2010) and (Vincent et al. 2011), differentiate (5.7) ݇ times (in space), 

݄௝

2 		

߲

ݐ߲

௝ݑ௞߲

൅	 ௞ߦ߲

߲௞ାଵ݂௝

߲ߦ௞ାଵ െ ሾ݂ሿ௅

߲௞ሺߜ ൅ ௅߳௅ሻ

߲ߦ ൅ ௞ሾ݂ሿோ

߲௞ሺߜோ ൅ ߳ோሻ

(5.18) 0. ൌ௞ 	 ߦ߲

Since ݂௝ is of degree ݇, 

߲௞ାଵ݂௝

߲ߦ௞ାଵ ൌ 0. (5.19)

Multiplying (5.18) by డೖ௨ೕ

డకೖ , 

݄௝

4 		

߲

௝ݑ௞߲ቆ	ݐ߲

ቇ௞ ߦ߲

ଶ

	െ	ሾ݂ሿ௅

௝ݑ௞߲

௞ߦ߲

߲௞ሺߜ ൅ ௅߳௅ሻ

߲ߦ ൅ ௞ሾ݂ሿோ

௝ݑ௞߲

௞ߦ߲

߲௞ሺߜோ ൅ ߳ோሻ

(5.20) 0. ൌ௞ 	 ߦ߲

Integrating on ܫ, we obtain 

݄௝

4

݀

݀ݐ න ቆ

௝ݑ௞߲

ቇ௞ ߦ߲

ଶ ଵ

ିଵ

(5.21) ෨ܴ ൌ ߦ݀

where 

ܴ෨ ൌ 		ሾ݂ሿ௅ න ߲௞ݑ௝

௞ߦ߲

ଵ

ିଵ

߲௞ሺߜ ൅ ௅߳௅ሻ

߲ߦ ௞݀ߦ െ ሾ݂ሿோ න ߲௞ݑ௝

௞ߦ߲

ଵ

ିଵ

߲௞ሺߜோ ൅ ߳ோሻ

(5.22) .ߦ݀ ௞ߦ߲

For the left hand side of (5.21), denote 

௞ܮ௞߲ ൌ௞ ݉

. ௞ߦ߲

Then, ݉௞ is a constant. In fact, ݉௞ ൌ ݇! ܽ௞ ൌ ሺ2݇ െ 1ሻ!/ሺ2௞ିଵሺ݇െ1ሻ!ሻ where ܽ௞ is the leading 

coefficient of ܮ ;௞however, this expression is not essential. Concerning the left hand side of (5.21), 
NASA/TM—2014-218135 17 

௝ݑ௞߲

௞ ,௝ݑ ൌ௞ ߦ߲

௞ܮ௞߲

(5.23) .௞݉ ௞ ,௝ݑ ൌ௞ ߦ߲

Consequently, recalling that ݑ௝, ௞ ൌ ݑ௝,௞ ሺݐሻ, the left hand side of (5.21) equals ௛ೕ

ଶ

ௗ

ௗ௧ 	൫݉௞

௞ ,௝ݑ	ଶ

ଶ ൯ where we 

use the fact that ׬ 1	݀ߦ ଵ

ିଵ ൌ 2. That is, (5.21) can be written as 

݄௝

4 		

݀

݀ݐ න ቆ

௝ݑ௞߲

ቇ௞ ߦ߲

ଶ ଵ

ିଵ

ൌ ߦ݀

݄௝

2

݀

௞݉൫ ݐ݀

௞ ,௝ݑ ଶ

ଶ ൯ൌܴ෨. (5.24)

We now relate ܴ෨ defined by (5.22) with ܴ of (5.17). The key fact employed is that ܴ෨ involves only the 

highest mode. First, by (3.8a) for ߜ ௅and (5.8a) for ߳௅, 

߲௞ሺߜ ൅ ௅߳௅ሻ

߲ߦ ௞ൌ ሺെ1ሻ௞ ൬

2݇ ൅ 1

2 ൅ ߳௞൰ ݉௞,

and by (3.8b) for ߜோ and (5.8b) for ߳ோ, 

߲௞ሺߜோ ൅ ߳ோሻ

߲ߦ ௞ൌ ൬

2݇ ൅ 1

2 ൅ ߳௞൰ ݉௞.

Due to the above two expressions and (5.23), Eq. (5.22) implies 

ܴ෨ ൌ 	 ሺ2݇ ൅ 1 ൅ 2߳௞ሻ ݉௞

ଶ ݑ௝, ௞ ሼ ሺെ1ሻ௞ሾ݂ሿ௅ െ ሾ݂ሿோ ሽ (5.25)

Using the above and the fact that ‖ܮ‖௞ଶ ൌ 2/ሺ2݇ ൅ 1ሻ, Eq. (5.17) implies 

ܴ ൌ

2߳௞

ሺ2݇ ൅ 1ሻሺ2݇ ൅ 1 ൅ 2߳௞ሻ݉௞

ଶ ܴ෨. (5.26)

Thus, by (5.24), 

ܴ ൌ	

߳௞

ሺ2݇ ൅ 1ሻሺ2݇ ൅ 1 ൅ 2߳௞ሻ ݄௝

௞ ,௝ݑሺ݀

ଶ ሻ

(5.27) . ݐ݀

What is crucial is that ܴ of (5.14), which involves no time derivative, can be expressed as above, which 

involves the time derivative and the highest component ݑ௝,௞. 

By the above, the left hand side of (5.15) yields 

݄௝

4 		

݀

௝ݑ න ݐ݀

ଶ

ଵ

ିଵ

	 ൌ	 ܴ	– 	ߦ݀

݄௝

4

݀

௝ݑ න ቊ ݐ݀

ଶ

ଵ

ିଵ

െ ߦ݀

4߳௞

ሺ2݇ ൅ 1ሻሺ2݇ ൅ 1 ൅ 2߳௞ሻ

௞ ,௝ݑ	

ଶ ቋ	. (5.28)

Using (5.16), the ݇-th mode of the above right hand side equals 

݄௝

4 		

݀

െ 1	൜	൤	ݐ݀

2߳௞

ሺ2݇ ൅ 1 ൅ 2߳௞ሻ

൨ 2

௞ ,௝ݑ 1 ൅ 2݇

ଶ ൠ . (5.29)

For the above to result in a norm, it is required that 
NASA/TM—2014-218135 18 

1 െ

2߳௞

ሺ2݇ ൅ 1 ൅ 2߳௞ሻ ൐ 0 (5.30)

or, ሺ2݇ ൅ 1 ൅ 2߳௞ െ 2߳௞ሻ/ሺ2݇ ൅ 1 ൅ 2߳௞ሻ ൐ 0 or ሺ2݇ ൅ 1ሻ/ሺ2݇ ൅ 1 ൅ 2߳௞ሻ ൐ 0 or ሺ2݇ ൅ 1 ൅ 2߳௞ሻ ൐

0, or 

߳௞ ൐ െ

2݇ ൅ 1

2 . (5.31)

Thus, instead of the standard norm (5.16), we employ a modified energy norm on ܧൌܧ ௝defined by 

ฮݑ௝ฮெ

ଶ ൌ ݄௝

2 ൥൬1	 െ 		

2߳௞

ሺ2݇ ൅ 1 ൅ 2߳௞ሻ

൰ 2

௞ ,௝ݑ 1 ൅ 2݇

ଶ ൅ ෍ 2

௜ ,௝ݑ 1 ൅ 2݅

ଶ

௞ିଵ

௜ୀ଴

	൩	. (5.32)

Or, equivalently, by (5.23), 

ฮݑ௝ฮெ

ଶ

ൌ		ൌ

݄௝

௝ݑ න ൥2

ଶ

ଵ

ିଵ

െ ߦ݀

2߳௞

ሺ2݇ ൅ 1 ൅ 2߳௞ሻሺ2݇ ൅ 1ሻ݉௞

ଶ න ቆ

௝ݑ௞߲

ቇ௞ ߦ߲

ଶ ଵ

ିଵ ݀ߦ൩. (5.33)

Thus, using the above norm, the family of schemes via the approximate Dirac delta function of ߜ൅߳, 

where the highest mode of ߜோ ൅ ߳ோ ൌ ߙ	௞ܮ ൅ ௞ߜோ, ௞ିଵ	 has a strictly positive coefficient ߙ௞, is energy 

stable. This completes the energy-stability proof. 

The following three remarks are in order. 

1. The energy-stability proof of Vincent et al. (2011) was extended to nonlinear conservation laws in 

Jameson, Vincent, and Castonguay (2012). There, it was also shown that similar to DG, the FR methods 

may encounter aliasing driven instability if the flux function is non-linear. They demonstrated that for 

nodal methods, the location of the solution points plays an important role in determining the extent of 

aliasing driven instabilities, and the strong quadrature points such as Gauss points are optimal in 

minimizing such instabilities. 

2. The energy-stability proof here also extends to nonlinear conservation laws provided that (a) the 

strong form S2 is employed and (b) the interface flux is an E-flux (Osher 1984, Jiang and Shu 1994); this 

extension is similar to that by Jameson, Vincent, and Castonguay (2012) mentioned above. 

3. It remains an open problem to prove or disprove that for the strong form S1, in the case of nonlinear 

conservation laws, the family of schemes discussed here is stable. The above proof does not extend for 

the following reason. A key ingredient of the proof is condition (5.19) or ߲௞ାଵ݂௝/߲ߦ௞ାଵ ൌ 0. For the 

strong form S1 of (3.13), the quantity ߲௞ାଵ݂௝/߲ߦ௞ାଵ is replaced by ߲௞ሺ࣪௞ሾሺ݂ሺݑ௛ሻሻకሿሻ/߲ߦ௞, which may 

not vanish. Also note that for nonlinear conservation laws and the strong form S1, whereas the stability of 

the family of schemes is not known, the stability for DG using the weak form (and thus the strong form 

S1) is well established as discussed in (Jiang and Shu 1994). 

6. Flux Reconstruction/Correction Procedure via Reconstruction 

(FR/CPR) Schemes 

The FR methods, which are obtained by different criteria of approximation for the (Dirac) delta 

function in the strong form S1, or different definitions for the correction functions in constructing ܨ ௝via 

(4.18), are discussed in this section. Their accuracy and stability properties by Fourier analyses are also 
NASA/TM—2014-218135 19 

mentioned. The current presentation of these methods provides additional relations not available in the 

literature. 

Focusing on the right boundary, recall that in the case of the DG scheme, the correction function 

݃ோ, ஽ீ is of degree ݇൅1 and defined by the ݇൅2 conditions of ݃ோ, ୈୋሺെ1ሻ ൌ 0, ݃ோ, ୈୋሺ1ሻ ൌ 1, and 

࣪௞ିଵ൫݃ோ, ୈୋ൯ൌ0 of (4.12). Its derivative is the approximate delta function ߜோ ൌ ߜோ, ௞ given in (3.8b). 

To obtain a family of schemes, the correction function ݃ோ of degree ݇൅1 is required to satisfy 

݃ோሺെ1ሻ ൌ 0, ݃ோሺ1ሻ ൌ 1, (6.1a,b)

together with ݇ additional conditions. As discussed before (4.18), when reconstructing the flux, the 

condition ݃ோሺ1ሻ ൌ 1 serves the purpose of recovering the common flux value ݂௝ାଵ/ଶ

ூ at the right interface 

whereas condition ݃ோሺെ1ሻ ൌ 0 leaves the value at the left interface unchanged. The ݇ additional 

conditions correspond to the fact that ݃ோ approximates the zero function on ሾെ1, 1ሻ (consistent with the 

remark after (4.11) that ݃ோ approximates the unit step function) in some sense. These approximation 

criteria must be chosen in a manner that the corresponding scheme is stable and accurate. 

6.1. Radau and Lobatto polynomials and their derivatives. The following polynomials play an 

important role in the definitions of various correction functions and their derivatives. 

Radau polynomials. Let the left Radau polynomial of degree ݇൅1 be defined by 

ܴ௅, ௞ାଵ ൌ 1

2 ሺܮ௞ାଵ ൅ ܮ௞ሻ. (6.2)

Whereas ܮ௞ାଵ vanishes at the ݇൅1 Gauss points, ܴ௅, ௞ାଵ vanishes at the ݇൅1 left Radau points. At the 

boundaries, 

ܴ௅, ௞ାଵሺെ1ሻ ൌ 0 and ܴ௅, ௞ାଵሺ1ሻ ൌ 1. (6.3a,b)

Since both ܮ ௞and ܮ௞ାଵ are orthogonal to ࡼ௞ିଵ, ܴ௅, ௞ାଵ also possesses this property. Equivalently, 

࣪௞ିଵ൫ܴ௅, ௞ାଵ൯ ൌ 0. (6.4)

Equations (6.3) and (6.4) above imply that ܴ௅, ௞ାଵ is the correction function ݃ோ, ୈୋ. Thus, by (4.10) and 

(3.8), 

ܴ௅, ௞ାଵᇱ ൌ ߜோ, ௞ ൌ ෍ 2݅ ൅ 1

௜ܮ 2

௞

௠ୀ଴

. (6.5)

In a similar manner, the right Radau polynomial of degree ݇൅1 is defined by 

ܴோ, ௞ାଵ ൌ ሺെ1ሻ௞ାଵ

2 ሺܮ௞ାଵ െ ܮ௞ሻ. (6.6)

Then ܴோ, ௞ାଵ vanishes at the ݇൅1 right Radau points. In addition, 

ܴோ, ௞ାଵሺെ1ሻ ൌ 1, ܴோ, ௞ାଵሺ1ሻ ൌ 0, (6.7a,b)

and 
NASA/TM—2014-218135 20 

࣪௞ିଵ൫ܴோ, ௞ାଵ൯ ൌ 0. (6.8)

Properties (6.7) and (6.8) above imply that ܴோ, ௞ାଵ is the correction function ݃௅, ୈୋ. By (4.15) and (3.8), 

ܴோ, ௞ାଵᇱ ൌ െߜ௅, ௞ ൌെ෍ 2݅ ൅ 1

2 ሺെ1ሻ௜ ܮ௜

௞

௠ୀ଴

. (6.9)

The plots of the first few left and right Radau polynomials and their derivatives are shown in 

Figures 6.1 and 6.2. 

(a) (b) 

Figure 6.1.—(a) Graphs of the left Radau polynomials ܴ௅, ௞, ݇൑4 (the dots are the 4 left Radau points) and 

(b) their derivatives. The ݇െ1 zeros of ܴ௅, ௞ᇱ

 together with ݔൌ1 form the ݇ right Radau points.

(a) (b) 

Figure 6.2.—Graphs of (a) the right Radau polynomials ܴோ, ௞, ݇൑4 (the dots are the 4 right Radau points) and 

(b) their derivatives. The ݇െ1 zeros of ܴோ, ௞ᇱ

 together with ݔ ൌ െ1 form the ݇ left Radau points. 

Lobatto polynomials. For any integer ݇൒2, let the Lobatto polynomial of degree ݇ be defined by 
NASA/TM—2014-218135 21 

Lo௞ ൌ ܮ ௞െ ܮ௞ିଶ. (6.10)

Then 

Lo௞ሺെ1ሻ ൌ Lo௞ሺ1ሻ ൌ 0. (6.11)

and, by definition (6.4), Lo௞ is orthogonal to ࡼ௞ିଷ. The zeros of the Lobatto polynomial of degree k are 

the k Lobatto points; they include the two boundaries 1. The derivative of Lo௞ାଵ, which has not been 

employed by FR authors and play an important role here, is given by 

Lo௞ାଵᇱ ൌ ሺ2݇ ൅ 1ሻܮ௞.) 6.12(

The plots of the Lobatto polynomials Lo௞, ݇ ൌ 1, ... , 4, and their derivatives are shown in Figure 6.3. 

6.2. FR schemes. Due to symmetry, we focus only on ݃ோ (݃௅ involves an extra factor (–1)k). For the 

DG scheme, as discussed after (6.14), 

݃ோ, ୈୋ ൌ ܴ௅, ௞ାଵ. (6.13)

Note that the correction function for the right boundary is the left Radau polynomial. 

The correction functions for the ESFR methods of Section 5 can now be derived. First, for simplicity, 

denote 

݃ൌ݃ோ.

(a) (b) 

Figure 6.3.—Graphs of (a) the Lobatto polynomials Lok, k = 1,...,4 and (b) their derivatives. 

 
NASA/TM—2014-218135 22 

Next, ݃ is assumed to satisfy the two conditions (6.1). For the ݇ additional conditions, instead of 

requiring ݃ to be orthogonal to ࡼ௞ିଵ as in the case of the DG scheme, we require it to be orthogonal to 

ࡼ௞ିଶ. It was verified via Fourier analysis (Huynh 2007, 2009) that if ݃ is orthogonal to ࡼ௞ିଶ the 

resulting scheme is stable (the converse is not true, however). 

Since both ܴ௅, ௞ାଵ and ܴ௅, ௞ are orthogonal to ࡼ௞ିଶ, if ݃ orthogonal to ࡼ௞ିଶ, it can be written as 

݃ ൌ ߙܴ௅,௞ ାଵ ൅ ሺ1 െ ߙሻܴ௅, ௞ (6.14)

where ߙ remains to be specified (note that the parameter ߙ here is completely independent from the ߙ for 

the location of the approximate delta function in Section 3). The current author considered only 0൑ߙ൑

1 in (Huynh 2009) for the obvious reason that such a correction function is a weighted average of the two 

Radau polynomials (correction functions). For ߙ൐1, the correction function becomes steeper than ݃஽ீ

and, in conjunction with an explicit time stepping method, the largest possible time step size or CFL limit 

becomes smaller than that for DG, a disadvantage. For ߙ൏0 the scheme is unstable. Concerning only 

energy stability and accuracy, however, ߙ can be larger than 1. 

The above implies, 

݃ ൌ ߙ൫ܴ௅, ௞ାଵ െ ܴ௅, ௞൯൅ܴ௅, ௞ ൌ ߙ

2 Lo௞ାଵ ൅ ܴோ, ௞. (6.15)

The second equality is a result of (6.10) and (6.2). Note that ߙൌ1 recovers the DG correction function. 

By (6.12), the above implies 

݃ᇱ ൌ ߙ

2 ሺ2݇ ൅ 1ሻܮ ൅ ௞ܴ௅,௞ 

ᇱ ൌ ߙ

2 ሺ2݇ ൅ 1ሻܮ ൅ ௞ߜோ, ௞ିଵ. (6.16)

Thus, by (5.10), we obtain the ESFR methods with 

ߙ ൌ௞ ߙ

2 ሺ2݇ ൅ 1ሻ, (6.17)

and the method is energy stable if ߙ൐0. 

Concerning accuracy, it was found by experiments using Fourier analyses (Huynh 2007) that if ݃ is 

orthogonal to ࡼ௠, then resulting scheme is accurate to order ݇൅݉൅2. Using this criteria, the DG 

scheme is accurate to order 2݇ ൅ 1, and the schemes of the ESFR family, order 2݇. The expected order of 

accuracy is ݇൅1. Due to this higher than expected accuracy, these methods are said to be super-accurate 

or to possess the property of super-convergence. 

Next, noteworthy members of the one parameter family of correction functions (6.14) are discussed. 

The first choice for ݃ is (6.13), which results in the DG scheme. 

The second choice for ݃, denoted by ݃ଶ or ݃୐୳୫୮, ୐୭ (for ‘lumping for Lobatto points’), is defined as 

follows. Since a steeper correction function tends to result in a scheme with a smaller CFL limit, to make 

݃ less steep, the extra condition is obtained by pushing one of the zeros to the left boundary, i.e., ߦ ൌ െ1

is a zero of multiplicity two. After some algebra (Huynh 2007), 

		݃ଶ ൌ ݃୐୳୫୮, ୐୭ ൌ ݇

2݇ ൅ 1 ܴ௅, ௞ାଵ ൅

݇൅1

2݇ ൅ 1 ܴ௅, ௞.
NASA/TM—2014-218135 23 

The function ݃ଶ has the following remarkable property. Among the ݇൅1 Lobatto points (and we can 

use them as solution points), ݃ଶ

ᇱ vanishes at ݇ of them; the exception is the right boundary where the 

right flux jump locates. 

The final choice here for ݃ requires that g vanishes at the ݇ Gauss points (Huynh2007), 

	݃ୋୟ ൌ ݇൅1

2݇ ൅ 1 ܴ௅, ௞ାଵ ൅

݇

2݇ ൅ 1 ܴ௅, ௞.

Whereas the staggered-grid (Kopriva and Kolias 1996) and SD schemes are mildly unstable, the above 

provides a modification using the ݇ Gauss points as flux points. The resulting scheme is both Fourier and 

energy stable for all ݇. 

Also note that the steepest slope of ݃, which takes place at the right boundary, relates to the time step 

size or CFL limit. For example, for DG, ݃ᇱ ൌ ܴ௅, ௞ାଵ

ᇱ ሺ1ሻ ൌ ሺ݇ ൅ 1ሻଶ/2; when an explicit Runge-Kutta 

method is employed, it is well known that the CFL limit for the DG scheme of degree ݇ is roughly 

proportional to 1/݇ଶ. 

As an example of correction functions, Fig. 6.4(a) shows the plots of the functions ݃ୈୋ, ݃ୋୟ, and ݃ଶ

for ݇ൌ3. Figure 6.4(b) shows the (top half of) spectra of the corresponding schemes. Note that all three 

methods are stable since the spectra lie in the left half of the complex plane. The spectra intersect the real 

axis at ݔୈୋୀ െ 19.2, ݔୋୟୀ െ 12.3, and ݔଶୀ െ 9.6. If the RK4 method is employed for time stepping, then 

the CFL limits are approximately .145, .227, and .289, respectively. 

Table 1 tabulates the orders of accuracy and errors of these semi-discrete schemes for ݇ൌ3. Here, the 

coarse mesh error corresponds to ݓ ൌ ߨ/4, and fine mesh error,	ݓ ൌ ߨ/8 so that errors are away from 

machine zero, and the calculation of order of accuracy is valid. 

Note that for each ݇, all errors have negative real parts—a fact consistent with the stability of the three 

schemes. In addition, the dominant part of the error for the DG scheme is real (dissipation type) whereas 

those for the other two schemes are imaginary (dispersion type). 

(a) Correction functions (b) Spectra 

Figure 6.4.—(a) Correction functions for k = 3, and (b) their spectra. 

 

g2

gGa

gDG
NASA/TM—2014-218135 24 

TABLE 1.—ORDERS OF ACCURACY AND ERRORS OF SCHEMES FOR ݇ൌ3

Scheme Order of 

accuracy 

Coarse mesh error, w  / 4 Fine mesh error, w /8

DG 7 i 7 8 1. 10 1. 10       i 10 11 4. 10 2. 10      

g Ga 6 i 7 6 3.1 10 1.3 10       i 9 8 1.2 10 1.1 10      

2 g 6 i 7 6 5.4 10 2.3 10       i 9 8 2.2 10 1.9 10      

A brief review of Fourier analysis for these type of methods can be found in the Appendix. 

Concerning open problems, it remains to be shown whether Fourier stability is equivalent to energy 

stability. 

Also note that there exists correction functions that are not orthogonal to any ࡼ௠, but the resulting 

scheme is Fourier stable. As an example, the correction function ݃ோ that lumps all corrections to the right 

boundary when the solution points are the Chebyshev-Lobatto points (Huynh 2007) possesses this 

property. It remains an open problem to identify all correction functions that yields stable schemes. 

7. Conclusions and Discussion 

In conclusion, a formulation for the discontinuous Galerkin (DG) method that leads to solutions using 

the differential form of the equation was presented. The formulation includes (a) a derivative calculation 

that involves only data within each cell with no data interaction among cells, and (b) for each cell, 

corrections to this derivative that deal with the jumps in fluxes at the cell boundaries and allow data 

across cells to interact. The derivative with no interaction is obtained by projection or interpolation, and 

the corrections are derived using the approximate (Dirac) delta functions. The formulation results in a 

family of schemes: different approximate delta functions give rise to different methods. It was shown that 

the current formulation is equivalent to the flux reconstruction (FR) formulation. Due to the use of 

approximate delta functions, an energy stability proof simpler than that of Vincent, Castonguay, and 

Jameson (2011) for a family of schemes was derived. Accuracy and stability of resulting schemes were 

examined via Fourier analyses. Some open problems were discussed. The current and FR formulations 

appear to be promising toward the goal of faster and more accurate solutions of conservation laws for 

practical applications. 
NASA/TM—2014-218135 25 

Appendix—A Brief Review of Fourier Stability and Accuracy Analysis 

Accuracy and stability of a numerical method are critical and can be derived via Fourier analyses. On 

ሺെ∞, ∞ሻ, consider the advection equation 

ݑ௧ ൅ ݑ௫ ൌ 0. (A.1)

The initial condition at ݐൌ0 is periodic: ݑ୧୬୧୲ ൌ ݁ூ௪௫ where ݓ is a real number between െߨ and ߨ called 

a wave number and, in this section, ܫ ൌ √െ1. Low frequency data corresponds to ݓ of small magnitude, 

high frequency, to ݓ near േߨ. The exact solution is ݑୣ୶ୟୡ୲ሺݔ, ݐሻ ൌ ݁ூ௪ሺ௫ି௧ሻ. At ݔൌ0 and ݐൌ0, 

ሺݑୣ୶ୟୡ୲ሻ௧ሺ0, 0ሻ ൌ െݓܫ.) A.2)

The cells are ܧ ௝ൌ ሾ݆ െ 1/2, ݆ ൅ 1/2ሿ. At time ݐൌ0, the data ݑ୧୬୧୲ on ܧ ௝is ݁ூ௪ሺ௝ାక/ଶሻ where |ߦ ൑ |1. 

Thus, the exact data on ܧ௝ିଵ is 

݁ିூ௪݁ூ௪ሺ௝ାక/ଶሻ.

In the local description, ݑ௝ሺߦሻ ൌ ∑ ݑ௝,	௜ ܮ௜

௞

௜ୀ଴ ሺߦሻ. For the calculation of eigenvalues below, we assume the 

data share the above property with the exact solution: ݑ௝ିଵሺߦሻ ൌ ݁ିூ௪ݑ௝ሺߦሻ, or, for all ݅, 0൑݅൑݇, 

ݑ௝ିଵ, ௜ ൌ ݁ିூ௪ݑ௝,௜ .) A.3)

To calculate the eigenvalues, the ݇൅1 values ݑ௝, ௜ are grouped together as a vector: with superscript T

denoting the transpose, set 

൯ ௞,௝ݑ,...,଴ ,௝ݑ൫ ൌ௝ ࢛

்

. (A.4)

The solution of (A.1) by the methods of the previous sections on each cell ܧ ௝can be expressed as 

݀

(.5A௝ (࢛଴࡯ ൅ ଵି௝࢛ଵି࡯ ൌ௝ ࢛ ݐ݀

where ࡯ିଵ and ࡯ ଴are ሺ݇൅1ሻ ൈ ሺ݇ ൅ 1ሻ matrices. Using (A.3), we replace ࢛௝ିଵ by ݁ିூ௪࢛௝. The spatial 

discretization results in 

݀

(.6A௝ (࢛ࡿ ൌ௝ ࢛ ݐ݀

where the ሺ݇൅1ሻ ൈ ሺ݇ ൅ 1ሻ matrix	ࡿ) for ‘space’ or ‘semi-discrete’) is given by 

ࡿൌ݁ିூ௪࡯ିଵ ൅ ࡯଴.) A.7)

Equation (A.6) is similar to the differential equation ݀ݑ/݀ݐ ൌ ݑߣ whose solution is ܿ݁ఒ௧ where c is an 

arbitrary constant. If Reሺߣሻ ൑ 0, the solution is stable. Here, the eigenvalues of ࡿ take the place of ߣ. For 

all schemes discussed here, ࡿ has ݇൅1 eigenvalues. Among these, the principal eigenvalue, denoted by 

ܵሺݓሻ, is the one approximating the exact value – ݓܫ of (A.2). All eigenvalues must lie in the left half of 

the complex plane for the semi-discretization to be stable. The collection of all eigenvalues forms the 

spectrum of the scheme. 
NASA/TM—2014-218135 26 

To find the order of accuracy, note that with a uniform mesh of width ݄, the principal eigenvalue ܵሺݓሻ

‘approximates’ െ݄߲/߲ݔ. A scheme is accurate to order m if S ‘approximates’ െ݄߲/߲ݔ to ܱሺ݄௠ାଵሻ; more 

precisely, for small w, 

ܵሺݓሻ ൌ െݓܫ ൅ ܱሺݓ௠ାଵሻ. (A.8)

In practice, we obtain the order of accuracy of a scheme by the following procedure. First, set ݓ to be, 

say, ߨ/4. We can estimate the error 

ܧሺݓሻ ൌ ܵሺݓሻ ൅ ݓܫ) A.9)

By halving the wave number ݓ) which is equivalent to doubling the number of mesh points), the error 

corresponding to ݓ/2) ൌ ߨ/8 (is 

ܧሺݓ/2ሻ ൌ ܵሺݓ/2ሻ ൅ ݓܫ/2.) A.10)

Since ܱሺሺݓ/2ሻ௠ାଵሻ ൎ ሺ1/2ሻ௠ାଵܱሺሺݓሻ௠ାଵሻ, for a scheme to be m-th order accurate, the following 

condition must hold: ܧሺݓሻ/ܧሺݓ/2ሻ ൎ 2௠ାଵ. That is, the order of accuracy is given by 

݉ ൎ Log ቆ ܧሺݓሻ

ሻ/2ݓሺܧ

ቇ Log ൗ ሺ2ሻ െ 1. (A.11)
NASA/TM—2014-218135 27 

References 

Y. Allaneau and A. Jameson (2011), Connections between the filtered discontinuous Galerkin method 

and the flux reconstruction approach to high order discretizations, Comput. Methods Appl. Mech. Engrg., 

200 (2011) 3628–3636 

F. Bassi and S. Rebay (1997a), A high-order accurate discontinuous finite element method for the 

numerical solution of the compressible Navier-Stokes equations, J. Comput. Phys., 131 (1997), pp. 267-279. 

F. Bassi and S. Rebay (1997b), High-order accurate discontinuous finite element solution for the 2D 

Euler equations, J. Comput. Phys., 138, pp. 251-285. 

F. Bassi and S. Rebay (2000), A high order discontinuous Galerkin method for compressible turbulent 

flows, in “Discontinuous Galerkin methods: Theory, Computation, and Application”, B. Cockburn, G. 

Karniadakis, and C.-W. Shu, editors, Lecture Notes in Computational Science and Engineering, Springer, 

pp. 77-88. 

P. Castonguay, P. E. Vincent, A. Jameson (2012), A New Class of High-Order Energy Stable Flux 

Reconstruction Schemes for Triangular Elements, Journal of Scientific Computing, Volume 51, Number 1, 

Pages 224-256. 

P. Castonguay, D. M. Williams, P. E. Vincent, A. Jameson (2013), Energy Stable Flux Reconstruction 

Schemes for Advection-Diffusion Problems. P. Computer Methods in Applied Mechanics and Engineering. 

Under Review. 

B. Cockburn, G. Karniadakis, and C.-W. Shu, Eds., (2000), Discontinuous Galerkin methods: Theory, 

Computation, and Application (Springer). 

B. Cockburn and C.-W. Shu (1998), The local discontinuous Galerkin methods for time-dependent 

convection diffusion systems, SIAM J. Numer. Anal., 35, pp. 2440-2463 

B. Cockburn and C.-W. Shu (2005), Foreword for the special issue on discontinuous Galerkin method, 

Journal of Scientific Computing, 22-23, pp. 1–3. 

B. Cockburn and C.-W. Shu (2009), Foreword for the special issue on discontinuous Galerkin method, 

Journal of Scientific Computing, 40, pp. 1–3. 

H. Gao and Z.J. Wang (2009), A High-Order Lifting Collocation Penalty Formulation for the NavierStokes Equations on 2-D Mixed Grids, AIAA-2009-3784. 

H. Gao and Z.J. Wang (2013), A Conservative Correction Procedure via Reconstruction Formulation 

with the Chain-Rule Divergence Evaluation”, J. Computational Physics 232, 7–13. 

H. Gao, Z.J. Wang, and H. T. Huynh (2013), Differential Formulation of Discontinuous Galerkin and 

Related Methods for the Navier-Stokes Equations, Commun. Comput. Phys, Vol. 13, No. 4, pp. 1013-1044. 

T. Haga, H. Gao, Z.J. Wang (2010), A High-Order Unifying Discontinuous Formulation for 3-D Mixed 

Grids, AIAA Paper 2010-540. 

T. Haga, H. Gao and Z. J. Wang (2011), A High-Order Unifying Discontinuous Formulation for the 

Navier-Stokes Equations on 3D Mixed Grids, Math. Model. Nat. Phenom., Vol. 6 (3), 28-56. 

J.S. Hesthaven and Tim Warburton (2008), Nodal Discontinuous Galerkin Methods (Springer). 

F.B. Hildebrand (1987), Introduction to Numerical Analysis, Second Edition, Dover Books on Advanced 

Mathematics. 

H.T. Huynh (2007), A flux reconstruction approach to high-order schemes including discontinuous 

Galerkin methods, AIAA Paper 2007-4079. 

H.T. Huynh (2009), A Reconstruction Approach to High-Order Schemes Including Discontinuous 

Galerkin for Diffusion, AIAA Paper 2009-403. 

H.T. Huynh, Z. J. Wang, and P. E. Vincent (2014), High-order methods for computational fluid 

dynamics: A brief review of compact differential formulations on unstructured grids, Computers and Fluids, 

to appear. 

A. Jameson (2010), A proof of the stability of the spectral difference method for all orders of accuracy, J. 

Sci. Comput. 45(1–3), 348–358. 

A. Jameson, P. E. Vincent, P. Castonguay (2012), On the Non-Linear Stability of Flux Reconstruction 

Schemes, Journal of Scientific Computing, Volume 50, Number 2, Pages 434-445 
NASA/TM—2014-218135 28 

G. Jiang and C.-W. Shu (1994), On cell entropy inequality for discontinuous Galerkin methods, Math. 

Comp., Vol. 62, Number 206, pp. 531-538. 

D.A. Kopriva and J.H. Kolias (1996), A conservative staggered-grid Chebyshev multidomain method for 

compressible flows, J. Comput. Phys. 125, 244. 

P. LaSaint and P.A. Raviart (1974), On a finite element method for solving the neutron transport 

equation, in Mathematical aspects of finite elements in partial differential equations, C. de Boor (ed.), 

Academic Press, pp. 89–145. 

C. Liang, A. Jameson, Z.J. Wang (2009a), Spectral difference method for compressible flow on 

unstructured grids with mixed elements. J. Comput. Phys. 228, 2847. 

C. Liang, S. Premasuthan, A. Jameson (2009b), High-order accurate simulation of low-Mach laminar 

flow past two side-by-side cylinders using spectral difference method. Comput. Struct. 87, 812. 

Y. Liu, M. Vinokur, and Z.J. Wang (2006), Discontinuous Spectral Difference Method for Conservation 

Laws on Unstructured Grids, J. Comput. Phys., 216, 780–801. 

P.M. Morse and H. Feshbach, Methods of Theoretical Physics, Vol. 1, McGraw-Hill, New York, 1953. 

S. Osher (1984), Riemann solvers, the entropy condition, and difference approximations, SIAM J. 

Numer. Anal. Vol. 21, No 2, 217-235. 

W.H. Reed and T.R. Hill (1973), Triangular mesh methods for the neutron transport equation, Los 

Alamos Scientific Laboratory Report, LA-UR-73-479. 

C.-W. Shu (2012), Discontinuous Galerkin method for time dependent problems: survey and recent 

developments. Preprint. 

P.E. Vincent, P. Castonguay, A. Jameson (2011a), A New Class of High-Order Energy Stable Flux 

Reconstruction Schemes, Journal of Scientific Computing, Volume 47, Number 1, Pages 50-72. 

P.E. Vincent, P. Castonguay, A. Jameson (2011b), Insights from von Neumann Analysis of High-Order 

Flux Reconstruction Schemes, J. Comput. Phys., Vol. 230, Issue 22, pp. 8134-8154. 

Z.J. Wang (2007), High-order methods for the Euler and Navier-Stokes equations on unstructured grids, 

Journal of Progress in Aerospace Sciences, Vol. 43, pp. 1–41. 

Z.J. Wang (2011), Adaptive High-order methods in Computational Fluid Dynamics, Advances in 

Computational Fluid Dynamics Vol. 2, World Scientific Publishing Co. 

Z.J. Wang et al. (2013), High-Order CFD Methods: Current Status and Perspective, Int. J. Numer. Meth. 

Fluids, 72, No. 8, pp. 811–845. 

Z.J. Wang and H. Gao (2009), A Unifying Lifting Collocation Penalty formulation including the 

discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids, J. 

Comput. Phys., 228, No. 2, pp. 8161–8186. 

Z.J. Wang, H. Gao and T. Haga (2011). A Unifying Discontinuous Formulation for Hybrid Meshes, in 

Adaptive High-Order Methods in Computational Fluid Dynamics, World Scientific, edited by Z.J. Wang. 

Z.J. Wang, Y. Liu, G. May, A. Jameson (2007), Spectral difference method for unstructured grids II: 

extension to the Euler equations. J. Sci. Comput. 32, 45 

Z.J. Wang, L. Zhang and Y. Liu (2004), Spectral (finite) volume method for conservation laws on 

unstructured grids IV: extension to two-dimensional Euler equations, J. Comput. Phys., 194, No. 2, pp. 716–

741. 

D.M. Williams, P. Castonguay, P. E. Vincent, A. Jameson (2011), An Extension of Energy Stable Flux 

Reconstruction to Unsteady, Non-linear, Viscous Problems on Mixed Grids, AIAA Paper 2011-3405. 

D.M. Williams, P. Castonguay, P. E. Vincent, A. Jameson (2013), Energy Stable Flux Reconstruction 

Schemes for Advection-Diffusion Problems on Triangles, Journal of Computational Physics, Under Review. 

M. Yu, Z. J. Wang, and H. Hu (2012), High-Fidelity Flapping-Wing Aerodynamics Simulations with a 

Dynamic Unstructured Grid Based Spectral Difference Method, in Proceedings of the 7th International 

Conference on Computational Fluid Dynamics, ICCFD7-4104. 

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-06-2014

2. REPORT TYPE

Technical Memorandum

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

On Formulations of Discontinuous Galerkin and Related Methods for Conservation

Laws

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Huynh, H., T.

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 794072.02.03.02.03

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-18938

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-2014-218135

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited

Subject Category: 59

Available electronically at http://www.sti.nasa.gov

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

A formulation for the discontinuous Galerkin (DG) method that leads to solutions using the differential form of the equation (as opposed to

the standard integral form) is presented. The formulation includes (a) a derivative calculation that involves only data within each cell with

no data interaction among cells, and (b) for each cell, corrections to this derivative that deal with the jumps in fluxes at the cell boundaries

and allow data across cells to interact. The derivative with no interaction is obtained by a projection, but for nodal-type methods, evaluating

this derivative by interpolation at the nodal points is more economical. The corrections are derived using the approximate (Dirac) delta

functions. The formulation results in a family of schemes: different approximate delta functions give rise to different methods. It is shown

that the current formulation is essentially equivalent to the flux reconstruction (FR) formulation. Due to the use of approximate delta

functions, an energy stability proof simpler than that of Vincent, Castonguay, and Jameson (2011) for a family of schemes is derived.

Accuracy and stability of resulting schemes are discussed via Fourier analyses. Similar to FR, the current formulation provides a unifying

framework for high-order methods by recovering the DG, spectral difference (SD), and spectral volume (SV) schemes. It also yields stable,

accurate, and economical methods.

15. SUBJECT TERMS

Discontinuous Galerkin; High-order methods; Conservation laws; Flux reconstruction

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU

18. NUMBER

 OF

 PAGES

36

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