High-Order Numerical Methods: When to Use Spectral, DG, and hp-FEM

TL;DR

High-order numerical methods deliver significantly higher accuracy per degree of freedom than traditional low-order approaches (FDM, FEM, FVM) by using higher-degree polynomial basis functions. There’s no single “best” method — the right choice depends on your problem’s geometry, solution smoothness, and conservation requirements. Spectral methods excel on simple domains with smooth solutions, discontinuous Galerkin (DG) unifies FEM flexibility with FVM-style conservation, spectral element methods combine spectral accuracy with complex geometry handling, and hp-FEM adapts both mesh size and polynomial order for optimal convergence across mixed-smoothness problems.

What Are High-Order Numerical Methods?

If you’ve spent time with finite difference (FDM), finite element (FEM), or finite volume (FVM) methods, you’re working with low-order discretizations — typically linear or quadratic polynomials per element. These methods are straightforward, robust, and work well on a wide range of problems. But they require fine meshes to achieve high accuracy, which means more degrees of freedom, more memory, and more computational cost.

High-order numerical methods take a different approach. Instead of refining the mesh (adding more cells), they increase the polynomial order within each cell. A third-degree polynomial already captures curvature that a linear polynomial misses, and a tenth-degree polynomial can resolve fine features with far fewer cells than a low-order mesh.

The result? Exponential convergence for smooth solutions instead of algebraic convergence, lower numerical dissipation, and — crucially — a method that scales much better on modern GPU hardware. The tradeoff is greater algorithmic complexity, more sophisticated stabilization techniques, and often more opaque error estimation.

This article picks up where our “FEM vs FVM vs FDM” guide leaves off. That piece covers the three traditional low-order methods. Here, we go deeper into the methods that modern computational scientists use when low-order approaches become too expensive or too inaccurate.

Spectral Methods: Exponential Accuracy on Simple Domains

Spectral methods are the oldest high-order approach, and they deliver the most dramatic convergence rates of any method we’ll discuss. Instead of piecewise polynomials on small cells, spectral methods use global basis functions — Fourier series, Chebyshev polynomials, or Legendre polynomials — defined over the entire domain.

For problems with smooth solutions on simple geometries, a spectral method with just 30–50 modes can match the accuracy of a low-order FEM mesh with hundreds of thousands of elements. The convergence is exponential: doubling the polynomial order doesn’t just halve the error — it can reduce the error by orders of magnitude.

This is why spectral methods dominate in fields where the geometry is simple and the physics is smooth: direct numerical simulation of turbulence, atmospheric modeling, and wave propagation in homogeneous media. The NUMATH group at Universidad Politécnica de Madrid has produced excellent visual comparisons showing how p-refinement (increasing polynomial order) reduces dissipation far more efficiently than h-refinement (adding cells) [NUMATH].

The catch: global basis functions require simple, structured meshes. You cannot easily wrap a Fourier series around a complicated aerodynamic shape. Spectral methods also struggle with discontinuities — a single jump in the solution ripples across the entire domain through the Gibbs phenomenon.

For complex geometries, you need the spectral element method, which we cover next. If you want a deeper dive into spectral method theory and implementation, see our guide to “Spectral Methods for PDE Simulation.”

Discontinuous Galerkin Methods: Flexibility Meets Conservation

Discontinuous Galerkin (DG) methods represent one of the most flexible frameworks in modern numerical methods. They combine two desirable properties that were previously hard to obtain together:

1. FEM-style geometric flexibility — unstructured meshes, complex boundaries, adaptive refinement
2. FVM-style conservation — exact local conservation through numerical fluxes at cell interfaces

The key insight is simple but powerful: DG drops inter-element continuity. Traditional FEM enforces continuity across element boundaries; DG replaces it with numerical fluxes, much like finite volume methods do. This discontinuity allows the method to handle discontinuous solutions (shock waves, material interfaces) without special treatment, and it enables explicit time stepping that’s cheap to parallelize.

The Wikipedia entry on DG methods provides a thorough mathematical formulation, including the scalar hyperbolic conservation law and the direct DG formulation for diffusion equations [Wikipedia DG].

For shock-capturing applications, hybrid limiters are essential. You blend a high-order reconstruction (accurate in smooth regions) with a low-order fallback (robust near discontinuities). This keeps DG from oscillating near shocks while preserving its high-order accuracy elsewhere.

If you’re working with hyperbolic conservation laws — fluid dynamics, gas dynamics, multiphase flow — DG is one of the most mature and well-tested approaches. See our “Discontinuous Galerkin Method Tutorial” for the derivation and Python implementation.

Spectral Element Methods: Best of Both Worlds

Spectral element methods (SEM) solve the geometry problem that holds spectral methods back. Instead of a single global basis over the whole domain, SEM breaks the domain into elements, each with its own local high-order polynomial expansion. Within each element, the basis functions are high-degree (often 3rd–15th order); across elements, the solution is only C⁰ continuous — just like standard FEM.

This gives you spectral accuracy on complex geometries. The GLL (Gauss–Legendre–Lobatto) node distribution is the standard choice for spectral element discretizations. GLL nodes cluster near element boundaries, which reduces the Gibbs phenomenon and provides excellent conditioning for the mass matrix.

The Nektar++ implementation paper bridges traditionally low-order FEM and spectral methods, providing essential guidance on matrix-free algorithms and node distributions [VoShKi-09].

Production-level solvers that use SEM:

• Nek5000 (Argonne) — highly scalable spectral element solver for thermal-hydraulics and fluid flows, used in nuclear engineering and combustion
• Nektar++ — open-source library for spectral/DG methods, widely used in academia
• Firedrake — high-level Python framework built on FEniCS, supporting SEM with automatic differentiation and adaptive solvers
• FEniCS — general-purpose FEM/DG solver, supports spectral elements alongside standard formulations

The Haubold 2023 thesis demonstrates how to implement both finite and spectral/hp element methods for optimal performance, covering mesh handling, node distributions, and matrix assembly strategies [Haubold 2023].

High-Order Finite Volume: Conservation With Spectral Resolution

Traditional finite volume methods are built on cell-centered or vertex-based reconstructions with linear or quadratic interpolants. High-order finite volume methods extend this idea to higher-degree polynomial reconstructions within each cell.

The advantage is immediate: you get spectral-like resolution while preserving the local conservation properties that make FVM the default choice for conservation laws. This is the practical hybrid that bridges traditional FVM with high-order accuracy.

High-order FVM methods are particularly useful when you need exact mass/energy/momentum conservation (from the finite volume formulation) but want to reduce numerical dissipation (from the high-order reconstruction). They’re less commonly discussed than DG or spectral elements, but they’re increasingly popular in CFD applications where conservation is non-negotiable but accuracy is equally important.

hp-FEM and p-Adaptivity: Choosing Your Convergence Strategy

This is where things get really interesting. The convergence behavior of high-order methods depends critically on which refinement strategy you use.

h-Refinement (Mesh Sizing)

Reduce element size, increase the number of elements. This is what low-order FEM practitioners do. Convergence is algebraic: error scales as O(h^p) where p is the polynomial degree. You get steady improvement — double the resolution, roughly halve the error.

p-Refinement (Polynomial Order)

Keep the mesh fixed, increase the polynomial degree within each element. For smooth solutions, convergence is exponential — the error drops like O(e^(-βp)) where β depends on solution smoothness. This is dramatically more efficient than h-refinement when the solution is smooth.

But here’s the critical caveat: p-refinement fails catastrophically around strong localized singularities. If your problem has a point singularity, a crack tip, or a discontinuity, increasing polynomial order won’t help — the solution isn’t smooth enough. You need h-refinement in those regions.

hp-Adaptivity (Both Strategies)

hp-FEM combines both approaches: it refines the mesh (h) where the solution is rough, and increases polynomial order (p) where the solution is smooth. Bringmann 2024 provided a rigorous proof showing that goal-oriented adaptive hp-FEM achieves full linear convergence with optimal complexity — the exponential rate of p-adaptivity combined with the flexibility of h-refinement [Bringmann 2024].

The NIST documentation on hp-FEM provides practical examples showing why hp-adaptivity achieves exponential convergence with respect to the number of degrees of freedom [NIST hp-FEM].

What I’d Choose

If your problem is smooth everywhere (heat equation, Laplace, wave propagation in homogeneous media), use p-refinement — exponential convergence with minimal mesh changes. If your problem has mixed smoothness (multiphase flow with material interfaces, shock-capturing, problems with singularities), use hp-adaptivity. And if your problem has strong localized singularities, h-refinement is mandatory near those features regardless of what you do elsewhere.

Practical Implementation: What You Actually Need to Build

The theory is elegant. The implementation is where most practitioners stumble. Here are the three practical details that matter most:

1. GLL Node Distributions

For spectral element discretizations, the choice of nodes determines both accuracy and conditioning. GLL nodes — Gauss–Legendre–Lobatto quadrature points — cluster near element boundaries, which reduces the Gibbs phenomenon and provides good mass matrix conditioning. They’re the default in most production codes.

2. Matrix-Free Algorithms

This is where modern high-order implementations diverge sharply from older approaches. Matrix-free algorithms (on-the-fly matrix-vector products) compute the action of the system matrix without ever assembling it. Instead, you evaluate the weak form directly at each iteration step.

Why does this matter? Because the full stiffness matrix for a high-order method with polynomial degree p on a single element is a p × p dense block. For p = 10, that’s 100 entries per element — and you multiply by the number of elements. Matrix-free evaluation avoids storing and solving with massive dense matrices, and it dramatically improves GPU/CPU cache locality. The Nektar++ implementation paper demonstrates this pattern extensively [VoShKi-09].

If you want to understand the GPU implementation side, see our “GPU Kernel Programming for Custom Physics Simulation.”

3. Spectral Vanishing Viscosity (SVV) Stabilization

Without stabilization, spectral methods produce spurious high-frequency oscillations — the same Gibbs phenomenon we discussed earlier. The standard stabilization technique for spectral/hp methods is spectral vanishing viscosity (SVV).

SVV selectively damps spurious high-frequency modes without degrading the accuracy of the low-frequency physical solution. The mechanism works by adding a controlled artificial viscosity that vanishes at low frequencies (preserving accuracy) and grows at high frequencies (damping oscillations).

Kirby & Sherwin’s 2006 paper is the standard reference for SVV stabilization, cited 151 times and considered the foundational guide for spectral/hp method stability [Kirby & Sherwin 2006].

When to Choose Each Method

This is the practical comparison. The table below maps the five high-order methods against the three traditional methods across the dimensions that matter most for method selection.

Method Comparison Table

Key takeaways from the table:

• Conservation matters: If you need exact local conservation (fluid dynamics, multiphase flow), DG or high-order FVM are your only high-order options. Spectral methods and spectral elements lack strict local conservation.
• Geometry matters: Spectral methods are off-limits for complex geometries. DG, spectral elements, and hp-FEM handle unstructured meshes.
• Smoothness matters: p-adaptivity delivers exponential convergence but fails near discontinuities. hp-adaptivity handles mixed smoothness. h-refinement handles singularities.
• Hardware matters: DG and spectral elements benefit most from GPU acceleration because their block-diagonal mass matrices yield local computations. Low-order methods don’t get the same hardware advantage.

My Recommendations

For smooth problems on simple domains: Spectral methods

If your geometry is a rectangle, circle, or sphere, and your solution is smooth everywhere, spectral methods will destroy low-order approaches in accuracy per degree of freedom. Just don’t try to force them onto complex meshes — that’s the #1 reason practitioners give up on spectral methods.

For conservation laws with complex geometry: Discontinuous Galerkin

DG is the most versatile high-order method for practical engineering problems. It handles unstructured meshes, enforces exact local conservation, and deals gracefully with discontinuities. If you’re doing CFD or gas dynamics, this is my go-to.

For complex geometry with high accuracy: Spectral element methods

SEM gives you spectral accuracy without sacrificing geometry flexibility. The GLL node distribution, matrix-free algorithms, and SVV stabilization make it production-ready. Nek5000, Nektar++, and Firedrake are all battle-tested SEM solvers.

For mixed smoothness problems: hp-FEM

Problems with smooth regions and localized singularities (crack propagation, multiphase interfaces, shock-capturing) benefit from hp-adaptivity. The Bringmann 2024 convergence proof confirms the theoretical foundation — exponential rates with optimal complexity.

Why High-Order Methods Matter Now

The computational landscape has shifted dramatically. GPU-dominant architectures are now common, and high-order methods are uniquely suited for this hardware because their block-diagonal mass matrices (DG, SEM) yield local computations that map beautifully to GPU threads. Low-order methods don’t get the same acceleration benefit.

Machine learning is also entering the space. Recent work (2023–2024) explores using deep neural networks as “experts” for quasi-optimal hp-refinement strategies, bypassing expensive heuristic searches while preserving exponential convergence [Hafeez et al. 2023].

If you’re still using linear FEM on a GPU without realizing how much more accurate your simulation could be with polynomial degree 5–10, you’re leaving accuracy on the table. The convergence rates alone justify the implementation effort.

Summary + Next Steps

High-order numerical methods aren’t a panacea — they require careful method selection, sophisticated stabilization, and often more complex implementation. But when the problem calls for it, they deliver accuracy per degree of freedom that low-order methods simply cannot match.

To decide which method serves your problem best:

1. Check conservation requirements — if exact local conservation is needed, choose DG or high-order FVM
2. Check geometry complexity — if complex geometry, avoid pure spectral methods
3. Check solution smoothness — if smooth, p-adaptivity is efficient; if mixed, hp-adaptivity; if singularities, h-refinement is mandatory near those features
4. Check hardware — if GPU-accelerated, DG and SEM get the biggest performance boost

Related Guides on MatForge

• FEM vs FVM vs FDM: A Practical Comparison for Beginners — foundational low-order method comparison
• Spectral Methods for PDE Simulation: Theory, Implementation, and When to Use Them — deep dive into spectral methods
• Discontinuous Galerkin Method Tutorial: Derivation, Fluxes, and Python Implementation — DG method walkthrough
• Mesh Quality and Convergence Studies: A Hands-On Guide — convergence testing practices
• Verification vs Validation in Scientific Simulations: A Practical Guide — V&V framework for code verification

What to Read Next

If you’re ready to move from understanding to implementation, the next steps are concrete: pick a solver (Nektar++ for research, Nek5000 for production), choose your polynomial order based on the smoothness analysis above, and implement matrix-free evaluation to avoid the dense-matrix bottleneck. The mesh-free approach to assembling the weak form is what separates modern high-order codes from legacy implementations.

If you want to dig deeper into verification methodology — truncation error estimation, the method of manufactured solutions, and code verification workflows — check our upcoming coverage of error estimation and verification methods.

This article complements our existing FEM vs FVM vs FDM guide by exploring the methods that computational scientists use when low-order approaches become too expensive or inaccurate. High-order methods aren’t always necessary, but when they are, understanding when and why to use each variant is the difference between an efficient simulation and one that wastes compute cycles.