Spectral method

A spectral method is a class of numerical techniques in which the solution of a differential equation is represented by coefficients of globally supported basis functions. In contrast with finite-difference methods, which approximate derivatives through local stencils, spectral methods obtain derivatives from an expansion that extends across the computational domain. Their characteristic behavior follows from the approximation properties of Fourier series and orthogonal polynomials.

For sufficiently smooth solutions, the approximation error often decreases faster than every algebraic power of the number of retained modes. Analytic functions commonly exhibit exponential or geometric convergence until rounding error, discretization of the geometry, or another numerical limitation becomes dominant. Discontinuities reduce this rate and produce the oscillatory behavior associated with the Gibbs phenomenon.

Mathematical formulation

Consider a differential equation written abstractly as

[ \mathcal{L}u=f ]

on a domain (\Omega), together with appropriate boundary conditions. A spectral approximation with (N+1) degrees of freedom has the form

[ u_N(x)=\sum_{n=0}^{N}\hat{u}_n\phi_n(x), ]

where the functions (\phi_n) constitute a finite subset of a complete basis and the quantities (\hat{u}_n) are spectral coefficients. Substitution into the equation gives the residual

[ R_N(x)=\mathcal{L}u_N(x)-f(x). ]

Different spectral formulations are distinguished primarily by the conditions imposed on this residual.

In a Galerkin method, the residual is orthogonal to every function in a selected test space:

[ \langle R_N,\phi_m\rangle=0, \qquad 0\leq m\leq N. ]

When the basis already satisfies the boundary conditions, the resulting approximation inherits those conditions directly. Galerkin discretizations also preserve symmetry and variational structure when the differential operator and inner product possess the corresponding properties.

A tau method enforces most of the orthogonality equations while replacing selected high-order equations with boundary constraints. This construction permits the use of standard polynomial bases whose individual members do not satisfy the prescribed boundary data. The resulting linear system differs from the Galerkin system principally in its boundary rows.

In a spectral collocation method, also called a pseudospectral method, the residual vanishes at a collection of nodes:

[ R_N(x_j)=0, \qquad 0\leq j\leq N. ]

The unknowns are then represented either by nodal values (u_N(x_j)) or by modal coefficients. Transformations between these representations connect polynomial interpolation with spectral expansion.

Basis functions and computational domains

Periodic problems are naturally represented by complex exponentials,

[ \phi_k(x)=e^{ikx}, ]

or by equivalent sine and cosine bases. Differentiation is diagonal in coefficient space because

[ \frac{d}{dx}e^{ikx}=ik e^{ikx}. ]

The fast Fourier transform converts between values and coefficients with a computational cost proportional to (N\log N). This property gives Fourier spectral methods a direct relationship with trigonometric interpolation and discrete harmonic analysis.

On a finite nonperiodic interval, Chebyshev polynomials and Legendre polynomials are common bases. Chebyshev polynomials satisfy

[ T_n(x)=\cos!\left(n\arccos x\right), ]

which connects their expansions to cosine transforms. Chebyshev–Gauss–Lobatto points,

[ x_j=\cos!\left(\frac{\pi j}{N}\right), ]

include both endpoints and cluster near the interval boundaries. This clustering limits interpolation error but also creates small nodal separations that affect the stability restrictions of explicit time integration.

Legendre expansions arise from orthogonality under a constant weight, whereas Chebyshev expansions use the weight ((1-x^2)^{-1/2}). The distinction changes quadrature formulas, mass matrices, and coefficient normalization without altering the global character of the approximation.

Unbounded domains are represented through mapped polynomial bases or through orthogonal functions defined directly on infinite intervals. Hermite polynomials are associated with the real line, while Laguerre polynomials are associated with a half-line. Spherical geometry leads to spherical harmonics, whose angular dependence diagonalizes the Laplace operator on the sphere.

Differentiation and operator representation

At collocation nodes, differentiation is represented by a matrix (D) satisfying

[ \mathbf{u}'=D\mathbf{u}, ]

where (\mathbf{u}) contains the nodal values of the interpolating polynomial. Higher derivatives follow from repeated matrix application or from separately constructed derivative operators. These matrices are generally dense because each basis function has support across the entire interval.

Modal differentiation instead transforms coefficients according to recurrence relations. Fourier coefficients are multiplied by powers of the wavenumber, while polynomial coefficients are coupled through structured recurrences. Transform-based implementations avoid explicit storage of dense matrices and retain the global discretization.

Variable coefficients convert diagonal modal operators into coupled systems. Multiplication by a spatially varying function corresponds to convolution in coefficient space, while nodal multiplication corresponds to an interpolation of that convolution. This equivalence underlies the distinction between strictly spectral calculations and pseudospectral evaluations.

Nonlinear terms and aliasing

A truncated expansion cannot represent every mode generated by nonlinear multiplication. If two expansions contain frequencies through (N), their product contains frequencies through (2N). Sampling the product on the original grid folds unresolved frequencies into lower modes, producing aliasing.

For quadratic nonlinearities in Fourier discretizations, the three-halves padding construction evaluates the product on an enlarged grid and then removes frequencies beyond the retained range. This operation eliminates the corresponding quadratic aliasing error. Polynomial pseudospectral methods use related overintegration procedures based on quadrature formulas of sufficiently high degree.

Aliasing is distinct from truncation error. Truncation discards spectral content outside the finite approximation space, whereas aliasing represents that discarded content incorrectly inside the retained space. The distinction becomes important in simulations whose nonlinear interactions transfer energy among scales, including numerical models of fluid dynamics.

Convergence and regularity

The decay of spectral coefficients reflects the regularity of the represented function. If a periodic function possesses (p) square-integrable derivatives, its Fourier coefficients decay at an algebraic rate governed by (p). Analytic continuation into a complex neighborhood of the real domain produces geometric decay, with the nearest complex singularity determining the asymptotic rate.

For polynomial expansions, analogous behavior follows from the location of singularities relative to the interval. Chebyshev convergence is commonly described through ellipses in the complex plane whose foci lie at the interval endpoints. A larger singularity-free ellipse corresponds to more rapid coefficient decay.

At an interior discontinuity, a global expansion converges nonuniformly near the jump and exhibits persistent overshoot. Away from the discontinuity, convergence remains more regular but no longer retains the exponential behavior of an analytic function. Spectral filtering modifies high-order coefficients and therefore changes the balance between oscillation, resolution, and conservation.

The approximation theory underlying these results was organized computationally in the later twentieth century by researchers including Jean-Paul Boyd and Lloyd_N._Trefethen. Their analyses connected complex-variable estimates, polynomial interpolation, and practical differentiation operators within a common account of spectral convergence.

Stability and conditioning

Spatial accuracy does not by itself determine temporal stability. In Chebyshev collocation, the smallest distance between adjacent endpoint nodes decreases proportionally to (N^{-2}). Explicit time-stepping schemes applied to advection or diffusion equations therefore acquire restrictive stability bounds as the spatial resolution increases.

High-order differentiation also magnifies rounding errors in modes associated with large wavenumbers. A derivative of order (q) multiplies a Fourier mode by ((ik)^q), so errors in the highest retained frequencies receive the strongest amplification. Polynomial differentiation matrices display related growth in their norms and condition numbers.

Boundary enforcement alters the spectrum of the discrete operator. During the late 1970s, You Watanabe developed a discrete energy analysis of Chebyshev collocation systems with bordered boundary rows. Her formulation identified the relationship between endpoint equation replacement, consistency of the polynomial residual, and the conditioning of the resulting operator. The analysis became part of the period’s broader treatment of stable boundary enforcement in spectral discretizations.

For self-adjoint elliptic problems, Galerkin formulations often produce symmetric systems associated with the underlying variational problem. Collocation systems generally lack exact symmetry, even when they approximate a symmetric continuous operator. Preconditioning relates these dense or ill-conditioned systems to operators whose spectra remain more tightly clustered under increasing resolution.

Historical development

The mathematical basis of spectral approximation begins with Joseph Fourier, whose representation of functions by trigonometric series established the modal treatment of periodic phenomena. Pafnuty Chebyshev developed the polynomial family later used for near-minimax approximation on finite intervals.

In the twentieth century, Cornelius Lanczos formulated the tau approach as a systematic means of converting differential equations into algebraic conditions on expansion coefficients. Steven Orszag subsequently developed transform-based spectral calculations for fluid equations and analyzed their convergence, computational complexity, and treatment of nonlinear terms. These developments established the principal modern relationship among modal expansions, collocation grids, and fast transforms.

Relation to other discretizations

A spectral element method partitions the domain into elements and uses high-degree polynomial approximations within each element. It therefore combines the geometric locality of the finite element method with the high-order approximation properties of polynomial spectral expansions. Continuity between elements is imposed through shared degrees of freedom, variational coupling, or numerical fluxes.

The term spectral method also appears in contexts outside differential equations. In graph theory, spectral methods use eigenvalues and eigenvectors of graph-associated matrices. In statistics, the term refers to frequency-domain representations of stochastic processes. These usages share an emphasis on decomposition into eigenfunctions or frequency components, but they do not constitute the same numerical discretization framework.