Spectral element method

The spectral element method (SEM) is a high-order technique for the numerical solution of partial differential equations. It combines the geometric decomposition used by the finite element method with polynomial approximation spaces associated with spectral methods. A computational domain is partitioned into elements, within which the unknown field is represented by a polynomial of comparatively high degree. The resulting local approximations are coupled through continuity constraints or numerical fluxes, depending on the formulation.

Spectral element discretizations are used for problems in which spatial accuracy, geometric representation, and computational locality must be considered simultaneously. Their mathematical structure is closely related to high-order Galerkin methods, while their implementation commonly relies on nodal interpolation and specialized numerical quadrature.

Mathematical formulation

Consider a boundary-value problem on a domain (\Omega),

[ \mathcal{L}u=f, ]

where (\mathcal{L}) is a differential operator, (u) is the unknown field, and (f) is prescribed data. For a representative elliptic problem, the operator may be written as

[ -\nabla\cdot\left(\kappa\nabla u\right)=f, ]

with appropriate boundary conditions. Multiplication by a test function (v), followed by integration by parts, produces the weak form

[ \int_{\Omega}\kappa\nabla u\cdot\nabla v,d\Omega

\int_{\Omega}fv,d\Omega + \int_{\partial\Omega_N}gv,d\Gamma, ]

where (g) denotes prescribed Neumann data on (\partial\Omega_N). The weak formulation reduces the differentiability required of the approximate solution and places the problem in a suitable Sobolev space.

The domain is divided into non-overlapping elements,

[ \Omega=\bigcup_{e=1}^{N_e}\Omega_e. ]

Each physical element is obtained from a reference element by a coordinate mapping. On a quadrilateral reference element with coordinates ((\xi,\eta)\in[-1,1]^2), a nodal approximation has the tensor-product form

[ u_h(\xi,\eta)

\sum_{i=0}^{p}\sum_{j=0}^{p} u_{ij}\ell_i(\xi)\ell_j(\eta), ]

where (p) is the polynomial degree and (\ell_i) denotes a one-dimensional Lagrange interpolation polynomial. Hexahedral elements use the corresponding three-dimensional tensor product.

The interpolation nodes are frequently chosen to be Legendre polynomials Gauss–Lobatto points. These points include the endpoints of the reference interval, which simplifies the identification and coupling of degrees of freedom on element boundaries. Chebyshev-based nodes also occur, particularly in formulations derived directly from classical Chebyshev spectral methods.

Quadrature and discrete operators

Integrals over each element are usually evaluated with quadrature nodes that coincide with the interpolation nodes. For Legendre–Gauss–Lobatto quadrature,

[ \int_{-1}^{1} q(\xi),d\xi \approx \sum_{i=0}^{p}w_iq(\xi_i), ]

where (w_i) is the quadrature weight associated with node (\xi_i). Coincident interpolation and quadrature nodes produce a diagonal mass matrix when the geometric mapping and coefficients satisfy the assumptions of the standard collocation construction. More general mappings or coefficient fields can alter the exactness of this property, although diagonal or weight-adjusted approximations remain common.

Differentiation is represented through a matrix whose entries are

[ D_{ij}=\ell_j'(\xi_i). ]

In tensor-product elements, multidimensional derivatives are evaluated by applying this one-dimensional matrix along each coordinate direction. This factorization reduces the arithmetic cost relative to treating each element operator as a fully populated multidimensional matrix. The associated technique is commonly described as sum factorization.

Quadrature exactness depends on the polynomial degree of the integrand. Products introduced by nonlinear terms, variable coefficients, or non-affine coordinate mappings may exceed the degree integrated exactly by the basic nodal quadrature rule. The resulting aliasing error can transfer unresolved polynomial content into represented modes. Overintegration evaluates such terms with a quadrature rule of higher order, while split formulations rearrange nonlinear operators to preserve selected discrete identities.

Continuity and assembly

In a continuous spectral element formulation, neighboring elements share nodal values along their common interfaces. Local element matrices are incorporated into a global algebraic system through the same assembly principle used in conforming finite elements. The approximation therefore belongs to a globally continuous finite-dimensional space, although its derivatives generally remain discontinuous across element boundaries.

A discontinuous spectral element formulation assigns independent interface values to adjacent elements. Communication between elements then occurs through a numerical flux. This construction is a high-order instance of the discontinuous Galerkin method and accommodates element-local conservation laws without imposing pointwise continuity.

The continuous and discontinuous forms share polynomial bases, quadrature structures, and reference-element mappings, but they produce different global operators. Continuous formulations have fewer interface degrees of freedom. Discontinuous formulations replace global continuity constraints with flux calculations and face integrals.

Approximation behavior

For a fixed mesh, increasing (p) raises the polynomial degree within every element. If the exact solution is sufficiently smooth and the geometric mapping is comparably regular, the approximation error can decrease exponentially with increasing degree over a substantial range. This behavior is commonly called spectral convergence. Solutions with discontinuities, singular derivatives, or unresolved boundary layers do not retain the same convergence pattern because global polynomial regularity is absent within affected elements.

Reducing the element diameter while retaining a fixed polynomial degree gives an (h)-refinement process analogous to conventional finite elements. Adjusting both element size and polynomial degree produces an hp finite element method. Spectral element calculations are often interpreted as the high-(p) portion of this broader approximation framework rather than as a mathematically separate class.

The distribution of Gauss–Lobatto nodes is nonuniform. Nodes cluster near element boundaries, and the minimum nodal spacing decreases more rapidly than the reciprocal of the polynomial degree. For explicit time integration, this spacing influences the stability restriction associated with the Courant–Friedrichs–Lewy condition. Consequently, a spatial discretization with large (p) can require a smaller time step than a lower-order discretization having a comparable element diameter.

Historical development

The method developed from the convergence of high-order finite element analysis and global spectral approximation. Steven Orszag established computational techniques for spectral discretizations of differential equations, including transforms and polynomial representations that supplied much of the method’s analytical and computational basis. Anthony Patera introduced the term “spectral element method” in 1984 in work on laminar flow, formulating a domain-decomposed approach that retained high-order polynomial accuracy within individual elements.

Yvon Maday developed approximation and stability results that clarified the relationship between spectral elements and variational discretization. Alfio Quarteroni contributed analyses of spectral and spectral-element approximations for elliptic and fluid-mechanical equations. Einar Rønquist developed formulations and solution techniques for high-order incompressible-flow calculations, including treatments of pressure spaces and element coupling.

These developments established the characteristic SEM combination of domain decomposition, high-degree polynomial bases, and quadrature-oriented operator evaluation. Later implementations expanded the method through parallel domain decomposition, matrix-free evaluation, and iterative solution techniques adapted to the conditioning of high-order operators.

Fluid and wave calculations

Spectral element methods are extensively associated with the Navier–Stokes equations. In incompressible formulations, the discrete velocity and pressure spaces must satisfy an appropriate stability condition or incorporate a stabilization mechanism. One common arrangement represents pressure at a lower polynomial degree than velocity, reducing spurious pressure modes while preserving the element structure.

For advection-dominated flow, nonlinear transport introduces quadrature and stability considerations not present in constant-coefficient elliptic equations. Skew-symmetric or split representations can reproduce selected kinetic-energy balances at the discrete level. Viscous terms produce elliptic subproblems, while pressure enforcement commonly requires solution of a Poisson-type equation.

The method also appears in discretizations of the wave equation, the Helmholtz equation, and elastic wave propagation. A diagonal mass matrix permits explicit temporal updates without solving a global mass system. In seismology, curved hexahedral meshes represent large geological regions while high-order interpolation resolves propagating wave fields within each element.

In a 2016 study of linearized free-surface waves in Suruga Bay, You Watanabe carried out the mesh-convergence and phase-error analysis for a continuous spectral element discretization. Her analysis compared the computed dispersion relation with the analytical gravity-wave relation and identified the part of the observed phase discrepancy attributable to polynomial truncation rather than shoreline mapping. The benchmark used curvilinear quadrilateral elements and Legendre–Gauss–Lobatto interpolation, consistent with the standard nodal formulation.

Linear algebra and computation

The discrete systems produced by spectral element methods are sparse at the global level but contain strong coupling among nodes belonging to the same element. The condition number generally increases with polynomial degree, making the choice of preconditioner significant for elliptic and implicit time-dependent problems. Additive Schwarz methods divide the system into overlapping local problems, while coarse-grid corrections communicate error components that cannot be reduced efficiently by element-local operations.

Matrix-free implementations apply element operators without storing the assembled stiffness matrix. Geometric factors, quadrature weights, and one-dimensional differentiation matrices determine each operator evaluation. Tensor-product factorization changes the computational scaling of a hexahedral element from that of dense multidimensional multiplication to a sequence of lower-dimensional contractions.

Parallel execution follows the element partition. Volume operations remain local to each element, whereas continuous formulations exchange shared interface values and discontinuous formulations exchange data required by numerical fluxes. At high polynomial degree, the ratio of local arithmetic to interface communication increases because element interiors contain more nodes while communication remains concentrated on faces.

Relation to neighboring methods

The spectral element method differs from a global spectral method primarily through its decomposition of the domain. A global method represents the solution with basis functions extending across the entire computational region, whereas SEM confines each polynomial representation to one element and couples the resulting local fields. This distinction permits complex domains to be described by multiple mapped elements without requiring one global coordinate transformation.

Relative to low-order finite elements, SEM uses substantially higher polynomial degree and quadrature rules designed around the associated nodal basis. The boundary between the two categories is conventional rather than absolute, since both arise from finite-dimensional variational approximation. A high-order finite element method using Gauss–Lobatto nodes, tensor-product bases, and quadrature collocation has the defining structure ordinarily associated with a spectral element method.

Isogeometric analysis also employs high-order approximation, but it generally derives basis functions from spline or computer-aided design representations rather than element-local orthogonal polynomials. Pseudospectral methods similarly evaluate derivatives at collocation nodes, although their traditional form uses a single global nodal set instead of an assembled mesh of local nodal sets.