Consider the following partial differential equation \begin{align} \frac{\partial u}{\partial t}+\frac{\partial f}{\partial x} &= g(x,t), \ \ x\in \Omega = [x_{L},x_{R}] \\ u(x,0) &= u_{0}(x) \\ u(x_{L},t) &= g_{1}(t) \end{align} where the flux function $f(u,x)$ is nonlinear. Let $\{x_{i}\}_{i=1}^{N+1}$ be a partition of $\Omega$ with $N$ elements such that \begin{equation} x_{L} = x_{1}<x_{2}<\cdots<x_{N}<x_{N+1} = x_{R} \end{equation} and let $D_{k} = [x_{k},x_{k+1}]$. Consider the weak formulation of the problem \begin{equation} \int_{D_{k}} \Bigg[ \frac{\partial u}{\partial t}\phi_{i}^{k} -f(u,x)\frac{d \phi_{i}^{k}}{d x}\Bigg] = \int_{D_{k}} g(x,t)\phi_{i}^{k}-\Big[f^{*}\phi_{i}^{k} \Big]_{x_{k}}^{x_{k+1}} \end{equation} where $f^{*}$ is called the numerical flux. If we consider approximations with polynomials of degree $p$, the approximation $u_{h}^{k}$ of $u$ in the element $D_{k}$ is given by \begin{equation} u_{h}^{k} = \sum_{j=1}^{p+1} \alpha_{j}^{k}\phi_{j}^{k} \end{equation} With the contribution of all elements, we get the following ODE \begin{equation} M \dot{\alpha}-r = b-l \end{equation} where the entries of the vectors $b$ and $r$ in the element $D_{k}$ are given by \begin{align} [b_{k}]_{i} &= \int_{D_{k}} g(x,t)\phi_{i}^{k} dx \\ [r_{k}]_{i} &= \int_{D_{k}} f(u,x)\frac{d \phi_{i}^{k}}{d x} dx \end{align} The Lax–Friedrichs numerical flux is given by \begin{equation} f^{*} = \{f(u)\}+\frac{C}{2}[\hspace{-0.6mm}[u]\hspace{-0.6mm}] \end{equation} where $C = \max{ |f'(u)| }$, $\{\cdot\}$ is the average and $[\hspace{-0.6mm}[\cdot]\hspace{-0.6mm}]$ is the jump.
The numerical flux vector $l$ in the element $k$ has entries \begin{align} \Big[f^{*}\phi_{i}^{k} \Big]_{x_{k}}^{x_{k+1}} =& \Bigg[ \frac{1}{2}f(u_{h}^{k},x_{k+1})+ \frac{1}{2}f(u_{h}^{k+1},x_{k+1})+\frac{C}{2}u_{h}^{k}-\frac{C}{2}u_{h}^{k+1} \Bigg] \phi_{i}^{k} (x_{k+1}) \\&-\Bigg[ \frac{1}{2}f(u_{h}^{k-1},x_{k})+ \frac{1}{2}f(u_{h}^{k},x_{k})+\frac{C}{2}u_{h}^{k-1}-\frac{C}{2}u_{h}^{k} \Bigg] \phi_{i}^{k} (x_{k}) \end{align} For $k=1$, we use the Dirichlet datum, i.e, $u_{h}^{0} = g_{1}(t)$. For $2\leq k \leq N-1$ we just follow the equation.

But I have problem when $k=N$, the reason is that I don't have $u_{h}^{N+1}$ in the following equation \begin{align} \Big[f^{*}\phi_{i}^{N} \Big]_{x_{N}}^{x_{N+1}} =& \Bigg[ \frac{1}{2}f(u_{h}^{N},x_{N+1})+ \frac{1}{2}f(u_{h}^{N+1},x_{N+1})+\frac{C}{2}u_{h}^{k}-\frac{C}{2}u_{h}^{N+1} \Bigg] \phi_{i}^{N} (x_{N+1}) \\&-\Bigg[ \frac{1}{2}f(u_{h}^{N-1},x_{k})+ \frac{1}{2}f(u_{h}^{N},x_{N})+\frac{C}{2}u_{h}^{N-1}-\frac{C}{2}u_{h}^{N} \Bigg] \phi_{i}^{N} (x_{N}) \end{align}
I have tried many things; first I take $\{f(u)\}\big|_{x_{N+1}} = f(u_{h}^{N})$ and $[\hspace{-0.6mm}[u]\hspace{-0.6mm}]\big|_{x_{N+1}} = u_{h}^{N}$ in the equation \begin{equation} f^{*} = \{f(u)\}+\frac{C}{2}[\hspace{-0.6mm}[u]\hspace{-0.6mm}] \end{equation} but it does not work.

Second, I know that this is not allowed but if I take $ u_{h}^{N+1} = g_{2}(t)$, i.e, if I consider a second Dirichlet datum in $x_{R}$, my implementation works fine. Clearly this only works if I know the solution $u$. Because in the case that I don't know the solution $u$, I just know $g_{1}(t)$.

My questions are:

  1. Do you know books or papers that explain how to solve this PDE with the DG method?
  2. Do you know how should I take $\{f(u)\}$, $[\hspace{-0.6mm}[u]\hspace{-0.6mm}]$ in the element $k=N$?

1 Answer 1


The Lax-Friedrichs flux roughly approximates a situation where the propagation of information can occur in both directions. This is why on the right face at $x_{N+1}$ you need to also know $u^{N+1}_h$ in order to compute this numerical flux.

As an alternative, you could use an upwinding flux (either just at the boundaries or in your entire domain).

These fluxes have the property that they only need information from one side of the interface. For example, take $$ f^* n^- = \begin{cases} f^- n^- & a^* n^- \ge 0\\ f^+ n^- & a^* n^- < 0 \end{cases}\\ a^* = \begin{cases} \frac{[[f]]}{[[u]]} & [[u]] \neq 0\\ \partial_u f^- & [[u]] = 0 \end{cases} $$ where $n^-$ is the outwards pointing normal, $f^-$ is $f(u_h)$ evaluated inside the element $k$ at the face, and $f^+$ is $f(u_h)$ evaluated outside the element the element $k$ at the face.

Notice that this is almost what you tried but didn't work, except it removes the jump penalty term entirely and decides whether to use $f^-$ or $f^+$ depending on if the wave is moving in or out.

If you have some insight into what $f(u,x)$ actually is, you might be able to determine the direction of the upwind flux even without needing to calculate $a^*$.

Upwind fluxes only work if you know at a face that waves only propagate in a single direction. I think this will always be true for your PDE, but there are other PDE systems where this is not true (ex.: the wave equation or Maxwell's equations).

Regardless of the choice of numerical flux you use, you might still need to specify boundary conditions on both sides if the wave speed on both domain boundaries are directed inwards.

For example, suppose you have $$f(u,x) = u \cos\left(\frac{\pi(x - x_L)}{x_R-x_L}\right)$$ You would need to know $u^{N+1}_h$ to solve this problem because you are advecting in $u$ from the right of your domain, and would also need to know $u^{0}_h$ since you are also advecting in $u$ from the left of your domain.

Useful resources:


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.