Discontinuous Galerkin Methods - From Finite Differences to High-Order Methods

Last Update: 2026-10-13

This training follows Finite Difference Methods I & II. Study the lecture notes by Shu [1] and complete the core exercises below. Additional references are listed in [2–5].

Purpose: Learn how discontinuous Galerkin methods extend finite difference ideas—upwinding, numerical dissipation, stability, and CFL restrictions—to element-wise polynomial approximation. The training combines independent implementation, numerical verification, and basic mathematical analysis.

The first-order model problem and the hyperbolic equation form the core training. The parabolic and dispersive equations are optional extensions, selected according to subsequent research needs.

References

1. Shu, Chi-Wang. Discontinuous Galerkin Methods: General Approach and Stability. Lecture notes.
2. Shu, Chi-Wang. High-order WENO and DG methods for time-dependent convection-dominated PDEs: a brief survey of several recent developments. Journal of Computational Physics 316 (2016), 598–613.
3. Di Pietro, Daniele Antonio; Ern, Alexandre. Mathematical Aspects of Discontinuous Galerkin Methods. Springer, 2012.
4. Ciarlet, Philippe G. The Finite Element Method for Elliptic Problems. North-Holland, 1978.
5. Brenner, Susanne C.; Scott, L. Ridgway. The Mathematical Theory of Finite Element Methods. Third edition. Springer, 2008.

General Implementation and Reporting Requirements

Use a uniform partition of $[0,1]$, with elements \(I_j=[x_{j-\frac12},x_{j+\frac12}],\) and the discontinuous approximation space \(V_h^k=\{v:\ v|_{I_j}\in\mathbb P^k(I_j)\}.\)

At an interface, $u^{-}$ and $u^{+}$ denote the traces from the left and right. Use a consistent jump convention throughout the report.

Some practical hints:

• Pay careful attention to your choice of polynomial basis on the reference element. Two common choices are modal basis (orthogonal Legendre polynomials) and nodal basis (Lagrange interpolating polynomials). Modal Legendre bases yield diagonal mass matrices, while Lagrange nodal bases produce dense mass matrices on each element. This difference strongly affects implementation, computational cost, and matrix inversion; it is a frequent source of beginner pitfalls.
• When performing convergence studies, always diagnose whether observed convergence rates represent true spatial‑order behaviour. Apparent order reduction can be artificially caused by insufficient quadrature accuracy, non‑negligible temporal integration error, or implementation bugs, rather than intrinsic properties of the DG scheme. Isolate these effects before drawing conclusions about the method.

For each implementation:

1. State the weak formulation, numerical flux, polynomial basis, and reference-to-physical element mapping.
2. Explain how the mass matrix, volume terms, and interface terms are assembled.
3. State the quadrature rules and ensure that quadrature errors do not dominate the reported convergence results.
4. For time-dependent problems, initialize the solution using the element-wise $L^2$ projection unless another choice is explicitly justified.
5. Report the time integrator, time-step rule, boundary conditions, and final time.
6. Submit independently written, reproducible code and a short report explaining the results.

For convergence studies, use \(h=\frac1{10},\frac1{20},\frac1{40},\frac1{80},\frac1{160},\ldots\) and compute \(\operatorname{order} = \frac{\log(E_h/E_{h/2})}{\log 2}.\)

Evaluate the $L^{2}$ error using sufficiently accurate quadrature. If the $L^\infty$ error is estimated by sampling, state the sampling procedure and check that further sampling refinement does not significantly change the result.

For time-dependent convergence tests, check the effect of temporal error by repeating selected runs with a smaller time step.

First-Order ODE Model Problem

Consider \(\begin{aligned} u_x&=f(x), &&0<x<1,\\ u(0)&=a. \end{aligned}\) Use the upwind numerical trace, with the prescribed value $a$ at the inflow boundary.

1. Formulation and comparison with finite differences

Derive the element-wise DG formulation.

Explain:

• What are the unknowns in the finite difference and DG methods?
• How does the boundary value enter the DG equations?
• How is information transferred between neighboring elements?
• What scheme is obtained when $k=0$?
• Why can the upwind DG system be solved sequentially from left to right?

2. Implementation and convergence

Implement the method for \(f(x)=\cos x,\qquad a=0,\) whose exact solution is \(u(x)=\sin x.\) Use $k=0,1,2$. Produce the $L^2$ and $L^\infty$ error tables and plot representative numerical solutions.

Check the expected accuracy of each polynomial degree.

3. Basic analysis

For the upwind DG method:

1. Prove existence and uniqueness of the discrete solution.
2. For $a=0$, prove an estimate of the form \(\|u_h\|_{L^2(0,1)} \leq C\|f\|_{L^2(0,1)},\) with $C$ independent of $h$.
3. Explain how the estimate changes when $a\neq0$.

Optional observation: Examine the error at the downwind element endpoints. Compare it with the global error and explain any special accuracy observed.

Linear Hyperbolic Equation

Consider \(u_t+u_x=0,\qquad 0<x<1,\) with periodic boundary conditions. The exact solution is \(u(x,t)=u_0(x-t),\) where the argument is interpreted periodically.

1. RKDG implementation

Use the upwind flux \(\widehat u=u^{-}.\)

Implement:

Use a time step proportional to $h$, with a CFL number appropriate to the polynomial degree and time integrator. State and test your choice.

Run the following initial conditions:

Smooth solution \(u_0(x)=\sin(2\pi x).\)

Discontinuous solution \(u_0(x)= \begin{cases} 1,&\frac14\leq x\leq\frac34,\\ 0,&\text{otherwise}. \end{cases}\)

For the smooth problem, report errors at $T=0.5$ and verify the convergence order. For both problems, plot solutions at $T=0.5$ and $T=1$.

Also check:

• conservation of the spatial integral;
• numerical dissipation;
• oscillations near discontinuities;
• sensitivity to the CFL number.

Do not expect the smooth-solution convergence order for the discontinuous initial condition.

2. Biased and central fluxes

Replace the upwind flux by \(\widehat u = \theta u^-+(1-\theta)u^+.\)

Test representative values, for example \(\theta=1,\quad \frac34,\quad \frac12,\quad \frac14.\)

1. Derive the semidiscrete energy identity \(\frac12\frac{d}{dt}\|u_h\|_{L^2}^2 + \left(\theta-\frac12\right)\sum_{\text{interfaces}}[u_h]^2=0.\)
2. Explain the distinction between dissipative, energy-conserving, and antidissipative fluxes.
3. Compare the numerical energy, solution profiles, and errors.
4. Distinguish semidiscrete stability from the stability of the chosen fully discrete method.

3. Basic theoretical analysis

For the semidiscrete, constant-coefficient, upwind DG method:

1. Prove $L^2$ stability.
2. Reproduce the optimal $L^2$ error estimate for sufficiently smooth solutions, following [1].
3. Identify the role of the projection used in the proof and explain why it is useful.

Present the argument in your own notation and make the cancellation of interface terms explicit.

4. Variable-coefficient extension

Consider \(u_t+(a(x)u)_x=0,\) where $a(x)$ is smooth and periodic, with flux \(\widehat f= \begin{cases} a(x)u^-,&a(x)\geq0,\\ a(x)u^+,&a(x)<0. \end{cases}\)

1. Derive a semidiscrete energy estimate, showing how $a_x$ enters the bound.
2. Explain why the $L^2$ norm need not be nonincreasing.
3. As an optional extension, reproduce an error estimate from the reference, stating its assumptions on $a$, the exact solution, and the projection.

TVD and TVB Limiters

Use the periodic linear transport problem from previous exercise.

1. Implementation

Study the limiter construction in [1] and implement a TVD limiter and a TVB limiter for $\mathbb P^1$. Extend to $\mathbb P^2$ after the $\mathbb P^1$ implementation has been verified.

State:

• the limiter formula and normalization;
• which coefficients or reconstructed quantities are modified;
• how the cell average is preserved;
• how limiting is incorporated into the Runge–Kutta stages.

2. Discontinuous solution

For the square-wave initial condition, compare:

• no limiter;
• TVD limiting;
• TVB limiting with $M=5,10,20,40$.

Report solution plots, overshoots and undershoots, conservation, and numerical smearing.

3. Smooth solution and the TVB parameter

For \(u_0(x)=\sin(2\pi x),\) repeat the convergence study with and without limiting.

1. Explain why a TVD limiter can reduce accuracy near smooth extrema.
2. Explain the role of the $Mh^2$ tolerance in the TVB modification.
3. Record how frequently and where the limiter is activated.
4. Investigate a suitable $M$ that preserves the expected accuracy for this test.

Do not interpret the resulting threshold as a universal constant: it depends on the test function and the precise limiter definition.

Parabolic Equation — Optional

Consider the forward heat equation \(u_t-u_{xx}=0,\qquad 0<x<1,\) with periodic boundary conditions and \(u(x,0)=\sin(2\pi x).\)

The exact solution is \(u(x,t)=e^{-4\pi^2t}\sin(2\pi x).\)

Use, for example, $T=0.1$. Choose the time step so that temporal errors do not dominate the spatial convergence study.

1. A naive DG discretization

Investigate the scheme \(\begin{aligned} \int_{I_j}(u_h)_t v\,dx + \int_{I_j}(u_h)_xv_x\,dx - \widehat{u_x}_{j+\frac12}v^{-}_{j+\frac12} +\widehat{u_x}_{j-\frac12}v^+_{j-\frac12} =0, \end{aligned}\) with \(\widehat{u_x} = \frac12\left((u_h)_x^-+(u_h)_x^+\right).\)

For $k=1,2$:

1. Implement the scheme and examine its behavior.
2. Report exact-solution errors where meaningful.
3. Compare solutions on meshes $h$ and $2h$, evaluating their difference on a common physical representation.
4. Explain why agreement between successive meshes is not sufficient to establish convergence to the correct solution.
5. Examine the discrete operator or energy relation to identify the scheme’s deficiencies.

2. Consistent diffusion discretizations

Implement one of the following first; additional comparisons are optional:

• LDG with central fluxes;
• LDG with alternating fluxes;
• symmetric interior penalty DG.

For the selected method:

1. Derive the weak formulation and specify all numerical fluxes.
2. For SIPG, state the penalty parameter and its scaling with $h$ and $k$.
3. Produce convergence tables for $k=1,2$.
4. Examine energy decay and time-step sensitivity.
5. Compare with the naive discretization.

Dispersive Wave Equation — Optional

Consider \(u_t+u_{xxx}=0,\qquad 0<x<1,\) with periodic boundary conditions and \(u(x,0)=\sin(2\pi x).\)

The exact solution is \(u(x,t)=\sin\left(2\pi x+(2\pi)^3t\right).\)

Implement either:

• an LDG formulation;
• an ultra-weak DG formulation.

For $k=1,2$:

1. Derive the formulation and specify the numerical fluxes.
2. Verify that the implemented flux choice has the intended stability properties.
3. Produce error tables and solution plots.
4. Examine amplitude and phase errors separately.
5. Check the effect of the time step, noting the restrictive scale associated with explicit treatment of a third derivative.

A comparison between the two formulations is an additional extension, not a prerequisite for completing the training.

Completion of the Core Training

The core training is complete when you can:

• derive and implement the DG formulation for a simple first-order problem;
• implement and verify RKDG for periodic linear transport;
• explain the effects of polynomial degree, numerical flux, CFL number, and time integration;
• prove the basic upwind energy estimate and follow the standard error analysis;
• implement TVD/TVB limiting and explain its effect on smooth and discontinuous solutions;
• produce reproducible numerical results and distinguish implementation errors from genuine numerical behavior.

Submit:

1. reproducible source code with a short README;
2. the required error tables and plots;
3. the core mathematical derivations;
4. a concise report discussing the main findings and any unresolved discrepancies.