Gauss-Lobatto quadrature and nodal points for FEM

Question

By using the Legendre-Gauss-Lobatto (LGL) quadrature formula (QF) and LGL nodal points one achives a diagonal mass-matrix for finite element problems. (More specifically, the spectral element method.)

Setting the $n$ quadrature points as nodal points, the polynomials of the ansatz and test space are of degree $n-1$.

The LGL QF is exact up to degree $2n-3$, but the mass-matrix should contain integrals of polynomials with degree $2n-2$.

Is it normal for this kind of method, that each integral is inexact?

Is it possible at all to combine a DG-method with this ansatz of the spectral element method (SEM)?

Any tips, pointers and literature will be appreciated.

Example:

Setting $n = 3$. Then the LGL QF is exactly Simpson's rule, which is exact up to degree $2n-3 = 3$. By using Langrange shape-functions on the (nodal =) quadrature points $\{x_1=-1,x_2=0,x_3=1\}$, we obtain a basis $\{l_1,l_2,l_3\}$ of our polynomial space, each basis-element is of degree $2$. The basis satisfies $$l_i(x_k) = \delta_{ik}.$$ The calculation of an $L^2$ scalar-product then needs $$\int_{-1}^1l_i(x)l_j(x) \mathrm{d}\,x \approx \sum_{k=1}^3w_k l_i(x_k)l_j(x_k) = w_i\delta_{ij},$$ which is an integral of a polynomial of degree 4, hence the LGL QF (or Simpon's rule for $n=3$) is only an approximation to the exact value.

Answer 1

Ok, here comes the answer promised in the comment section. Let's start the other way round, going from a general grid to Gaussian grids and further constructions such as spectral elements.

In grid methods, one basically selects a number of $N+1$ gridpoints $\{x_k\}_{k=0}^{N}$. As basis functions, one can use Lagrange polynomials constructed over these nodes, each one having polynomial order $N$. For this setup, using Newton-Cotes-like quadrature weights $w_k$ (see e.g. here), one obtains a quadrature method that is able to integrate any polynomial of order up to $N$. For a general gridpoint distribution, this is the best you can get.

Now, if you want to evaluate integrals of the mass matrix, $M_{ij} =\int_{-1}^1 w(x) L_i(x) L_j(x)\,dx$, one basically has two options:

1. Use the implied quadrature method to evaluate $M_{ij} = \sum_{k=0}^N w_k L_i(x_k)L_j(x_k) = w_{i}\delta_{ij}$ and thus obtain a diagonal mass matrix. As can be easily seen by the involved polynomial orders, the quadrature is only an approximation (the integrand has order $2N$, but the quadrature rule is exact only for order $N$).
2. Use any other quadrature method, usually one that is accurate for the involved polynomial orders. This, however, in general won't lead to a diagonal mass matrix.

You see, it's not optimal in either way. The general construction using cardinal basis functions defined over a grid together with the corresponding quadrature rule is either orthogonal or accurate (in polynomial sense).

Enter Gaussian integration: The idea is, that by using special grids and weights, one obtains an integration accurate to order $2N+1$. This is great because these special grids basically solve the conflict shown before for a general grid. By this, one obtains the benefits of both points 1. and 2. above at the same time, i.e. an orthogonal and accurate basis at the same time. (One even gets one integration order more as required for the mass matrix, so one can also integrate the spatial operator "$x$" exactly. This fact stems from the three-term recurrence formula).

Now, Gaussian grids are not always appropriate for modelling real world problems. For boundary problems, for example, it is convenient to prescribe function values at the boundary of the grid. This leads to the Gauss-Lobatto construction, where one fixes the points $x_0=-1$ and $x_N=1$. By this, of course, one looses integration accuracy, and the result of Lobatto is a grid distribution where the loss is as minimal as possible, namely two orders, so that one obtains an integration method accurate up to order $2N-1$. Now you run into the same problems as shown above, as already the mass matrix is not evaluated exactly. As already mentioned in the comments, a common approach is to pick alternative 1. and willingly accept the approximation (That is only one polynomial order short from "exact").

The Gauss-Lobatto grid now also offers a convenient way to construct spectral elements, which is basically to stick together several GL grids and also connect the basis functions correponding to the boundary nodes in a continouus way. The benefits here are, that one gets banded and better conditioned derivative matrices, and that the integration accuracy at least inside elements is quite high (but not accross elements). As already mentioned in the comments, the FE-DVR is an instance of this idea used heavily in quantum mechanics, but the idea itself is older, and dates back at least to a paper of the scicomp-user @L.Young from 1977.