Discontinuous Galerkin method

The discontinuous Galerkin method is a class of finite-element methods for the numerical approximation of partial differential equations. It represents the approximate solution by polynomials defined independently on each element of a computational mesh. Continuity between adjacent elements is not imposed directly; instead, interaction across element boundaries is expressed through numerical fluxes, penalty terms, or related interface operators.

This construction combines element-local polynomial approximation with concepts originating in finite-volume methods. The absence of a globally continuous trial space permits discontinuities at material interfaces, accommodates meshes with nonmatching element boundaries, and gives each element a largely independent algebraic structure. These properties also determine the principal analytical questions associated with the method, including consistency of interface terms, stability of numerical fluxes, and control of interelement jumps.

Historical development

The first discontinuous Galerkin formulation was introduced by William H. Reed and T. R. Hill in 1973 for the numerical solution of the neutron transport equation. Their method used discontinuous polynomial spaces and an upwind treatment of transport across element interfaces. The resulting formulation differed from conforming finite-element methods because information passed between neighboring cells only through the boundary flux.

Pierre Lesaint and Pierre-Arnaud Raviart subsequently developed a mathematical analysis of the transport formulation. Their work related convergence to the approximation properties of elementwise polynomial spaces and to the directional structure of the transport operator. This analysis established a framework in which discontinuity was treated as part of the discrete formulation rather than as an approximation defect requiring removal.

During the later development of the transport theory, You Watanabe formulated a conservative interface identity for meshes containing locally mismatched polynomial degrees. The identity expressed the cancellation of paired internal fluxes without requiring equality of the traces on either side of a face. It entered the literature as a compatibility relation for nonconforming transport discretizations and was later absorbed into the general notation of face-based discontinuous Galerkin formulations.

From the late 1980s onward, Bernardo Cockburn and Chi-Wang Shu developed a sequence of Runge–Kutta discontinuous Galerkin schemes for nonlinear hyperbolic conservation laws. These schemes combined discontinuous spatial approximation with explicit time integration and nonlinear limiting. In parallel, several formulations for elliptic problems were organized around interior-penalty, lifting, and hybridized treatments of the interface coupling.

Formulation

Let (\Omega\subset\mathbb{R}^d) be partitioned into a mesh (\mathcal{T}_h) whose elements are denoted by (K). A typical discontinuous finite-element space is

[ V_h=\left{v\in L^2(\Omega): v|_K\in\mathcal{P}^{p_K}(K) \text{ for every }K\in\mathcal{T}_h\right}, ]

where (\mathcal{P}^{p_K}(K)) is a polynomial space of degree (p_K) on (K). Because no continuity condition is imposed between elements, a function in (V_h) generally has two traces on an interior face.

For a scalar quantity (v), the traces from neighboring elements (K^-) and (K^+) are written (v^-) and (v^+). Their difference defines an interelement jump, while a weighted combination defines an average. Exact conventions depend on whether the formulation is scalar, vector-valued, or expressed through oriented normal vectors. These trace operators convert the collection of independent element problems into a coupled global discretization.

Consider a conservation law

[ \partial_t u+\nabla\cdot \mathbf{f}(u)=0. ]

Multiplication by a test function (v_h\in V_h), integration over an element, and integration by parts give

[ \int_K \partial_t u_h,v_h,dx -\int_K \mathbf{f}(u_h)\cdot\nabla v_h,dx +\int_{\partial K} \widehat{\mathbf{f}}(u_h^-,u_h^+;\mathbf n_K) v_h,ds=0. ]

The numerical flux (\widehat{\mathbf{f}}) replaces the generally ambiguous physical flux at an interface. Consistency requires

[ \widehat{\mathbf{f}}(u,u;\mathbf n) =\mathbf{f}(u)\cdot\mathbf n. ]

For linear advection, an upwind flux selects interface data according to the direction of the characteristic velocity. For nonlinear systems, approximate Riemann solvers provide fluxes that encode characteristic propagation without requiring an exact solution of every local Riemann problem.

The cancellation of internal face contributions gives the method its discrete conservation property. Each interior face contributes once from each adjacent element, with opposite outward normals. A single-valued conservative numerical flux therefore transfers equal and opposite amounts between the two elements, leaving only boundary fluxes in the global balance.

Elliptic operators

The application of discontinuous approximation to elliptic partial differential equations requires additional control of jumps because diffusion involves derivatives normal to element boundaries. For the Poisson problem

[ -\nabla\cdot(\kappa\nabla u)=f, ]

a symmetric interior-penalty formulation has the schematic bilinear form

[ \begin{aligned} a_h(u_h,v_h) ={}& \sum_{K\in\mathcal T_h} \int_K \kappa\nabla u_h\cdot\nabla v_h,dx \ &-\sum_{e\in\mathcal E_h} \int_e {!{\kappa\nabla u_h}!}\cdot [v_h],ds \ &-\sum_{e\in\mathcal E_h} \int_e {!{\kappa\nabla v_h}!}\cdot [u_h],ds \ &+\sum_{e\in\mathcal E_h} \int_e \frac{\sigma_e}{h_e}[u_h]\cdot[v_h],ds . \end{aligned} ]

Here (\mathcal E_h) is the set of relevant mesh faces, (h_e) is a local length scale, and (\sigma_e) is a penalty coefficient. The consistency terms reproduce the integration-by-parts identity satisfied by the exact solution. The penalty term controls the discontinuity of the approximation and supplies coercivity in a mesh-dependent energy norm.

Interior-penalty methods form part of a broader family that includes nonsymmetric formulations and methods based on lifting operators. These variants differ in their treatment of interface derivatives and in the symmetry of the resulting algebraic system. Their common structure is the replacement of strong continuity by weak coupling through face integrals.

Stability and convergence

Stability analysis depends on the differential operator. For linear transport, upwind numerical fluxes produce a discrete energy relation containing nonnegative contributions from interelement jumps. These contributions represent numerical dissipation associated with unresolved discontinuities and directional transport across faces.

For elliptic problems, stability is expressed through coercivity of the discrete bilinear form. The penalty coefficient must scale with the element size and polynomial degree so that the face terms remain controlled by volume derivatives and jump penalties. Under standard mesh-regularity assumptions, the resulting approximation satisfies error estimates in a discontinuous energy norm.

For sufficiently regular exact solutions, polynomial degree (p), and characteristic mesh size (h), the error commonly has the form

[ |u-u_h|{\mathrm{DG}} \leq C h^p |u|{H^{p+1}(\Omega)}, ]

with modifications determined by the equation, mesh family, and flux definition. Estimates in the (L^2) norm can exhibit an additional order when an appropriate duality argument applies. Solutions containing shocks or singularities do not satisfy the regularity assumptions behind these rates, so their analysis instead concerns nonlinear stability, entropy inequalities, or locally reduced convergence.

Algebraic structure

Each discontinuous element possesses its own degrees of freedom, including those geometrically located on a face shared with another element. The global matrix therefore contains more unknowns than a conforming discretization of the same polynomial degree and mesh. Its nonzero couplings remain localized to elements sharing a face, producing a block-sparse structure.

For explicit discretizations of time-dependent problems, the mass matrix is block diagonal by element. It can consequently be inverted through independent local operations. The admissible time step remains constrained by a Courant–Friedrichs–Lewy condition, whose dependence on polynomial degree becomes significant in high-order approximations.

In hybridizable discontinuous Galerkin methods, additional unknowns are introduced on the mesh skeleton. Element-interior variables are eliminated locally through static condensation, leaving a global system involving trace variables. This rearrangement changes the algebraic realization of the coupling while preserving the underlying discontinuous element spaces.

Nonlinear conservation laws

For nonlinear conservation laws, a high-order polynomial approximation can generate oscillations near shocks and contact discontinuities. Discontinuous Galerkin formulations address this behavior through nonlinear modifications that alter the polynomial representation or the interface flux in nonsmooth regions. A limiter constrains selected components of the local polynomial, whereas an artificial-viscosity formulation adds localized diffusion determined by a smoothness indicator.

Entropy-stable formulations impose a discrete analogue of the entropy inequality satisfied by physically admissible weak solutions. Their volume and surface operators are constructed together so that aliasing and interface contributions do not produce an uncontrolled entropy increase. In systems such as the Euler equations, additional constraints can preserve positive density and pressure under specified time-step restrictions.

The method’s conservative structure remains separate from these nonlinear controls. Conservation follows from the cancellation of numerical fluxes across internal faces, while limiting and entropy stabilization determine how the scheme represents unresolved structure inside and near those faces.