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.


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.

  • $\begingroup$ Regarding your two questions: which integrals do you mean are inexact? The LGL quadrature integrates polynomials up to the given degree, and provides an approximation to all other integrands. Next, there actually is a spectral element version of this idea, which is heavily used at least in physics under the name FE-DVR. $\endgroup$ – davidhigh Dec 23 '20 at 16:53
  • $\begingroup$ Thanks for pointing out FE-DVR. I found some papers which lead me to this documentation of a python-package and there I found the following: The DVR representation not variational because both the potential and overlap matrix elements are approximated by the quadrature and are not exact, although with moderately dense grids that property barely noticeable. $\endgroup$ – Dagon Dec 23 '20 at 17:54
  • $\begingroup$ The inexactness is also called mass lumping and is quite standard in the DG community. $\endgroup$ – ConvexHull Dec 23 '20 at 21:58
  • 1
    $\begingroup$ The reason why LGL falls short is because essentially you are adding a $(1-x^2)$ as weight function to the integral. By this you effectively lose two orders, and therefore your example integrates only polynomials up to three (and not five). $\endgroup$ – davidhigh Dec 23 '20 at 22:06
  • $\begingroup$ One thing that you should keep in mind is that, in general, you never integrate "stiffness" matrices exactly since those are not polynomial functions but rational functions. If your elements aren't distorted then you indeed have polynomial functions. $\endgroup$ – nicoguaro Dec 24 '20 at 16:32

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.

  • $\begingroup$ A small suggestion: I would appreciate it if you describe the connection/differences between the quadrature rules in more detail. For example $M_{ij} =\int_{-1}^1 L_i(x) L_j(x)\,dx\approx \sum_{k=0}^N w_k L_i(x_k)L_j(x_k) = w_{i}\delta_{ij}$ on Gauss-Legendre nodes is actually the well known Gauss quadrature rule with weighting function $1$ and therefore exact for polynomials of degree $2N+1$. $\endgroup$ – ConvexHull Dec 29 '20 at 17:10
  • $\begingroup$ @ConvexHull: thanks for the suggestion, but could you be a bit more specific? Basically, your formula says it already. You can do the quadrature on arbitrary grids, or you can use Gaussian grids, in which case everytjing gurns out quite nicely and you get "spectral accuracy". By the way, I'm missing a weight function in the mass integral (because one is not tied to Gauss-Legendre)... I'll edit that in. But again, please specify your question. $\endgroup$ – davidhigh Dec 29 '20 at 17:48
  • $\begingroup$ I am only a little bit confused. The way you describe it, the reader assumes that the implied quadrature method on general nodes to evaluate $M_{ij} =\int_{-1}^1 w(x) L_i(x) L_j(x)\,dx$ is only exact for order N. To derive the Legendre-Gauss or Lobatto quadrature weights you do not have to do anything special. Choosing the right nodes gives you directly the implied weights. $\endgroup$ – ConvexHull Dec 29 '20 at 18:19
  • $\begingroup$ @ConvexHull: ok, maybe that was is a bit confusing. The line of explanation goes from a general grid to Gaussian grids to Lobatto grids. For a general grid, and also for equally spaced grids, the order of integration is N (resp. N+1 for an equally spaced grid with odd numer of nodes). A general grid is suited to explain the discrepancy between accuracy and convenient formulation (diagonality) that is shared also by Lobatto grids and cured only for Gaussian (and Radau) grids. $\endgroup$ – davidhigh Dec 29 '20 at 20:35

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.