Fast and Well-Conditioned Spectral Method

This repository contains a Python reimplementation of the ultraspherical spectral method for solving boundary value problems (BVPs), originally proposed by Olver & Townsend (2013), as part of the following White Paper!

Overview

Classical spectral methods solve differential equations by expanding the unknown solution in global basis functions (e.g., Chebyshev polynomials), achieving super-algebraic convergence for smooth solutions. However, they often suffer from poor conditioning and dense linear systems at high resolution.

The ultraspherical spectral method resolves this tension by:

• Expressing differentiation sparsely via basis changes between ultraspherical polynomial families
• Producing almost-banded linear systems (banded interior + dense boundary rows)
• Using QR factorization with Givens rotations that exploits the filled-in structure
• Applying a diagonal preconditioner that yields bounded condition numbers

We implement the full discretization framework, the fast almost-banded QR algorithm, and replicate key numerical experiments from the original paper.

Features

• Chebyshev and ultraspherical polynomial coefficient-space operators
• Sparse differentiation operators $\mathcal{D}_\lambda$
• Banded conversion operators $\mathcal{S}_\lambda$ between ultraspherical bases
• Multiplication operators $\mathcal{M}_\lambda[a]$ for variable coefficients
• Almost-banded QR factorization via Givens rotations
• Diagonal preconditioning for bounded condition numbers
• Adaptive truncation with residual monitoring
• Data caching for efficient re-plotting

Numerical Experiments

We replicate Figures 3.1–3.3 from Olver & Townsend (2013):

1. Airy equation ($\varepsilon = 10^{-9}$): Highly oscillatory solution with spectral convergence
2. Condition number study: Demonstrates bounded conditioning with diagonal preconditioning
3. Boundary layer problem ($\varepsilon = 10^{-7}$): Interior layers with sharp transitions

Results

Our implementation achieves:

• Spectral (super-algebraic) convergence: Cauchy error drops from $10^5$ to $10^{-15}$
• Bounded condition number $\kappa_2(A_n R_n) = O(1)$ with preconditioning
• Clean resolution of stiff boundary layers at extreme $\varepsilon$ values

Setup

Install dependencies:

uv pip install numpy scipy matplotlib

Run experiments:

python ultraspherical_experiments.py

Figures are saved to figs/. Computed data is cached in data/ for fast re-plotting.

References

• S. Olver and A. Townsend, "A Fast and Well-Conditioned Spectral Method," SIAM Review 55(3), pp. 462–489, 2013. DOI: 10.1137/120865458
• C. Canuto et al., Spectral Methods: Fundamentals in Single Domains, Springer, 2006.
• N. J. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, 2002.