Spectral Methods

Pdes

Spectral methods represent the unknown function in a GLOBAL basis of smooth functions — typically Fourier modes on periodic domains, Chebyshev polynomials on intervals, spherical harmonics on the sphere — rather than as values on a mesh. Derivatives are computed by differentiating the basis functions exactly; the discrete problem is a linear system in the coefficient space. The defining property of spectral methods: for SMOOTH solutions, the error converges FASTER than any polynomial rate in the number of basis functions — SPECTRAL or EXPONENTIAL convergence. A 32-point Chebyshev spectral solve can match the accuracy of a several-thousand-point finite-difference solve when the underlying solution is smooth.

The canonical reference is J. P. Boyd's Chebyshev and Fourier Spectral Methods (Dover, 2001) — by the same author whose "Devil's Invention" review anchors the hyperasymptotics page. Boyd's book is the standard practical guide and is freely available online from his website.

Why spectral convergence

For a function $u \in C^{k}$ on $[- 1, 1]$ , the Chebyshev-coefficient expansion $u (x) = \sum_{n} a_{n} T_{n} (x)$ has coefficients decaying like $a_{n} \sim 1/ n^{k}$ — fast for high-regularity $u$ . For a $C^{\infty}$ function the coefficients decay FASTER than any polynomial, and for an ANALYTIC function (extension to a strip around $[- 1, 1]$ in the complex plane) they decay EXPONENTIALLY: $∣ a_{n} ∣ \leq C ρ^{- n}$ for some $ρ > 1$ .

Truncating the series after $N$ terms gives error $\sim ρ^{- N}$ for analytic functions — exponential decay in $N$ . Finite differences give algebraic decay $\sim h^{p} = N^{- p}$ for some fixed integer $p$ . Algebraic vs exponential is the same gap as polynomial vs exponential, and it is the entire reason spectral methods are unbeatable on smooth problems.

The catch: for non-smooth solutions (shocks, kinks, jumps), spectral convergence DEGRADES. Discontinuities produce $\sim 1/ n$ coefficient decay (Gibbs phenomenon), giving the same $O (N^{- 1})$ convergence as a first-order finite-difference scheme — or worse, oscillations near the discontinuity. Spectral methods shine for smooth problems and require careful handling for problems with shocks.

Chebyshev collocation

Among the spectral families, Chebyshev polynomials $T_{n} (x) = cos (n arccos x)$ are the workhorse for non-periodic problems on a bounded interval. Two reasons. (1) The Chebyshev grid $x_{k} = cos (kπ / N)$ for $k = 0, 1, \dots, N$ clusters near the endpoints ( $\sim 1/ N^{2}$ spacing) — exactly what's needed to avoid the Runge phenomenon (polynomial interpolation at equispaced nodes gives wild oscillations near the endpoints). (2) Differentiation in this basis can be computed by a DIFFERENTIATION MATRIX: a dense $(N + 1) \times (N + 1)$ matrix $D$ such that $(D u)_{j}$ equals $u^{'} (x_{j})$ for any polynomial of degree $\leq N$ sampled at the Chebyshev nodes.

The differentiation matrix has explicit entries (Canuto-Hussaini-Quarteroni-Zang; Trefethen's textbook Spectral Methods in MATLAB):

D_{ij} = ⎩ ⎨ ⎧ \frac{c _{i} ( - 1 ) ^{i + j}}{c _{j} ( x _{i} - x _{j} )}, - \sum_{k \neq = i} D_{ik}, i \neq = j i = j

with $c_{0} = c_{N} = 2$ and $c_{k} = 1$ for $0 < k < N$ . The diagonal entries are determined by the requirement that $D$ annihilate constants ( $D 1 = 0$ ). This is the negative-sum trick — the cleanest way to compute diagonals stably.

Higher derivatives are obtained by repeated application: $(D^{2} u)_{j} = u^{''} (x_{j})$ for polynomials of degree $\leq N$ . Dirichlet boundary conditions $u (\pm 1) = 0$ are imposed by STRIPPING the first and last rows and columns of $D^{2}$ (which correspond to the boundary points) — the resulting $(N - 1) \times (N - 1)$ linear system is solved for the interior unknowns. Other boundary conditions can be enforced by modifying the boundary rows directly.

Code: 1D Poisson with spectral convergence

A concrete benchmark. Solve $- u^{''} (x) = f (x)$ on $[- 1, 1]$ with $u (\pm 1) = 0$ , choosing $u (x) = (1 - x^{2}) cos (5 x)$ as the exact solution (a smooth, oscillatory profile that satisfies the boundary conditions automatically), and check the error at several discretization sizes.

# Chebyshev spectral collocation for 1D Poisson on [-1, 1].
#
#   -u''(x) = f(x),   u(-1) = u(1) = 0.
#
# Discretize on the Chebyshev nodes x_k = cos(k pi / N) and approximate
# differentiation by the (N+1)x(N+1) Chebyshev differentiation matrix D.
# Then D² approximates d²/dx². Apply BCs by stripping the boundary rows
# and columns. Solve the resulting (N-1)x(N-1) linear system.
#
# Choose a smooth test problem:  exact u(x) = (1 - x²) cos(5 x).
# Compute f(x) = -u''(x) analytically, run the solver at several N,
# observe the EXPONENTIAL convergence in N (until round-off floors).

import numpy as np

def cheb_diff_matrix(N):
    """Trefethen's cheb.m. Returns (N+1)x(N+1) differentiation matrix
    and the Chebyshev nodes x_k = cos(k pi / N), k = 0, ..., N."""
    if N == 0:
        return np.zeros((1, 1)), np.array([1.0])
    x = np.cos(np.pi * np.arange(N + 1) / N)
    c = np.array([2.0] + [1.0]*(N - 1) + [2.0]) * (-1)**np.arange(N + 1)
    X = np.tile(x, (N + 1, 1)).T
    dX = X - X.T
    D = np.outer(c, 1.0/c) / (dX + np.eye(N + 1))
    D -= np.diag(D.sum(axis=1))
    return D, x

def u_exact(x):
    return (1 - x**2) * np.cos(5 * x)

def f_rhs(x):
    """f = -u''  for  u = (1 - x²) cos(5 x)."""
    # u'  = -2 x cos(5x) - 5 (1 - x²) sin(5x)
    # u'' = -2 cos(5x) + 20 x sin(5x) - 25 (1 - x²) cos(5x)
    return 2*np.cos(5*x) - 20*x*np.sin(5*x) + 25*(1 - x**2)*np.cos(5*x)

print("Chebyshev solve: -u'' = f, exact u(x) = (1-x²) cos(5x), u(±1) = 0")
print(f"  {'N':>5s}  {'||err||_inf':>14s}")
for N in [8, 12, 16, 20, 24, 32, 48, 64]:
    D, x = cheb_diff_matrix(N)
    D2 = D @ D
    # Strip the boundary rows and columns (Dirichlet u(±1) = 0)
    D2_int = D2[1:-1, 1:-1]
    u_int  = np.linalg.solve(-D2_int, f_rhs(x[1:-1]))
    u = np.zeros(N + 1)
    u[1:-1] = u_int
    err = np.max(np.abs(u - u_exact(x)))
    print(f"  {N:>5d}  {err:>14.4e}")
print("(Error drops by orders of magnitude per refinement — spectral convergence.)")

Output:

Chebyshev solve: -u'' = f, exact u(x) = (1-x²) cos(5x), u(±1) = 0
      N     ||err||_inf
      8      2.9938e-02
     12      8.5614e-05
     16      3.7458e-08
     20      8.1330e-12
     24      6.1062e-14
     32      2.3315e-15
     48      2.2538e-14
     64      6.3949e-14
(Error drops by orders of magnitude per refinement — spectral convergence.)

Three orders of magnitude per $N$ increment of 4. Reaching machine precision (~ $1 0^{- 14}$ ) at $N = 24$ — twenty-five degrees of freedom for fourteen digits of accuracy. Compare to finite differences on the same problem: hitting $1 0^{- 14}$ with central second-order FD would need $N \sim 1 0^{7}$ grid points (since FD error $\sim h^{2} \sim 1/ N^{2}$ , and $(1 0^{7})^{- 2} = 1 0^{- 14}$ ). Six orders of magnitude fewer degrees of freedom for the same accuracy. The error floors at $N \geq 32$ at the round-off limit — the conditioning of $D^{2}$ scales like $N^{4}$ , so accuracy beyond $ϵ_{mach} N^{4} \sim 1 0^{- 16} \cdot 1 0^{6} = 1 0^{- 10}$ is gradually lost; with the negative-sum trick we still squeeze out $1 0^{- 14}$ .

Fourier spectral methods on periodic domains

For periodic boundary conditions, the natural basis is plain Fourier modes $e^{ik x}$ . Differentiation in this basis is MULTIPLICATION BY $ik$ in coefficient space:

u (x) = k \sum \overset{u}{^}_{k} e^{ik x} ⟹ u^{(m)} (x) = k \sum (ik)^{m} \overset{u}{^}_{k} e^{ik x} .

With the FFT, you can transform between physical and Fourier representations in $O (N lo g N)$ — so a typical spectral time-step on a periodic domain costs $O (N lo g N)$ , beating any local-stencil method's $O (N)$ cost by only a log factor while delivering exponential accuracy in $N$ . The operator splitting page uses this for the kinetic part of the time-dependent Schrödinger equation.

What spectral methods are good for

Smooth solutions on simple geometries: rectangle, annulus, sphere, periodic box, etc. Climate codes, turbulence simulation in periodic boxes, plasma physics, geophysical fluid dynamics — all use spectral methods at the highest-accuracy end.
Eigenvalue problems: spectral discretization of $\hat{L} u = λ u$ gives a dense small eigenvalue problem with exponentially-accurate eigenpairs. Routine in hydrodynamic stability analysis and quantum mechanics.
Periodic Schrödinger / quantum dynamics: split-step Fourier methods (see TDSE) treat the kinetic operator $- (1/2) \nabla^{2}$ exactly in Fourier space, achieving spectral accuracy for the spatial discretization.
Pseudo-spectral methods for nonlinear PDEs: evaluate nonlinearities (e.g. $u^{2}$ ) in physical space, evaluate derivatives in Fourier space; the workhorse for incompressible Navier-Stokes simulations in periodic geometries.

What they're bad for

Complex geometry: spectral methods need the geometry to match the basis. Use finite elements, spectral elements (high-order FE), or DG methods on unstructured meshes.
Shocks and discontinuities: Gibbs oscillations destroy the spectral convergence and pollute the solution far from the discontinuity. Limiter techniques and filtered spectral methods help but rarely match the elegance of WENO finite-volume methods for shock-dominated problems.
Very stiff problems: the differentiation matrix has condition number $\sim N^{4}$ for the second derivative; ill-conditioning can be a problem at very high $N$ .
Adaptive mesh refinement: spectral methods are global; refining "locally" is awkward. Spectral elements are the standard fix.

Separation of variables — spectral methods are the numerical embodiment of the same eigenfunction-expansion philosophy.
Finite difference methods — the algebraic-convergence alternative.
Finite elements — global vs local basis functions; spectral elements bridge the two.
Time-dependent Schrödinger — split-step Fourier method is a spectral time integrator.
Hyperasymptotics — Boyd's other major contribution; both topics share his pedagogical style and connect through the analytical-vs-numerical theme.