In my earlier post on the finite element method, I described approximation by piecewise polynomials. For conforming degree-$p$ elements applied to a sufficiently regular elliptic solution on a shape-regular mesh, the energy-norm error is typically $O(h^p)$. Other norms and reduced regularity give different rates.
This is algebraic convergence under fixed polynomial degree. On a quasi-uniform mesh in $d$ dimensions, replacing $h$ by $h/2^8$ increases the number of degrees of freedom by approximately $2^{8d}$ and reduces an $O(h^p)$ error by $2^{8p}$. Adaptive meshes and higher-order finite elements change this comparison.
Spectral methods offer a different convergence regime. Analytic solutions can yield exponential coefficient decay and errors of the form $O(e^{-cN})$ in one-dimensional truncation order $N$. Finitely differentiable solutions give algebraic rates, while infinitely differentiable nonanalytic solutions generally give superalgebraic rather than necessarily exponential convergence. The comparison with finite elements depends on dimension, regularity, geometry, and whether high-order or spectral elements are used.
The Core Idea
A spectral method approximates the solution $u$ of a PDE by a linear combination of globally supported basis functions:
\[u_N(x) = \sum_{k=0}^{N} \hat{u}_k \phi_k(x),\]where ${\phi_k}$ is a sequence of smooth, globally defined functions (Fourier modes, Chebyshev polynomials, Legendre polynomials), and $\hat{u}_k$ are the expansion coefficients to be determined.
The contrast with finite elements is immediate: FEM uses piecewise polynomial basis functions that are zero outside a small patch of elements. Spectral basis functions are nonzero everywhere in the domain. Each spectral basis function “sees” the entire solution, which is why spectral methods can extract more information from a given number of degrees of freedom.
Fourier Methods on Periodic Domains
On the interval $[0, 2\pi]$ with periodic boundary conditions, the natural spectral basis is the Fourier series:
\[u_N(x) = \sum_{k=-N/2}^{N/2} \hat{u}_k e^{ikx}.\]The coefficients $\hat{u}_k$ are the Fourier coefficients of $u$. The quality of the approximation is controlled by how fast these coefficients decay. For a function with $p$ continuous derivatives, $\lvert\hat{u}_k\rvert = O(\lvert k\rvert^{-p})$ as $\lvert k\rvert \to \infty$, giving algebraic decay. For a real analytic function (one that extends to a holomorphic function in a strip around the real axis), the coefficients decay exponentially: $\lvert\hat{u}_k\rvert \leq C e^{-\sigma \lvert k\rvert}$ for some $\sigma > 0$. This exponential decay of coefficients is the origin of spectral convergence.
Differentiation. For a periodic Fourier series, the $k$-th coefficient of $u’$ is $ik\hat{u}_k$. This identity is exact for the truncated representation. Truncation, aliasing, roundoff, and the conditioning of discrete differentiation still affect a numerical method.
The FFT. Evaluating a Fourier series at $N$ equally spaced points, and computing the Fourier coefficients from $N$ function values, costs $O(N \log N)$ via the Fast Fourier Transform. This makes Fourier spectral methods exceptionally efficient: the cost per degree of freedom is nearly optimal.
Gibbs phenomenon. Near an isolated jump, Fourier partial sums exhibit a limiting overshoot of about 9% of the jump height. For a piecewise smooth function with a nonzero jump, the coefficients generally decay as $O(\lvert k\rvert^{-1})$. The convergence rate depends on the norm: the global $L^2$ truncation error is typically of order $N^{-1/2}$, while pointwise convergence fails at the jump and is nonuniform nearby.
Chebyshev and Legendre Methods on Bounded Intervals
For non-periodic problems on $[-1, 1]$, the natural spectral bases are orthogonal polynomials. The two most widely used are:
Chebyshev polynomials $T_k(x) = \cos(k \arccos x)$. They are defined by the recurrence $T_0 = 1$, $T_1 = x$, $T_{k+1}(x) = 2xT_k(x) - T_{k-1}(x)$. The Chebyshev nodes $x_j = \cos(\pi j / N)$ cluster near the endpoints $\pm 1$, which is essential for controlling the Runge phenomenon.
Legendre polynomials $P_k(x)$, defined by the recurrence $P_0 = 1$, $P_1 = x$, $(k+1)P_{k+1} = (2k+1)xP_k - kP_{k-1}$. They are $L^2$-orthogonal with respect to the uniform weight on $[-1, 1]$.
Both families achieve spectral convergence for smooth functions: if $u$ is analytic on $[-1, 1]$, the expansion coefficients decay exponentially and the $N$-term approximation satisfies
\[\lVert u - u_N \rVert_{L^\infty} \leq C e^{-\sigma N}\]for constants $C, \sigma > 0$ depending on the analytic continuation of $u$.
The Chebyshev expansion has the practical advantage that it reduces to a cosine transform: evaluating the expansion at Chebyshev nodes costs $O(N \log N)$ via the DCT (discrete cosine transform). Chebyshev methods are also closely related to Fourier methods through the substitution $x = \cos\theta$.
Solving PDEs: The Galerkin and Collocation Approaches
Given a basis ${\phi_k}$, there are two standard ways to enforce the PDE.
Galerkin spectral method. Require the residual to be orthogonal to the span of the first $N$ basis functions:
\[\langle F(u_N), \phi_j \rangle = 0 \quad \text{for } j = 0, \ldots, N.\]For a linear PDE $Lu = f$, this gives the system $\langle L u_N, \phi_j \rangle = \langle f, \phi_j \rangle$, which in matrix form is a dense linear system for the coefficients $\hat{u}_k$. The resulting matrix is often well-conditioned for elliptic problems when the right basis is used.
Spectral collocation (pseudospectral) method. Require the PDE to be satisfied exactly at a set of collocation points ${x_j}_{j=0}^N$:
\[(Lu_N)(x_j) = f(x_j) \quad \text{for all } j.\]The collocation points are typically the Gauss quadrature nodes (roots of $\phi_{N+1}$) or the Gauss-Lobatto nodes (roots of $(1-x^2)\phi_N’$). Collocation is often simpler to implement than Galerkin for variable-coefficient or nonlinear problems because it avoids computing inner products analytically.
The two approaches are closely related and have similar convergence properties. In practice, spectral collocation is more widely used in fluid dynamics and wave equations; Galerkin is more common in problems where variational structure is important.
Where Spectral Methods Excel
Regular solutions on simple geometries. Fourier methods fit periodic rectangular domains, while polynomial spectral methods fit intervals and tensor-product domains. Analytic data and solutions can produce exponential convergence, but the number of modes needed for a given tolerance depends on singularities in the complex plane, boundary conditions, conditioning, and dimension.
Fluid dynamics. Spectral methods have dominated direct numerical simulation of turbulence since the 1970s. For the Navier-Stokes equations in a box with periodic boundary conditions, Fourier methods are the tool of choice. Many of the benchmark results in turbulence research were computed spectrally.
Wave propagation. High-order accuracy matters for wave equations because low-order methods introduce numerical dispersion: the computed wave speed differs from the true wave speed, and the error accumulates over long propagation distances. Spectral methods, with their exponential accuracy for smooth solutions, suffer far less numerical dispersion.
Where Spectral Methods Struggle
Complex geometry. A Fourier method lives on a rectangle, a Chebyshev method on an interval. Handling L-shaped domains, holes, or patient-specific geometry requires either domain decomposition (spectral element methods) or coordinate transformations that can spoil the exponential convergence. Finite elements handle complex geometry naturally; spectral methods do not.
Discontinuous solutions. As noted, Gibbs oscillations contaminate the solution near any discontinuity. Remedies exist (filtering, WENO reconstruction, Gegenbauer post-processing) but they add complexity and do not fully restore spectral convergence.
Preconditioning. The dense matrices arising in spectral methods can be ill-conditioned, especially for high-order problems. Designing efficient preconditioners is more difficult than for FEM, where the sparse, banded structure facilitates many standard preconditioner constructions.
Why I Know This
Most of my direct numerical work uses finite element methods, because EIT involves complex boundary shapes and potentially discontinuous conductivities, exactly the settings where FEM excels over spectral methods. But spectral methods appear indirectly in my work in two ways.
First, Fourier approximation and spectral regularization both analyze expansions in basis or singular vectors, but their convergence theories are not identical. Spectral cutoff for an inverse problem is governed by singular-value decay, noise, and source conditions, whereas Fourier approximation is governed primarily by function regularity.
Second, the spectral element method (SEM), which combines spectral accuracy within each element with FEM flexibility in geometry, is increasingly used for forward solvers in wave-based inverse problems. Understanding spectral methods makes SEM more transparent.
Spectral methods make the role of regularity in approximation explicit. Finite smoothness gives algebraic decay, additional smoothness can give superalgebraic decay, and analyticity can give exponential decay. This relationship appears throughout approximation theory and operator equations.