# Discontinuous Galerkin method

In applied mathematics, discontinuous Galerkin methods (DG methods) form a class of numerical methods for solving differential equations. They combine features of the finite element and the finite volume framework and have been successfully applied to hyperbolic, elliptic, parabolic and mixed form problems arising from a wide range of applications. DG methods have in particular received considerable interest for problems with a dominant first-order part, e.g. in electrodynamics, fluid mechanics and plasma physics.

Discontinuous Galerkin methods were first proposed and analyzed in the early 1970s as a technique to numerically solve partial differential equations. In 1973 Reed and Hill introduced a DG method to solve the hyperbolic neutron transport equation.

The origin of the DG method for elliptic problems cannot be traced back to a single publication as features such as jump penalization in the modern sense were developed gradually. However, among the early influential contributors were Babuška, J.-L. Lions, Joachim Nitsche and Miloš Zlámal. DG methods for elliptic problems were already developed in a paper by Garth Baker in the setting of 4th order equations in 1977. A more complete account of the historical development and an introduction to DG methods for elliptic problems is given in a publication by Arnold, Brezzi, Cockburn and Marini. A number of research directions and challenges on DG methods are collected in the proceedings volume edited by Cockburn, Karniadakis and Shu.

## Overview

Much like the continuous Galerkin (CG) method, the discontinuous Galerkin (DG) method is a finite element method formulated relative to a weak formulation of a particular model system. Unlike traditional CG methods that are conforming, the DG method works over a trial space of functions that are only piecewise continuous, and thus often comprise more inclusive function spaces than the finite-dimensional inner product subspaces utilized in conforming methods.

As an example, consider the continuity equation for a scalar unknown ${\displaystyle \rho }$ in a spatial domain ${\displaystyle \Omega }$ without "sources" or "sinks" :

${\displaystyle {\frac {\partial \rho }{\partial t}}+\nabla \cdot \mathbf {J} =0,}$

where ${\displaystyle \mathbf {J} }$ is the flux of ${\displaystyle \rho }$.

Now consider the finite-dimensional space of discontinuous piecewise polynomial functions over the spatial domain ${\displaystyle \Omega }$ restricted to a discrete triangulation ${\displaystyle \Omega _{h}}$, written as

${\displaystyle S_{h}^{p}(\Omega _{h})=\{v_{|\Omega _{e_{i}}}\in P^{p}(\Omega _{e_{i}}),\ \ \forall \Omega _{e_{i}}\in \Omega _{h}\}}$

for ${\displaystyle P^{p}(\Omega _{e_{i}})}$ the space of polynomials with degrees less than or equal to ${\displaystyle p}$ over element ${\displaystyle \Omega _{e_{i}}}$ indexed by ${\displaystyle i}$. Then for finite element shape functions ${\displaystyle N_{j}\in P^{p}}$ the solution is represented by

${\displaystyle \rho _{h}=\sum _{j=1}^{\text{dofs}}\rho _{j}^{i}(t)N_{j}^{i}({\boldsymbol {x}}),\quad \forall {\boldsymbol {x}}\in \Omega _{e_{i}}.}$

Then similarly choosing a test function

${\displaystyle \varphi _{h}({\boldsymbol {x}})=\sum _{j=1}^{\text{dofs}}\varphi _{j}^{i}N_{j}^{i}({\boldsymbol {x}}),\quad \forall {\boldsymbol {x}}\in \Omega _{e_{i}},}$

multiplying the continuity equation by ${\displaystyle \varphi _{h}}$ and integrating by parts in space, the semidiscrete DG formulation becomes:

${\displaystyle {\frac {d}{dt}}\int _{\Omega _{e_{i}}}\rho _{h}\varphi _{h}\,d{\boldsymbol {x}}+\int _{\partial \Omega _{e_{i}}}\varphi _{h}\mathbf {J} _{h}\cdot {\boldsymbol {n}}\,d{\boldsymbol {x}}=\int _{\Omega _{e_{i}}}\mathbf {J} _{h}\cdot \nabla \varphi _{h}\,d{\boldsymbol {x}}.}$

## Scalar hyperbolic conservation law

A scalar hyperbolic conservation law is of the form

{\displaystyle {\begin{aligned}\partial _{t}u+\partial _{x}f(u)&=0\quad {\text{for}}\quad t>0,\,x\in \mathbb {R} \\u(0,x)&=u_{0}(x)\,,\end{aligned}}}

where one tries to solve for the unknown scalar function ${\displaystyle u\equiv u(t,x)}$, and the functions ${\displaystyle f,u_{0}}$ are typically given.

### Space discretization

The ${\displaystyle x}$-space will be discretized as

${\displaystyle \mathbb {R} =\bigcup _{k}I_{k}\,,\quad I_{k}:=\left(x_{k},x_{k+1}\right)\quad {\text{for}}\quad x_{k}

Furthermore, we need the following definitions

${\displaystyle h_{k}:=|I_{k}|\,,\quad h:=\sup _{k}h_{k}\,,\quad {\hat {x}}_{k}:=x_{k}+{\frac {h_{k}}{2}}\,.}$

### Basis for function space

We derive the basis representation for the function space of our solution ${\displaystyle u}$. The function space is defined as

${\displaystyle S_{h}^{p}:=\left\lbrace v\in L^{2}(\mathbb {R} ):v{\Big |}_{I_{k}}\in \Pi _{p}\right\rbrace \quad {\text{for}}\quad p\in \mathbb {N} _{0}\,,}$

where ${\displaystyle {v|}_{I_{k}}}$ denotes the restriction of ${\displaystyle v}$ onto the interval ${\displaystyle I_{k}}$, and ${\displaystyle \Pi _{p}}$ denotes the space of polynomials of maximal degree ${\displaystyle p}$. The index ${\displaystyle h}$ should show the relation to an underlying discretization given by ${\displaystyle \left(x_{k}\right)_{k}}$. Note here that ${\displaystyle v}$ is not uniquely defined at the intersection points ${\displaystyle (x_{k})_{k}}$.

At first we make use of a specific polynomial basis on the interval ${\displaystyle [-1,1]}$, the Legendre polynomials ${\displaystyle (P_{n})_{n\in \mathbb {N} _{0}}}$, i.e.,

${\displaystyle P_{0}(x)=1\,,\quad P_{1}(x)=x\,,\quad P_{2}(x)={\frac {1}{2}}(3x^{2}-1)\,,\quad \dots }$

Note especially the orthogonality relations

${\displaystyle \left\langle P_{i},P_{j}\right\rangle _{L^{2}([-1,1])}={\frac {2}{2i+1}}\delta _{ij}\quad \forall \,i,j\in \mathbb {N} _{0}\,.}$

Transformation onto the interval ${\displaystyle [0,1]}$, and normalization is achieved by functions ${\displaystyle (\varphi _{i})_{i}}$

${\displaystyle \varphi _{i}(x):={\sqrt {2i+1}}P_{i}(2x-1)\quad {\text{for}}\quad x\in [0,1]\,,}$

which fulfill the orthonormality relation

${\displaystyle \left\langle \varphi _{i},\varphi _{j}\right\rangle _{L^{2}([0,1])}=\delta _{ij}\quad \forall \,i,j\in \mathbb {N} _{0}\,.}$

Transformation onto an interval ${\displaystyle I_{k}}$ is given by ${\displaystyle \left({\bar {\varphi }}_{ki}\right)_{i}}$

${\displaystyle {\bar {\varphi }}_{ki}:={\frac {1}{\sqrt {h_{k}}}}\varphi _{i}\left({\frac {x-x_{k}}{h_{k}}}\right)\quad {\text{for}}\quad x\in I_{k}\,,}$

which fulfill

${\displaystyle \left\langle {\bar {\varphi }}_{ki},{\bar {\varphi }}_{kj}\right\rangle _{L^{2}(I_{k})}=\delta _{ij}\quad \forall \,i,j\in \mathbb {N} _{0}\forall \,k\,.}$

For ${\displaystyle L^{\infty }}$-normalization we define ${\displaystyle \varphi _{ki}:={\sqrt {h_{k}}}{\bar {\varphi }}_{ki}}$, and for ${\displaystyle L^{1}}$-normalization we define ${\displaystyle {\tilde {\varphi }}_{ki}:={\frac {1}{\sqrt {h_{k}}}}{\bar {\varphi }}_{ki}}$, s.t.

${\displaystyle \|\varphi _{ki}\|_{L^{\infty }(I_{k})}=\|\varphi _{i}\|_{L^{\infty }([0,1])}=:c_{i,\infty }\quad {\text{and}}\quad \|{\tilde {\varphi }}_{ki}\|_{L^{1}(I_{k})}=\|\varphi _{i}\|_{L^{1}([0,1])}=:c_{i,1}\,.}$

Finally, we can define the basis representation of our solutions ${\displaystyle u_{h}}$

{\displaystyle {\begin{aligned}u_{h}(t,x):=&\sum _{i=0}^{p}u_{ki}(t)\varphi _{ki}(x)\quad {\text{for}}\quad x\in (x_{k},x_{k+1})\\u_{ki}(t)=&\left\langle u_{h}(t,\cdot ),{\tilde {\varphi }}_{ki}\right\rangle _{L^{2}(I_{k})}\,.\end{aligned}}}

Note here, that ${\displaystyle u_{h}}$ is not defined at the interface positions.

Besides, prism bases are employed for planar-like structures, and are capable for 2-D/3-D hybridation.

### DG-scheme

The conservation law is transformed into its weak form by multiplying with test functions, and integration over test intervals

{\displaystyle {\begin{aligned}\partial _{t}u+\partial _{x}f(u)&=0\\\Rightarrow \quad \left\langle \partial _{t}u,v\right\rangle _{L^{2}(I_{k})}+\left\langle \partial _{x}f(u),v\right\rangle _{L^{2}(I_{k})}&=0\quad {\text{for}}\quad \forall \,v\in S_{h}^{p}\\\Leftrightarrow \quad \left\langle \partial _{t}u,{\tilde {\varphi }}_{ki}\right\rangle _{L^{2}(I_{k})}+\left\langle \partial _{x}f(u),{\tilde {\varphi }}_{ki}\right\rangle _{L^{2}(I_{k})}&=0\quad {\text{for}}\quad \forall \,k\;\forall \,i\leq p\,.\end{aligned}}}

By using partial integration one is left with

{\displaystyle {\begin{aligned}{\frac {\mathrm {d} }{\mathrm {d} t}}u_{ki}(t)+f(u(t,x_{k+1})){\tilde {\varphi }}_{ki}(x_{k+1})-f(u(t,x_{k})){\tilde {\varphi }}_{ki}(x_{k})-\left\langle f(u(t,\,\cdot \,)),{\tilde {\varphi }}_{ki}'\right\rangle _{L^{2}(I_{k})}=0\quad {\text{for}}\quad \forall \,k\;\forall \,i\leq p\,.\end{aligned}}}

The fluxes at the interfaces are approximated by numerical fluxes ${\displaystyle g}$ with

${\displaystyle g_{k}:=g(u_{k}^{-},u_{k}^{+})\,,\quad u_{k}^{\pm }:=u(t,x_{k}^{\pm })\,,}$

where ${\displaystyle u_{k}^{\pm }}$ denotes the left- and right-hand sided limits. Finally, the DG-Scheme can be written as

{\displaystyle {\begin{aligned}{\frac {\mathrm {d} }{\mathrm {d} t}}u_{ki}(t)+g_{k+1}{\tilde {\varphi }}_{ki}(x_{k+1})-g_{k}{\tilde {\varphi }}_{ki}(x_{k})-\left\langle f(u(t,\,\cdot \,)),{\tilde {\varphi }}_{ki}'\right\rangle _{L^{2}(I_{k})}=0\quad {\text{for}}\quad \forall \,k\;\forall \,i\leq p\,.\end{aligned}}}

## Scalar elliptic equation

A scalar elliptic equation is of the form

{\displaystyle {\begin{aligned}-\partial _{xx}u&=f(x)\quad {\text{for}}\quad x\in (a,b)\\u(x)&=g(x)\,\quad {\text{for}}\,\quad x=a,b\end{aligned}}}

This equation is the steady-state heat equation, where ${\displaystyle u}$ is the temperature. Space discretization is the same as above. We recall that the interval ${\displaystyle (a,b)}$ is partitioned into ${\displaystyle N+1}$ intervals of length ${\displaystyle h}$.

We introduce jump ${\displaystyle [{}\cdot {}]}$ and average ${\displaystyle \{{}\cdot {}\}}$ of functions at the node ${\displaystyle x_{k}}$:

${\displaystyle [v]{\Big |}_{x_{k}}=v(x_{k}^{+})-v(x_{k}^{-}),\quad \{v\}{\Big |}_{x_{k}}=0.5(v(x_{k}^{+})+v(x_{k}^{-}))}$

The interior penalty discontinuous Galerkin (IPDG) method is: find ${\displaystyle u_{h}}$ satisfying

${\displaystyle A(u_{h},v_{h})+A_{\partial }(u_{h},v_{h})=\ell (v_{h})+\ell _{\partial }(v_{h})}$

where the bilinear forms ${\displaystyle A}$ and ${\displaystyle A_{\partial }}$ are

${\displaystyle A(u_{h},v_{h})=\sum _{k=1}^{N+1}\int _{x_{k-1}}^{x_{k}}\partial _{x}u_{h}\partial _{x}v_{h}-\sum _{k=1}^{N}\{\partial _{x}u_{h}\}_{x_{k}}[v_{h}]_{x_{k}}+\varepsilon \sum _{k=1}^{N}\{\partial _{x}v_{h}\}_{x_{k}}[u_{h}]_{x_{k}}+{\frac {\sigma }{h}}\sum _{k=1}^{N}[u_{h}]_{x_{k}}[v_{h}]_{x_{k}}}$

and

${\displaystyle A_{\partial }(u_{h},v_{h})=\partial _{x}u_{h}(a)v_{h}(a)-\partial _{x}u_{h}(b)v_{h}(b)-\varepsilon \partial _{x}v_{h}(a)u_{h}(a)+\varepsilon \partial _{x}v_{h}(b)u_{h}(b)+{\frac {\sigma }{h}}{\big (}u_{h}(a)v_{h}(a)+u_{h}(b)v_{h}(b){\big )}}$

The linear forms ${\displaystyle \ell }$ and ${\displaystyle \ell _{\partial }}$ are

${\displaystyle \ell (v_{h})=\int _{a}^{b}fv_{h}}$

and

${\displaystyle \ell _{\partial }(v_{h})=-\varepsilon \partial _{x}v_{h}(a)g(a)+\varepsilon \partial _{x}v_{h}(b)g(b)+{\frac {\sigma }{h}}{\big (}g(a)v_{h}(a)+g(b)v_{h}(b){\big )}}$

The penalty parameter ${\displaystyle \sigma }$ is a positive constant. Increasing its value will reduce the jumps in the discontinuous solution. The term ${\displaystyle \varepsilon }$ is chosen to be equal to ${\displaystyle -1}$ for the symmetric interior penalty Galerkin method; it is equal to ${\displaystyle +1}$ for the non-symmetric interior penalty Galerkin method.