Spectral Methods and High-Order Numerical Methods 2026

 

Spectral Methods and High-Order Numerical Methods

Proposed 2026 English Course Outline

 

Part I — Approximation and spectral representations

Week 1 — Foundations of high-order discretization

• Model PDE classes: elliptic, parabolic, hyperbolic, dispersive

• Weighted-residual framework

• Galerkin, Petrov–Galerkin, tau, collocation and least-squares formulations

• Approximation error, consistency, stability and convergence

• Sobolev regularity versus spectral convergence

• Modal versus nodal representations

• Reproducible numerical experiments in Julia, Python or MATLAB

Week 2 — Fourier approximation and the FFT

• Fourier series and transforms

• Truncation, interpolation and projection

• Parseval identities

• Trigonometric interpolation

• Discrete Fourier transform and FFT

• Spectral differentiation and integration

• Aliasing, convolution and the two-thirds/three-halves de-aliasing rules

• Gibbs phenomena and filtering

Week 3 — Fourier methods for periodic PDEs

• Fourier–Galerkin and Fourier collocation methods

• Poisson, heat, wave and Helmholtz equations

• Burgers, Korteweg–de Vries, Kuramoto–Sivashinsky and Allen–Cahn equations

• Conservation laws and invariant drift

• Semi-discrete stability

• Explicit, implicit, IMEX and exponential time integrators

Part II — Orthogonal polynomials and fast transforms

Week 4 — General orthogonal-polynomial systems

• Weighted inner products

• Three-term recurrence and Jacobi matrices

• Zeros and interlacing

• Gaussian, Radau and Lobatto quadrature

• Christoffel–Darboux kernels

• Reproducing and projection kernels

• Modified weights and Christoffel transformations

• Stable polynomial evaluation

Week 5 — Jacobi, Legendre and Chebyshev systems

• Jacobi Sturm–Liouville operator

• Self-adjointness and weighted integration by parts

• Parameter-shift and derivative identities

• Legendre and Chebyshev specializations

• Connection coefficients between polynomial families

• Endpoint behaviour and singular Jacobi weights

• Polynomial interpolation and projection estimates

This retains the Jacobi material from the original course but makes the kernel, Sturm–Liouville and parameter-transport structure explicit.

Week 6 — Fast algorithms for orthogonal polynomials

• Fast synthesis and analysis

• Discrete cosine and sine transforms

• Fast Chebyshev–Legendre transforms

• Fast connection-coefficient transforms

• Clenshaw and barycentric algorithms

• Structured multiplication and differentiation operators

• Almost-banded matrices

• Adaptive coefficient-space computation

• Multivariate orthogonal polynomials and Koornwinder constructions

A 2025 graduate spectral-methods syllabus already places fast orthogonal-polynomial transforms, connection coefficients, multivariate constructions and exponential integration near the centre of the subject rather than treating them as optional implementation details.

Part III — Boundary-value and eigenvalue problems

Week 7 — Second-order boundary-value problems

• Dirichlet, Neumann, Robin and mixed conditions

• Weak formulations and coercivity

• Legendre- and Chebyshev–Galerkin methods

• Boundary-adapted modal bases

• Tau and collocation formulations

• Equivalence and non-equivalence of formulations

• Sparse stiffness and mass matrices

• Error estimates

Week 8 — Modern sparse spectral solvers

• Ultraspherical and coefficient-space methods

• Banded differentiation and conversion operators

• Integral reformulations

• Spectral integration

• Conditioning of differentiation matrices

• Diagonal and operator preconditioning

• Adaptive truncation

• Residual and coefficient-tail error indicators

• Matrix-free Krylov solution

Modern courses now routinely pair classical collocation with ultraspherical methods, fast structured matrices and conditioning analysis.

Week 9 — Higher-order and constrained problems

• Fourth- and higher-order equations

• Generalized Jacobi bases

• Dual Petrov–Galerkin methods

• Biharmonic and Cahn–Hilliard equations

• Mixed formulations

• Differential-algebraic boundary constraints

• Exact enforcement of multiple boundary conditions

• Nullspaces, compatibility and gauge conditions

Week 10 — Spectral eigenvalue problems

• Sturm–Liouville eigenproblems

• Generalized matrix eigenvalue formulations

• Spurious eigenvalues and spectral pollution

• Non-self-adjoint operators

• Pseudospectra and resolvent growth

• Polynomial and rational filtering

• Contour and shift-invert methods

• Validation by residual and backward error

Part IV — Time evolution, singularity and nonlocality

Week 11 — Time-dependent PDEs

• Method of lines

• CFL restrictions

• Implicit–explicit schemes

• Exponential and integrating-factor methods

• Operator splitting

• Symplectic and energy-preserving integration

• Long-time error and invariant preservation

• Adaptive temporal and spectral resolution

Week 12 — Integral, fractional and nonlocal equations

• Volterra and Fredholm integral equations

• Spectral Nyström methods

• Weakly singular kernels

• Jacobi methods for endpoint singularities

• Fractional derivatives and fractional Sturm–Liouville systems

• Nonlocal diffusion

• Delay differential equations

• Singular quadrature and kernel compression

Week 13 — Unbounded and semi-infinite domains

• Hermite polynomials and functions

• Laguerre polynomials and functions

• Rational Chebyshev bases

• Mapped spectral methods

• Scaling and translation parameters

• Oscillatory and localized solutions

• Transparent and absorbing boundary treatments

• Domain truncation versus native unbounded-domain bases

Part V — Multiple dimensions and complex geometry

Week 14 — Tensor-product and multidimensional methods

• Rectangles, cuboids and periodic boxes

• Tensor-product bases

• Kronecker structure

• Fast diagonalization

• Alternating-direction solvers

• Disks, cylinders, balls and spherical shells

• Spherical harmonics

• Coordinate singularities

• Sparse grids and low-rank tensor representations

Week 15 — Spectral and hp-element methods

• Reference-to-physical element maps

• High-order nodal and modal elements

• Gauss–Lobatto quadrature

• Conforming and discontinuous formulations

• Static condensation

• Sum factorization

• Curved elements and geometric aliasing

• h-, p- and hp-adaptivity

• Mortar and interface coupling

• Parallel and GPU-oriented implementation

Spectral-element methods remain a core bridge between global spectral accuracy and geometric flexibility; established courses continue to organize the subject around Fourier, polynomial and spectral-element stages. (Department of Mathematics)

Week 16 — Stability on curved and nonlinear systems

• Metric identities and discrete geometric conservation laws

• Summation-by-parts structure

• Split forms and entropy stability

• De-aliasing on mapped elements

• Shock detection and spectral viscosity

• Positivity and realizability

• Incompressible Navier–Stokes projection methods

• Transition, turbulence and under-resolved computation

Part VI — Modern extensions and verification

Week 17 — Scientific machine learning and spectral operators

• Solution operators between function spaces

• Fourier and spectral neural operators

• Learning in coefficient space

• Spectral residual losses

• Parseval-based training objectives

• Aliasing and resolution-transfer failure

• Hybrid solver–learner architectures

• Training-data and discretization dependence

• Comparison with classical spectral solvers at equal accuracy

• Stability, conservation and posterior certification

Neural spectral methods now use orthogonal bases and coefficient-space residuals directly, but they should enter after the approximation, stability and solver machinery has been established. They are an extension of the numerical carrier, not a substitute for discretization analysis.

Week 18 — Error certification and computational audit

• A priori versus a posteriori error

• Residual, truncation and quadrature decomposition

• Roundoff and conditioning

• Manufactured solutions

• Grid and polynomial-order refinement

• Cross-formulation verification

• Conservation and stability audits

• Reproducibility and benchmark design

• Cost versus accuracy

• When spectral convergence fails

Week 19 — Final projects and research presentations

Suggested projects:

• Navier–Stokes or magnetohydrodynamics on periodic or curved domains

• High-order Schrödinger eigenvalue computations

• Fractional or nonlocal PDEs

• Spectral-element wave propagation

• Non-normal operator pseudospectra

• Adaptive ultraspherical solvers

• Fast Jacobi connection transforms

• Spectral methods for kinetic equations

• Structure-preserving Cahn–Hilliard or Hamiltonian systems

• Classical spectral solver versus neural operator under equal error and cost constraints

2026 core reading structure

1. Shen, Tang and Wang, Spectral Methods: Algorithms, Analysis and Applications — foundational analysis and classical algorithms.

2. Trefethen, Spectral Methods in MATLAB and Approximation Theory and Approximation Practice — computational approximation and collocation.

3. Olver, Slevinsky and Townsend, “Fast Algorithms Using Orthogonal Polynomials” — modern transform and structured-operator layer. A 2025 graduate course lists this as a primary reference.

4. Kopriva, Implementing Spectral Methods — spectral-element implementation.

5. Karniadakis and Sherwin, Spectral/hp Element Methods for Computational Fluid Dynamics — complex geometry and hp discretization.

6. Selected recent papers on structure-preserving schemes, fractional spectral methods, fast tensor algorithms and neural spectral operators.

The principal 2026 change is not removal of the classical material. It is a reordering:

orthogonal basis
→ fast transform
→ sparse operator
→ stable PDE discretization
→ complex geometry
→ certified computation
→ learned operator extension.

This prevents machine learning, dense collocation matrices or software demonstrations from replacing the approximation and operator structure on which the methods depend.

Spectral Methods and High-Order Numerical Methods — 2026

Throughout, P_N denotes the space of algebraic polynomials of degree at most N; N denotes modal, nodal, or Fourier resolution; h denotes a physical mesh scale where applicable; and ||.||_X denotes the norm of a stated function space X.

Part I — Approximation and spectral representations

Week 1 — Foundations of high-order discretization

Model PDE classes: elliptic, parabolic, hyperbolic, dispersive

A high-order method must first be matched to the analytical class of the differential operator. Elliptic equations, such as

-div(a grad u) + c u = f,

are global boundary-value problems. Their numerical stability is normally inherited from coercivity, ellipticity, or an inf-sup condition. Parabolic equations, such as

u_t - div(a grad u) = f,

combine spatial smoothing with stiff temporal evolution. A spectral spatial discretization can approximate the solution exponentially while simultaneously generating eigenvalues of size O(N^2) or larger, so temporal stability becomes a separate constraint. Hyperbolic equations, such as

u_t + a · grad u = 0,

transport information along characteristics. Their numerical treatment is controlled by boundary fluxes, dispersion, dissipation, and the discrete propagation speed. Dispersive equations, such as

u_t + u_xxx = 0,

do not smooth in the parabolic sense; instead, different Fourier modes travel with different phase velocities. Small phase errors can therefore accumulate into large long-time solution errors.

A useful abstract representation is

M u_t + A(u) = f,

where M is a mass operator and A may be coercive, dissipative, skew-adjoint, Hamiltonian, or nonlinear. The continuous energy identity should guide the discretization. If

dE(u)/dt = D(u) + B(u),

where D is interior production or dissipation and B is a boundary flux, then a structurally faithful semidiscretization should possess a corresponding discrete identity. Spectral approximation alone does not guarantee this. The operator class determines the trial space, test space, boundary treatment, time integrator, and admissible stabilization.

Weighted-residual framework

Consider

L u = f in Omega,
B u = g on partial Omega.

Choose a finite-dimensional trial space

V_N = span{phi_0,...,phi_N}

and represent the approximate solution as

u_N = sum_{j=0}^N c_j phi_j.

The differential residual is

R_N = L u_N - f.

A weighted-residual method imposes

<R_N,w_i> = 0, for i = 0,...,N,

where the functions w_i span a test space W_N. This one construction includes Galerkin, Petrov-Galerkin, tau, collocation, least-squares, discontinuous Galerkin, and many spectral-element formulations.

The test space determines which part of the residual is suppressed. Orthogonality in L2 does not imply smallness in a stronger norm, and a small residual does not automatically imply a small solution error. If the continuous inverse satisfies

||v||_X <= C_L ||L v||_Y,

then

||u-u_N||_X <= C_L ||R_N||_Y

provided the boundary conditions and domains of the operators are correctly matched. The complete numerical object is therefore the triple

trial representation × residual norm × inverse stability.

Integration by parts is often used to lower the differentiation order imposed on the trial functions:

integral (L u_N) v_N = a(u_N,v_N) + boundary terms.

Those boundary terms cannot be discarded without justification. They must be eliminated by strong boundary conditions, incorporated naturally into the weak formulation, or replaced by consistent numerical fluxes.

Galerkin, Petrov-Galerkin, tau, collocation and least-squares formulations

A Galerkin method uses the same trial and test space. For

a(u,v) = l(v),

the approximation satisfies

a(u_N,v_N) = l(v_N) for every v_N in V_N.

If a is symmetric and coercive, the resulting matrix often inherits symmetry and positive definiteness. A Petrov-Galerkin method uses V_N != W_N. This allows the test functions to absorb derivatives, enforce directional stability, or produce sparse matrices. Its central stability condition is

inf_{u_N != 0} sup_{v_N != 0} a(u_N,v_N) / (||u_N||_V ||v_N||_W) >= beta_N.

Uniform stability requires beta_N not to collapse as N grows.

A tau method enforces the differential equation in the lower modes and replaces a number of highest-mode equations by boundary conditions. Collocation instead enforces

L u_N(x_i) = f(x_i)

at selected nodes. It is computationally direct but must be interpreted through interpolation and quadrature theory; pointwise enforcement is not automatically stable in the continuum norm. Least-squares methods minimize

J_N(v) = ||L v-f||_Y^2 + ||Bv-g||_Z^2.

They provide a direct residual measure and normally yield positive-definite systems, but the normal-equation operator behaves like L*L, which can square the condition number. The formulations may use the same polynomial space yet produce different discrete spectra, conditioning, conservation laws, and boundary errors.

Approximation error, consistency, stability and convergence

The approximation error measures what the finite-dimensional space could achieve if the best approximation were known:

E_N(u) = inf_{v_N in V_N} ||u-v_N||.

Consistency measures whether the exact solution satisfies the discrete equations. Stability measures how residual or data perturbations are amplified by the discrete inverse. These components combine schematically as

||u-u_N|| <= C_stab [E_N(u) + consistency defect + data defect].

For a conforming coercive Galerkin method, Cea’s estimate is

||u-u_N||_V <= (M/alpha) inf_{v_N in V_N} ||u-v_N||_V,

where

|a(u,v)| <= M ||u||_V ||v||_V

and

a(v,v) >= alpha ||v||_V^2.

For quadrature-based or nonconforming methods, a Strang estimate adds the defect between the exact and discrete bilinear forms.

Spectral convergence is an approximation property, not a replacement for stability analysis. An exponentially accurate basis combined with an unstable boundary closure still produces divergence. Conversely, a stable discretization cannot converge faster than the regularity of the exact solution permits. Computed errors should therefore be separated into projection error, quadrature error, aliasing error, time-discretization error, linear-solver error, and floating-point error.

Sobolev regularity versus spectral convergence

For a function u in H^s(-1,1), polynomial projection commonly satisfies an algebraic estimate of the form

||u-P_N u||_{H^m} <= C N^(m-s) ||u||_{H^s},

for 0 <= m <= s. Periodic Fourier projection has the corresponding estimate

||u-P_N u||_{H^m} <= C N^(m-s) ||u||_{H^s}.

The coefficient decay reflects this regularity. If u has only s derivatives, the coefficients decay algebraically. If u extends analytically to a complex neighborhood, they decay geometrically:

|u_hat_n| <= C rho^(-n), with rho > 1.

For polynomial approximation on [-1,1], rho is related to the largest Bernstein ellipse in which the analytic continuation remains regular.

Analyticity alone does not guarantee rapid convergence at moderate N. A boundary layer such as

u(x) = exp(-(1+x)/epsilon)

is entire, but its variation occurs on scale epsilon; the resolution required before geometric convergence becomes visible grows as epsilon decreases. Coordinate mappings can also change regularity. A smooth physical solution may appear less regular in reference coordinates if the geometric map has nearly singular derivatives. Regularity must therefore be assessed relative to both the basis and the coordinate carrier.

Modal versus nodal representations

A modal representation stores expansion coefficients:

u_N(x) = sum_{k=0}^N u_hat_k phi_k(x).

This form exposes coefficient decay, sparsity of differential operators, filtering, and regularity. A nodal representation stores values

u_j = u_N(x_j)

at interpolation nodes. It is convenient for boundary data, nonlinear products, physical-space diagnostics, and coupling to external models. The transform between the two representations is

u_j = sum_k V_jk u_hat_k, where V_jk = phi_k(x_j).

Direct inversion of a general Vandermonde matrix is numerically undesirable. Fourier, Chebyshev, and some related transforms use FFT-based algorithms. Legendre and Jacobi transforms require structured recurrences or fast connection algorithms.

Linear differential operators are often sparse in coefficient space, while nonlinear operations are inexpensive in nodal space. A typical pseudospectral step is

modal coefficients -> nodal values -> pointwise nonlinear operation -> modal projection.

This cycle introduces aliasing if the nodal grid cannot represent the product. Modal and nodal descriptions are not competing discretizations; they are complementary coordinate systems on the same finite-dimensional state, and every transform has an associated cost and rounding error.

Reproducible numerical experiments in Julia, Python or MATLAB

A reproducible computation specifies the PDE, domain, coefficients, boundary conditions, basis normalization, node set, quadrature rule, time integrator, nonlinear solver, tolerances, precision, and software environment. The approximation should be tested in more than one norm. A small pointwise error does not certify a small residual, and a small residual may be misleading if the operator inverse is ill-conditioned.

For a stationary problem, report

solution error,
differential residual,
boundary residual,
condition estimate,
coefficient tail.

For a time-dependent problem, spatial and temporal convergence should be separated. First choose a time step small enough that spatial refinement dominates; then choose a spatial resolution large enough that temporal refinement dominates. Manufactured solutions provide exact convergence tests, while invariant and energy checks expose structural defects. Randomized experiments require fixed seeds and distributions over runs rather than a selected realization. The implementation should also verify exact discrete properties, such as matrix symmetry, summation by parts, mass conservation, or zero-mode preservation.

Week 2 — Fourier approximation and the FFT

Fourier series and transforms

For a 2pi-periodic function,

u(x) = sum_{k in Z} u_hat_k exp(i k x),

where

u_hat_k = (1/2pi) integral_0^(2pi) u(x) exp(-i k x) dx.

The Fourier basis diagonalizes constant-coefficient differential operators:

d^m/dx^m exp(i k x) = (i k)^m exp(i k x).

Thus differentiation becomes multiplication of each coefficient by (i k)^m. Translation also becomes diagonal:

u(x-a) has coefficient exp(-i k a) u_hat_k.

Parseval’s identity gives

||u||_L2^2 = 2pi sum_k |u_hat_k|^2.

Sobolev regularity is represented by weighted coefficients:

||u||_Hs^2 ~ sum_k (1+k^2)^s |u_hat_k|^2.

The Fourier transform on the real line replaces the integer mode k by a continuous frequency xi. It converts differentiation into multiplication by i xi and convolution into multiplication. Periodic spectral methods use the discrete mode structure, while unbounded-domain analysis often uses the continuous transform or mapped bases.

Truncation, interpolation and projection

The Fourier projection onto modes |k| <= K is

P_K u = sum_{|k|<=K} u_hat_k exp(i k x).

It is the best L2 approximation in the retained trigonometric space. Fourier interpolation instead constructs a trigonometric polynomial matching u at equispaced nodes. Its coefficients are discrete Fourier coefficients rather than exact integrals. Projection and interpolation coincide only when the unresolved modes contribute no aliases to the retained frequencies.

The truncation error is

u-P_K u = sum_{|k|>K} u_hat_k exp(i k x).

Hence

||u-P_K u||_L2^2 = 2pi sum_{|k|>K} |u_hat_k|^2.

Interpolation error contains both the discarded tail and its aliases. For sufficiently smooth periodic functions the two have comparable asymptotic accuracy, but nonlinear calculations can make their distinction operationally important.

Parseval identities

For periodic functions,

(u,v)_L2 = 2pi sum_k u_hat_k conjugate(v_hat_k).

Consequently,

||u_x||_L2^2 = 2pi sum_k k^2 |u_hat_k|^2.

This converts differential energy estimates into coefficient estimates. For the heat equation,

u_t = nu u_xx,

one obtains

(1/2) d/dt ||u||_L2^2 = -nu ||u_x||_L2^2.

In coefficient form,

d/dt |u_hat_k|^2 = -2 nu k^2 |u_hat_k|^2.

Parseval also permits nonlinear diagnostics. The total resolved energy is the sum of modal energies, while the energy in a band measures the distribution across scales. In discrete computations, the normalization of the DFT must be stated because factors of N or 2pi can otherwise corrupt physical energy calculations.

Trigonometric interpolation

For N equispaced nodes

x_j = 2pi j/N,

the interpolant has the form

I_N u(x) = sum_{k in K_N} u_tilde_k exp(i k x),

where K_N is the retained frequency set. The coefficients are obtained from the DFT. The interpolation nodes are optimal for periodic translation-invariant problems because they permit exact discrete orthogonality:

sum_j exp(i(k-l)x_j) = N delta_(k,l mod N).

The derivative of the interpolant is computed by multiplying each mode by i k and applying the inverse transform. This is usually preferable to explicitly forming the dense Fourier differentiation matrix. Point evaluation away from the grid can be performed through the trigonometric barycentric formula or a nonuniform FFT when many off-grid evaluations are required.

Discrete Fourier transform and FFT

The DFT is

u_tilde_k = sum_{j=0}^{N-1} u_j exp(-2pi i j k/N),

with an inverse normalization determined by convention. Direct evaluation costs O(N^2). The FFT factors N into smaller radices and evaluates the same transform in O(N log N) operations. The FFT is exact as an algebraic transform up to floating-point rounding; it is not an approximation separate from the DFT.

Practical issues include the ordering of negative frequencies, treatment of the Nyquist mode for even N, normalization, and real-valued symmetry:

u_hat_(-k) = conjugate(u_hat_k).

Real FFT routines exploit this symmetry. Multidimensional transforms are formed by successive one-dimensional FFTs along each coordinate. The arithmetic complexity is nearly optimal, but memory layout and data movement often dominate on modern architectures.

Spectral differentiation and integration

For a Fourier interpolant,

u_x = sum_k i k u_hat_k exp(i k x),
u_xx = sum_k -k^2 u_hat_k exp(i k x).

Higher derivatives amplify high-frequency coefficients by powers of k, so roundoff and unresolved noise are magnified. Spectral integration divides by i k for nonzero modes:

v_hat_k = u_hat_k/(i k).

The zero mode is an integration constant and must be fixed separately. For Poisson’s equation,

-u_xx = f,

the solution is

u_hat_k = f_hat_k/k^2, for k != 0.

The zero mode gives the compatibility condition f_hat_0=0 and the nullspace of additive constants.

Dense differentiation matrices are useful for analysis or small problems but should normally not be used for large periodic computations. FFT differentiation costs O(N log N) and preserves the exact modal structure.

Aliasing, convolution and the two-thirds/three-halves de-aliasing rules

If

u(x) = sum_p u_hat_p exp(i p x)
and
v(x) = sum_q v_hat_q exp(i q x),

then

(uv)_hat_k = sum_(p+q=k) u_hat_p v_hat_q.

This convolution doubles the maximum degree. On an N-point grid, frequencies differing by multiples of N are indistinguishable:

exp(i(k+mN)x_j) = exp(i k x_j).

Unresolved high-frequency products therefore fold into low modes. This is aliasing.

For quadratic nonlinearities, the three-halves rule pads a retained N-mode spectrum to approximately 3N/2 points, transforms to physical space, forms the product, transforms back, and truncates to the original modes. The equivalent two-thirds rule retains only the lowest two-thirds of the grid-supported modes. Exact constants depend on the indexing convention. Cubic and higher nonlinearities require more padding or staged products.

De-aliasing is not only a matter of pointwise accuracy. Aliased terms can violate discrete conservation and inject energy into incorrect modes. Conservative, advective, and skew-symmetric formulations may react differently to the same aliasing.

Gibbs phenomena and filtering

If u has a jump discontinuity, Fourier coefficients decay only as 1/|k|. The truncated series develops an overshoot near the jump that approaches a nonzero fraction of the jump magnitude. Increasing N narrows the oscillatory region but does not eliminate the peak overshoot. Away from the discontinuity, convergence remains rapid.

A modal filter replaces

u_hat_k -> sigma(|k|/K) u_hat_k.

An exponential filter may use

sigma(eta) = exp(-alpha eta^p),

where p controls the filter order. Low modes are nearly unchanged while high modes are damped. Spectral viscosity instead adds a high-mode dissipative operator. Filtering stabilizes under-resolved calculations but modifies the PDE, lowers effective resolution, and can alter phase speeds. The filter must therefore be treated as part of the numerical model and documented explicitly.

Week 3 — Fourier methods for periodic PDEs

Fourier-Galerkin and Fourier collocation methods

Fourier-Galerkin methods project the residual onto the retained Fourier modes. Fourier collocation methods enforce the differential equation at equispaced nodes. For constant-coefficient linear equations they are algebraically equivalent after transformation. For variable coefficients or nonlinear terms, collocation forms interpolated products and introduces aliasing, whereas exact Galerkin projection requires convolution over all contributing modes.

Consider

u_t = L u + N(u).

In coefficient space, a constant-coefficient linear operator is diagonal:

d u_hat_k/dt = lambda_k u_hat_k + N_hat_k.

This exact modal diagonalization is the principal efficiency advantage of Fourier methods. The nonlinear term is commonly formed in physical space and de-aliased. Stability should be checked in the discrete norm induced by the trapezoidal rule, which is exact for appropriate products of retained trigonometric polynomials.

Poisson, heat, wave and Helmholtz equations

For periodic Poisson,

-u_xx = f,

the modal equations are

k^2 u_hat_k = f_hat_k.

The zero mode imposes

f_hat_0 = 0.

The mean of u is free and is fixed by a gauge such as u_hat_0=0.

For the heat equation,

u_t = nu u_xx,

each mode satisfies

u_hat_k(t) = exp(-nu k^2 t) u_hat_k(0).

This yields exact spatial evolution but also eigenvalues of magnitude O(N^2), which constrain explicit time steps.

For the wave equation,

u_tt = c^2 u_xx,

each mode is an oscillator:

u_hat_k'' + c^2 k^2 u_hat_k = 0.

Long-time accuracy depends strongly on phase error.

For Helmholtz,

-u_xx - kappa^2 u = f,

the modal denominator is k^2-kappa^2. Near resonance the solution is sensitive to perturbations, and exact resonance requires compatibility with the corresponding eigenspace.

Burgers, Korteweg-de Vries, Kuramoto-Sivashinsky and Allen-Cahn equations

Viscous Burgers is

u_t + u u_x = nu u_xx.

The nonlinear term can be written as u u_x, (u^2/2)_x, or a split combination. These are analytically equivalent but not discretely equivalent under aliasing. The split form can improve energy behavior.

Korteweg-de Vries is

u_t + 6u u_x + u_xxx = 0.

The third derivative is diagonal in Fourier space and purely imaginary, so it creates dispersive stiffness rather than dissipative stiffness. Accurate long-time computation requires control of phase and invariant drift.

Kuramoto-Sivashinsky is

u_t + u u_x + u_xx + u_xxxx = 0.

Low modes are linearly unstable, while high modes are strongly damped by u_xxxx. IMEX and exponential schemes are natural because the fourth derivative creates severe explicit restrictions.

Allen-Cahn is

u_t = epsilon^2 u_xx - F'(u).

It is a gradient flow with energy

E(u) = integral [epsilon^2 |u_x|^2/2 + F(u)] dx.

A stable discretization should reproduce monotone energy decay or a controlled discrete analogue.

Conservation laws and invariant drift

Suppose a PDE preserves I(u). A semidiscrete method preserves it only if the discrete operator reproduces the continuous cancellation. For periodic Burgers,

integral u^2 u_x dx = (1/3) integral (u^3)_x dx = 0.

A collocation product may violate this because the cubic expression is not represented exactly on the grid. De-aliasing or a split formulation can recover the cancellation.

Mass is often tied to the zero mode and can be preserved exactly. Momentum, energy, enstrophy, or Hamiltonians may require a specific spatial form and time integrator. Invariant drift is a structural diagnostic: a solution can look visually accurate while accumulating incorrect phase, energy, or amplitude over long times. Reported computations should therefore include invariant histories, not only final snapshots.

Semi-discrete stability

After spatial discretization,

M_N du_N/dt = F_N(u_N).

For linear systems,

du_N/dt = A_N u_N,

stability depends on the spectrum, numerical range, and non-normal transient growth of A_N. A normal matrix is largely characterized by its eigenvalues; a non-normal matrix is not. For nonlinear systems, derive a discrete energy estimate:

dE_N/dt <= C E_N + G_N.

The constant should remain controlled as N grows. If the estimate acquires large powers of N absent from the continuous problem, the formulation may be using an inappropriate norm or boundary closure.

For constant-coefficient periodic problems, Fourier differentiation preserves skew-adjointness and self-adjointness exactly in the discrete Fourier inner product. Variable coefficients and nonlinear products require additional product-rule or split-form analysis.

Explicit, implicit, IMEX and exponential time integrators

Explicit Runge-Kutta methods are limited by the spectral radius of the spatial operator. Diffusion produces

dt = O(N^-2),

while fourth-order dissipation produces

dt = O(N^-4).

Fully implicit methods remove this linear restriction but require linear or nonlinear solves. An IMEX method treats stiff linear terms implicitly and nonlinear terms explicitly:

(I-dt L)u^(n+1) = u^n + dt N(u^n).

Exponential integrators use the variation-of-constants formula:

u(t+dt) = exp(dt L)u(t) + integral_0^dt exp((dt-s)L)N(u(t+s)) ds.

The integral is approximated using phi-functions. Integrating-factor methods transform

v(t)=exp(-tL)u(t),

while ETD Runge-Kutta methods retain high temporal order without resolving the linear stiffness explicitly. The method should be chosen according to stiffness, non-normality, conservation requirements, and the cost of solving the implicit or exponential action.

Part II — Orthogonal polynomials and fast transforms

Week 4 — General orthogonal-polynomial systems

Weighted inner products

Let w(x)>0 on I=(a,b) with w in L1(I). Define

(u,v)_w = integral_a^b u(x)v(x)w(x) dx.

A polynomial sequence {p_n} is orthogonal if

(p_n,p_m)_w = gamma_n delta_nm.

The weight determines both approximation geometry and endpoint sensitivity. Jacobi weights

w_(alpha,beta)(x) = (1-x)^alpha (1+x)^beta

are integrable for alpha,beta>-1. Positive parameters suppress endpoints; negative parameters emphasize them.

Normalization affects formulas but not orthogonality. Monic polynomials have leading coefficient one. Orthonormal polynomials satisfy gamma_n=1. Classical Jacobi, Legendre, and Chebyshev normalizations simplify recurrence and endpoint expressions. Differential operators may map one weighted space into another, so every derivative identity must specify the source and target weights.

Three-term recurrence and Jacobi matrices

Every orthogonal-polynomial family for a positive measure satisfies

x p_n = a_(n+1) p_(n+1) + b_n p_n + a_n p_(n-1),

with a_n>0 in orthonormal normalization. Multiplication by x is therefore represented by the symmetric tridiagonal Jacobi matrix

J = tridiag(a_n,b_n,a_(n+1)).

The eigenvalues of the truncated Jacobi matrix are Gaussian quadrature nodes. The first components of normalized eigenvectors determine the quadrature weights. This is the Golub-Welsch construction.

The recurrence evaluates the basis in O(N) operations at a point. Forward recurrence is usually stable inside the orthogonality interval but can become unstable outside it. Clenshaw’s backward recurrence evaluates a complete expansion without explicitly forming all basis values and is generally more stable.

Zeros and interlacing

A degree-n orthogonal polynomial for a positive measure has exactly n simple zeros in the interior of the support. The zeros of p_n and p_(n+1) interlace. This follows from orthogonality: insufficient sign changes would permit multiplication by a lower-degree polynomial that contradicts orthogonality.

Zeros can be computed as eigenvalues of the Jacobi matrix, by Newton iteration using derivative recurrences, or through asymptotic initial guesses followed by refinement. Interlacing supplies safe bracketing intervals. As n grows, the zeros cluster according to the equilibrium measure of the interval; for [-1,1], the limiting density has endpoint clustering proportional to (1-x^2)^(-1/2).

Gaussian, Radau and Lobatto quadrature

Gaussian quadrature with N+1 nodes is exact for every polynomial of degree at most 2N+1. Its nodes are the zeros of p_(N+1). Gauss-Radau quadrature fixes one endpoint and is exact through degree 2N; Gauss-Lobatto fixes both endpoints and is exact through degree 2N-1.

Positive quadrature weights provide a discrete positive norm:

||v||_N^2 = sum_j w_j |v(x_j)|^2.

For polynomials within the exactness range, this norm equals or is uniformly equivalent to the continuous weighted norm. Variable coefficients, curved Jacobians, and nonlinear products can exceed the exactness degree, producing quadrature error and aliasing. Over-integration or exact coefficient projection is then required.

Christoffel-Darboux kernels

The polynomial projection kernel is

K_N(x,y) = sum_(j=0)^N p_j(x)p_j(y)/gamma_j.

It reproduces every q in P_N:

q(x) = integral_I K_N(x,t) q(t) w(t) dt.

The Christoffel-Darboux identity collapses this quadratic modal sum to two adjacent basis functions:

K_N(x,y) = c_N [p_(N+1)(x)p_N(y)-p_N(x)p_(N+1)(y)]/(x-y).

On the diagonal,

K_N(x,x) = c_N [p_(N+1)'(x)p_N(x)-p_N'(x)p_(N+1)(x)].

This reduction is fundamental for quadrature weights, interpolation, endpoint asymptotics, and fast evaluation.

Reproducing and projection kernels

The operator

P_N u(x) = integral_I K_N(x,y)u(y)w(y)dy

is the orthogonal projector onto P_N. The diagonal K_N(x,x) measures the local concentration of the polynomial space. Its reciprocal,

lambda_N(x) = 1/K_N(x,x),

is the Christoffel function. It gives the minimum weighted norm of a polynomial constrained to equal one at x and controls point evaluation:

|q(x)|^2 <= K_N(x,x)||q||_w^2.

Near endpoints, K_N(x,x) often grows rapidly, reflecting the high cost of pointwise control. This explains endpoint sensitivity and motivates sampling densities proportional to K_N(x,x)w(x) in stable randomized least-squares approximation.

Modified weights and Christoffel transformations

Multiplying a measure by a linear factor, such as (x-a)w(x), generates a new orthogonal family related to the original family by a Christoffel transformation. The endpoint-restricted kernel itself becomes an orthogonal polynomial for the modified weight. For example, if a is an endpoint,

K_N(x,a)

is proportional to a degree-N orthogonal polynomial associated with (x-a)w(x).

The inverse modification, division by a linear factor plus a mass correction, is a Geronimus transformation. These transformations explain parameter shifts in Jacobi families, modifications induced by boundary conditions, and rank-one updates of Jacobi matrices. They also provide structured connections between basis systems rather than requiring dense coefficient conversion.

Stable polynomial evaluation

Evaluating

p(x) = sum_(k=0)^N c_k phi_k(x)

by forming monomials is unstable for large N. Clenshaw’s algorithm uses the three-term recurrence backward and evaluates the expansion in O(N) operations with controlled error. Near interpolation nodes, barycentric interpolation is preferred:

p(x) = [sum_j lambda_j f_j/(x-x_j)] / [sum_j lambda_j/(x-x_j)].

At a node, the exact stored value is returned directly to avoid cancellation.

Derivative evaluation should use differentiated recurrences, coefficient differentiation, or barycentric differentiation formulas. Repeated application of dense differentiation matrices can be ill-conditioned; coefficient-space operators are often preferable.

Week 5 — Jacobi, Legendre and Chebyshev systems

Jacobi Sturm-Liouville operator

Jacobi polynomials J_n^(alpha,beta) are orthogonal under

w_(alpha,beta)(x) = (1-x)^alpha(1+x)^beta.

They satisfy

L_(alpha,beta) J_n^(alpha,beta) = lambda_n J_n^(alpha,beta),

where

L_(alpha,beta)u = -w_(alpha,beta)^(-1) d/dx [w_(alpha+1,beta+1)u']

and

lambda_n = n(n+alpha+beta+1).

Expanded,

L u = (x^2-1)u'' + [alpha-beta+(alpha+beta+2)x]u'.

The operator is singular at the endpoints because the leading coefficient vanishes there. The weight and flux factor are precisely matched so that the endpoint boundary terms vanish on the natural operator domain.

Self-adjointness and weighted integration by parts

Let

p(x)=w_(alpha+1,beta+1)(x).

Then

L u = -w^(-1)(p u')'.

Weighted integration by parts gives

(L u,v)_w = integral_-1^1 p u'v' dx - [p u'v]_-1^1.

On the natural Jacobi domain, the boundary term vanishes, so

(L u,v)_w = (u,L v)_w.

This self-adjointness yields real eigenvalues, orthogonality of distinct modes, and positive semidefinite energy:

(L u,u)_w = ||u'||_(w_(alpha+1,beta+1))^2.

The endpoint term must be checked for the chosen parameters and functions. Formal integration by parts is insufficient when weights are singular or when trial functions fall outside the natural domain.

Parameter-shift and derivative identities

Jacobi differentiation lowers degree and raises both parameters:

d/dx J_n^(alpha,beta) = (n+alpha+beta+1)/2 J_(n-1)^(alpha+1,beta+1).

More generally,

d^m/dx^m J_n^(alpha,beta) = c_(n,m) J_(n-m)^(alpha+m,beta+m),

where

c_(n,m) = Gamma(n+m+alpha+beta+1) / [2^m Gamma(n+alpha+beta+1)].

Differentiation is therefore sparse if the target coefficients are represented in a shifted Jacobi space. Multiplication by 1-x, 1+x, or 1-x^2 generates reverse parameter-shift identities. These structured maps are the basis of ultraspherical and generalized Jacobi spectral methods.

Legendre and Chebyshev specializations

Legendre polynomials correspond to alpha=beta=0 and weight w=1. They satisfy

-d/dx[(1-x^2)L_n'] = n(n+1)L_n.

Their natural unweighted L2 orthogonality is convenient for Galerkin formulations.

Chebyshev polynomials of the first kind satisfy

T_n(cos theta)=cos(n theta)

and are orthogonal under

w(x)=(1-x^2)^(-1/2).

Chebyshev coefficients are cosine-series coefficients, enabling FFT-based transforms. Chebyshev nodes cluster strongly near the endpoints and support stable barycentric interpolation. Legendre methods often have cleaner variational structure; Chebyshev methods often have faster transforms and simpler nodal implementations.

Connection coefficients between polynomial families

Given two bases,

p_n^(A)(x) = sum_(k=0)^n C_(k,n) p_k^(B)(x).

The connection matrix is upper triangular. Direct dense conversion costs O(N^2), but adjacent Jacobi parameter shifts are banded or admit stable recurrences. Large parameter changes can be composed from elementary shifts or accelerated by divide-and-conquer and asymptotic transforms.

Connection coefficients permit the operator to be represented in the basis where it is sparse while the solution is stored in another basis. They are also essential for multidomain coupling, differentiation, boundary-adapted bases, and comparison of modal tails across polynomial families.

Endpoint behaviour and singular Jacobi weights

The endpoint values satisfy

J_n^(alpha,beta)(1) ~ n^alpha,
|J_n^(alpha,beta)(-1)| ~ n^beta.

Thus pointwise growth depends on the parameters even when the weighted norm remains controlled. When alpha or beta is negative, the weight is singular but integrable for parameters greater than -1. Such weights are useful when the solution has known endpoint singularities.

A basis should incorporate the expected asymptotic form. If

u(x) ~ (1-x)^mu v(x)

with smooth v, factoring out (1-x)^mu can restore rapid coefficient decay. Ignoring the singularity and using ordinary Legendre or Chebyshev polynomials usually produces algebraic convergence.

Polynomial interpolation and projection estimates

Let P_N u be the weighted orthogonal projection and I_N u a polynomial interpolant at Gauss-type nodes. Projection is best in the weighted L2 norm:

||u-P_Nu||_w = inf_(v in P_N) ||u-v||_w.

Interpolation must additionally control the Lebesgue constant and point-evaluation stability. Chebyshev-Lobatto interpolation has a Lebesgue constant growing only logarithmically:

Lambda_N = O(log N).

Endpoint clustering prevents the exponential instability of equally spaced polynomial interpolation.

For analytic functions,

||u-I_Nu||_infinity <= C rho^(-N).

For finite regularity, the convergence is algebraic. Derivative errors lose powers of N; a small function error does not imply an equally small high-order derivative error.

Week 6 — Fast algorithms for orthogonal polynomials

Fast synthesis and analysis

Synthesis maps coefficients to values:

u_hat -> {u(x_j)}.

Analysis maps values or weighted samples back to coefficients. A dense transform costs O(N^2). Fourier and Chebyshev transforms use FFTs and cost O(N log N). Legendre and Jacobi transforms can be accelerated using asymptotic oscillatory representations, low-rank block compression, divide-and-conquer recurrences, or conversion through Chebyshev space.

The transform must be stable over the full degree range, including endpoint-dominated modes. An asymptotic formula valid only for high degrees is normally combined with direct recurrence in the low-degree block.

Discrete cosine and sine transforms

For Chebyshev points

x_j = cos(j pi/N),

the identity

T_k(x_j)=cos(k j pi/N)

turns Chebyshev synthesis into a discrete cosine transform. Different endpoint and parity conventions produce DCT types I through IV. Sine transforms similarly represent bases that vanish at endpoints.

These transforms diagonalize or simplify many one-dimensional operators and allow O(N log N) coefficient-value conversion. Correct normalization and endpoint half-weights are essential. A single missing factor of two corrupts quadrature, energy, and inverse transforms.

Fast Chebyshev-Legendre transforms

Chebyshev and Legendre coefficients represent the same polynomial in bases with different orthogonality measures. Direct conversion is triangular and dense. Fast algorithms exploit the asymptotic relation between both families, where high-degree polynomials resemble oscillatory cosine modes with slowly varying amplitudes.

A practical transform divides the connection matrix into blocks that are numerically low rank, applies compressed matrix multiplication, and handles low modes directly. The result approaches O(N log^2 N) or similar complexity, depending on the algorithm. Stability must be verified for both forward and inverse conversion because the two bases have different endpoint scaling.

Fast connection-coefficient transforms

For general Jacobi parameters, connection matrices can be represented through recurrence relations, hierarchical low-rank structure, or compositions of elementary parameter shifts. Integer shifts often yield banded operators. Noninteger shifts are denser but still structured.

Fast conversion is important when differentiation moves coefficients from (alpha,beta) to (alpha+1,beta+1) and subsequent terms must be combined in a common space. Without fast conversion, a mathematically sparse operator can become computationally dense.

Clenshaw and barycentric algorithms

Clenshaw’s algorithm evaluates a modal expansion from the recurrence relation. If

phi_(k+1) = (a_k x+b_k)phi_k-c_k phi_(k-1),

the coefficient sum is accumulated backward, avoiding explicit formation of all high-degree basis values.

Barycentric interpolation evaluates a nodal polynomial by

p(x) = [sum_j lambda_j f_j/(x-x_j)]/[sum_j lambda_j/(x-x_j)].

The common factors in the interpolation product cancel analytically, producing excellent numerical stability. Differentiation matrices can be derived from the barycentric weights, but their repeated use remains more ill-conditioned than coefficient-space differentiation.

Structured multiplication and differentiation operators

Multiplication by x is tridiagonal in an orthogonal-polynomial basis. Multiplication by a degree-m polynomial is banded with bandwidth proportional to m. A smooth nonpolynomial coefficient has rapidly decaying modal coefficients, so multiplication can often be approximated by a banded operator after truncating its expansion.

Jacobi differentiation is diagonal between parameter-shifted spaces. Conversion back to a common basis is banded for adjacent shifts. This produces sparse operator factorizations:

differentiate -> convert -> multiply -> combine.

Preserving this factorization avoids dense differentiation matrices and improves conditioning.

Almost-banded matrices

Boundary-value operators in coefficient space are often banded except for a small number of dense boundary rows. Such matrices are called almost banded. A second-order differential equation may produce a banded interior operator plus two dense rows enforcing endpoint conditions.

Adaptive QR or structured elimination can solve an almost-banded system in nearly linear complexity. The dense rows must be retained symbolically rather than spread through the matrix by premature elimination. This is one of the central algorithmic advantages of coefficient-space spectral methods.

Adaptive coefficient-space computation

An adaptive solver does not fix N in advance. It expands the system until the coefficient tail falls below a tolerance and the residual is small. A robust criterion examines a block of trailing coefficients rather than the final coefficient alone:

max_(N-m<=k<=N) |u_hat_k| <= tol * scale.

The solver must distinguish genuine convergence from temporary coefficient cancellation. It should also monitor residual coefficients, boundary defects, and conditioning. When the coefficient tail decays slowly, the method may increase N, change basis parameters, split the domain, or introduce a coordinate map.

Multivariate orthogonal polynomials and Koornwinder constructions

Tensor-product bases are natural on rectangles but inefficient or geometrically inappropriate on simplices and curved domains. Koornwinder-type polynomials construct orthogonal bases on triangles and tetrahedra from nested Jacobi factors. A typical triangular basis has the structure

P_p^(alpha,beta)(xi) (1-eta)^p P_q^(gamma,delta_p)(eta).

The parameter dependence ensures orthogonality under the simplex Jacobian. Differentiation and multiplication retain structured sparsity, though less simply than in one dimension. Stable evaluation requires hierarchical recurrences rather than direct monomial expansion.

Part III — Boundary-value and eigenvalue problems

Week 7 — Second-order boundary-value problems

Dirichlet, Neumann, Robin and mixed conditions

Consider

-epsilon u'' + p(x)u' + q(x)u = f, for x in (-1,1).

General separated conditions are

a_- u(-1)+b_- u'(-1)=c_-,
a_+ u(1)+b_+ u'(1)=c_+.

Dirichlet conditions prescribe values, Neumann conditions prescribe outward derivatives, and Robin conditions prescribe a linear combination. In a strong spectral formulation all conditions are imposed directly. In a weak formulation, Dirichlet conditions are usually essential and built into the trial space, while Neumann or Robin conditions arise naturally from integration by parts.

Inhomogeneous conditions can be homogenized by writing

u = u_0 + u_b,

where u_b satisfies the boundary conditions and u_0 satisfies homogeneous conditions. This prevents boundary data from contaminating every modal equation.

Weak formulations and coercivity

For the symmetric model

-u''+alpha u=f,

multiply by v and integrate:

integral u'v' dx + alpha integral uv dx - [u'v] = integral fv dx.

After incorporating natural boundary terms, define

a(u,v)=integral u'v' + alpha uv dx + boundary Robin terms.

Coercivity means

a(v,v) >= c ||v||_H1^2

on the admissible space, possibly after imposing a gauge for a Neumann nullspace. For nonsymmetric advection-diffusion, coercivity may hold after integrating the advection term:

integral p u' u = [p u^2/2] - integral p' u^2/2.

If coercivity fails, an inf-sup estimate or stabilization is required.

Legendre- and Chebyshev-Galerkin methods

A Legendre-Galerkin approximation expands the solution in a basis satisfying the boundary conditions. Orthogonality under the unweighted inner product produces sparse mass and stiffness matrices. A Chebyshev-Galerkin method uses the Chebyshev weight, permitting FFT-based transforms but modifying the weighted weak form.

The basis should be chosen so that differentiation and boundary enforcement are sparse. Generic polynomial bases lead to dense matrices even when the operator is simple. The formulation and basis are therefore designed together.

Boundary-adapted modal bases

For a second-order problem, a useful basis combines three neighboring polynomials:

phi_k = P_k + a_k P_(k+1)+b_k P_(k+2).

The coefficients are selected so that

a_- phi_k(-1)+b_- phi_k'(-1)=0,
a_+ phi_k(1)+b_+ phi_k'(1)=0.

For homogeneous Dirichlet Legendre problems,

phi_k = L_k-L_(k+2).

These functions satisfy the endpoint constraints exactly. Their derivatives have simple orthogonality properties, often making the stiffness matrix diagonal and the mass matrix pentadiagonal. Boundary adaptation converts a globally constrained polynomial space into a sparse modal system.

Tau and collocation formulations

In a tau method, the residual coefficients vanish in all but the highest modes. Those highest equations are replaced by boundary constraints. For a second-order equation, two tau equations are normally replaced. The placement and scaling of these rows affect conditioning and possible spurious modes.

Collocation chooses nodes {x_j} and imposes

L u_N(x_j)=f(x_j)

at interior nodes plus boundary conditions at the endpoints. The method is easy to formulate with differentiation matrices but produces dense systems. Stability follows from discrete norm equivalence and interpolation estimates, not simply from pointwise satisfaction.

Equivalence and non-equivalence of formulations

For constant-coefficient linear problems and compatible quadrature, Galerkin, tau, and collocation formulations may be related by invertible changes of basis. This equivalence can fail for variable coefficients, nonlinearities, approximate quadrature, or different boundary treatments.

Two methods may return the same polynomial in exact arithmetic but have very different conditioning. Conversely, visually similar matrices may correspond to different weak problems. Equivalence should be established at the operator level:

A_2 = S A_1 T

with controlled transformations S and T, rather than assumed from matching dimensions.

Sparse stiffness and mass matrices

With basis functions {phi_k}, define

S_jk = integral phi_k' phi_j' dx,
M_jk = integral phi_k phi_j dx.

Boundary-adapted Legendre bases can make S diagonal and M banded. The system

(S+alpha M)u=f

can then be solved in O(N) operations for fixed bandwidth.

Variable coefficients broaden the bandwidth according to their polynomial degree. Smooth coefficients can be truncated spectrally, producing controllably almost-banded matrices. Matrix sparsity is a consequence of basis and operator compatibility, not merely of polynomial orthogonality.

Error estimates

For a coercive Galerkin scheme,

||u-u_N||_H1 <= C inf_(v_N in V_N) ||u-v_N||_H1.

Duality can improve the L2 estimate by one power when the adjoint problem has additional regularity:

||u-u_N||_L2 <= C N^-1 ||u-u_N||_H1

in representative smooth cases.

Collocation estimates include interpolation and quadrature defects. Boundary errors must be included explicitly if conditions are not built into the space. For analytic solutions the error decays geometrically until conditioning or floating-point error dominates. For endpoint singularities the rate is algebraic unless the basis or domain decomposition resolves the singular structure.

Week 8 — Modern sparse spectral solvers

Ultraspherical and coefficient-space methods

Chebyshev coefficients are convenient for representing functions, but repeated derivatives are dense when all terms are forced back into the Chebyshev basis. Ultraspherical methods retain derivatives in shifted Gegenbauer spaces. The derivative map is sparse:

d^m/dx^m : C^(lambda) -> C^(lambda+m).

Each lower-order term is converted to the highest parameter space using sparse conversion operators. A differential operator becomes a banded combination of differentiation, multiplication, and conversion matrices, plus dense boundary rows.

This avoids the severe conditioning of dense collocation differentiation matrices and supports adaptive infinite-dimensional QR.

Banded differentiation and conversion operators

In ultraspherical coefficients,

d/dx T_n = n C_(n-1)^(1),

and subsequent derivatives remain sparse in higher-parameter spaces. Conversion from parameter lambda to lambda+1 is banded because adjacent ultraspherical families satisfy short recurrence relations.

For

L u = a_m(x)u^(m)+...+a_0(x)u,

each term is mapped into the common lambda=m space. If the coefficient functions have short or truncated expansions, the complete operator remains banded or almost banded.

Integral reformulations

Differential equations can be preconditioned by integrating rather than differentiating. For

u''=f,

write

u = I^2 f + c_0+c_1 x.

Boundary conditions determine c_0,c_1. Integral operators smooth coefficient sequences and are compact or nearly compact, producing systems of the form

(I+K)u=g.

Second-kind equations have condition numbers that remain bounded or grow mildly with resolution. Integral reformulation is therefore both an analytical and numerical preconditioner.

Spectral integration

Coefficient integration uses sparse recurrence relations. In Chebyshev form,

integral T_n dx

is a combination of T_(n+1) and T_(n-1) for n>=2. Repeated spectral integration is stable because it divides high-frequency coefficients rather than multiplying them.

Integration constants carry the nullspace of differentiation. They must be retained as explicit unknowns and fixed by boundary or compatibility conditions. Discarding them is a common source of rank deficiency and incorrect boundary closure.

Conditioning of differentiation matrices

A first-order polynomial collocation differentiation matrix has norm growing roughly as O(N^2) because of endpoint clustering. A second derivative can have norm or condition growth of order O(N^4) or worse, depending on scaling and boundary enforcement. This growth reflects both the unbounded differential operator and the coordinate representation.

The condition number of a matrix is not itself the whole stability result. One should distinguish coefficient conditioning, nodal conditioning, and backward error. Nevertheless, dense high-order differentiation matrices become increasingly difficult to solve accurately without preconditioning.

Diagonal and operator preconditioning

Diagonal preconditioners rescale rows or columns to equalize modal magnitudes. Operator preconditioners approximate the inverse of the highest derivative or principal part. In coefficient space, integration provides a natural right preconditioner:

L D^(-m) = I + compact/lower-order terms.

For variable-coefficient elliptic problems, preconditioners may use the constant or separable principal part. A valid preconditioner should produce resolution-independent clustering of eigenvalues or singular values while remaining inexpensive to apply.

Adaptive truncation

An adaptive solver extends the coefficient vector until the tail is resolved. A typical stopping rule compares the envelope of the final coefficients with the norm of the full coefficient vector:

max tail <= tol * max(1,||u_hat||).

The residual should also be computed in a space appropriate to the differential operator. A small coefficient tail can be misleading if the basis poorly represents a localized feature or if the operator amplifies the tail strongly.

Adaptive truncation can be combined with automatic domain splitting, basis-parameter changes, or precision increases. The resolution becomes an output of the solver rather than an input.

Residual and coefficient-tail error indicators

The coefficient tail estimates unresolved approximation content. The residual estimates violation of the equation. They should be used together. For a boundedly invertible operator,

||u-u_N||_X <= C ||f-Lu_N||_Y.

The constant C may be estimated analytically, through a discrete resolvent bound, or by a posteriori validation. Boundary residuals require separate measurement.

A plateauing coefficient tail may indicate roundoff, a singularity, insufficient precision, or a basis mismatch. A plateauing residual with a decaying tail may indicate solver tolerance or operator assembly error.

Matrix-free Krylov solution

Large multidimensional spectral systems are often applied as operators rather than assembled. A matrix-free Krylov method requires products

v -> A v

computed through transforms, derivatives, coefficient multiplications, and geometric mappings. GMRES is appropriate for general nonsymmetric systems; conjugate gradients require a symmetric positive-definite operator in the chosen inner product.

The preconditioner determines practical convergence. Matrix-free execution reduces storage but does not automatically reduce iteration count. Residual norms reported by Krylov algorithms should be related to the physical residual, especially when left or right preconditioning changes the computed norm.

Week 9 — Higher-order and constrained problems

Fourth- and higher-order equations

A fourth-order equation requires four boundary conditions and greater regularity in a direct weak formulation. For the biharmonic model

u''''=f,

a conforming Galerkin formulation based on two integrations by parts requires trial functions in H^2. Spectral bases can satisfy this globally, but their endpoint constraints must include both value and derivative conditions when clamped boundaries are imposed.

Higher derivatives increase stiffness and condition growth. Coefficient-space formulations remain sparse if derivatives are represented through parameter shifts or integration preconditioning.

Generalized Jacobi bases

Generalized Jacobi functions incorporate endpoint factors:

phi_n(x) = (1-x)^r (1+x)^s J_n^(alpha,beta)(x).

The exponents enforce repeated homogeneous boundary conditions. For a 2m-order problem, choosing factors with sufficient endpoint multiplicity creates a basis in the correct constrained space.

These bases retain derivative identities and weighted orthogonality. They are especially effective when the operator naturally maps one generalized Jacobi family into another.

Dual Petrov-Galerkin methods

Odd-order equations and non-self-adjoint operators often benefit from distinct trial and test bases. A dual Petrov-Galerkin design chooses spaces so that

(L phi_k,psi_j)

is diagonal or narrowly banded. The trial basis satisfies one set of boundary conditions, while the dual test basis satisfies the adjoint conditions.

Stability requires an inf-sup estimate, not merely a sparse matrix. The duality should reflect the continuous adjoint operator and boundary form. Properly designed, it converts odd-order differential operators into well-conditioned sparse systems.

Biharmonic and Cahn-Hilliard equations

The biharmonic equation is

Delta^2 u=f.

In one dimension, clamped conditions are u=u'=0; simply supported conditions differ. Mixed formulations introduce v=-Delta u, reducing the problem to two second-order equations.

Cahn-Hilliard is

u_t = Delta mu,
mu = F'(u)-epsilon^2 Delta u.

It conserves mass and dissipates free energy:

dE/dt = -||grad mu||^2.

A spectral method should preserve the zero mode and use a time integrator consistent with energy decay. Direct fourth-order formulations are compact, while mixed formulations offer simpler boundary handling and preconditioning.

Mixed formulations

A mixed formulation introduces auxiliary variables to lower differential order. For biharmonic,

v=-Delta u,
-Delta v=f.

For incompressible flow, pressure is a Lagrange multiplier enforcing div u=0. Mixed methods require compatibility between spaces, typically expressed by an inf-sup condition.

Spectral spaces can satisfy compatibility exactly through divergence-free bases or carefully paired velocity-pressure polynomial spaces. An arbitrary pair may produce spurious pressure modes or rank deficiency.

Differential-algebraic boundary constraints

Boundary conditions, conservation laws, or gauge constraints can produce a differential-algebraic system:

M u_dot = F(u)+C^T lambda,
C u = d.

The multiplier lambda enforces the constraint. The index of the DAE determines how many times the constraint must be differentiated to expose an evolution equation. Time integrators must preserve the constraint or apply projection.

In spectral discretizations, algebraic constraints may involve dense endpoint evaluation rows. Their scaling is important because poor scaling can dominate the condition number.

Exact enforcement of multiple boundary conditions

A differential equation of order m requires m independent boundary conditions for a regular two-point problem. They can be imposed by basis construction, tau rows, bordering, nullspace projection, or Lagrange multipliers.

Basis construction gives exact satisfaction for every coefficient vector but may complicate variable coefficients. Tau and bordering are more flexible but introduce dense rows. Nullspace projection constructs a basis for the kernel of the boundary matrix. Exact independence must be verified: redundant or incompatible conditions produce singular systems.

Nullspaces, compatibility and gauge conditions

Neumann Poisson has a constant nullspace:

-u''=f,
u'(-1)=u'(1)=0.

Solvability requires

integral f dx = 0.

The solution is unique only after fixing a gauge such as

integral u dx=0.

Similar nullspaces occur in pressure, electromagnetic potentials, periodic antiderivatives, and rigid-body modes. The discrete nullspace should match the continuous nullspace. Artificial removal of a mode can destroy conservation, while failure to fix the gauge leaves the linear system singular.

Week 10 — Spectral eigenvalue problems

Sturm-Liouville eigenproblems

A regular Sturm-Liouville problem is

-(p u')'+q u=lambda w u

with self-adjoint boundary conditions. The weak form is

a(u,v)=lambda m(u,v),

where

a(u,v)=integral [p u'v'+q uv],
m(u,v)=integral w uv.

A conforming Galerkin discretization yields

A c = lambda M c.

For positive p,w, A is symmetric and M positive definite. Eigenvalues obey variational principles, and the low spectrum is approximated particularly accurately when the eigenspaces are smooth.

Generalized matrix eigenvalue formulations

The matrix pencil

A x=lambda B x

should be treated according to the properties of A and B. If B is positive definite, a Cholesky reduction converts the problem to standard symmetric form. If B is singular, infinite eigenvalues or algebraic constraints may be present.

Boundary enforcement can introduce rows that do not belong to the mass operator. Tau formulations may therefore produce singular pencils. The physical finite eigenvalues must be separated from constraint-associated eigenvalues.

Spurious eigenvalues and spectral pollution

Spectral pollution occurs when discrete eigenvalues converge to points outside the true spectrum or when nonphysical modes persist under refinement. Causes include nonconforming spaces, inconsistent boundary conditions, singular matrix pencils, truncated unbounded domains, and non-self-adjoint sensitivity.

Validation requires residuals, convergence of eigenvectors, comparison across formulations, and variational bounds where available. Agreement of eigenvalues at two nearby resolutions is insufficient if both discretizations share the same defect.

Non-self-adjoint operators

For A != A*, eigenvectors need not be orthogonal and may be severely ill-conditioned. A small perturbation can cause a large eigenvalue displacement. The eigenvalue condition number is related to left and right eigenvectors:

kappa(lambda)=||x||||y||/|y* x|.

When y*x is small, the eigenvalue is sensitive. Spectral discretizations of advection, open systems, delay equations, and absorbing-boundary problems commonly produce non-normal matrices. Their behavior cannot be inferred from eigenvalues alone.

Pseudospectra and resolvent growth

The epsilon-pseudospectrum is

Lambda_epsilon(A) = {z : ||(zI-A)^(-1)|| > 1/epsilon}

together with the spectrum. Equivalently, it contains eigenvalues of A+E for perturbations ||E||<epsilon.

Large resolvent norm away from eigenvalues indicates transient amplification and spectral sensitivity. For a semigroup,

exp(tA),

pseudospectral geometry can predict short-time growth even when every eigenvalue lies in the stable half-plane. Numerical computation uses smallest singular values of zI-A, contour algorithms, or resolvent actions.

Polynomial and rational filtering

To isolate a spectral region, apply a polynomial p(A) that is large on desired eigenvalues and small elsewhere. Chebyshev polynomials provide near-minimax filters on intervals. Rational filters approximate spectral projectors using shifted inverses:

r(A)=sum_j omega_j (A-z_j I)^(-1).

Filtering accelerates Krylov eigensolvers and separates interior eigenvalues. Its effectiveness depends on the spectral geometry and conditioning of the shifted solves.

Contour and shift-invert methods

Shift-invert transforms

A x=lambda x

into

(A-sigma I)^(-1)x = (lambda-sigma)^(-1)x.

Eigenvalues near sigma become dominant. Contour methods use the spectral projector

P = (1/2pi i) integral_Gamma (zI-A)^(-1) dz.

Quadrature replaces the contour integral by a sum of shifted linear solves. These methods parallelize naturally over quadrature nodes and can compute all eigenvalues in a region.

Validation by residual and backward error

For an approximate eigenpair (lambda_hat,x_hat), the residual is

r=A x_hat-lambda_hat B x_hat.

A normwise backward error asks for the smallest perturbations Delta A,Delta B making the pair exact. A representative measure is

eta = ||r|| / [(||A||+|lambda_hat|||B||)||x_hat||].

For non-normal problems, a small residual does not imply a small forward eigenvalue error unless the eigenvalue condition number is also controlled. Report residual, conditioning, and resolution dependence together.

Part IV — Time evolution, singularity and nonlocality

Week 11 — Time-dependent PDEs

Method of lines

Spatial discretization produces an ODE or DAE:

M du/dt = F(u,t).

The mass matrix may be diagonal, sparse, dense, or singular. The method of lines separates spatial and temporal analysis but does not eliminate their interaction. Spectral spatial discretization often creates a wide range of eigenvalues, making the ODE stiff.

Before selecting a time method, inspect the eigenvalues, numerical range, nonlinear Lipschitz scale, invariants, and mass-matrix structure. Explicit treatment may be adequate for transport but inefficient for diffusion or high-order dissipation.

CFL restrictions

For explicit methods, stability requires dt lambda_j to lie inside the method’s stability region. Fourier advection gives

dt=O(N^-1),

while polynomial collocation can yield a more restrictive endpoint-driven scale near O(N^-2) for first derivatives. Diffusion commonly gives O(N^-2) in Fourier space and potentially stronger restrictions in unpreconditioned polynomial nodal form. Fourth derivatives produce O(N^-4) scales.

The relevant quantity is the actual discrete spectral radius or numerical range, not a memorized power law. Mappings and element sizes introduce factors of h^-m.

Implicit-explicit schemes

Split

u_t=L u+N(u),

where L is stiff and inexpensive to invert, while N is nonstiff or nonlinear. A first-order IMEX step is

(I-dt L)u^(n+1)=u^n+dt N(u^n).

Higher-order IMEX Runge-Kutta and multistep schemes require coupled explicit and implicit order conditions. Stability can degrade when the split operators do not commute or when the explicit part is itself stiff. The implicit solve should exploit spectral diagonalization, banded coefficient matrices, or separable preconditioners.

Exponential and integrating-factor methods

The exact linear evolution is incorporated through

u(t+dt)=exp(dt L)u(t)+integral_0^dt exp((dt-s)L)N(u(t+s))ds.

ETD methods approximate the integral using

phi_1(z)=(exp(z)-1)/z,

and higher phi-functions. In Fourier space, diagonal L makes these functions inexpensive. For general matrices, Krylov, contour, or rational approximations compute their action.

Integrating-factor methods define

v=exp(-tL)u,

so

v_t=exp(-tL)N(exp(tL)v).

This removes linear stiffness but can create rapidly varying transformed nonlinearities for strongly non-normal or dissipative operators.

Operator splitting

If

u_t=(A+B)u,

Lie splitting is

exp(dt A)exp(dt B),

with first-order accuracy. Strang splitting is

exp(dt A/2)exp(dt B)exp(dt A/2),

with second-order accuracy. The leading error depends on commutators such as [A,B].

Splitting is attractive when each subproblem is exactly solvable or cheaply implicit. Boundary conditions can cause order reduction because the domains of A and B may not share compatible traces. Nonlinear splitting requires careful interpretation of each flow.

Symplectic and energy-preserving integration

Hamiltonian systems satisfy

z_t=J grad H(z),

with skew matrix J. Symplectic integrators preserve the symplectic form and normally exhibit bounded long-time energy error rather than exact energy preservation. Discrete-gradient and average-vector-field methods can preserve an energy exactly.

For Hamiltonian PDEs, the spatial discretization must first preserve a finite-dimensional Poisson or symplectic structure. A symplectic time integrator cannot repair a spatial discretization that has already destroyed it.

Long-time error and invariant preservation

Classical local truncation error does not fully predict long-time behavior. Phase error, modified Hamiltonians, invariant drift, and weakly damped modes can dominate. Backward-error analysis interprets a geometric integrator as the exact flow of a nearby modified system over long intervals.

Long-time validation should report phase, invariants, modal energy, and statistical quantities where appropriate. A small short-time norm error may conceal secular drift.

Adaptive temporal and spectral resolution

Temporal adaptation estimates local time error and changes dt. Spectral adaptation monitors coefficient tails, residuals, or localized indicators and changes N. The two controls interact: increasing N may introduce stiffer eigenvalues and force a smaller explicit time step.

A robust controller prevents rapid oscillation of N or dt, transfers states conservatively between resolutions, and separates transient coefficient growth from persistent under-resolution. In multidomain methods, local p-adaptation may be preferable to global resolution increases.

Week 12 — Integral, fractional and nonlocal equations

Volterra and Fredholm integral equations

A Fredholm equation has the form

u(x)-lambda integral_a^b K(x,t)u(t)dt=f(x).

A Volterra equation has an upper limit depending on x:

u(x)-integral_a^x K(x,t)u(t)dt=f(x).

Fredholm operators are global and may have nontrivial spectra. Volterra operators are triangular in physical causality and are often uniquely solvable by iteration. Smooth kernels admit rapid spectral approximation; singular or weakly singular kernels require adapted quadrature and bases.

Spectral Nystrom methods

A Nystrom method approximates the integral directly by quadrature:

integral K(x,t)u(t)dt ~ sum_j w_j K(x,t_j)u(t_j).

Collocating at the same nodes gives a matrix system for the nodal values. For smooth kernels and Gaussian quadrature, convergence can be spectral. The matrix is generally dense, so fast multipole, hierarchical, low-rank, or FFT structure may be needed at large scale.

Quadrature error must be uniform in the target variable x. Near singular diagonals, ordinary Gaussian rules lose accuracy.

Weakly singular kernels

A kernel such as

K(x,t)=|x-t|^-mu K_0(x,t), with 0<mu<1,

is integrable but not smooth along the diagonal. Polynomial approximation of the unmodified integrand converges slowly. Techniques include singularity subtraction, product integration, graded meshes, coordinate transformations, and fractional Jacobi bases.

The singular part should be integrated analytically or with a weight-adapted quadrature. Treating it as an ordinary smooth kernel may produce apparently convergent values with incorrect asymptotic rates.

Jacobi methods for endpoint singularities

If a Volterra or fractional problem produces

u(x)=(1+x)^mu v(x)

with smooth v, choose a weighted Jacobi basis that includes (1+x)^mu. The remaining coefficient sequence can then decay rapidly.

Fractional integration maps weighted Jacobi functions to parameter-shifted weighted Jacobi functions, producing sparse or diagonal operational matrices. The basis is selected from the operator’s singularity structure rather than from generic polynomial convenience.

Fractional derivatives and fractional Sturm-Liouville systems

Fractional derivatives are nonlocal. Riemann-Liouville and Caputo derivatives differ in how endpoint data enter. For 0<alpha<1,

D_RL^alpha u = d/dx I^(1-alpha)u,

while the Caputo derivative applies the fractional integral after an ordinary derivative.

Fractional Sturm-Liouville systems combine left and right fractional derivatives with weighted inner products. Their eigenfunctions are often generalized Jacobi functions rather than classical polynomials. Self-adjointness requires a fractional integration-by-parts identity and appropriate endpoint spaces.

Nonlocal diffusion

A representative nonlocal diffusion operator is

L_delta u(x)=integral_[|y-x|<delta] gamma_delta(x,y)[u(y)-u(x)]dy.

As delta->0, suitable kernels recover local diffusion. Spectral approximation can diagonalize translation-invariant kernels on periodic domains. On bounded domains, the treatment of interactions outside the domain determines the nonlocal boundary condition.

Nonlocal operators often produce dense matrices, but convolution, separability, low-rank compression, or fast summation can reduce cost.

Delay differential equations

A delay equation has

u'(t)=F(u(t),u(t-tau)).

The state at time t is the history segment on [t-tau,t], so the natural phase space is infinite-dimensional. Spectral methods approximate the history by polynomials or collocation values.

State-dependent delays are more difficult because the evaluation point depends on the solution. Stability analysis leads to transcendental characteristic equations, and spectral collocation can approximate their roots. History interpolation must preserve the requested temporal order.

Singular quadrature and kernel compression

Singular quadrature isolates a known singular factor:

integral s(x,t)g(t)dt.

Product-integration rules interpolate g while integrating s exactly against the interpolation basis. Logarithmic, algebraic, and Cauchy singularities each require different formulas.

Kernel compression approximates a smooth off-diagonal block by low rank:

K(x,t) ~ sum_(r=1)^R a_r(x)b_r(t).

Hierarchical matrices, adaptive cross approximation, butterfly methods, and fast multipole schemes exploit this structure. Near-diagonal singular blocks remain uncompressed or receive specialized treatment.

Week 13 — Unbounded and semi-infinite domains

Hermite polynomials and functions

Hermite polynomials are orthogonal on the real line under exp(-x^2). Hermite functions incorporate the Gaussian factor and form an orthonormal basis of unweighted L2(R). They diagonalize the quantum harmonic oscillator and are natural for Gaussian-localized solutions.

Differentiation and multiplication by x are tridiagonal through raising and lowering relations. The global Gaussian envelope can be inappropriate for algebraic tails or translated structures, so scaling and centering are essential.

Laguerre polynomials and functions

Laguerre polynomials are orthogonal on [0,infinity) under x^alpha exp(-x). Laguerre functions absorb part of the weight and are used for semi-infinite domains. Their recurrence and derivative identities produce sparse operators.

They are effective for exponentially decaying solutions. Slowly decaying or oscillatory solutions may require scaling, rational bases, domain decomposition, or complex contour methods.

Rational Chebyshev bases

Map an unbounded interval to a finite interval, for example

x = L(1+y)/(1-y), with y in [-1,1).

Then expand the mapped function in Chebyshev polynomials of y. The resulting rational Chebyshev functions cluster points near the finite endpoint and distribute them algebraically toward infinity.

The scale L controls where resolution is concentrated. The map introduces variable metric factors into derivatives and may create endpoint singularities in reference space. Its suitability depends on the physical decay and length scale.

Mapped spectral methods

A map x=g(xi) converts

d/dx = [g'(xi)]^-1 d/dxi.

Higher derivatives contain derivatives of g. A good map regularizes the solution and distributes points according to its variation. A poor map introduces large metric factors and worsens conditioning.

Maps can resolve boundary layers, semi-infinite domains, coordinate singularities, and localized structures. The transformed PDE, boundary conditions, quadrature Jacobian, and energy norm must all be derived consistently.

Scaling and translation parameters

Hermite and Laguerre bases can be scaled:

phi_n((x-x_0)/L).

The center x_0 tracks localization, while L matches the width. An unsuitable scale causes slow coefficient decay even for a smooth function.

Parameters may be chosen from moments, minimization of tail coefficients, residual reduction, or physical length scales. Adaptive scaling changes the basis over time, requiring stable coefficient transport and accounting for the induced basis-motion terms.

Oscillatory and localized solutions

Highly oscillatory solutions require phase-aware bases. Standard polynomial approximation may need O(k) modes to resolve exp(i k x). Alternatives include WKB enrichment, plane-wave bases, mapped coordinates, and oscillatory quadrature.

Localized pulses on the real line may be represented by Hermite functions, rational bases, or moving-domain decompositions. The basis should match both envelope and phase. Resolving only the envelope while misrepresenting the phase produces severe dispersion error.

Transparent and absorbing boundary treatments

When an unbounded problem is truncated, artificial boundaries can reflect outgoing waves. A transparent boundary condition represents the exact exterior Dirichlet-to-Neumann map, usually as a nonlocal operator in time or along the boundary.

Approximate absorbing conditions use local asymptotic expansions. Perfectly matched layers introduce a complex coordinate stretch that damps outgoing waves without reflection in the continuous model. Discrete reflection depends on profile smoothness, layer thickness, and spectral resolution.

Domain truncation versus native unbounded-domain bases

Truncation permits standard finite-domain methods but introduces boundary-model error. Native Hermite, Laguerre, or rational bases avoid a finite artificial boundary but may poorly match the tail or produce dense variable-coefficient operators.

The decision should compare three errors:

tail truncation,
basis approximation,
conditioning.

No method is universally superior. Exponential decay favors Hermite or Laguerre systems; algebraic decay often favors rational mappings; localized wave propagation may favor truncation with a PML.

Part V — Multiple dimensions and complex geometry

Week 14 — Tensor-product and multidimensional methods

Rectangles, cuboids and periodic boxes

On Cartesian product domains, multidimensional bases are formed from one-dimensional bases:

phi_(i,j)(x,y)=phi_i(x)psi_j(y).

Fourier bases are natural for periodic directions; Chebyshev or Legendre bases are natural for bounded directions. Mixed domains can use Fourier in one coordinate and Jacobi in another.

Constant-coefficient separable operators inherit tensor structure. Boundary conditions can be imposed independently by direction when the geometry and operator permit separation.

Tensor-product bases

If

u_N(x,y)=sum_i sum_j u_ij phi_i(x)psi_j(y),

then differentiation in x acts only on the first index and differentiation in y only on the second. The coefficient array should be treated as a tensor rather than flattened prematurely.

Tensor-product approximation gives exponential convergence for analytic functions but suffers from dimensional growth:

degrees of freedom ~ (N+1)^d.

Sparse grids, low-rank tensors, and anisotropic resolution attempt to reduce this cost.

Kronecker structure

A separable two-dimensional operator produces matrices such as

A = A_x kron M_y + M_x kron A_y.

Matrix-vector products can be performed by reshaping the vector into a matrix:

A_x U M_y^T + M_x U A_y^T.

This avoids assembling the full Kronecker matrix. Eigenvalues and preconditioners can often be constructed from one-dimensional factors.

Fast diagonalization

For

A_x U M_y^T + M_x U A_y^T = F,

solve generalized one-dimensional eigensystems

A_x X=M_x X Lambda_x,
A_y Y=M_y Y Lambda_y.

Transforming with X and Y reduces the problem to scalar divisions by

lambda_i^x+lambda_j^y.

Fast diagonalization gives a direct solver or powerful preconditioner for separable elliptic operators. Variable coefficients and curved mappings break exact separability but may be treated as perturbations.

Alternating-direction solvers

ADI methods split multidimensional operators into directional parts. For a stationary Sylvester equation, successive shifted one-dimensional solves reduce the error. For time evolution, alternating implicit directions avoid a full multidimensional implicit solve.

The convergence depends on shift selection and spectral intervals. Noncommuting variable-coefficient operators introduce splitting error. ADI is particularly effective when each directional solve is banded or diagonal in a spectral basis.

Disks, cylinders, balls and spherical shells

Curvilinear domains require bases adapted to regularity at coordinate singularities. On a disk, Fourier modes in angle combine with radial Jacobi or Zernike polynomials. A mode exp(i m theta) must have radial behavior proportional to r^|m| near the origin.

Cylinders combine angular Fourier, radial orthogonal, and axial polynomial or Fourier bases. Balls and shells use spherical harmonics with radial Jacobi-type functions. Shells avoid the origin and are easier than full balls with respect to regularity.

Spherical harmonics

Spherical harmonics satisfy

-Delta_S2 Y_l^m = l(l+1)Y_l^m.

They form an orthonormal basis on the sphere. A scalar field is expanded as

u(theta,phi)=sum_l sum_(m=-l)^l u_lm Y_l^m.

The spherical Laplacian is diagonal. Vector and tensor fields require vector spherical harmonics or spin-weighted harmonics to preserve geometric meaning. Fast spherical harmonic transforms combine FFTs in longitude with associated Legendre transforms in latitude.

Coordinate singularities

Polar and spherical coordinates contain apparent singular factors such as 1/r and 1/sin theta. A smooth Cartesian field satisfies compatibility relations that cancel these singularities. A naive tensor polynomial basis may violate them and generate spurious modes.

Regularity should be enforced through mode-dependent radial factors, parity constraints, or bases derived from Cartesian polynomials. Coordinate singularities are representation failures, not physical singularities.

Sparse grids and low-rank tensor representations

Sparse grids select tensor basis indices satisfying a constraint such as

i_1+...+i_d <= N.

For mixed-regularity functions, this reduces degrees of freedom relative to a full tensor product. The convergence rate depends on anisotropic derivative bounds.

Low-rank representations approximate

U(i_1,...,i_d) ~ sum_(r=1)^R a_r^(1)(i_1)...a_r^(d)(i_d).

Tensor trains and hierarchical Tucker formats provide more stable high-dimensional variants. Differential operators with tensor structure can be applied without reconstructing the full array, but nonlinear operations may increase rank and require recompression.

Week 15 — Spectral and hp-element methods

Reference-to-physical element maps

Each physical element K is obtained from a reference element K_hat by

x=F_K(xi).

Gradients transform as

grad_x u = J_K^-T grad_xi u,

and integrals include |det J_K|. Curved maps make the metric coefficients variable even when the physical PDE has constant coefficients.

The map must remain invertible and shape regular. Poor Jacobian conditioning amplifies derivative errors and can destroy approximation quality.

High-order nodal and modal elements

Nodal elements store values at interpolation points and simplify fluxes, nonlinearities, and boundary coupling. Modal elements store coefficients in hierarchical bases and simplify p-adaptivity, filtering, and static condensation.

The two forms are related by local transforms. A well-designed code may evaluate nonlinear terms nodally while solving or preconditioning in a modal basis. Local polynomial degree can vary by element.

Gauss-Lobatto quadrature

Gauss-Lobatto points include element endpoints, making them convenient for conforming assembly and numerical fluxes. With N+1 points, the rule is exact through degree 2N-1.

Using the same nodes for interpolation and quadrature produces a diagonal mass matrix in collocated formulations, but the quadrature is not exact for all products of degree-N polynomials when curved geometry or nonlinear coefficients are present. This under-integration creates aliasing.

Conforming and discontinuous formulations

A conforming spectral element method enforces continuity across elements and approximates a global H1 space. A discontinuous Galerkin method permits independent element traces and couples elements through numerical fluxes.

DG methods handle nonconforming meshes, local conservation, and hyperbolic fluxes naturally. Their stability depends on penalty parameters or upwind fluxes. Conforming methods have fewer interface unknowns but are less flexible under local refinement and complex coupling.

Static condensation

Element unknowns are divided into interior and boundary components:

[A_ii A_ib; A_bi A_bb][u_i;u_b]=[f_i;f_b].

Eliminate

u_i=A_ii^-1(f_i-A_ib u_b).

The global system then contains only interface unknowns. Interior solutions are recovered element by element. At high order, this substantially reduces the global solve and enables efficient local factorizations.

Sum factorization

Naive evaluation of a d-dimensional tensor polynomial costs approximately O(N^(2d)). Sum factorization applies one-dimensional transforms successively, reducing cost toward O(d N^(d+1)).

It is essential for high-order finite and spectral elements. Modern performance is often limited by memory traffic rather than floating-point count, so data layout and batched tensor contractions matter.

Curved elements and geometric aliasing

On curved elements, metric terms and Jacobians are nonpolynomial or higher-degree polynomial functions. Their product with the solution can exceed quadrature exactness. Aliased geometric terms can violate free-stream preservation, symmetry, and conservation.

Remedies include over-integration, polynomial projection of metric terms, compatible curl-form metric identities, and split formulations. Geometry is part of the discrete operator and must be approximated with the same structural care as the solution.

h-, p- and hp-adaptivity

h-adaptivity subdivides elements. p-adaptivity increases polynomial degree. hp-adaptivity combines both. Smooth regions favor increasing p; singular or nonsmooth regions favor decreasing h.

Modal coefficient decay can distinguish smooth from nonsmooth behavior. Exponential decay suggests p-refinement; algebraic decay or coefficient stagnation suggests h-refinement. Reliable adaptivity also requires error localization and conformity management.

Mortar and interface coupling

Mortar methods couple nonmatching traces through an intermediate interface space. Continuity or flux balance is imposed weakly by projection. The mortar space must be chosen to satisfy stability and avoid loss of conservation.

Interfaces may join different polynomial degrees, element types, or physical models. Projection operators should preserve constants and relevant flux moments. Poorly designed interpolation can inject or remove energy.

Parallel and GPU-oriented implementation

High-order methods have high arithmetic intensity because each degree of freedom participates in many local tensor operations. Matrix-free sum-factorized kernels are therefore well suited to GPUs and vector processors.

Efficient execution requires contiguous data, batched element kernels, minimized global synchronization, and fused operator stages. Global communication arises from interface exchange and Krylov reductions. Performance should be measured by time to a specified error, not merely throughput per degree of freedom.

Week 16 — Stability on curved and nonlinear systems

Metric identities and discrete geometric conservation laws

A transformed conservation law contains contravariant fluxes and metric factors. Continuous metric identities ensure that a constant state remains constant under the mapped divergence. The discrete metric terms must satisfy analogous identities.

Failure produces free-stream errors even when the physical flux is constant. Metric factors constructed through compatible discrete curls or projections can enforce the identities. Moving meshes additionally require a discrete geometric conservation law relating Jacobian evolution to grid velocity.

Summation-by-parts structure

A discrete derivative matrix D with mass matrix H satisfies

H D + D^T H = B,

where B represents endpoint evaluation. This is the discrete integration-by-parts identity. It permits energy estimates that mirror the continuum derivation.

Gauss-Lobatto collocation with compatible quadrature provides an SBP structure. Boundary and interface terms are then controlled by simultaneous approximation terms or numerical fluxes. SBP is an algebraic conservation carrier, not merely a differentiation formula.

Split forms and entropy stability

Nonlinear conservation laws satisfy an entropy identity involving an entropy variable v=partial U/partial u. Direct nodal products generally violate the continuous product rule. Split forms rewrite nonlinear terms symmetrically to recover a discrete energy or entropy balance.

Flux differencing uses a two-point entropy-conservative flux and adds dissipation to obtain entropy stability. The resulting semidiscrete inequality is

d/dt sum_j H_j U(u_j) <= boundary contribution.

The proof depends on symmetry, consistency, and an exact discrete telescoping relation.

De-aliasing on mapped elements

Nonlinear fluxes, material coefficients, and geometric factors all increase polynomial degree. Over-integration evaluates the weak form with a quadrature capable of representing more of the product. Polynomial dealiasing projects each product into the approximation space before differentiation.

Mapped elements require geometric and physical dealiasing together. Over-integrating only the physical flux does not fix aliased metric products. The required quadrature order depends on the polynomial degree of the solution, geometry, and nonlinearity.

Shock detection and spectral viscosity

Smooth spectral approximations oscillate near shocks. A shock detector measures modal decay, local variation, entropy residual, or interelement jumps. Stabilization is activated only where the smoothness indicator signals under-resolution.

Spectral viscosity damps high modes while preserving low-mode accuracy. In element methods, artificial viscosity may vary in space and polynomial degree. The viscosity profile must remain positive, sufficiently smooth, and scaled to the local wave speed and element size.

Positivity and realizability

Physical states may require density, pressure, concentration, or distribution functions to remain nonnegative. High-order polynomials can violate positivity between nodes even when nodal values are positive.

Limiters rescale the high-order correction toward a positive cell average:

u_new = u_bar + theta(u_high-u_bar), with 0<=theta<=1.

For systems, positivity must be checked in the physical admissible set. Moment methods additionally require realizability: the moments must correspond to a nonnegative underlying distribution.

Incompressible Navier-Stokes projection methods

The equations are

u_t+u·grad u=-grad p+nu Delta u,
div u=0.

Projection methods first compute an intermediate velocity and then solve a pressure Poisson equation to remove its divergent component. Boundary conditions for the pressure correction determine splitting accuracy and may create boundary layers.

Spectral methods can use divergence-free bases, pressure-velocity formulations, or Fourier projection in periodic directions. Nonlinear de-aliasing and kinetic-energy-preserving forms are central to under-resolved stability.

Transition, turbulence and under-resolved computation

In transitional and turbulent flows, the smallest dynamically relevant scales may not be resolved. Spectral accuracy on the retained modes does not imply physical accuracy if energy transfer to unresolved modes is misrepresented.

Under-resolved methods require controlled dissipation, explicit subgrid models, implicit LES behavior, or entropy-stable fluxes. Validation should examine spectra, transfer functions, dissipation rates, coherent structures, and statistical convergence. A visually plausible turbulent field is not a sufficient benchmark.

Part VI — Modern extensions and verification

Week 17 — Scientific machine learning and spectral operators

Solution operators between function spaces

A PDE solution operator maps input functions to output functions:

G : a(.) -> u(.).

The input may be a coefficient, forcing term, initial condition, geometry, or boundary data. Unlike finite-dimensional regression, operator learning must remain meaningful under changes of discretization.

A spectral representation approximates both input and output in coefficient spaces. The learned map should specify the function-space norms, admissible input class, and target resolution. Without this, apparent generalization may be only interpolation on a fixed grid.

Fourier and spectral neural operators

A Fourier neural operator layer has the schematic form

v_(l+1)(x)=sigma(W v_l(x)+F^-1[R_l(k)F(v_l)(k)]).

Only a finite set of modes is transformed by learned matrices R_l(k). The global convolution is efficient through FFTs. Spectral neural operators replace Fourier modes with polynomial, spherical-harmonic, wavelet, or problem-adapted modes.

The learned spectral multiplier is generally nonlinear through repeated activation and mixing. It should not be confused with the exact diagonalization of a linear PDE operator. Boundary conditions and geometry require explicit treatment.

Learning in coefficient space

Instead of training on nodal values, represent

a(x)=sum a_k phi_k(x),
u(x)=sum u_k phi_k(x),

and learn the map {a_k}->{u_k}. This separates approximation error from regression error and permits mode-dependent weighting. Low modes may dominate physical energy, while high modes encode fine-scale corrections.

Coefficient normalization is essential because modal magnitudes can span many orders. The model should preserve known symmetries, parity, reality constraints, and conservation modes where possible.

Spectral residual losses

Given a predicted u_theta, define the PDE residual

R_theta=L u_theta-f.

A spectral residual loss is

L_res=sum_k w_k |R_hat_k|^2.

The weights determine the residual norm. Unweighted coefficients correspond to an L2-type norm; derivative-weighted coefficients correspond to Sobolev norms. Boundary residuals and initial conditions require separate terms.

A small residual does not prove a small solution error without a stability estimate. Training can also minimize residuals in modes that are easy for the model while ignoring difficult localized defects. Independent physical-space and boundary validation is required.

Parseval-based training objectives

For an orthonormal basis,

||u-u_ref||_L2^2=sum_k |u_k-u_ref,k|^2.

Parseval therefore permits exact coefficient-space loss evaluation. Sobolev losses use

sum_k (1+lambda_k)^s |error_k|^2

when the basis diagonalizes an appropriate elliptic operator.

Mode weighting should reflect the intended application. Relative weighting by the reference amplitude can overemphasize tiny noisy coefficients. Absolute, relative, energy-based, and task-based losses should be distinguished.

Aliasing and resolution-transfer failure

Neural spectral layers form nonlinear products or activations in physical space, which generate unresolved modes. If these are transformed back without padding, the network contains aliasing just like a pseudospectral PDE solver.

A model trained on one grid may exploit the grid’s aliasing pattern and fail when evaluated at another resolution. Resolution invariance requires de-aliasing, consistent projection, or a continuous operator definition whose discretizations commute approximately with refinement.

Hybrid solver-learner architectures

A learned component can approximate a closure, preconditioner, coarse correction, constitutive law, or initial guess while a classical solver enforces the governing equation. For example,

u_(n+1)=u_n+M_theta(r_n)

uses a learned residual correction inside an iterative method. The outer solver supplies consistency and a stopping criterion.

Hybrid designs are safer when the learned component accelerates rather than replaces certification. Failure can be detected through residual stagnation, invariant violation, or out-of-distribution diagnostics.

Training-data and discretization dependence

Training data inherit the numerical errors of the generating solver. A network may learn dispersion, artificial diffusion, boundary artifacts, or unresolved turbulence rather than the continuum operator.

Datasets should record solver type, resolution, tolerance, and estimated error. Training and testing on data generated by the same discretization can conceal shared bias. Cross-discretization tests reveal whether the learned object approximates the PDE operator or merely a particular numerical map.

Comparison with classical spectral solvers at equal accuracy

A fair comparison fixes the target error and measures total cost, including data generation, training, inference, setup, and repeated solves. A neural operator may be expensive to train but cheap per query; a direct spectral solver may be superior for a small number of queries.

Compare residual, boundary error, conservation, robustness, memory, and extrapolation, not only runtime. Both methods must use comparable precision and hardware. Report where the learned model fails and whether a classical solver can detect or repair the failure.

Stability, conservation and posterior certification

Architectural constraints can preserve mass, symmetry, equivariance, or positivity. Stability may be imposed through contractive layers, energy-bounded multipliers, or dissipative parameterizations. These constraints are useful but do not automatically certify the complete nonlinear operator.

Posterior certification evaluates the PDE residual, boundary residual, invariant defects, and an error estimator after inference. A learned prediction can be accepted, corrected by a solver, or rejected. Certification converts a statistical approximation into a controlled numerical component.

Week 18 — Error certification and computational audit

A priori versus a posteriori error

An a priori estimate predicts error from regularity and resolution before computation:

||u-u_N|| <= C N^(m-s)||u||_Hs.

An a posteriori estimate uses the computed solution and residual:

||u-u_N||_X <= C_stab ||f-Lu_N||_Y + boundary terms.

A priori estimates guide method design; a posteriori estimates guide adaptation and acceptance. Both depend on constants whose size matters. An asymptotic rate without a usable constant may not certify a finite computation.

Residual, truncation and quadrature decomposition

Let the assembled discrete solution satisfy a quadrature-based operator L_N^Q. Then

f-Lu_N = [f-f_N] + [L_N^Q u_N-Lu_N] + [f_N-L_N^Q u_N].

These terms represent data approximation, operator/quadrature error, and algebraic residual. A coefficient tail measures truncation but not necessarily quadrature error. Each component should be evaluated separately when diagnosing a convergence plateau.

Roundoff and conditioning

Floating-point arithmetic satisfies

fl(a op b)=(a op b)(1+delta), with |delta|<=u_machine

in a standard model. Forward error is controlled by backward error multiplied by condition number. Spectral differentiation amplifies high-mode roundoff, while interpolation at stable nodes is usually backward stable.

When exponential convergence reaches machine precision, increasing N can worsen the result. Extended precision, better scaling, orthonormal bases, and preconditioning can move the roundoff floor.

Manufactured solutions

Choose an exact function u_exact, apply the differential operator analytically to define f, and derive compatible boundary data. The code should recover the predicted convergence rate.

Manufactured solutions should test variable coefficients, mixed boundaries, curved geometry, nonlinear terms, and all implementation branches. Choosing a polynomial exactly representable by the basis tests assembly but not convergence. Nonpolynomial analytic and finite-regularity examples are both needed.

Grid and polynomial-order refinement

For h-refinement, estimate the observed rate from

p_obs = log(e_h/e_(h/2))/log 2.

For spectral p-refinement, inspect semilog error versus N for geometric decay and log-log error versus N for algebraic decay. Refinement should continue until the asymptotic regime and the error floor are visible.

Refining only the nodal grid while keeping the polynomial space or quadrature inconsistent may not test the intended method. Every coupled resolution parameter should be stated.

Cross-formulation verification

Solve the same problem with two structurally different formulations, such as Galerkin and collocation, coefficient and nodal methods, or direct and mixed formulations. Agreement within independent error estimates is stronger evidence than refinement of one code path.

Shared libraries, quadrature, or boundary routines can create common-mode errors, so genuinely independent implementations provide better verification.

Conservation and stability audits

For each invariant or energy law, compute the discrete balance defect. If

dE/dt = -D+B,

evaluate

defect = dE_N/dt + D_N-B_N.

For conservative problems, track mass, momentum, energy, and boundary fluxes. For dissipative problems, verify sign and rate. Stability audits should include transient growth, not only final boundedness.

Reproducibility and benchmark design

A benchmark defines equations, parameters, geometry, boundary data, reference values, norms, tolerances, and hardware conventions. It should include enough information to reproduce the calculation without reverse engineering.

Reference solutions require their own validation. A high-resolution run is not automatically exact. Benchmarks should cover smooth, singular, stiff, non-normal, nonlinear, and geometry-sensitive regimes rather than a single favorable problem.

Cost versus accuracy

The relevant performance curve is error versus wall-clock time, memory, or energy use. High-order methods have larger per-degree cost but can reach small errors with far fewer unknowns. Low-order methods may dominate at low accuracy or for nonsmooth solutions.

Setup cost, factorization reuse, transform cost, parallel efficiency, and number of right-hand sides all matter. Complexity stated only as a function of degrees of freedom can conceal high constants and communication costs.

When spectral convergence fails

Failure occurs when the solution is nonsmooth, the basis mismatches endpoint behavior, the geometry map is irregular, the PDE develops shocks, quadrature aliases nonlinear products, boundary conditions are inconsistent, or conditioning reaches the floating-point floor.

The response should identify the mechanism. Remedies include domain decomposition, singularity factoring, coordinate mapping, weighted bases, filtering, viscosity, over-integration, higher precision, or a different formulation. Increasing N without diagnosing the failure can increase cost while reducing accuracy.

Week 19 — Final projects and research presentations

Navier-Stokes or magnetohydrodynamics on periodic or curved domains

A complete project should specify the spatial basis, pressure or divergence constraint, nonlinear form, de-aliasing, time integration, and energy diagnostics. Curved domains require metric identities and consistent geometric terms. Magnetohydrodynamics additionally requires control of div B=0, magnetic energy, and coupling between velocity and field.

High-order Schrodinger eigenvalue computations

The project should formulate the self-adjoint or non-self-adjoint eigenproblem, choose a basis adapted to the potential and domain, and validate eigenpairs by residual and convergence. Singular Coulomb potentials require coordinate adaptation, basis enrichment, or domain decomposition. Interior eigenvalues may require shift-invert or contour methods.

Fractional or nonlocal PDEs

A project should define the precise fractional derivative or nonlocal kernel, boundary interpretation, singularity structure, and fast application strategy. Validation should include known asymptotics, convergence under refinement, and comparison with an independent quadrature or transform method.

Spectral-element wave propagation

The implementation should examine dispersion, numerical fluxes, CFL limits, curved geometry, and absorbing boundaries. Measure phase and amplitude error over many wavelengths rather than only short-time norms. Sum factorization and matrix-free execution should be included in the performance study.

Non-normal operator pseudospectra

Construct a spectral discretization of a non-self-adjoint operator, compute eigenvalues, resolvent norms, and pseudospectral contours, and connect them to transient growth. Validate that observed pseudospectral features persist under resolution and are not discretization artifacts.

Adaptive ultraspherical solvers

Build a coefficient-space boundary-value solver using sparse differentiation and conversion operators, dense boundary rows, and adaptive QR. Test smooth, boundary-layer, and singular problems. Compare conditioning and cost with Chebyshev collocation.

Fast Jacobi connection transforms

Implement direct and fast transforms between selected Jacobi families. Measure complexity, stability, endpoint accuracy, and inverse consistency. Demonstrate a differential operator whose sparse representation depends on rapid parameter conversion.

Spectral methods for kinetic equations

Possible systems include Vlasov-Poisson, Fokker-Planck, or Boltzmann-type models. The project must address phase-space dimensionality, positivity, conservation, velocity truncation, and filamentation. Low-rank or sparse tensor strategies may be needed.

Structure-preserving Cahn-Hilliard or Hamiltonian systems

For Cahn-Hilliard, preserve mass and verify energy decay. For Hamiltonian systems, preserve symplectic or Poisson structure and examine long-time invariant behavior. Compare a structure-preserving time method against a conventional high-order integrator at equal short-time accuracy.

Classical spectral solver versus neural operator under equal error and cost constraints

Train a neural operator and construct a classical spectral solver for the same parameterized PDE. Establish independent reference solutions. Compare total training and inference cost, repeated-solve break-even point, residuals, boundary error, conservation, resolution transfer, and out-of-distribution behavior. Posterior correction by the classical solver should be tested.

Core reading structure

The foundational layer is supplied by classical spectral approximation, orthogonal-polynomial analysis, Fourier methods, Galerkin formulations, and spectral-element discretization. The modern computational layer consists of fast transforms, ultraspherical and almost-banded solvers, tensor-product operator application, geometric conservation, entropy-stable split forms, adaptive coefficient methods, and matrix-free high-order implementation. The verification layer consists of residual estimates, backward error, invariant audits, cross-formulation comparison, and reproducible benchmarks. Scientific machine learning belongs after these layers because its predictions must be represented, differentiated, stabilized, compared, and certified within the same function-space and operator framework.

The complete 2026 progression is therefore

operator class
-> approximation space
-> modal/nodal representation
-> sparse or transform-based operator
-> stable boundary and time treatment
-> multidimensional geometry
-> nonlinear structural preservation
-> adaptive error control
-> learned operator extension
-> posterior certification.

The course is mathematically strong but architecturally incomplete. In xSCD terms, it is an extensive payload collection rather than a fully functioning instructional system.

The global course state can be typed as

C_course = ⟨function spaces, operators, bases, discretizations, algorithms, geometries, solvers, certificates⟩.

Its intended FSB triangle is:

M — what the learner retains: operator structure, approximation structure, and numerical mechanism.

Q — what counts as the same method: not merely the same basis or matrix, but the same trial/test carrier, boundary transport, residual norm, and stability class.

Ψ — what distinguishes methods: approximation rate, stability, conditioning, conservation, computational cost, geometric compatibility, and liftback to the continuum problem.

The course contains all three components, but it does not yet enforce their interaction systematically.

1. Primitive failure

The principal failure is not lack of content. It is that the course is organized primarily as a sequence of topics:

Fourier
→ orthogonal polynomials
→ boundary problems
→ time evolution
→ multidimensional geometry
→ machine learning
→ verification.

That is a taxonomic progression. xSCD requires a functioning-system progression:

carrier
→ failure of current carrier
→ new transport
→ new residue/debt
→ certificate
→ liftback.

The course frequently explains individual mechanisms well, but the transitions between mechanisms are not always causally prosecuted.

For example:

Fourier methods
→ Jacobi methods

is presented because the domain changes from periodic to nonperiodic. The stronger xSCD transition is:

periodic translation-invariant carrier
→ boundary trace cannot be represented without Gibbs debt
→ counterkernel: Fourier periodic identification is invalid
→ weighted Sturm–Liouville carrier
→ boundary-compatible polynomial transport
→ new debts: endpoint growth, quadrature, conditioning.

That distinction matters. Without it, students learn a catalogue of methods. With it, they learn when and why a carrier must be replaced.

2. N-clock evaluation: local refinement within each fixed carrier

The N-clock asks whether each topic is developed deeply enough while its underlying carrier remains fixed.

Part I — Approximation and spectral representations

This is the strongest part of the course. It establishes operator classes, weighted residuals, approximation, stability, and modal/nodal representations.

Its local weakness is overcompression. Week 1 introduces too many authority-level ideas simultaneously:

PDE classification
× weighted residual
× Galerkin
× stability
× convergence
× regularity
× reproducibility.

The concepts are correct, but the local state is not allowed to stabilize before the next distinction is introduced.

The missing local kernel is:

operator
→ trial carrier
→ test carrier
→ residual
→ certificate.

Every method should be instantiated through that five-object spine before comparison.

N-clock status:

MATHEMATICAL CONTENT: CERT

LOCAL INSTRUCTIONAL CLOSURE: FRONTIER_PAYLOAD

Part II — Orthogonal polynomials and fast transforms

This part has a coherent mathematical carrier:

positive measure
→ orthogonal family
→ three-term recurrence
→ Jacobi matrix
→ projection kernel
→ structured operators
→ fast transforms.

The strongest feature is the transition from abstract orthogonality to computational structure. The Christoffel–Darboux kernel, parameter shifts, recurrence relations, and sparse differentiation all belong to one native system.

The main local defect is that three different notions of “kernel” appear without a formal disambiguation protocol:

reproducing kernel

integral-equation kernel

learned operator kernel.

xSCD would require a type declaration each time:

kernel := carrier × domain × codomain × composition law.

Without this, lexical identity can be mistaken for structural identity.

N-clock status:

ORTHOGONAL-POLYNOMIAL CARRIER: CERT

KERNEL TYPE DISCIPLINE: FRONTIER_PAYLOAD

Part III — Boundary-value and eigenvalue problems

This part correctly introduces boundary-adapted bases, tau methods, collocation, coefficient-space solvers, nullspaces, and spectral validation.

Its central unresolved debt is the continuum–matrix transport. The course explains how matrices arise, but the following map is not made first-class:

continuous operator domain
→ finite trial/test quotient
→ matrix or pencil
→ computed vector
→ reconstructed function
→ continuum residual
→ continuum claim.

The matrix is often treated as the operational object, although in xSCD it is only an intermediate carrier.

The most important missing distinction is:

small algebraic residual
≠
small differential residual
≠
small continuum solution error.

N-clock status:

DISCRETIZATION MECHANISMS: CERT

MATRIX-TO-CONTINUUM LIFTBACK: FRONTIER_PAYLOAD

Part IV — Time evolution, singularity, and nonlocality

This part contains three separate carrier changes:

static operator → evolution operator

local differential operator → singular/nonlocal operator

bounded domain → unbounded domain.

Each is substantial enough to require an explicit carrier-mutation event.

The course currently treats them as consecutive subject areas. xSCD would require the failure that forces each transition:

static discretization fails to encode temporal stiffness

local derivative fails to encode memory or long-range interaction

finite-domain basis creates artificial-boundary residue.

The content is correct, but these changes are not marked as L-clock mutations; they appear as N-clock continuation.

N-clock status within each topic:

CERT

Clock typing:

MISCLASSIFIED: several L-clock mutations presented as ordinary topic progression

Part V — Multidimensional geometry and high-order elements

This is the most important globalization layer. It transports one-dimensional spectral structure into multidimensional, mapped, and decomposed domains.

The course correctly introduces:

tensor products

Kronecker structure

curvilinear maps

metric identities

summation by parts

geometric conservation

entropy stability.

The central debt is that geometry first appears as a domain complication and only later becomes an operator component. In xSCD, the geometric map is part of the carrier from the moment differentiation is transported:

D_x = J^-T D_ξ.

Once that map is introduced, every later identity depends on its regularity, invertibility, metric consistency, and quadrature representation.

The missing persistent object is a geometry debt ledger:

Δ_geo = ⟨Jacobian error, metric identity defect, mapping regularity, quadrature mismatch, interface mismatch⟩.

N-clock status:

LOCAL GEOMETRIC METHODS: CERT

GEOMETRIC DEBT TRANSPORT: FRONTIER_PAYLOAD

Part VI — Machine learning and verification

The course correctly places scientific machine learning after the classical numerical machinery. That sequencing is sound.

However, Week 17 and Week 18 are architecturally reversed in authority. Verification should not first appear as a late topic. It should be active from Week 1.

Week 18 contains the actual liftback machinery:

residual
× backward error
× conditioning
× conservation
× cross-formulation verification
× reproducibility.

These are not terminal topics. They are authority conditions for every previous topic.

Machine learning is correctly described as an extension rather than a replacement, but the admissibility test should be sharper:

learned operator admission
requires

discretization independence
× residual control
× stability
× conservation
× resolution transport
× posterior correction or rejection.

N-clock status:

ML CONTENT: CERT AS OPTIONAL PAYLOAD

CERTIFICATION PLACEMENT: ARCHITECTURALLY LATE

3. G-clock evaluation: globalization across the course

The G-clock asks whether locally valid results descend and compose into a global course-level system.

The course has many local closures:

Fourier diagonalization

Jacobi sparse differentiation

Galerkin coercivity

ultraspherical bandedness

tensor-product separability

SBP energy estimates

residual certification.

What is missing is an explicit descent map connecting them.

The global transport should be:

continuous law
→ basis representation
→ discrete operator
→ solver
→ computed state
→ diagnostic
→ liftback.

At present, each chapter often restarts this process implicitly. The learner must reconstruct the common architecture.

The course therefore needs a persistent cross-week invariant table:

Carrier

Identity relation

Observation family

Transport map

Boundary residue

Aliasing residue

Conditioning debt

Geometry debt

Certificate

Liftback.

Without this, local closures remain isolated.

G-clock failure:

LOCAL CERTIFICATES DO NOT AUTOMATICALLY DESCEND TO A COURSE-WIDE CERTIFICATE.

Examples:

orthogonal basis
does not descend automatically to
stable PDE solver.

sparse matrix
does not descend automatically to
well-conditioned algorithm.

small coefficient tail
does not descend automatically to
small continuum error.

energy-stable semidiscretization
does not descend automatically to
energy-stable fully discrete method.

accurate training loss
does not descend automatically to
accurate solution operator.

G-clock status:

GLOBALIZATION: FRONTIER_PAYLOAD

4. L-clock evaluation: carrier mutations

The course contains at least eight genuine carrier mutations:

1. periodic functions → nonperiodic weighted polynomial spaces

2. smooth functions → endpoint-singular weighted spaces

3. dense nodal differentiation → sparse coefficient operators

4. static problems → time-evolution systems

5. local operators → nonlocal/fractional operators

6. bounded domains → unbounded-domain bases or mapped carriers

7. single domains → multidomain spectral elements

8. deterministic solvers → learned solution operators

xSCD requires each mutation to be justified by an exact counterkernel. The course currently provides the rationale in prose, but not as an admission rule.

The required mutation schema is:

old carrier C

CK(C) = exact failure

missing payload

new carrier C'

transport C → C'

new residue generated

liftback test.

Example:

C = global polynomial basis on one interval

CK = localized nonsmoothness destroys global coefficient decay

missing payload = locality

C' = spectral-element decomposition

transport = restriction + reference mapping + interface coupling

new residue = interface flux/mortar debt

liftback = global conservation and convergence.

L-clock status:

MUTATIONS PRESENT

MUTATION CERTIFICATES NOT EXPLICIT

5. Course-level counterkernels

The following counterkernels should be active throughout the course.

CK_BASIS_IS_METHOD

A basis does not determine a numerical method. Trial space, test space, quadrature, boundary treatment, solver, and norm are also required.

CK_SPECTRAL_CONVERGENCE_LAUNDERING

Analytic approximation does not imply stable or accurate computation.

CK_BOUNDARY_ERASURE

Integration by parts exports boundary terms. They cannot disappear unless a trace condition, natural boundary law, or numerical flux discharges them.

CK_QUADRATURE_IDENTITY

Quadrature exactness on a polynomial class does not imply exactness for nonlinear, variable-coefficient, or mapped products.

CK_MATRIX_CONTINUUM_COLLAPSE

A matrix eigenvalue, residual, or condition number is not identical to the corresponding continuum object.

CK_CONDITION_NUMBER_MONISM

A large matrix condition number does not by itself prove instability, and a small condition number does not certify continuum accuracy.

CK_LOCAL_ENERGY_GLOBAL_STABILITY

A stable spatial discretization does not imply stability after time discretization.

CK_GEOMETRY_AS_DATA

The geometry map is an operator component, not passive input.

CK_RESOLUTION_INVARIANCE_ASSUMED

A method working at one N, mesh, or grid does not thereby define a discretization-independent operator.

CK_LEARNED_SOLVER_REPLACEMENT

A learned approximation without residual, stability, and liftback is not a numerical solver certificate.

6. Debt and residue ledger

The course discusses most of the relevant residues, but does not carry them forward as persistent state.

The course-level ledger should be

Δ_course =

Δ_boundary

+ Δ_aliasing

+ Δ_quadrature

+ Δ_conditioning

+ Δ_time

+ Δ_geometry

+ Δ_interface

+ Δ_model

+ Δ_floating_point

+ Δ_liftback.

Each subsection should either:

discharge a debt

transform a debt

introduce a new debt

or

prove that a debt is absent.

Examples:

Fourier padding discharges part of Δ_aliasing.

Boundary-adapted bases discharge Δ_boundary but may introduce basis-conditioning debt.

Mapped elements introduce Δ_geometry.

IMEX methods reduce linear stiffness debt but leave splitting and explicit-nonlinearity debt.

Neural operators introduce Δ_model and Δ_resolution.

Week 18 should not introduce this ledger; it should close the ledger accumulated from Weeks 1–17.

7. Liftback audit

The course’s strongest conceptual claim is that numerical structure must be traced back to the continuum problem. This is stated repeatedly but not operationalized uniformly.

Every method should end with the same liftback test:

computed coefficients or nodal values

→ reconstructed function u_N

→ boundary trace check

→ continuum residual R_N

→ stability/inverse estimate

→ error claim

→ cost and reproducibility record.

The course currently permits weaker closure:

matrix solved
→ plot produced
→ method accepted.

That must be prohibited.

The universal certificate should be:

CERT_method

iff

carrier typed

∧ transport defined

∧ boundary residue discharged

∧ stability shown

∧ approximation error bounded

∧ algebraic solve verified

∧ continuum residual measured

∧ liftback justified.

8. Missing payload

The exact missing payload is not more subject matter. It is a uniform operating protocol embedded in every subsection.

For each topic, add six mandatory paragraphs:

1. Native carrier

What objects are represented, and in which function space?

2. Primitive failure

Why is the previous method or carrier insufficient?

3. Transport

How is the continuous object moved into the finite representation?

4. Debt and residue

What is lost, approximated, aliased, truncated, or exported to a boundary?

5. Certificate

Which theorem or numerical test licenses the result?

6. Liftback

How is the computed object returned to a continuum statement?

This overlay does not alter the TOC. It changes the course from a compendium into a functioning discovery and certification system.

9. Terminal adjudication

Local content coverage:

CERT

Mathematical breadth:

CERT

Classical-to-modern sequencing:

CERT

Local method explanations:

CERT WITH UNEVEN DEPTH

Course-wide carrier continuity:

FRONTIER_PAYLOAD

Debt propagation:

FRONTIER_PAYLOAD

Mutation certification:

FRONTIER_PAYLOAD

Continuum liftback:

FRONTIER_PAYLOAD

Overall terminal:

FRONTIER_PAYLOAD

The course is not deficient in mathematics. It is deficient in explicit system transport.

10. Next executable prosecution step

Do not add another topic.

Rewrite Week 1, subsection by subsection, using the invariant xSCD template:

carrier
→ primitive failure
→ transport
→ debt
→ counterkernel
→ certificate
→ liftback
→ next admissible mutation.

Then propagate that template unchanged through all 19 weeks.

The first executable unit is:

Week 1.1 Model PDE classes

For each of the four PDE classes, produce:

native state carrier

energy identity

boundary residue

spectral transport

discrete energy identity

stability certificate

failure counterkernel

time-discretization debt

continuum liftback.

That is the point at which the course becomes xSCD-complete rather than merely comprehensive.