GLL quadrature and the diagonal mass matrix

Spectral Methods

The mass matrix in FEM is tridiagonal — a small price to pay for first-order elements. In SEM with high $N$ it would naturally be dense. SEM avoids that fate by a single structural choice that this chapter unpacks. The result is a mass matrix that's diagonal, no matter how high the polynomial order. That single fact is most of why SEM is fast.

What we're trying to compute

For nodal Lagrange basis functions $ℓ_{i}$ on a single reference element $[- 1, 1]$ :

M_{ij} = \int_{- 1}^{1} ℓ_{i} (x) ℓ_{j} (x) d x .

The integrand is a polynomial of degree $2 N$ (product of two degree- $N$ polynomials). For most basis choices this gives a fully dense $(N + 1) \times (N + 1)$ matrix.

GLL quadrature

Gauss-Lobatto-Legendre quadrature on $[- 1, 1]$ uses the GLL nodes themselves as integration points, with closed-form weights:

\int_{- 1}^{1} f (x) d x \approx k = 0 \sum N w_{k} f (x_{k}), w_{k} = \frac{2}{N ( N + 1 ) [ P _{N} ( x _{k} ) ] ^{2}} .

GLL quadrature is exact for polynomials up to degree $2 N - 1$ . Compare with Gauss-Legendre (without the Lobatto endpoints) which is exact up to $2 N + 1$ — two degrees more. GLL gives up two degrees of accuracy in exchange for landing the integration points on the endpoints, which lets neighboring elements share boundary nodes cleanly. The "Lobatto" in the name means "with the endpoints."

def gll_weights(nodes: np.ndarray) -> np.ndarray:
    """GLL quadrature weights at the GLL nodes.

    The closed form is  w_k = 2 / (N (N+1) [P_N(x_k)]²),
    where N = len(nodes) - 1 is the polynomial order.
    """
    N = len(nodes) - 1
    w = np.empty(N + 1)
    for k in range(N + 1):
        p, _, _ = _legendre(N, nodes[k])
        w[k] = 2 / (N * (N + 1) * p * p)
    return w

The diagonalization

Now the trick. Apply GLL quadrature to the mass-matrix integral, using the same GLL nodes the basis was built on:

M_{ij}^{GLL} = k = 0 \sum N w_{k} ℓ_{i} (x_{k}) ℓ_{j} (x_{k}) .

The Kronecker-delta property strikes: $ℓ_{i} (x_{k}) = δ_{ik}$ . Every term in the sum is zero unless $i = j = k$ . So:

M_{ij}^{GLL} = w_{i} δ_{ij}

Diagonal. The off-diagonal entries — which would have been nonzero in the exact mass matrix — vanish under GLL quadrature.

What's the catch?

The integrand $ℓ_{i} ℓ_{j}$ has degree $2 N$ , but GLL quadrature is only exact through degree $2 N - 1$ . So the diagonal mass matrix is wrong by one quadrature degree — an under-integration error of order $O (h^{N})$ on a smooth problem.

For SEM that's a great trade. Mass matrix entries on the diagonal are correct to within a tiny relative error (the diagonal entries equal $w_{i}$ exactly when the basis is orthogonal to itself in the discrete inner product — which it is, here). The "missing degree" only shows up in off-diagonals that we threw away anyway. And in exchange:

# In any production SEM code, you'd never compute the off-diagonal
# entries of M. The mass matrix simply IS the vector of GLL weights —
# stored as a 1D array, applied as elementwise multiplication.
M = gll_weights(nodes)               # length N+1, not (N+1)^2

def apply_mass(u, M): return M * u            # O(N), not O(N²)
def apply_inverse_mass(u, M): return u / M    # also O(N)

Solving $M c = b$ goes from a linear system ( $O (N^{2})$ per element) to a single elementwise division ( $O (N)$ ). For time-stepping schemes that need $M^{- 1}$ at every step, the saving is the difference between "viable" and "not viable" at high $N$ .

The two mass matrices, side by side

The visualizer below shows both: the exact mass matrix computed via high-order Gauss-Legendre (used as ground truth), and the GLL-quadrature mass matrix from the same basis. Same nodes. Same basis. Different quadrature. The visible difference: a dense matrix vs. a diagonal one.

def mass_matrix_exact(nodes: np.ndarray) -> np.ndarray:
    """Mass matrix computed with high-order Gauss-Legendre quadrature.
    Integrates ℓ_i ℓ_j (degree 2N) to machine precision."""
    N = len(nodes) - 1
    Q = 2 * N + 4                           # plenty for degree-2N polynomials
    g_pts, g_wts = gauss_legendre(Q)
    M = np.zeros((N + 1, N + 1))
    for i in range(N + 1):
        for j in range(N + 1):
            M[i, j] = sum(
                g_wts[q] * lagrange(i, x, nodes) * lagrange(j, x, nodes)
                for q, x in enumerate(g_pts)
            )
    return M


def mass_matrix_gll(nodes: np.ndarray) -> np.ndarray:
    """Mass matrix computed with GLL quadrature *on the basis nodes*.
    The Kronecker-delta property collapses the integral to a diagonal."""
    return np.diag(gll_weights(nodes))

order N

view

off-diagonal max — exact: 0.0353off-diagonal max — GLL: 0.0e+0

GLL weights: [0.048, 0.277, 0.432, 0.488, 0.432, 0.277, 0.048]

Things to notice:

The exact mass matrix has visible structure around the diagonal but with non-zero bands extending several rows out. That structure encodes "neighboring basis functions overlap" — true, just not the whole story we want.
The GLL mass matrix is purely diagonal at every $N$ . The off-diagonal-max readout confirms it: bit-pattern zero, modulo floating-point noise.
The diagonal entries of the GLL matrix are the GLL weights — the same ones you'd use for regular GLL quadrature. The mass matrix, viewed as a vector, is the weight vector.
Bumping $N$ grows the matrix size but doesn't change the diagonal-ness. SEM scales to high orders without the mass-matrix story getting harder.

What just got built

On top of gll_points and lagrange from chapter 1, we now have gll_weights and the diagonal mass matrix. Three lines of Python encode the entire "mass matrix" of an SEM element:

nodes = gll_points(N)
weights = gll_weights(nodes)
M = weights         # done. this is the mass matrix.

Chapter 3 does the same job for the stiffness matrix via the spectral differentiation matrix $D_{ij} = ℓ_{j}^{'} (x_{i})$ . The differentiation matrix won't be diagonal — the gradient of one basis function isn't zero at the other nodes — but it has its own closed-form structure that makes the resulting stiffness matrix dense but cheap to apply.