Spectral Methods

Цели урока

• Explain exponential convergence and why spectral methods beat finite differences for smooth problems
• Implement Chebyshev collocation for boundary value problems
• Apply Galerkin spectral methods and understand spectral element decomposition
• Handle nonlinear terms via pseudospectral approach and 3/2-rule de-aliasing

Предварительные знания

• nm-22
• nm-14

• nm-22
• nm-14

Spectral methods deliver solutions of unprecedented accuracy for smooth problems - weather forecast models, plasma simulations, fluid turbulence, and quantum chemistry all rely on spectral solvers to achieve results that finite difference codes would need orders-of-magnitude more grid points to match.

• ECMWF global weather model uses spherical harmonics (a spectral method) for dynamical core - this design, unchanged for 40 years, remains the gold standard for 10-day forecasts
• Direct numerical simulation (DNS) of turbulence: pseudospectral Fourier codes like HIT3D resolve Kolmogorov scales on petascale machines, informing turbulence models used in aircraft design
• Quantum chemistry codes (DALTON, Gaussian) use Gaussian-type orbital bases that share the exponential-convergence philosophy of spectral methods
• Laser pulse propagation in nonlinear optics: split-step Fourier method (a pseudospectral technique) simulates femtosecond pulses with mm-scale spatial resolution

History of Spectral Methods

1966: Orszag and Patterson introduce the pseudospectral method, showing FFT makes nonlinear PDE solving practical. 1969: Orszag's landmark turbulence DNS paper demonstrates spectral DNS on 32^3 grid - first rigorous turbulence simulation. 1975: ECMWF adopts spectral model for global weather prediction; operational forecast accuracy improves dramatically

Exponential Convergence: Why Spectral Methods Win

Finite difference and finite element methods converge algebraically: halving the mesh spacing h roughly divides the error by 2^p for some fixed order p. Spectral methods break this barrier entirely. For analytic (infinitely smooth) functions the error decreases exponentially with the number of modes N, often written as O(exp(-cN)). Practitioners call this spectral accuracy or infinite-order convergence.

If a function belongs to the Sobolev space H^k for all k (i.e., has derivatives of every order), then the N-th Fourier coefficient decays faster than any polynomial in 1/N. The truncated series therefore converges faster than any polynomial rate.

The trade-off is global support: each basis function spans the entire domain, so a discontinuity or sharp gradient causes the Gibbs phenomenon, polluting all modes. Spectral methods therefore shine on smooth, periodic (or carefully treated non-periodic) problems.

Chebyshev Collocation for Boundary Value Problems

Periodic Fourier methods apply naturally to problems on circles or periodic boxes. Non-periodic problems on an interval [-1, 1] call for Chebyshev polynomials T_n(x) = cos(n arccos(x)). Collocation at Chebyshev-Gauss-Lobatto (CGL) points clusters nodes near the endpoints, counteracting the Runge phenomenon that plagues equispaced interpolation.

The Chebyshev D matrix has condition number O(N^2) and D^2 has O(N^4). For large N, iterative solvers or preconditioning are preferred over direct LU factorization.

Galerkin Method and Spectral Elements

The Galerkin approach projects the PDE onto a finite-dimensional function space: write u_N as a sum of basis functions and require the residual to be orthogonal to the same basis. For Fourier or Chebyshev bases this yields diagonal or near-diagonal mass matrices, making time-stepping efficient.

Spectral element methods (SEM) inherit the best of both worlds. The domain is partitioned into elements (like FEM), and within each element a high-order Chebyshev or Legendre expansion is used. Continuity is enforced at element boundaries. This achieves exponential convergence in the polynomial degree p while remaining geometrically flexible - complex shapes handled by the mesh, accuracy handled by the spectral expansion.

SEM supports both h-refinement (adding more elements) and p-refinement (raising polynomial degree within an element). Exponential convergence reappears when p is increased with fixed smooth elements.

Pseudospectral Methods and De-aliasing

Galerkin methods handle linear terms elegantly, but nonlinear terms such as u*(du/dx) create a complication: multiplication in physical space is convolution in spectral space, which couples high-frequency modes and produces aliasing errors - energy artificially appearing at low wavenumbers from modes above the cutoff.

The pseudospectral approach embraces physical-space multiplication and combats aliasing with the 3/2-rule or 2/3-rule (zeroing the top 1/3 of modes before transform). This reduces the effective resolution slightly but removes the aliasing instability entirely. All major spectral solvers - nek5000, Dedalus, SpectralDNS - use this strategy.

Without dealiasing, aliased energy accumulates at high wavenumbers, leading to a numerical instability called spectral blocking. In turbulence simulations this ruins energy spectra within a few eddy turnover times.

Related Topics

Chebyshev collocation matrices appear as building blocks in spectral element codes for fluid dynamics and structural mechanics

• Pseudospectral Turbulence DNS — Related topic
• Split-Step Fourier Method — Related topic
• Spectral Element Methods — Related topic

Key Takeaways

• Spectral methods achieve exponential (infinite-order) convergence for smooth, analytic functions - error drops geometrically with the number of modes
• Chebyshev collocation solves non-periodic BVPs to near-machine-precision with remarkably few grid points; CGL node clustering prevents Runge oscillations
• Galerkin spectral formulations decouple linear PDEs into independent mode ODEs, enabling exact exponential integrators with no stiffness
• Pseudospectral methods handle nonlinearities in physical space; the 3/2-rule dealiasing eliminates aliasing instability by zero-padding before multiplication

Вопросы для размышления

• A function has a jump discontinuity at x=0 on [-pi, pi]. How would the L2 error of a truncated Fourier series behave as N grows, and what numerical technique could recover exponential convergence?
• The ECMWF spectral model uses triangular truncation at wavenumber T_co = 1279, corresponding to roughly 16 km horizontal resolution. Estimate how many grid points a 4th-order finite difference scheme would need to match the same accuracy for a typical synoptic-scale wave.
• In a pseudospectral Navier-Stokes simulation with N=512 Fourier modes, what is the minimum padded size for aliasing-free quadratic nonlinearities, and what fraction of compute time would padding add?

Связанные уроки

• nm-22
• nm-14
• nm-24-parallel-linalg
• nm-21