Spectral methods

Deep dive

How it works

In: N grid values of a smooth field (periodic, or on [−1, 1] at Chebyshev points) · Out: its derivative or the PDE solution at the same points · State: none for a derivative; the field and time step for a PDE run · Cost: O(N log N) per derivative (two FFTs)

Fourier pseudo-spectral on [0, 2π), N points x_j = 2πj/N (N even):

1. FFT grid values u_j → û_k, k = 0..N−1, interpreted as wavenumbers 0, 1, …, N/2−1, −N/2 (Nyquist), …, −1.
2. Derivative: multiply û_k by i·k; zero the Nyquist coefficient (its i·k derivative is not real-valued on the grid).
3. Inverse FFT, divide by N.
4. Nonlinear term u·u_x: compute both on the grid, multiply pointwise, transform back. The product of two band-limited fields has modes up to 2·k_max, which fold back (alias) onto low modes. Zeroing |k| > N/3 in both factors (2/3 rule) keeps the aliased part outside the kept range.
5. Time stepping: method of lines with RK4 or an IMEX scheme treating the linear stiff term (diffusion −νk²) exactly via integrating factor e^(−νk²t).

Chebyshev (non-periodic [−1, 1]): sample at x_j = cos(πj/N), take a DCT to get a_k, use the recurrence for derivative coefficients, impose boundary conditions by replacing rows (collocation) or by basis recombination (Galerkin). Points cluster near the ends with spacing O(1/N²), which is where the O(N⁻⁴) diffusion limit comes from.

Code

#[derive(Clone, Copy, Debug, Default)]
pub struct C { pub re: f64, pub im: f64 }
impl C {
    fn mul(self, o: C) -> C { C { re: self.re * o.re - self.im * o.im, im: self.re * o.im + self.im * o.re } }
}

/// In-place iterative radix-2 FFT; `inverse` uses e^{+i}, without the 1/N factor.
pub fn fft(a: &mut [C], inverse: bool) {
    let n = a.len();
    assert!(n.is_power_of_two());
    let mut j = 0;
    for i in 1..n { // bit-reversal permutation
        let mut bit = n >> 1;
        while j & bit != 0 { j ^= bit; bit >>= 1; }
        j |= bit;
        if i < j { a.swap(i, j); }
    }
    let sign = if inverse { 1.0 } else { -1.0 };
    let mut len = 2;
    while len <= n {
        let ang = sign * 2.0 * std::f64::consts::PI / len as f64;
        for s in (0..n).step_by(len) {
            for k in 0..len / 2 {
                let w = C { re: (ang * k as f64).cos(), im: (ang * k as f64).sin() };
                let (u, v) = (a[s + k], a[s + k + len / 2].mul(w));
                a[s + k] = C { re: u.re + v.re, im: u.im + v.im };
                a[s + k + len / 2] = C { re: u.re - v.re, im: u.im - v.im };
            }
        }
        len <<= 1;
    }
}

/// Signed wavenumber of FFT bin k on a 2*pi-periodic grid of n points.
fn wavenumber(k: usize, n: usize) -> f64 { if k <= n / 2 { k as f64 } else { k as f64 - n as f64 } }

/// Spectral derivative of order `m` of periodic samples on [0, L).
pub fn spectral_derivative(u: &[f64], m: u32, length: f64) -> Vec<f64> {
    let n = u.len();
    let mut a: Vec<C> = u.iter().map(|&x| C { re: x, im: 0.0 }).collect();
    fft(&mut a, false);
    let scale = 2.0 * std::f64::consts::PI / length;
    for k in 0..n {
        let kk = wavenumber(k, n) * scale;
        if m % 2 == 1 && k == n / 2 { a[k] = C::default(); continue; } // Nyquist: no real odd derivative
        let mut f = C { re: 1.0, im: 0.0 };
        for _ in 0..m { f = f.mul(C { re: 0.0, im: kk }); } // (ik)^m
        a[k] = a[k].mul(f);
    }
    fft(&mut a, true);
    a.iter().map(|c| c.re / n as f64).collect()
}

/// Solve u'' = f on a periodic domain (mean of f must be ~0); returns the zero-mean solution.
pub fn periodic_poisson(f: &[f64], length: f64) -> Vec<f64> {
    let n = f.len();
    let mut a: Vec<C> = f.iter().map(|&x| C { re: x, im: 0.0 }).collect();
    fft(&mut a, false);
    let scale = 2.0 * std::f64::consts::PI / length;
    for k in 0..n {
        let kk = wavenumber(k, n) * scale;
        a[k] = if k == 0 { C::default() } else { C { re: -a[k].re / (kk * kk), im: -a[k].im / (kk * kk) } };
    }
    fft(&mut a, true);
    a.iter().map(|c| c.re / n as f64).collect()
}

/// Dealiased pseudo-spectral product u * v (2/3 rule): truncate |k| > n/3 before and after.
pub fn dealiased_product(u: &[f64], v: &[f64]) -> Vec<f64> {
    let n = u.len();
    let cut = |a: &mut [C]| for k in 0..n { if wavenumber(k, n).abs() > n as f64 / 3.0 { a[k] = C::default(); } };
    let to_spec = |x: &[f64]| { let mut a: Vec<C> = x.iter().map(|&r| C { re: r, im: 0.0 }).collect(); fft(&mut a, false); a };
    let (mut a, mut b) = (to_spec(u), to_spec(v));
    cut(&mut a); cut(&mut b);
    fft(&mut a, true); fft(&mut b, true);
    let mut p: Vec<C> = (0..n).map(|j| C { re: a[j].re * b[j].re / (n * n) as f64, im: 0.0 }).collect();
    fft(&mut p, false);
    cut(&mut p);
    fft(&mut p, true);
    p.iter().map(|c| c.re / n as f64).collect()
}

Complexity and constants

• One derivative = 2 FFTs of N points: about 2 × 5N log₂N flops. A 3D DNS step at 1024³ is dominated by FFTs and their all-to-all transposes across nodes.
• Measured with the code above, derivative of e^(sin x) on [0, 2π): max error 4.3·10⁻³ at N = 8, 1.8·10⁻⁷ at N = 16, 5.6·10⁻¹⁵ at N = 32 (machine precision); second-order central differences give 9.4·10⁻² and 2.6·10⁻² at N = 16 and 32. Squaring cos(5x) on 16 points without dealiasing puts a spurious 0.25 into mode 6 (cos 10x folded back); the 2/3-rule product returns exactly 0.5.
• Explicit stability (method of lines): Fourier advection u_t + c u_x = 0 needs c·Δt ≤ O(Δx) with a constant set by the time stepper's stability region on the imaginary axis (RK4: |λΔt| ≤ 2√2, λ_max = c·N/2); Fourier diffusion ν·Δt·(N/2)² ≤ 2.78 for RK4; Chebyshev diffusion scales with N⁴.
• Resolution: Fourier needs just over 2 points per wavelength, Chebyshev about π (Trefethen); low-order finite differences need several times more for the same phase accuracy.

Efficiency

Cost: one derivative is 2 FFTs, about 2 × 5N log₂N flops (order of), so a Fourier or Chebyshev time step costs a few FFTs per right-hand side evaluation. Bound: FFT arithmetic on one node, memory bandwidth and the all-to-all transposes in multi-node 3D runs (1024³ DNS).

Levers:

• Spectral instead of second-order differences — error 5.6·10⁻¹⁵ at N = 32 against 2.6·10⁻² for central differences (measured above) — smooth, analytic data.
• Real-input FFTs (r2c/c2r) — about half the work and memory — real fields.
• Integrating factor / ETDRK4 for the linear term — step set by the nonlinear term only (Kassam-Trefethen) — stiff diffusion.
• 2/3 truncation or 3/2 padding — exact for quadratic terms, removes aliasing — pseudo-spectral products.
• Implicit Chebyshev diffusion — avoids the O(N⁻⁴) step limit — non-periodic boxes.

Naive → standard → state of the art:

• Dense differentiation matrix · O(N²) per derivative
• FFT-based pseudo-spectral with dealiasing · O(N log N) per derivative
• Slab/pencil-decomposed distributed FFTs, GPU spectral elements (NekRS) · multi-node scaling

Implementation notes

• Real-input FFTs (r2c/c2r) halve the work and memory; store N/2+1 complex coefficients.
• Wavenumber ordering: bins N/2+1..N−1 are negative frequencies; getting this wrong gives a derivative that is right for low modes and wildly wrong for the upper half.
• Integrating factor / ETDRK4 for stiff linear parts (Kassam-Trefethen): exact treatment of −νk² allows steps limited by the nonlinear term only.
• Dealiasing: 2/3 truncation or 3/2 padding (transform to 3N/2 points for the product); both exact for quadratic terms.
• Filtering: exponential filters on the top modes stabilize long runs with marginal resolution, at the price of formal accuracy.
• Non-periodic: Chebyshev via DCT-I, or Legendre with Gauss-Lobatto quadrature; boundary conditions by tau method or basis recombination (e.g. φ_k = T_k − T_{k+2} for Dirichlet).

Pitfalls

• Forgetting the 1/N normalization (or applying it twice) in the inverse FFT.
• Keeping the Nyquist mode for first derivatives: produces a non-real or oscillating result for real input.
• Applying spectral methods to non-smooth data: Gibbs ringing does not shrink with N (overshoot near a jump stays about 9 %).
• Non-periodic data in a Fourier method: the implicit periodic extension has a jump at the boundary; accuracy collapses to first order.
• Explicit time steps chosen from the grid spacing alone: spectral differentiation has larger eigenvalues than FD (factor up to π for first derivatives), so the CFL constant is smaller.

Variants and alternatives

In production

Dedalus (general spectral PDE framework), Chebfun (MATLAB, Chebyshev technology), ECMWF IFS (spectral transform on the sphere), turbulence DNS codes (pseudo-spectral in periodic boxes), FFTW and cuFFT as the transforms underneath, Nek5000/NekRS (spectral elements).

Implementations

Sources

• Lloyd N. Trefethen, Spectral Methods in MATLAB, SIAM, 2000. book
• John P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd ed., Dover, 2001.
• Claudio Canuto, M. Yousuff Hussaini, Alfio Quarteroni, Thomas A. Zang, Spectral Methods: Fundamentals in Single Domains, Springer, 2006. book
• Steven A. Orszag, On the Elimination of Aliasing in Finite-Difference Schemes by Filtering High-Wavenumber Components, Journal of the Atmospheric Sciences 28(6), 1971. paper
• Aly-Khan Kassam, Lloyd N. Trefethen, Fourth-Order Time-Stepping for Stiff PDEs, SIAM Journal on Scientific Computing 26(4), 2005. paper
• Keaton J. Burns, Geoffrey M. Vasil, Jeffrey S. Oishi, Daniel Lecoanet, Benjamin P. Brown, Dedalus: A flexible framework for numerical simulations with spectral methods, Physical Review Research 2, 023068, 2020. paper · code