I am trying to implement a second order DG-method using a monomial basis and explicit Euler in time. I have written down some of the theory, which I present below:

Theory Consider the linear advection function over $\Omega\subset\mathbb{R}$ and a class of testfunctions $v\in C_0^1(\Omega)$ $$ \int_\Omega\partial_tuvd x+a\int_\Omega\partial_xuvd x=0\;.\hspace{70pt}(1) $$ Now, divide the computational domain $\Omega=\cup_{i=1}^N\Omega_i$ where $\Omega_i$ are disjoint, and let $h_i=x_{i+1}-x_{i}$ be the cell length. Integration by parts yield of $(1)$ yields $$ \int_{\Omega_i}\partial_tuv d x-a\int_{\Omega_i}uv' d x+\left[auv\right]_{x_i}^{x_{i+1}}=0\;.\hspace{30pt}(2) $$ The exact solution $u$ is now be approximated with $U_i$ over the cell $\Omega_i$ in a finite dimensional subset of $C_0^1(\Omega_i)$. Then the approximation of $u$ in cell $\Omega_i$ can be written $$ U_i=\sum_{j=0}^pc_{ij}\phi_j $$ where $c_{ij}$ are time-depentent coefficients. In order to get a Galerkin method the basis functions and the testfunctions should be in the same space. For simplicity let this be the space of polynomials of degree $p$, i.e. $U_i,v\in\Pi_p$. For $v=\phi_k$, with $k=0,\ldots,p$, inserting into $(2)$ yields $$ \int_{\Omega_i}\partial_t\sum_{j=0}^pc_{ij}\phi_j\phi_k dx-a\int_{\Omega_i}\sum_{j=0}^pc_{ij}\phi_j\phi_k' d x+\left[a\sum_{j=0}^pc_{ij}\phi_j\phi_k\right]_{x_i}^{x_{i+1}}=0\;.\hspace{10pt}(3) $$ Let us define a suitable inner product on $\Omega_i$ following the $L^2$-space inner product definition $$ (u,v)=\int_{x_L}^{x_R}uvd x\;. $$ Due to the smoothness properties of the testfunctions, the derivatives and the sums can be moved out from under the integral sign, and we can rewrite $(3)$ as $$ \sum_{j=0}^p\left[\partial_tc_{ij}(\phi_j,\phi_k)-ac_{ij}(\phi_j,\phi_k')\right]+\left[a\sum_{j=0}^pc_{ij}\phi_j\phi_k\right]_{x_i}^{x_{i+1}}=0\;.\hspace{30pt}(4) $$ Let us now make some simplifications. First let $h=h_i$ for all cells, i.e. an equidistant grid. Secondly, restrict $p=1$ with the monomial basis $$ \phi_0=1,\qquad\phi_1=\frac{x-x_L}{x_R-x_L} $$ giving $$ (\phi_0,\phi_0)=h,\quad(\phi_0,\phi_1)=\frac{h}{2},\quad(\phi_1,\phi_1)=\frac{h}{3},\\\quad(\phi_0,\phi_0')=0,\quad(\phi_0,\phi_1')=1, \quad(\phi_1,\phi_1')=\frac{1}{2}\;. $$

We are now ready to write this as a system $$ \partial_tMc_i-aM^fc_i+\left[M^s\overline{f}^N\right]_{x_i}^{x_{i+1}}=0 $$ where $M=(m_{jk})$ where $m_{jk}=(\phi_j,\phi_k)$ and $M^f=(m^f_{jk})$ where $m^f_{jk}=(\phi_j,\phi_k')$ and $$ M^s=\begin{pmatrix}\phi_0(x)&0\\0&\phi_1(x)\end{pmatrix},\qquad\overline{f}^N=\begin{pmatrix}f^N(x)\\f^N(x)\end{pmatrix},\qquad c_i=\begin{pmatrix}c_{i,0}(t)\\c_{i,1}(t)\end{pmatrix} $$ for some numerical flux function $f^N$. We will use the upwind flux, i.e. $$ f^N(U_{i-1},U_i)=aU_{i-1}(x_{i-1/2}) $$ The last step to complete the scheme is to apply the explicit Euler method in time. First rewrite as $$ \partial_t c_i=aM^{-1}M^fc_i-M^{-1}\left[M^s\overline{f}^N\right]_{x_i}^{x_{i+1}} $$ We get $$ c_i^{n+1}=c_i^n+\Delta t\left(aM^{-1}M^fc_i^n-M^{-1}\left[M^s\overline{f}^N\right]_{x_i}^{x_{i+1}}\right) $$

Implementaiton The following is in Matlab

Assume we have an initial pulse $f(x)=e^{-x^2}$ over $[-4,4]$. Then I extract the initial coefficient using $$ u(x_L)=c_0\qquad u(x_R)=c_0+c_1 $$ i.e.

V =[1 0;-1 1];       
u=zeros(2, Ncells);
cini = V*u;

where $x$ is a matrix containing the the end points of each cell. Then I run the above, but the solution diverge quickly.

Questions Does the theory seem okay, or are there any flaws? I can't figure out where it goes wrong. I can live with it if the implementation is a little off, but as long as I am not trying to program something that won't work it would be a relief.


Just glancing through it quickly, it seems that you are using a second-order discretization in space with forward Euler in time. That is unstable for any fixed CFL number (on fine enough grids). You should use, e.g., a higher-order Runge-Kutta method. It is possible to check the stability by computing the eigenvalues of your semi-discretization and using method of lines stability analysis.

  • $\begingroup$ Thank you for answering. I know that the central differences with explicit Euler in time is unconditionally unstable, but there are other methods like Lax-Wendroff, Beam-Warming that are second order and stable for some CFL condition. So, this must mean that the DG-method above is a central difference scheme essentially, but I don't find it very transparent. But, essentially changing to Heun's method would then be stable I assume. How can I check the stability via von Neumann analysis? I am used to FD-schemes and going to FV- or DG-schemes is not straight-forward. $\endgroup$ – user136475 Dec 13 '14 at 9:30
  • $\begingroup$ One of your comments is correct; I have accordingly changed my answer. $\endgroup$ – David Ketcheson Dec 13 '14 at 10:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.