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);
u(:)=f(x(:));
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.
