GLL points and the nodal Lagrange basis

Spectral Methods

Every SEM thing we'll do downstream — mass matrices, stiffness matrices, derivatives, quadrature, time stepping — sits on top of one decision: where to put the nodes. This chapter justifies the answer and builds the cardinal basis.

The basis on a single element

Inside one element (parameterized to $[- 1, 1]$ for convenience), pick $N + 1$ nodes and use the nodal Lagrange basis:

ℓ_{j} (x) = k \neq = j \prod \frac{x - x _{k}}{x _{j} - x _{k}}, j = 0, 1, \dots, N .

Each $ℓ_{j}$ is a degree- $N$ polynomial that equals 1 at node $x_{j}$ and 0 at every other node. Just like FEM hat functions, but smooth and high-order. The Kronecker-delta property carries through: a function expanded in this basis, $u (x) = \sum_{j} c_{j} ℓ_{j} (x)$ , has coefficients $c_{j} = u (x_{j})$ — values at nodes, no transformation.

def lagrange(j: int, x: float, nodes: np.ndarray) -> float:
    """The j-th cardinal Lagrange polynomial at point x."""
    xj = nodes[j]
    prod = 1.0
    for k, xk in enumerate(nodes):
        if k != j:
            prod *= (x - xk) / (xj - xk)
    return prod


# Property check: ℓ_j is 1 at its own node and 0 at every other node.
nodes = gll_points(8)
for j in range(len(nodes)):
    print(f"ℓ_{j} at nodes: {[round(lagrange(j, x, nodes), 4) for x in nodes]}")

ℓ_0 at nodes: [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
ℓ_1 at nodes: [0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
ℓ_2 at nodes: [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
...
ℓ_8 at nodes: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]

Why GLL specifically

Three node distributions to consider on $[- 1, 1]$ : equispaced (the obvious choice — and the one that will fail), Chebyshev ( $x_{k} = - cos (πk / N)$ , used in global spectral methods), and Gauss-Lobatto-Legendre.

Equispaced nodes look natural and behave terribly. Their Lebesgue constant — a quantitative measure of how badly interpolation can amplify input error — grows like $2^{N} / (N lo g N)$ . Exponentially bad. The visible consequence is Runge's phenomenon: interpolating a smooth function at equispaced nodes produces enormous oscillations near the endpoints, growing with $N$ . The visualizer below shows it directly.

Chebyshev and GLL nodes both cluster near the endpoints, which suppresses the Runge oscillations. Their Lebesgue constants grow only like $lo g N$ — exponentially better. Both give spectral interpolation accuracy on smooth functions.

What sets GLL apart from Chebyshev is the matched quadrature. Chebyshev nodes are the natural points for Chebyshev expansions; GLL nodes are the quadrature points of GLL quadrature. Using GLL nodes for both the basis and the quadrature gives a rare alignment that makes the mass matrix diagonal — chapter 2's punchline. Chebyshev doesn't have that property; that's why SEM picks GLL even though Chebyshev would interpolate equally well.

Computing the nodes

The N+1 GLL nodes are: $\pm 1$ plus the $N - 1$ interior roots of the derivative of the $N$ -th Legendre polynomial, $P_{N}^{'} (x)$ .

Newton iteration on $P_{N}^{'}$ , starting from the Chebyshev nodes $- cos (πk / N)$ as initial guesses, converges to machine precision in a few iterations. Compute $P_{N}$ via the three-term recurrence $(n + 1) P_{n + 1} = (2 n + 1) x P_{n} - n P_{n - 1}$ ; derive $P_{N}^{'}$ and $P_{N}^{''}$ from the standard Legendre ODE.

import numpy as np

def gll_points(N: int) -> np.ndarray:
    """Return the N+1 Gauss-Lobatto-Legendre nodes on [-1, 1].

    Endpoints are ±1; the N-1 interior nodes are roots of P'_N(x).
    Newton iteration on the Legendre polynomial derivative, starting
    from Chebyshev points (always close enough to converge fast).
    """
    if N == 1:
        return np.array([-1.0, 1.0])
    nodes = np.empty(N + 1)
    nodes[0], nodes[-1] = -1.0, 1.0

    for k in range(1, N):
        x = -np.cos(np.pi * k / N)             # initial guess
        for _ in range(60):
            p, dp, ddp = _legendre(N, x)
            dx = -dp / ddp
            x += dx
            if abs(dx) < 1e-14:
                break
        nodes[k] = x
    return nodes


def _legendre(N: int, x: float):
    """P_N(x), P'_N(x), P''_N(x) via three-term recurrence + ODE."""
    pNm2, pNm1 = 1.0, x
    for n in range(1, N):
        pN = ((2 * n + 1) * x * pNm1 - n * pNm2) / (n + 1)
        pNm2, pNm1 = pNm1, pN

    p, p_prev = (1.0, 0.0) if N == 0 else ((x, 1.0) if N == 1 else (pNm1, pNm2))
    dp = N * (p_prev - x * p) / (1 - x * x)
    ddp = (2 * x * dp - N * (N + 1) * p) / (1 - x * x)
    return p, dp, ddp

For $N \leq 30$ or so this is essentially instant — a few microseconds per call. Production codes precompute and cache the nodes; we recompute on the fly because the cost is negligible at our resolution.

See it

Top plot: the $N + 1$ nodes for the chosen distribution (GLL or equispaced) and the cardinal Lagrange basis functions living on them. Click any node, or use the prev/next buttons, to spotlight one basis function — notice it's exactly 1 at its own node and 0 at the others.

Bottom plot: Runge's function $1/ (1 + 25 x^{2})$ interpolated at the same nodes. The truth is the muted gray curve; the interpolant is the accent polyline through the orange nodes. Switch between GLL and equispaced and ramp $N$ up — at $N = 16$ or higher the equispaced interpolant has visible spikes near the endpoints; the GLL one tracks the truth even at the highest order shown.

nodes

order N

spotlightℓ₂ at x = -0.6772max |interp| = 1.00

Things to read off:

With equispaced nodes at $N = 16$ or $N = 20$ , the basis functions develop large excursions near $\pm 1$ . Those excursions are the same mechanism that makes equispaced interpolation fail: the basis isn't well-behaved, and Lagrange interpolation amplifies that ill-conditioning into output error.
With GLL nodes, every basis function stays bounded near 1 in magnitude — even at $N = 20$ . The clustering at the endpoints is the visible mechanism for the logarithmic Lebesgue constant.
The max |interp| readout is the maximum magnitude of the interpolated curve. At GLL it stays at $\approx 1$ (the truth peaks at 1). At equispaced it grows past 2, then 5, then more, as $N$ rises — Runge's phenomenon as a number.

What just got built

Two ingredients ready for chapter 2:

gll_points(N) — returns $N + 1$ nodes on $[- 1, 1]$ .
lagrange(j, x, nodes) — evaluates the $j$ -th cardinal basis function at $x$ .

Chapter 2 builds GLL quadrature on top of these and shows how integrating the basis against itself collapses the mass matrix to a diagonal. That's the moment SEM gets cheap.