Numerical Gauss-Lobatto quadrature in DG and instability

Question

I am trying to implement the DG method for the solution of non-linear hyperbolic problems. The aim is mainly educational. I want to make sure I understand exactly what is happening and not just using some ready-made toolbox and think that I know what is happening when in fact I do not. Therefore, I am writing the complete code from scratch.

I started with the linear advection equation and I managed to get a working solution for the nodal approximation for both Gauss-Legendre and Gauss-Lobatto nodes with Lagrange polynomials as the basis function. A naive transition to the non-linear Burgers' equation, where numerical quadratures were used for the stiffness matrix, was unsuccessful so I did the numerical quadrature version for the linear advection as well. It is not at all stable either. I would like to understand the reason for that and find out the solution for a stable numerical method.

I do understand that the Gauss-Lobato is only accurate for polynomials of order $2n-3$ and I do understand that the diagonal mass matrix resulting from such quadrature is not exact. The polynomial order is $k=2$, the number of GL points is $n=3$. The integrands in the mass matrix are of order $4$ and in the stiffness matrix of order $3$.

From that I would assume that the stiffness matrix should be evaluated exactly in case of the linear advection, but I am getting a different form.

The scheme is: $$ M \frac{\partial y}{\partial t} = K y - (f_r - f_l) / dx $$

together with the RK SSP43 solver in time. Numerical fluxes use the upwind flux scheme.

When using these matrices, derived in Mathematica by analytical integration, the scheme works well:

$$ M = \begin{pmatrix} 4 & 2 & -1\\ 2 & 16 & 2 \\ -1 & 2 & 4 \end{pmatrix}/15\\ K = \begin{pmatrix} -3 & -4 & 1\\ 4 & 0 & -4 \\ -1 & 4 & 3 \end{pmatrix} / 6\\ $$

With the numerically-derived stiffness matrix it produces very oscillatory solutions after a very short time, no matter which mass-matrix is used:

$$ M = \begin{pmatrix} 1 & 0 & 0\\ 0 & 3 & 0 \\ 0 & 0 & 1 \end{pmatrix}/4\\ K = \begin{pmatrix} -3 & -3 & 1\\ 4 & 0 & -4 \\ -1 & 3 & 3 \end{pmatrix} / 8\\ $$

---

My Mathematica code for the derivations is probably not very elegant, but I can share it too

In[4]:= Clear["Global`*"]

In[5]:= a = ({
   {a1},
   {a2},
   {a3}
  })

Out[5]= {{a1}, {a2}, {a3}}

In[7]:= nodes = {-1, 0, 1}

Out[7]= {-1, 0, 1}

In[8]:= \[Phi][\[Xi]_] = ({
    {InterpolatingPolynomial[{{nodes[[1]], 1}, {nodes[[2]], 
        0}, {nodes[[3]], 0}}, \[Xi]]},
    {InterpolatingPolynomial[{{nodes[[1]], 0}, {nodes[[2]], 
        1}, {nodes[[3]], 0}}, \[Xi]]},
    {InterpolatingPolynomial[{{nodes[[1]], 0}, {nodes[[2]], 
        0}, {nodes[[3]], 1}}, \[Xi]]}
   }) // FullSimplify

Out[8]= {{1/2 (-1 + \[Xi]) \[Xi]}, {1 - \[Xi]^2}, {1/
   2 \[Xi] (1 + \[Xi])}}

In[9]:= \[Phi][\[Xi]_] = {InterpolatingPolynomial[{{nodes[[1]], 
      1}, {nodes[[2]], 0}, {nodes[[3]], 0}}, \[Xi]], 
   InterpolatingPolynomial[{{nodes[[1]], 0}, {nodes[[2]], 
      1}, {nodes[[3]], 0}}, \[Xi]], 
   InterpolatingPolynomial[{{nodes[[1]], 0}, {nodes[[2]], 
      0}, {nodes[[3]], 1}}, \[Xi]]} // FullSimplify

Out[9]= {1/2 (-1 + \[Xi]) \[Xi], 1 - \[Xi]^2, 1/2 \[Xi] (1 + \[Xi])}

In[10]:= D\[Phi][\[Xi]_] = D[\[Phi][\[Xi]], \[Xi]] // FullSimplify

Out[10]= {-(1/2) + \[Xi], -2 \[Xi], 1/2 + \[Xi]}

In[11]:= M = Integrate[TensorProduct[\[Phi][x], \[Phi][x]], {x, -1, 1}]

Out[11]= {{4/15, 2/15, -(1/15)}, {2/15, 16/15, 2/15}, {-(1/15), 2/15, 
  4/15}}

In[12]:= S = 
 Integrate[TensorProduct[\[Phi][x], D\[Phi][x]], {x, -1, 1}] // 
  FullSimplify

Out[12]= {{-(1/2), 2/3, -(1/6)}, {-(2/3), 0, 2/3}, {1/6, -(2/3), 1/2}}

In[13]:= S2 = S * 6

Out[13]= {{-3, 4, -1}, {-4, 0, 4}, {1, -4, 3}}

In[14]:= K2 = Transpose[S2]

Out[14]= {{-3, -4, 1}, {4, 0, -4}, {-1, 4, 3}}

In[15]:= {{-3, -4, 1}, {4, 0, -4}, {-1, 4, 3}}


Out[15]= {{-3, -4, 1}, {4, 0, -4}, {-1, 4, 3}}

In[16]:= TPM = TensorProduct[\[Phi][x], \[Phi][x]]

Out[16]= {{1/4 (-1 + x)^2 x^2, 1/2 (-1 + x) x (1 - x^2), 
  1/4 (-1 + x) x^2 (1 + x)}, {1/2 (-1 + x) x (1 - x^2), (1 - x^2)^2, 
  1/2 x (1 + x) (1 - x^2)}, {1/4 (-1 + x) x^2 (1 + x), 
  1/2 x (1 + x) (1 - x^2), 1/4 x^2 (1 + x)^2}}

In[17]:= Mnum = 
 1/4 * (TPM /. x -> -1) + 3/4*(TPM /. x -> 0) + 1/4*(TPM /. x -> 1)

Out[17]= {{1/4, 0, 0}, {0, 3/4, 0}, {0, 0, 1/4}}

In[19]:= TP = TensorProduct[\[Phi][x], D\[Phi][x]]

Out[19]= {{1/2 (-1 + x) (-(1/2) + x) x, -((-1 + x) x^2), 
  1/2 (-1 + x) x (1/2 + x)}, {(-(1/2) + x) (1 - x^2), -2 x (1 - 
     x^2), (1/2 + x) (1 - x^2)}, {1/
   2 (-(1/2) + x) x (1 + x), -x^2 (1 + x), 1/2 x (1/2 + x) (1 + x)}}

In[20]:= Snum = 
 1/4 * (TP /. x -> -1) + 3/4*(TP /. x -> 0) + 1/4*(TP /. x -> 1)

Out[20]= {{-(3/8), 1/2, -(1/8)}, {-(3/8), 0, 3/8}, {1/8, -(1/2), 3/8}}

In[21]:= Snum2 = Snum*8

Out[21]= {{-3, 4, -1}, {-3, 0, 3}, {1, -4, 3}}

In[22]:= Knum = Transpose[Snum2]

Out[22]= {{-3, -3, 1}, {4, 0, -4}, {-1, 3, 3}}
```

Answer 1

About stability

A spatial semi-discretization of the linear advection equation $$ y_t + y_x = 0 $$ using an upwind flux might look something like $$ M \mathbf{y}_t + K \mathbf{y} = -(y_0 - y^*) \mathbf{e}_0. $$ Here, $\mathbf{y} = (y_0, \dots, y_n)^\top$ is the discrete solution, $M$ and $K$ are the so-called mass and stiffness matrices, and $\mathbf{e}_0 = (1,0,\dots,0)^\top \in \mathbb{R}^{n+1}$. This would be the discretization on a single element, and $y^*$ is the data entering this element from its neighbour to the left. Alternatively, it could be data from a Dirichlet condition at the left boundary.

The stability of such a discretization has nothing to do with finite difference, finite volume, finite element/DG methods. It is a consequence of the properties of $M$ and $K$ only. To see this, begin by exploring the stability properties inherent to the PDE: Multiply by $y$ and integrate over some domain $[a,b]$ to obtain $$ \int_a^b yy_t \text{d}x = -\int_a^b yy_x \text{d}x. $$ With the notation $\| y \| = \sqrt{\int_a^b y^2 \text{d}x}$ it should be clear that the left-hand side is the time-derivative of $\frac{1}{2}\| y \|^2$. Carrying out the integration (by parts) of the right-hand side thus gives $$ \| y \|^2_t = y^2(a) - y^2(b) \leq y^2(a). $$ To keep the solution bounded, some kind of condition must therefore be imposed on $y(a)$.

Returning to the discrete setting, consider momentarily what happens in the absence of the numerical flux: I will take the liberty of bluntly stating the properties we should impose on $M$ and $K$ and motivate them later. Thus, we demand the following:

• $M$ is symmetric positive definite,
• $K + K^\top = \mathbf{e}_n \mathbf{e}_n^\top - \mathbf{e}_0 \mathbf{e}_0^\top = \text{diag}(-1,0,\dots,0,1)$.

Consider $M \mathbf{y}_t + K \mathbf{y} = 0$ and left-multiply by $\mathbf{y}^\top$, then add the transpose of the result. This gives $$ \mathbf{y}^\top M \mathbf{y}_t + \mathbf{y}_t^\top M \mathbf{y} = - \mathbf{y}^\top (K + K^\top) \mathbf{y}. $$ With the notation $\| \mathbf{y} \|_M = \sqrt{\mathbf{y}^\top M \mathbf{y}}$ it should be clear that the left-hand side is the time-derivative of $\| \mathbf{y} \|_M^2$. Using the listed property of $K$, we find that $$ \frac{\text{d}}{\text{d}t} \| \mathbf{y} \|_M^2 = y_0^2 - y_n^2 \leq y_0^2. $$ This is exactly the same behaviour that we saw for the continuous problem. Again, some kind of condition is necessary on $y_0$, and this is where the upwind flux enters the picture. I leave it to you as an exercise to show that the same steps as outlined above will lead to a contribution from the upwind flux of the form $(y^*)^2 - (y_0 - y^*)^2$.

What did we do?

We demanded that the discretization should have the same properties as the PDE. Why we should consider the particular norm $\| y \|$ is a bit subtle. For this discussion, it suffices to understand that if $\| y \|$ is bounded, then so is $y$. This was our stability statement. To reach it, all we needed was integration by parts.

Demanding that $M$ be symmetric positive definite allowed us to define a discrete norm "similar" to the continuous one. In fact, we could introduce the matrix $D = M^{-1}K$ and note that $D$ is an approximation of the operator $\frac{\partial}{\partial x}$. The condition we imposed on $K$ can then be restated in terms of the inner product $(\mathbf{f}, \mathbf{g})_M = \mathbf{f}\top M \mathbf{g}$ as $$ (\mathbf{f}, D \mathbf{g})_M = f_n g_n - f_0 g_0 - (D \mathbf{f}, \mathbf{g})_M. $$ Compare this to integration by parts: Let $(f,g) = \int_a^b f g \text{d}x$. Then $$ (f, \partial g / \partial x) - f(b) g(b) - f(a) g(a) - (\partial f / \partial x, g). $$ In words, we selected a discretization that preserved integration by parts, hence the inherent stability of the continuous solution. From this alone it is possible to show that $M$ approximates integrals; it defines a quadrature rule.

What about numerical methods?

Numerical analysis is a branch of mathematics and mathematics is the study of mathematical objects. A mathematical object, roughly speaking, is something that satisfies certain well-defined mathematical properties. And what a mathematical property is, well ... The point is that

the finite difference/volume/element/... methods are not mathematical objects

While they take their starting points from certain mathematical considerations, there is no limit to the tweaks and tricks and shenanigans you can add to any of these methods and still call it a method. They have no inherent properties. They are not mathematically meaningful objects of study. So what do they do? They provide you with a recipe for assigning numbers to the elements in the matrices $M$ and $K$. There is absolutely no reason to expect that they will provide you matrices with desirable properties.

However, as we know from experience, sometimes they do. There are ways to be intelligent about this, and that is what numerical analysists try to be.

Discontinuous Galerkin Use Legendre-Gauss-Lobatto nodes, and the Gauss-Lobatto quadrature rule, and express the solution in terms of Lagrange interpolation. If you do your calculations right, you will end up with the matrices $$ M = \frac{b-a}{2}\begin{pmatrix} 1 & & \\ & 4 & \\ & & 1 \end{pmatrix}/3, \qquad K = \pmatrix{ -1/2 & 2/3 & -1/6 \\ -2/3 & 0 & 2/3 \\ 1/6 & -2/3 & 1/2 }. $$ They possess the desired properties. Note that $M$ must scale with the size of the element. It does this since it approximates an integral. $D$ approximates a derivative so it would scale as $2/(b-a)$. $K$ is scale invariant.

Finite differences Consider the trapezoidal rule and centered differences: $$ M = \frac{b-a}{n} \begin{pmatrix} \frac{1}{2} & & & & \\ & 1 & & & \\ & & \ddots & & \\ & & & 1 & \\ & & & & \frac{1}{2} \end{pmatrix}, \qquad K = \begin{pmatrix} -\frac{1}{2} & \frac{1}{2} & & & \\ -\frac{1}{2} & 0 & \frac{1}{2} & & \\ & & \ddots & & \\ & & -\frac{1}{2} & 0 & \frac{1}{2} \\ & & & -\frac{1}{2} & \frac{1}{2} \end{pmatrix}. $$ They have the same desired properties and are thus stable by exactly the same token as the DG method.

Traditionally finite difference methods have been considered "structured" methods that require some type of regular grid. This is a misunderstanding: All methods are structured/unstructured in the same way. We can comfortably use the $M$ and $K$ above within a multi-element setting and distribute those elements as we like. "Structure" is a demand on the grid points within an element. Indeed, in DG the Legendre-Gauss-Lobatto nodes are "structured". It's just that finite difference methods traditionally were used on a single element.

Finite volumes At the opposite end of the spectrum, finite volume methods use many elements with only a single grid point within those elements (i.e. $n=0$). Indeed, consider $$ M = (b-a), \qquad K = 0. $$ Here we have to demand that $\mathbf{e}_0 = \mathbf{e}_n = 1$ when $n=0$. Again, stability follows for the same reasons as before.

Summary of main points

• Study properties of your PDE, then mimic those properties in the discrete setting. How to do that is a mathematical question.
• Numerical methods are not mathematical objects. There is an infinite number of methods that will return completely meaningless approximations to any given problem. You can freely call any of them "finite ... method" if you like. Who's to tell you differently?
• Numerical methods provide suggestions for how to approximate the operators present in your problem (derivatives, integrals, etc.) Often this is done by assuming some representation of the solution, then tracing the consequences of that assumption (expansions in terms of basis functions lead to finite elements, piecewise constants lead to finite volumes, point values lead to finite differences, ...). It is not necessary to make those assumptions in order to obtain the operators. We didn't assume anything when we wrote down the example operators in this answer. But if you want to, I won't stop you. It may lead to some interesting insights.

Burger's equation

To limit this very long answer, let me just point out how to make the same principles work for Burger's equation, $$ u_t + u u_x = 0. $$ You will have to use the fact that $uu_x = (uu_x + (u^2)_x)/3$ and discretize as $$ M \mathbf{u}_t + \frac{1}{3} \left( \text{diag}(\mathbf{u}) K \mathbf{u} + K \text{diag}(\mathbf{u}) \mathbf{u} \right) = \mathbf{0}. $$ Here, it is necessary to additionally demand that $M$ is diagonal. A suitable numerical flux may be $(u_l^2 + u_l u_r + u_r^2)/6$. In fact, the splitting of $uu_x$ and this flux are intimately connected for reasons that is beyond the scope of this discussion.