Spectral Methods and High-Order Numerical Methods 2026 xCSD

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

 

Spectral Methods and High-Order Numerical Methods — 2026

xSCD-integrated expert course text

The course treats every numerical method as a typed functioning system rather than as a formula collection. For each method, the native carrier consists of the function space, operator domain, boundary law, geometry, basis, quadrature, and discrete state representation. Transport maps connect continuum functions to coefficients or nodal values and return computed states to continuum claims. Every approximation introduces an explicit debt ledger containing truncation, aliasing, quadrature, boundary, conditioning, temporal, geometric, interface, modeling, and roundoff residues. A result is accepted only when the relevant stability estimate converts measured residuals into a justified error statement. Refinement within a fixed carrier is an N-clock operation; assembly and multidimensional descent are G-clock operations; a change of basis, geometry, formulation, or model forced by a counterexample is an L-clock mutation.

Part I — Approximation and spectral representations

Week 1 — Foundations of high-order discretization

Model PDE classes: elliptic, parabolic, hyperbolic, dispersive

The classification of a partial differential equation determines which quantities propagate, which quantities are constrained globally, and which boundary data are admissible. For an elliptic equation such as -div(a grad u)+cu=f, uniform ellipticity requires a(x)xi·xi >= a0|xi|², and the natural certificate is a coercive weak form. The solution is determined globally by the forcing and boundary conditions; local marching is not meaningful. For a parabolic equation such as u_t-div(a grad u)=f, the state evolves while spatial diffusion smooths high frequencies. The energy relation is 1/2 d||u||²/dt + integral a grad u·grad u = (f,u)+boundary flux. Spatial spectral accuracy can therefore coexist with severe temporal stiffness because the largest resolved diffusion eigenvalue grows rapidly with resolution.

A hyperbolic equation such as u_t+a·grad u=0 transports information along characteristics. Its energy changes through inflow, outflow, and coefficient divergence, so an omitted inflow condition is not a small technical defect but an incomplete operator carrier. A dispersive equation such as u_t+beta u_xxx=0 preserves L2 energy on a periodic domain because the third derivative is skew-adjoint. Its principal numerical risk is phase error rather than artificial growth. The four classes should therefore be distinguished through their energy transport, boundary flux, resolvent behavior, and time-scale structure, not only through derivative order. In xSCD terms, each class has a different admissible identity relation and different failure observers: elliptic coercivity loss, parabolic stiffness, hyperbolic flux inconsistency, and dispersive phase drift.

Weighted-residual framework

Let L u=f with boundary operator B u=g. Choose a finite trial space V_N=span{phi_0,...,phi_N} and write u_N=sum_j c_j phi_j. The residual is R_N=L u_N-f. A weighted-residual method chooses test functions psi_i and imposes <R_N,psi_i>=0. This does not mean that R_N=0; it means only that the residual component visible to the chosen test space vanishes. If P_W denotes projection onto the test space, the discrete equations enforce P_W R_N=0, while (I-P_W)R_N remains unresolved residue.

The numerical method is not determined by the polynomial basis alone. It is determined by the trial carrier, test carrier, residual pairing, boundary transport, quadrature, and reconstruction map. A continuum error estimate requires an inverse or stability relation such as ||u-u_N||_X <= C ||L(u-u_N)||_Y, together with control of the boundary defect. The weighted-residual framework is therefore a transport system: the continuum equation is projected into a finite observation family, solved algebraically, reconstructed, and audited through a continuum residual. A small algebraic residual in the assembled matrix is not automatically a small differential residual because scaling, preconditioning, and basis transformations change the algebraic norm.

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

A Galerkin method uses the same trial and test space and solves a(u_N,v_N)=l(v_N) for every v_N in V_N. If a is bounded and coercive, the matrix inherits a stable variational interpretation. A Petrov–Galerkin method uses V_N != W_N; its stability is controlled by the discrete inf-sup constant beta_N = inf_u sup_v a(u,v)/(||u||_V||v||_W). Sparsity or directional bias in the test basis is useful only if beta_N remains bounded away from zero.

A tau method enforces the differential equation in a selected set of modal equations and replaces some highest-mode equations by boundary constraints. The omitted modal residual is a named residue rather than an invisible approximation. Collocation imposes L u_N(x_i)=f(x_i) at nodes and enforces boundary equations separately. Its certificate depends on interpolation stability, node distribution, and discrete norm equivalence. Least squares minimizes ||L v-f||²+||Bv-g||², yielding a positive semidefinite normal system, but the normal operator behaves like L*L and can square the condition number. These formulations may generate the same polynomial space and still differ in conditioning, conservation, eigenvalues, and off-grid residual. Their equivalence must be demonstrated by a bounded transport map, not inferred from matching dimensions.

Approximation error, consistency, stability and convergence

The total numerical error is not a single object. It is the composition of approximation debt, consistency debt, stability amplification, algebraic-solver error, and floating-point error. The best-approximation term is inf_(v_N in V_N)||u-v_N||. Consistency measures whether the exact solution satisfies the discrete equations after projection, quadrature, geometry approximation, and boundary transport. Stability measures how the inverse discrete operator amplifies these defects.

For a coercive conforming Galerkin method, Cea’s estimate gives ||u-u_N||_V <= (M/alpha) inf_(v_N in V_N)||u-v_N||_V. If the form is approximated numerically, a Strang estimate adds sup_w |a(v_N,w)-a_N(v_N,w)|/||w||. Spectral approximation can make the first term exponentially small while the quadrature or stability term remains dominant. The xSCD counterkernel is an unstable boundary closure built on an exponentially accurate basis: basis accuracy survives, but the method diverges. Convergence is therefore certified only when approximation, consistency, and stability are all transported through the same norms.

Sobolev regularity versus spectral convergence

Spectral convergence is controlled by the regularity of the represented function relative to the chosen basis and coordinates. If u in H^s, a typical polynomial or Fourier projection estimate is ||u-P_Nu||_(H^m) <= C N^(m-s)||u||_(H^s) for m<=s. If u extends analytically to a complex neighborhood, coefficients often satisfy |u_hat_n| <= C rho^(-n), giving geometric convergence until rounding error or conditioning dominates.

Analyticity alone does not guarantee rapid convergence at practical resolutions. A boundary layer u(x)=exp(-(1+x)/epsilon) is entire but contains a scale epsilon, so the onset of geometric decay moves to larger N as epsilon decreases. A mapped coordinate can either resolve or create apparent singular behavior. Coefficient decay is therefore an observer, not a complete diagnosis. Slow decay may indicate a true singularity, a basis mismatch, localized structure, an irregular geometry map, or insufficient precision. An L-clock mutation is triggered when N-clock refinement fails to reduce both coefficient tail and continuum residual.

Modal versus nodal representations

A modal representation stores coefficients in u_N(x)=sum_k u_hat_k phi_k(x). It exposes regularity through coefficient decay and often yields sparse differentiation, multiplication, filtering, and preconditioning. A nodal representation stores values u_j=u_N(x_j). It simplifies boundary imposition, pointwise nonlinearities, diagnostics, and coupling to data or geometry.

The two carriers are connected by u_j=sum_k V_jk u_hat_k, where V_jk=phi_k(x_j). The transform must be stable and invertible on the finite space. Equality of nodal vectors and equality of coefficient vectors are different identity relations; they become equivalent only after reconstruction. Nonlinear pseudospectral methods deliberately alternate carriers: modal coefficients are transformed to nodes, nonlinear products are formed, and the result is projected back. This cycle creates aliasing whenever the nodal carrier is too small for the product bandwidth. The transform, nonlinear evaluation, and projection must therefore be audited as one composite transport.

Reproducible numerical experiments in Julia, Python or MATLAB

A reproducible experiment specifies the continuum problem, discrete carrier, solver, precision, and validation procedure. The problem manifest includes equation, coefficients, domain, boundary and initial data, and known invariants. The discretization manifest includes basis normalization, nodes, quadrature, geometry map, polynomial degree, dealiasing rule, and time step. The solver manifest includes algorithm, preconditioner, stopping criterion, nonlinear tolerance, and arithmetic precision.

A valid output contains more than a plot. For stationary problems it should report the reconstructed continuum residual, boundary residual, coefficient tail, algebraic residual, condition estimate, and error against an independently certified reference. For time-dependent systems it should report spatial and temporal refinement separately, invariant drift, and accumulated phase or energy error. Randomized experiments require fixed seeds and distributions over repeated runs. Reproducibility is part of the certificate because a result that cannot be replayed lacks provenance and cannot support a numerical claim.

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(ikx), where u_hat_k=(1/2pi) integral_0^(2pi)u(x)exp(-ikx)dx. The carrier is the periodic function space L2(S1), not an interval function whose endpoint values merely happen to coincide. Fourier modes diagonalize constant-coefficient derivatives: D^m exp(ikx)=(ik)^m exp(ikx). Translation is also diagonal because the coefficient of u(x-a) is exp(-ika)u_hat_k.

Parseval’s relation ||u||²=2pi sum_k|u_hat_k|² transports continuous energy into coefficient space. Sobolev regularity becomes weighted square summability: ||u||_(H^s)² ~ sum_k(1+k²)^s|u_hat_k|². The Fourier carrier is therefore ideal when the operator and geometry respect periodic translation. Its L-clock mutation trigger is nonperiodic boundary data, localized singularity, irregular geometry, or strong coefficient variation that destroys diagonal structure.

Truncation, interpolation and projection

The Fourier projection P_Ku=sum_(|k|<=K)u_hat_k exp(ikx) is the best L2 approximation in the retained trigonometric space. Interpolation constructs a trigonometric polynomial that agrees with u at equispaced nodes. Its discrete coefficients satisfy an aliasing identity of the form u_tilde_k=sum_m u_hat_(k+mM), where M is the number of grid points.

Projection discards unresolved modes. Interpolation folds unresolved modes into visible frequencies. These are different residues, even though both approximations may converge at comparable rates for smooth data. A high-frequency mode and its low-frequency alias are identical on the sampling grid but distinct in the continuum carrier. This is a canonical xSCD identity failure: nodal equivalence modulo M is not continuum frequency equivalence. Projection and interpolation must therefore be distinguished whenever nonlinear products, stability, or spectral diagnostics are involved.

Parseval identities

Parseval converts inner products into coefficient sums: (u,v)=2pi sum_k u_hat_k conjugate(v_hat_k). Hence ||u_x||²=2pi sum_k k²|u_hat_k|². For the heat equation u_t=nu u_xx, one obtains 1/2 d||u||²/dt=-nu||u_x||², and modewise d|u_hat_k|²/dt=-2nu k²|u_hat_k|².

In discrete calculations, Parseval is exact only under a consistent transform normalization and mode convention. Real-valued fields satisfy u_hat_(-k)=conjugate(u_hat_k), and the Nyquist mode in an even grid requires special treatment. A numerical energy assembled with an inconsistent normalization is not an approximation to the intended continuum energy; it is a differently scaled observable. Parseval is therefore both an analytical identity and a carrier-typing rule for diagnostics.

Trigonometric interpolation

With nodes x_j=2pi j/M, the interpolation basis satisfies discrete orthogonality because sum_j exp(i(k-l)x_j)=M delta_(k,l mod M). The interpolant is evaluated through its DFT coefficients. Differentiation is performed by multiplying each mode by ik and transforming back.

The grid identifies frequencies modulo M: exp(i(k+mM)x_j)=exp(ikx_j). This identity makes the FFT possible but also creates aliasing. Interpolation is stable for smooth periodic functions because the nodes match the periodic geometry. It is unsuitable for nonperiodic functions because the periodic quotient forces an artificial endpoint identification, producing Gibbs oscillations and only algebraic global convergence.

Discrete Fourier transform and FFT

The DFT is u_tilde_k=sum_(j=0)^(M-1)u_j exp(-2pi ijk/M). Direct evaluation costs O(M²). The FFT factors the transform into smaller structured operations and computes the same finite mapping in O(M log M) work. The FFT changes the implementation carrier, not the mathematical finite transform.

Important implementation details include transform normalization, negative-frequency ordering, Nyquist treatment, real-valued symmetry, data layout, and multidimensional stride order. The algorithmic debt consists of floating-point roundoff and permutation error. A useful implementation certificate is the round-trip test ||F^-1(Fu)-u||/||u||, together with scaling versus M. Speed is a cost certificate and does not replace accuracy or inverse consistency.

Spectral differentiation and integration

For Fourier coefficients, differentiation is u_hat_k -> ik u_hat_k, and the mth derivative is (ik)^m u_hat_k. High derivatives amplify high-mode roundoff and unresolved noise. Integration is u_hat_k -> u_hat_k/(ik) for k!=0. The zero mode is an integration constant and must be retained explicitly.

For periodic Poisson, -u_xx=f gives u_hat_k=f_hat_k/k² for nonzero k. The zero mode produces the compatibility condition f_hat_0=0, while u_hat_0 is fixed by a gauge such as zero mean. Nullspace and compatibility are therefore part of the operator carrier. Dividing by k² without declaring the zero mode produces an untyped inverse.

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

If u and v have Fourier coefficients, the product satisfies (uv)_hat_k=sum_(p+q=k)u_hat_pv_hat_q. A product of two fields resolved to degree K contains modes up to 2K. Evaluating it on the original grid folds the excess modes into lower frequencies.

For quadratic nonlinearities, the three-halves rule pads the spectrum to a larger grid, transforms to physical space, forms the product, transforms back, and truncates to the original modal range. The equivalent two-thirds rule retains only a safe lower portion of the grid-supported spectrum. Dealiasing does not assert that discarded modes are physically zero; it ensures that their projection does not corrupt retained modes. Cubic or higher nonlinearities require larger padding or staged projection.

Gibbs phenomena and filtering

A discontinuity forces Fourier coefficients to decay algebraically, typically as 1/|k|. Truncation produces oscillations whose spatial width shrinks with N but whose peak overshoot remains finite. Filtering modifies coefficients through u_hat_k -> sigma(|k|/K)u_hat_k. An exponential filter might use sigma(eta)=exp(-alpha eta^p).

Filtering introduces model debt because the filtered state solves a modified discrete problem. Spectral viscosity similarly adds high-frequency dissipation. These tools can stabilize under-resolved calculations and suppress nonphysical oscillations, but they alter phase, conservation, and effective resolution. The filter parameters, affected modes, and induced energy loss must be recorded as part of the method.

Week 3 — Fourier methods for periodic PDEs

Fourier–Galerkin and Fourier collocation methods

Fourier–Galerkin enforces modal residual orthogonality, while Fourier collocation enforces the equation at grid points. For constant-coefficient linear operators and exact transforms, the two formulations descend to the same diagonal modal equations. For nonlinear or variable-coefficient terms, collocation forms interpolated products and therefore carries aliasing residue.

A genuine equivalence requires a commuting diagram between continuum projection, nodal interpolation, modal transformation, and discrete operator application. Visual agreement does not establish equivalence. In xSCD terms, Galerkin and collocation share a periodic state carrier but use different observation families and different quotient relations for the residual.

Poisson, heat, wave and Helmholtz equations

Periodic Poisson exposes compatibility and gauge. Heat exposes dissipative stiffness because eigenvalues behave as -nu k². Wave equations expose phase accuracy because each Fourier mode satisfies u_hat_k''+c²k²u_hat_k=0. Helmholtz exposes resonance and resolvent amplification because u_hat_k=f_hat_k/(k²-kappa²).

These equations demonstrate that diagonalization does not imply harmless inversion. Near resonance, the formula exists but small perturbations in forcing or kappa produce large changes in the solution. The relevant certificate includes denominator separation or resolvent bounds, not merely successful division.

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

Viscous Burgers, u_t+u u_x=nu u_xx, combines nonlinear transfer with diffusion. The advective, conservative, and split forms are analytically equivalent but can differ discretely under aliasing. Korteweg–de Vries, u_t+6u u_x+u_xxx=0, combines nonlinear steepening with dispersive phase transport and possesses invariants whose drift is a long-time diagnostic.

Kuramoto–Sivashinsky, u_t+u u_x+u_xx+u_xxxx=0, contains linearly unstable low modes and strongly damped high modes, producing a stiff mixed spectrum. Allen–Cahn, u_t=epsilon²u_xx-F'(u), is a gradient flow with energy E=integral[epsilon²|u_x|²/2+F(u)]dx. Each equation therefore tests a different structural certificate: energy balance, phase, instability–dissipation coexistence, or free-energy decay.

Conservation laws and invariant drift

A continuous invariant I(u) becomes a discrete quantity I_N(u_N). Its derivative generally has the form dI_N/dt=R_alias+R_quad+R_boundary+R_time. Exact preservation requires the discrete algebra to reproduce the continuum cancellation.

For periodic Burgers, integral u²u_x dx=0 follows from a product rule and periodic boundary cancellation. A collocation method without dealiasing may violate this identity. Invariant drift can reveal a structurally incorrect computation long before pointwise error becomes visually obvious. The full observer family includes mass, momentum, energy, spectral transfer, and boundary flux.

Semidiscrete stability

After spatial discretization, M u_dot=A_N u+N_N(u). A semidiscrete certificate may be an energy estimate, an eigenvalue bound, or a resolvent bound. For normal Fourier operators, eigenvalues often characterize linear stability. For non-normal variable-coefficient or coupled systems, transient growth can occur even when all eigenvalues lie in the stable half-plane.

A stable spatial discretization is only locally closed. Its composition with a time integrator requires a separate certificate. This prevents the invalid inference that energy stability of the semidiscrete operator automatically implies stability of the fully discrete scheme.

Explicit, implicit, IMEX and exponential time integrators

Explicit Runge–Kutta methods require the scaled eigenvalues dt lambda_j to lie inside the stability region. Diffusion therefore gives a restriction near dt=O(N^-2) in Fourier space, while fourth-order dissipation gives O(N^-4). Implicit methods remove this linear restriction but require solves.

IMEX schemes treat stiff linear terms implicitly and nonlinear terms explicitly. Exponential integrators use u(t+dt)=exp(dtL)u(t)+integral exp((dt-s)L)N(u(t+s))ds and approximate the integral through phi-functions. Integrating-factor methods transform away the linear part. Each approach introduces a different debt: nonlinear explicit stability, implicit-solver error, matrix-function approximation, or splitting commutators. The time method is part of the full numerical carrier rather than an interchangeable postprocessing step.

Part II — Orthogonal polynomials and fast transforms

Week 4 — General orthogonal-polynomial systems

Weighted inner products

Let w(x)>0 on (a,b) and define (u,v)_w=integral_a^b u v w dx. A family {p_n} is orthogonal if (p_n,p_m)_w=gamma_n delta_nm. The measure w(x)dx is part of the carrier; the same algebraic polynomial family under a different measure represents a different approximation system.

Weights redistribute approximation sensitivity. Jacobi weights (1-x)^alpha(1+x)^beta emphasize or suppress endpoints depending on the parameters. Weighted regularity and unweighted regularity are distinct. Boundary traces may remain large even when the weighted norm is small, so endpoint observers must be included explicitly.

Three-term recurrence and Jacobi matrices

Orthogonal polynomials satisfy x p_n=a_(n+1)p_(n+1)+b_n p_n+a_n p_(n-1) in orthonormal normalization. Multiplication by x is represented by a symmetric tridiagonal Jacobi matrix. This triadic closure is the algebraic mechanism behind stable evaluation, quadrature, zero interlacing, and many sparse operators.

The finite Jacobi matrix is the compression of the multiplication operator to the first N+1 modes. Its eigenvalues are Gaussian quadrature nodes. Positivity of the recurrence coefficients is the carrier certificate supplied by the underlying positive measure.

Zeros and interlacing

A degree-n orthogonal polynomial has n simple roots inside the support of the measure, and roots of consecutive degrees interlace. This follows from sign-change and orthogonality arguments. Interlacing provides safe brackets for root refinement and explains the stability of Gaussian nodes.

Asymptotically, roots cluster according to the equilibrium density of the interval. Endpoint clustering is not a numerical accident; it is encoded by the measure and polynomial family. A computed root set should be audited for reality, simplicity, interval inclusion, and interlacing.

Gaussian, Radau and Lobatto quadrature

Gaussian quadrature with N+1 nodes is exact for polynomials through degree 2N+1. Radau fixes one endpoint and is exact through degree 2N; Lobatto fixes both endpoints and is exact through degree 2N-1. The rule is typed by its measure, nodes, weights, and exactness space.

Exactness does not extend automatically to nonlinear products, variable coefficients, mapped Jacobians, or singular integrands. Under-integration becomes quadrature residue and can create aliasing or violate energy identities. Positive weights provide a useful discrete norm, but equivalence with the continuum norm is limited to the declared polynomial range.

Christoffel–Darboux kernels

The projection kernel is K_N(x,y)=sum_(j=0)^N p_j(x)p_j(y)/gamma_j. It reproduces every polynomial in P_N. The Christoffel–Darboux formula compresses the quadratic modal sum to an expression involving only p_N and p_(N+1) divided by x-y.

This is an exact finite-carrier compression, not an approximation. On the diagonal it produces a derivative formula. The term “kernel” is explicitly typed because projection kernels, Green functions, integral kernels, and learned convolution kernels have different domains, codomains, measures, and composition laws.

Reproducing and projection kernels

The projector is P_Nu(x)=integral K_N(x,y)u(y)w(y)dy. The diagonal K_N(x,x) quantifies the amplification of point evaluation, while lambda_N(x)=1/K_N(x,x) is the Christoffel function. The inequality |q(x)|²<=K_N(x,x)||q||_w² shows why endpoint evaluation can be ill-conditioned.

These quantities guide stable sampling, weighted least squares, and interpolation analysis. Large diagonal growth represents a genuine local concentration effect and must be carried into trace and conditioning estimates.

Modified weights and Christoffel transformations

Multiplying a measure by (x-a) or another low-degree factor creates a new orthogonal family connected to the original through structured formulas. Endpoint-restricted projection kernels become orthogonal polynomials for the modified weight. The inverse operation is related to Geronimus transformations.

A weight change is an L-clock mutation because it changes orthogonality, recurrence coefficients, projection, and norms. Connection coefficients are the transport map between the old and new modal carriers.

Stable polynomial evaluation

Direct monomial evaluation becomes unstable at high degree because powers are badly scaled and coefficients can cancel strongly. Clenshaw’s recurrence evaluates an orthogonal expansion backward in O(N) operations. Barycentric interpolation evaluates a nodal polynomial without explicitly forming the ill-conditioned Vandermonde inverse.

The implementation should be validated by round-trip transform tests, comparison with high precision, and endpoint stress cases. An algebraically correct formula can still be a numerically invalid carrier if it amplifies rounding errors excessively.

Week 5 — Jacobi, Legendre and Chebyshev systems

Jacobi Sturm–Liouville operator

Jacobi polynomials J_n^(alpha,beta) are orthogonal under w_(alpha,beta)=(1-x)^alpha(1+x)^beta. They satisfy L_(alpha,beta)J_n=lambda_nJ_n, where L u=-w^-1 d/dx[w_(alpha+1,beta+1)u'] and lambda_n=n(n+alpha+beta+1).

The operator domain includes endpoint flux conditions. The vanishing of the leading coefficient at x=±1 makes the differential expression singular, but the weighted operator is self-adjoint on its natural domain. The function space, measure, and endpoint behavior are therefore inseparable parts of the carrier.

Self-adjointness and weighted integration by parts

Weighted integration by parts gives (Lu,v)_w-(u,Lv)_w=boundary flux. When the domain conditions make the flux vanish, L is self-adjoint and (Lu,u)_w=||u'||_(w_(alpha+1,beta+1))².

This identity is the certificate for real eigenvalues, orthogonality, and positive energy. Formal differentiation without endpoint control is insufficient. The course therefore retains boundary flux explicitly until the domain discharges it.

Parameter-shift and derivative identities

Jacobi differentiation satisfies dJ_n^(alpha,beta)/dx = (n+alpha+beta+1)/2 J_(n-1)^(alpha+1,beta+1). Differentiation lowers degree and raises both weight parameters. Higher derivatives continue the same transport.

Sparse differentiation is obtained by allowing the target carrier to change. Forcing all derivatives back into one fixed basis creates unnecessary density. Parameter conversion operators then combine terms in a common space. This is a canonical xSCD case where controlled carrier mutation preserves structure.

Legendre and Chebyshev specializations

Legendre polynomials correspond to alpha=beta=0 and use the unweighted inner product. Their variational structure is convenient for Galerkin methods. Chebyshev polynomials satisfy T_n(cos theta)=cos(n theta) and are orthogonal under (1-x²)^(-1/2), enabling FFT-based transforms.

Although both span P_N, they define different projections, quadratures, endpoint scaling, and discrete norms. A basis selection should therefore be justified by the operator and computational transport: Legendre for unweighted weak forms, Chebyshev for fast cosine transforms and nodal interpolation.

Connection coefficients between polynomial families

If p_n^A=sum_(k<=n)C_(k,n)p_k^B, the triangular matrix C transports coefficients between bases. Adjacent Jacobi parameter shifts often yield banded maps; large or noninteger shifts are denser but structured.

A connection algorithm is certified by forward accuracy, inverse consistency, endpoint stability, and complexity. Connection matrices permit each operator term to be represented in its sparsest natural basis and then transported to a common target space.

Endpoint behaviour and singular Jacobi weights

Endpoint values satisfy asymptotic laws such as J_n^(alpha,beta)(1)~n^alpha and |J_n^(alpha,beta)(-1)|~n^beta. Negative but integrable parameters create singular weights that can represent endpoint singularities efficiently.

If a solution behaves like (1-x)^mu v(x), factoring the known singular term or selecting a compatible Jacobi carrier can restore rapid coefficient decay. Persistent algebraic decay in an ordinary Legendre expansion is a mutation signal rather than a reason for unlimited degree growth.

Polynomial interpolation and projection estimates

Orthogonal projection is best in its weighted L2 norm. Interpolation additionally depends on the Lebesgue constant. Chebyshev–Lobatto nodes give Lambda_N=O(log N), while equally spaced polynomial interpolation can be exponentially unstable.

For analytic functions, interpolation and projection errors typically decay geometrically. Derivative errors lose powers of N because differentiation amplifies high modes. A small function error therefore does not certify a small high-order derivative error.

Week 6 — Fast algorithms for orthogonal polynomials

Fast synthesis and analysis

Synthesis maps coefficients to values, while analysis maps sampled values to coefficients. Dense transforms cost O(N²). Fourier and Chebyshev transforms use FFTs; Legendre and Jacobi transforms exploit asymptotic oscillation, low-rank blocks, recurrence structure, or conversion through Chebyshev space.

The fast algorithm is represented as T_fast=T_exact+E_alg, where E_alg includes compression, asymptotic, and floating-point residues. Complexity must be reported together with accuracy and inverse consistency.

Discrete cosine and sine transforms

Chebyshev nodes x_j=cos(jpi/N) satisfy T_k(x_j)=cos(kjpi/N), reducing Chebyshev transforms to DCTs. Sine transforms represent parity or endpoint-vanishing bases. Different DCT and DST types encode different endpoint conventions and normalizations.

Using the wrong transform type changes the represented polynomial. The convention is therefore carrier metadata, not an implementation detail.

Fast Chebyshev–Legendre transforms

Chebyshev–Legendre conversion is dense in direct form. Fast algorithms exploit the fact that high-degree Legendre and Chebyshev modes share related oscillatory asymptotics. Blockwise low-rank approximation or divide-and-conquer reduces cost toward quasi-linear complexity.

Low modes and endpoint-sensitive regions are normally handled by direct stable recurrence. Both forward and inverse transforms must be validated because endpoint scaling differs between the bases.

Fast connection-coefficient transforms

General Jacobi connection transforms can be built from adjacent parameter shifts, hierarchical low-rank approximations, or recurrence-based algorithms. Integer shifts often preserve narrow bandwidth; fractional shifts generate structured dense matrices.

Their value is not only speed. They preserve the sparse operator graph differentiate -> shift parameters -> multiply -> convert -> combine, avoiding premature dense assembly.

Clenshaw and barycentric algorithms

Clenshaw evaluates modal series through backward recurrence. Barycentric formulas evaluate nodal interpolants using precomputed weights. Each algorithm belongs to a different native carrier and is stable because it respects the structure of that carrier.

Derivative evaluation can use differentiated recurrences or barycentric differentiation, but repeated dense nodal differentiation may amplify conditioning more than coefficient-space differentiation.

Structured multiplication and differentiation operators

Multiplication by x is tridiagonal in an orthogonal basis. Multiplication by a polynomial of degree m is banded with width proportional to m. Smooth coefficients can be spectrally truncated to obtain controlled almost-banded multiplication.

Differentiation is diagonal or narrowly banded between shifted Jacobi spaces. Operator assembly should retain this composition rather than flattening everything into a dense matrix, which erases sparsity, provenance, and often conditioning advantages.

Almost-banded matrices

Coefficient-space boundary-value operators are frequently banded except for a few dense boundary rows. The dense rows encode global trace conditions and should remain explicit. Adaptive QR or structured elimination solves such systems efficiently.

Absorbing the boundary rows into a dense matrix preserves algebraic equivalence but destroys the visible separation between local differential transport and global constraints.

Adaptive coefficient-space computation

An adaptive solver extends the coefficient vector until the trailing envelope, residual coefficients, and boundary defects fall below tolerance. A single small final coefficient is not a reliable criterion because coefficients may oscillate or cancel.

When increasing N fails, the method must diagnose whether to change the weight, map, domain partition, formulation, or precision. N-clock refinement becomes an L-clock mutation when the carrier no longer matches the solution structure.

Multivariate orthogonal polynomials and Koornwinder constructions

On triangles and tetrahedra, Koornwinder-type bases combine nested Jacobi polynomials with coordinate-dependent factors. The parameter shifts encode the simplex Jacobian and preserve orthogonality. These bases avoid unstable monomial representations and yield structured differentiation and mass operators.

The multidimensional carrier introduces index-set growth and coordinate-coupling debt. Orthogonality of one-dimensional factors does not automatically certify multidimensional conditioning or boundary behavior.

Part III — Boundary-value and eigenvalue problems

Week 7 — Second-order boundary-value problems

Dirichlet, Neumann, Robin and mixed conditions

For -epsilon u''+p(x)u'+q(x)u=f, separated boundary conditions are a_-u(-1)+b_-u'(-1)=c_- and a_+u(1)+b_+u'(1)=c_+. The differential expression and boundary law together define the operator. Changing boundary conditions changes its domain, nullspace, spectrum, and solvability.

Inhomogeneous conditions are homogenized by writing u=v+u_b, where u_b satisfies the traces. This isolates boundary transport from the interior solve. Neumann problems may require compatibility and a gauge.

Weak formulations and coercivity

For -u''+alpha u=f, integration by parts gives integral u'v'+alpha uv = integral fv + boundary contribution. Dirichlet conditions are usually built into the trial space; Neumann or Robin terms enter naturally.

Coercivity requires a(v,v)>=c||v||_(H1)² on the admissible space. For pure Neumann problems, constants form a nullspace and coercivity holds only on the zero-mean quotient. Boundary terms remain explicit until discharged by the domain.

Legendre- and Chebyshev–Galerkin methods

Legendre–Galerkin uses the unweighted inner product and boundary-adapted Legendre bases. Chebyshev–Galerkin uses the Chebyshev weight and benefits from fast transforms. Both can achieve spectral convergence, but their weak forms, mass matrices, and trace scaling differ.

The basis is selected together with the operator and quadrature. Reusing a Chebyshev basis inside an unweighted variational identity without transporting the weight is an untyped formulation.

Boundary-adapted modal bases

A basis phi_k=P_k+a_kP_(k+1)+b_kP_(k+2) can be chosen to satisfy two homogeneous boundary conditions. For Dirichlet Legendre problems, phi_k=L_k-L_(k+2). Such bases span the constrained polynomial space exactly.

The neighboring-mode combination exploits recurrence and endpoint formulas, producing diagonal or banded stiffness matrices and narrow mass matrices. Boundary debt is discharged by construction.

Tau and collocation formulations

Tau methods replace selected modal equations by boundary equations. The displaced high-mode residual must be recognized. Collocation imposes interior equations at nodes and boundary equations at endpoints, producing a dense differentiation system.

Their polynomial outputs may coincide in special constant-coefficient cases, but conditioning and residual distribution can differ. Off-grid continuum residuals are required for liftback.

Equivalence and non-equivalence of formulations

Two formulations are equivalent only when there exists an invertible, uniformly controlled transformation mapping one finite operator, data vector, and reconstruction rule to the other. Matching degree and node count is insufficient.

Variable coefficients, nonlinear products, approximate quadrature, and different trace enforcement can break equivalence. Cross-formulation comparison is therefore a diagnostic, not an assumed identity.

Sparse stiffness and mass matrices

Define S_jk=integral phi_k'phi_j' and M_jk=integral phi_kphi_j. A boundary-adapted Legendre basis can make S diagonal and M pentadiagonal. Fixed-band systems can be solved in linear time.

Variable coefficients broaden bandwidth according to their modal content. Truncating a smooth coefficient creates an explicitly controlled multiplication residue.

Error estimates

Galerkin quasi-optimality gives ||u-u_N||_(H1)<=C inf_(v_N)||u-v_N||_(H1). Duality can improve L2 error when the adjoint problem is regular. Collocation estimates add interpolation and quadrature terms.

The complete liftback is algebraic residual -> reconstructed polynomial -> differential residual -> boundary residual -> inverse estimate -> continuum error.

Week 8 — Modern sparse spectral solvers

Ultraspherical and coefficient-space methods

Repeated differentiation becomes dense if every derivative is immediately converted back into Chebyshev coefficients. Ultraspherical methods retain each derivative in its natural shifted Gegenbauer space. Differentiation is sparse, and conversion between adjacent spaces is banded.

A variable-coefficient differential operator becomes an almost-banded composition plus dense boundary rows. The method’s central xSCD move is permitting controlled carrier change under differentiation.

Banded differentiation and conversion operators

Chebyshev differentiation maps to ultraspherical parameter one; further derivatives map to higher parameters. Conversion operators align lower-order terms with the highest-order target space.

Every coefficient vector must therefore carry basis metadata. Adding vectors from different ultraspherical spaces without conversion is an invalid operation.

Integral reformulations

Integrating the differential equation transforms an unbounded differentiation operator into a second-kind equation resembling I+K. For u''=f, u=I²f+c_0+c_1x, with constants fixed by the boundary conditions.

Integration acts as an operator preconditioner because it suppresses high-frequency growth. The constants retain the differentiation nullspace and cannot be discarded.

Spectral integration

Chebyshev and Jacobi antiderivatives obey short coefficient recurrences. Integration divides high-mode amplitudes and is usually better conditioned than differentiation.

Repeated integration requires explicit storage of integration constants. Boundary conditions then determine them through a small dense system.

Conditioning of differentiation matrices

Polynomial collocation differentiation matrices have rapidly growing norms because endpoint nodes cluster and differentiation is unbounded. A first derivative can scale near N²; second derivatives can grow near N⁴, depending on representation and boundary rows.

Matrix condition numbers must be interpreted in the physical norm and basis scaling. Nevertheless, unpreconditioned dense differentiation eventually loses significant digits.

Diagonal and operator preconditioning

Diagonal scaling equalizes modal magnitudes. Operator preconditioning approximates the inverse principal differential operator, often through integration. The goal is a system whose spectral or singular-value distribution remains controlled as N grows.

A valid preconditioner must preserve boundary conditions, gauges, and nullspaces. It is not merely a matrix that reduces GMRES iterations.

Adaptive truncation

The solver increases N until a trailing coefficient block, residual, and boundary error all meet tolerance. The coefficient tail is interpreted relative to the basis and operator order.

Failure to converge triggers analysis of singularity, layer scale, geometry, precision, or formulation. Increasing N indefinitely is not an admissible universal response.

Residual and coefficient-tail error indicators

Coefficient tails estimate unresolved representation content. Residuals estimate equation violation. Boundary defects measure operator-domain violation. All three are necessary because each can be small while another remains large.

An a posteriori bound uses ||e||_X<=C_stab||R||_Y, with C_stab derived or estimated. A residual without a stability constant does not determine the solution error.

Matrix-free Krylov solution

In multidimensional problems, the operator is applied through transforms, tensor contractions, differentiation, coefficient multiplication, and geometry factors rather than assembled. GMRES handles general nonsymmetric systems; conjugate gradients require symmetry and positive definiteness in the chosen inner product.

Preconditioned algebraic residuals must be converted back to physical residuals before certification. Matrix-free execution reduces memory, not automatically iteration count.

Week 9 — Higher-order and constrained problems

Fourth- and higher-order equations

A fourth-order operator requires four boundary conditions and higher trial regularity in a direct weak form. The biharmonic model can be treated directly in H² or through a mixed reduction. High derivatives produce stronger stiffness and conditioning growth.

Coefficient-space parameter shifts and integration reformulations retain sparsity. Boundary constraints must match the operator order and remain independent.

Generalized Jacobi bases

Generalized Jacobi functions multiply Jacobi polynomials by endpoint factors such as (1-x)^r(1+x)^s. These factors enforce repeated value or derivative conditions and represent known endpoint singularities.

The basis is selected so the differential operator maps between related weighted spaces with narrow bandwidth. This combines boundary discharge and sparse transport.

Dual Petrov–Galerkin methods

Odd-order and non-self-adjoint operators often admit trial and test bases adapted to the operator and adjoint boundary conditions. The matrix a(phi_k,psi_j) can become diagonal or banded.

Sparsity alone is insufficient; the trial–test pairing requires a uniform inf-sup certificate. The adjoint boundary form determines the correct dual carrier.

Biharmonic and Cahn–Hilliard equations

The biharmonic equation Delta²u=f can use clamped or simply supported boundaries, which define different operator domains. Cahn–Hilliard is u_t=Delta mu, mu=F'(u)-epsilon²Delta u. It conserves mass and dissipates free energy: dE/dt=-||grad mu||².

A valid spectral discretization preserves the zero mode and reproduces a discrete energy law. Mixed and direct fourth-order formulations introduce different conditioning and boundary debts.

Mixed formulations

Introducing auxiliary variables lowers differential order: v=-Delta u, -Delta v=f. Mixed incompressible formulations use pressure as a multiplier enforcing divergence-free velocity.

The spaces must satisfy compatibility or an inf-sup condition. Arbitrary polynomial pairs can create spurious multiplier modes. The mixed carrier enlarges the state while simplifying local differentiation.

Differential-algebraic boundary constraints

Constraints produce systems M u_dot=F(u)+C^Tlambda, C u=d. The multiplier enforces a manifold condition. The DAE index measures how many differentiations expose explicit evolution.

Time integration must preserve or project onto the constraint. Dense trace rows and scaling strongly affect numerical conditioning.

Exact enforcement of multiple boundary conditions

Boundary constraints may be enforced through adapted bases, tau rows, bordering, nullspace projection, or multipliers. Basis construction satisfies them identically; bordering is more flexible; nullspace projection reduces the state to admissible coordinates.

The boundary matrix must have full rank. Redundant conditions produce singular systems, while incompatible conditions yield no solution.

Nullspaces, compatibility and gauge conditions

Neumann Poisson has a constant nullspace and requires integral f=0. Pressure in incompressible flow is defined up to a constant. Periodic antiderivatives also carry a zero mode.

The numerical carrier should reproduce the correct continuum nullspace. A gauge fixes a representative in the quotient. Removing a coefficient without declaring the gauge can destroy conservation or consistency.

Week 10 — Spectral eigenvalue problems

Sturm–Liouville eigenproblems

A Sturm–Liouville problem is -(p u')'+q u=lambda w u with self-adjoint boundary conditions. Its weak form yields A c=lambda M c, with symmetric A and positive-definite M under appropriate assumptions.

Variational principles give upper or lower eigenvalue bounds and explain accurate approximation of low modes. The finite pencil approximates the operator only after boundary and weight transport are correctly included.

Generalized matrix eigenvalue formulations

For A x=lambda B x, a positive-definite B permits reduction to a standard symmetric problem. Singular B introduces infinite eigenvalues or algebraic constraints. Tau boundary rows can produce singular pencils.

Physical finite eigenvalues must be distinguished from constraint-associated roots. Matrix-pencil structure is part of the carrier.

Spurious eigenvalues and spectral pollution

Spurious eigenvalues can arise from nonconforming spaces, inconsistent boundary conditions, singular pencils, unbounded-domain truncation, or non-self-adjoint sensitivity. Apparent convergence at two resolutions may still reflect a common structural defect.

Validation requires residuals, eigenspace convergence, cross-formulation comparison, and variational bounds where available. The reconstructed eigenfunction must satisfy the continuum boundary conditions and differential equation.

Non-self-adjoint operators

When A != A*, right eigenvectors need not be orthogonal. Left and right eigenvectors determine the eigenvalue condition number kappa(lambda)=||x||||y||/|y*x|. Small y*x indicates high sensitivity.

Eigenvalues alone do not determine transient growth. Non-normal systems require resolvent and pseudospectral observers.

Pseudospectra and resolvent growth

The epsilon-pseudospectrum consists of z where ||(zI-A)^-1||>1/epsilon. It also describes eigenvalues reachable by perturbations of norm below epsilon.

Large pseudospectral regions reveal transient amplification and spectral instability. Computation uses smallest singular values, contour methods, or resolvent actions. Resolution continuation is required to distinguish continuum features from discretization artifacts.

Polynomial and rational filtering

Polynomial filters amplify desired spectral intervals and damp others. Chebyshev filters are effective when the unwanted spectrum lies in a known real interval. Rational filters use shifted inverses and approximate spectral projectors.

Filter accuracy depends on approximation error and shifted-solve residuals. These errors propagate into eigenspace estimates and must be included.

Contour and shift-invert methods

Shift-invert maps eigenvalues near sigma to dominant reciprocals (lambda-sigma)^-1. Contour methods approximate P=(1/2pi i)integral_Gamma(zI-A)^-1dz through quadrature.

The solves at contour nodes can be parallelized. Their residuals and quadrature errors control projector accuracy. The contour must avoid poorly resolved or highly non-normal regions unless resolvent growth is accounted for.

Validation by residual and backward error

For an approximate pair, r=A x_hat-lambda_hat B x_hat. A normalized backward error is eta=||r||/[(||A||+|lambda_hat|||B||)||x_hat||].

Small backward error means the pair is exact for a nearby pencil. Forward eigenvalue error additionally depends on conditioning. Reporting residual without sensitivity is inadequate for non-normal problems.

Part IV — Time evolution, singularity and nonlocality

Week 11 — Time-dependent PDEs

Method of lines

Spatial discretization produces M u_dot=F(u,t). The mass matrix may be diagonal, sparse, dense, or singular. The method of lines separates spatial and temporal construction but not their stability.

The discrete spectrum, numerical range, invariants, and stiffness determine the time integrator. The fully discrete method is the composition of the spatial carrier and time map.

CFL restrictions

Explicit stability requires dt lambda_j to remain inside the integrator’s stability region. Fourier advection gives a scale near N^-1, while polynomial nodal methods may be limited more severely by endpoint clustering. Diffusion and fourth-order terms introduce N^-2 and N^-4 scales in representative settings.

Actual restrictions should be derived from the assembled operator or numerical range, including geometry and element size. Power laws are estimates, not certificates.

Implicit–explicit schemes

IMEX methods split u_t=L u+N(u), treating L implicitly and N explicitly. A first-order step is (I-dtL)u_(n+1)=u_n+dtN(u_n).

The split creates commutator and explicit-nonlinearity debt. An implicit linear part does not remove stiffness in the nonlinear component. High-order IMEX methods require coupled order and stability conditions.

Exponential and integrating-factor methods

Exponential integrators use the exact linear propagator and approximate the nonlinear variation-of-constants integral through phi-functions. Integrating-factor methods define v=exp(-tL)u.

When L is diagonal in Fourier space, these methods are inexpensive. For general matrices, Krylov or rational approximations introduce matrix-function error that must be certified separately.

Operator splitting

Lie splitting is first order; Strang splitting is second order when the operator domains and commutators behave appropriately. Boundary conditions can cause order reduction because suboperators may not preserve the same trace space.

The splitting residue is governed by commutators such as [A,B]. It is not captured by solving each subproblem accurately.

Symplectic and energy-preserving integration

For z_t=J grad H(z), symplectic methods preserve the symplectic form and generally control long-time energy error. Discrete-gradient or average-vector-field methods can preserve an energy exactly.

A geometric time method requires a spatial semidiscretization with a compatible finite Hamiltonian or Poisson structure. Time integration cannot restore structure already destroyed spatially.

Long-time error and invariant preservation

Long-time accuracy depends on phase, modified energy, invariant drift, and weakly damped modes. A method with small local error may accumulate secular phase error.

Validation should include invariant histories, spectra, statistics, and phase rather than only endpoint norms. Backward-error analysis interprets a geometric method as exact evolution of a nearby modified system over long times.

Adaptive temporal and spectral resolution

Temporal controllers estimate local time error; spectral controllers monitor coefficient tails and residuals. Increasing spatial resolution can increase stiffness and force smaller explicit steps, so the controllers must communicate.

Changing resolution requires stable projection or prolongation and preservation of conserved modes. Repeated oscillation of N or dt should be suppressed by hysteresis or predictive control.

Week 12 — Integral, fractional and nonlocal equations

Volterra and Fredholm integral equations

Fredholm equations couple the entire interval: u-lambda integral_a^b K(x,t)u(t)dt=f. Volterra equations use an upper limit x and retain causal triangular structure.

Although both discretize to dense matrices, their solvability mechanisms differ. Fredholm problems can encounter eigenvalue singularities; Volterra problems often admit stable forward continuation.

Spectral Nyström methods

Nyström methods approximate the integral through quadrature and collocate at nodes. For smooth kernels and Gaussian rules, convergence can be spectral.

The quadrature error depends on the target point x. Near diagonal singularities, ordinary smooth quadrature fails. Dense matrix cost may require hierarchical or low-rank acceleration.

Weakly singular kernels

Kernels such as |x-t|^-mu K_0(x,t) with 0<mu<1 are integrable but nonsmooth. Direct polynomial approximation converges slowly. Product integration, singularity subtraction, graded coordinates, or weighted bases separate the singular carrier from the smooth factor.

The singular part should be integrated analytically or with a matched weight. Treating it as smooth hides the dominant residue.

Jacobi methods for endpoint singularities

If u(x)=(1+x)^mu v(x) with smooth v, a weighted Jacobi basis can represent the singular factor exactly and expand only the smooth remainder. Fractional integrals often map one weighted Jacobi family to another.

This is an L-clock mutation from an ordinary polynomial carrier to a singularity-adapted carrier, justified by persistent algebraic coefficient decay.

Fractional derivatives and fractional Sturm–Liouville systems

Riemann–Liouville and Caputo derivatives differ in endpoint treatment. For 0<alpha<1, D_RL^alpha u=d/dx I^(1-alpha)u, while Caputo applies fractional integration to u'.

The boundary carrier and admissible initial data therefore differ. Fractional Sturm–Liouville systems use fractional integration by parts and generalized Jacobi functions. Classical boundary intuition cannot be transferred without retyping the operator.

Nonlocal diffusion

A nonlocal diffusion operator may be L_delta u(x)=integral gamma_delta(x,y)[u(y)-u(x)]dy. The horizon, kernel symmetry, and exterior interaction law define the operator.

Periodic translation-invariant kernels can be diagonalized by Fourier modes. Bounded-domain operators are dense and require fast convolution, low-rank approximation, or hierarchical summation. Truncating exterior interactions changes the model.

Delay differential equations

A delay equation u'(t)=F(u(t),u(t-tau)) has a history segment as its state. The natural carrier is a function space over [t-tau,t], not the instantaneous value.

Spectral collocation approximates the history. State-dependent delays complicate evaluation and stability. History interpolation must match the order of the time method.

Singular quadrature and kernel compression

Product-integration rules interpolate the smooth factor and integrate the singular weight exactly. Logarithmic, algebraic, and principal-value singularities require different formulas.

Away from the diagonal, smooth kernel blocks can be approximated by low rank. The operator-norm compression error must be smaller than the requested solution error after stability amplification. Near-diagonal blocks remain specialized.

Week 13 — Unbounded and semi-infinite domains

Hermite polynomials and functions

Hermite polynomials are orthogonal under a Gaussian weight; Hermite functions form an orthonormal basis of unweighted L2(R). They are natural for Gaussian-localized states and diagonalize the harmonic oscillator.

Differentiation and multiplication by x are tridiagonal through raising and lowering operators. Scaling and translation determine whether the basis matches the physical localization.

Laguerre polynomials and functions

Laguerre systems are defined on [0,infinity) with weight x^alpha exp(-x). Laguerre functions absorb part of the weight and are useful for exponentially decaying solutions.

They are less suitable for algebraic tails or highly oscillatory fields unless combined with scaling, mappings, or domain decomposition.

Rational Chebyshev bases

A map such as x=L(1+y)/(1-y) transforms [-1,1) to [0,infinity). Chebyshev expansion in y becomes a rational approximation in x.

The parameter L controls point distribution. The map introduces variable metric factors and may create reference-coordinate singularities. Its suitability depends on the tail scale.

Mapped spectral methods

For x=g(xi), d/dx=(g')^-1 d/dxi. Higher derivatives involve derivatives of g. The map must be included in the operator, quadrature, norm, and boundary transport.

A good map resolves layers or localization. A poor map produces large metric coefficients and worsens conditioning. Mapping only the nodes without transforming the analysis is invalid.

Scaling and translation parameters

A basis phi_n((x-x_0)/L) uses x_0 to follow location and L to match width. Parameters may be selected through moment estimates, residual minimization, or tail decay.

Changing them over time requires coefficient transport and may introduce basis-motion terms. Adaptive scaling is a carrier mutation, not a cosmetic resampling.

Oscillatory and localized solutions

Highly oscillatory functions require enough degrees of freedom to resolve phase. Standard polynomials may need O(k) modes for frequency k. WKB enrichment, plane waves, or phase-fitted mappings can reduce cost.

Localized pulses require a basis matching both envelope and position. Phase and amplitude errors must be audited separately.

Transparent and absorbing boundary treatments

Truncating an unbounded domain creates artificial reflection. Exact transparent conditions use a nonlocal Dirichlet-to-Neumann operator. Approximate absorbing conditions use asymptotic local operators. PMLs apply complex coordinate stretching.

Reflection error depends on layer thickness, profile smoothness, resolution, and incidence angle. It must be separated from interior discretization error.

Domain truncation versus native unbounded-domain bases

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

The decision balances tail approximation, conditioning, localization, and solver structure. No carrier is universally superior.

Part V — Multiple dimensions and complex geometry

Week 14 — Tensor-product and multidimensional methods

Rectangles, cuboids and periodic boxes

Cartesian product domains support tensor products of one-dimensional bases. Fourier modes are used in periodic directions; polynomial bases in bounded directions. Mixed Fourier–Jacobi carriers are common.

Separable operators inherit product structure, but boundary conditions and variable coefficients can couple directions. The multidimensional carrier is not merely a flattened vector space.

Tensor-product bases

A two-dimensional expansion is u=sum_ij u_ij phi_i(x)psi_j(y). Derivatives act on one index at a time. Treating the coefficients as a tensor preserves operator structure and reduces cost.

Full tensor resolution grows as (N+1)^d, producing dimensionality debt. Anisotropic degrees should reflect directional regularity.

Kronecker structure

Separable discretizations yield matrices such as A_x kron M_y + M_x kron A_y. Matrix-vector products are performed by reshaping coefficients and applying one-dimensional matrices on each side.

Avoiding explicit Kronecker assembly reduces memory and exposes preconditioners. Flattening is an implementation choice, not the native operator form.

Fast diagonalization

Generalized one-dimensional eigendecompositions diagonalize separable elliptic operators. In transformed coordinates, each coefficient is divided by a sum of directional eigenvalues.

Fast diagonalization acts as a direct solver or preconditioner. Variable coefficients and curved geometry create nonseparable residue handled iteratively.

Alternating-direction solvers

ADI alternates one-dimensional implicit solves with selected shifts. Its convergence depends on spectral intervals and shift quality.

Noncommuting operators introduce splitting debt. ADI is particularly effective when directional solves remain banded or diagonal.

Disks, cylinders, balls and spherical shells

Polar and spherical domains require regularity at coordinate singularities. On a disk, angular Fourier modes combine with radial factors behaving like r^|m|. On a ball, spherical harmonics combine with radial Jacobi functions.

Shells avoid the origin and simplify regularity. Full balls require mode-dependent origin constraints.

Spherical harmonics

Spherical harmonics satisfy -Delta_S2 Y_l^m=l(l+1)Y_l^m. They diagonalize the spherical Laplacian. Vector and tensor fields require vector or spin-weighted harmonics to preserve tangent geometry.

Componentwise scalar expansions can violate divergence or tangency constraints. The carrier type must match the field type.

Coordinate singularities

Terms such as 1/r or 1/sin theta are singular in coordinates but may act regularly on smooth physical fields. Compatible parity and mode behavior cancel them.

A naive tensor basis can violate these compatibility conditions and generate spurious modes. Regularity should be tested in intrinsic or Cartesian coordinates.

Sparse grids and low-rank tensor representations

Sparse grids retain index sets with bounded combined degree and work well for mixed regularity. Low-rank tensors approximate coefficient arrays through separated factors.

Both methods introduce truncation or rank debt. Nonlinear operations increase rank and require recompression. Their accuracy depends on anisotropic and separability structure in the solution.

Week 15 — Spectral and hp-element methods

Reference-to-physical element maps

Each physical element is obtained from a reference element by x=F_K(xi). Gradients transform through J^-T, and integrals include |det J|.

The map must remain invertible and sufficiently regular. Geometry approximation contributes directly to operator and quadrature error.

High-order nodal and modal elements

Nodal elements store values at interpolation points; modal elements use hierarchical coefficients. Nodal form simplifies fluxes and nonlinearities; modal form simplifies filtering, adaptivity, and static condensation.

Stable local transforms connect them. A mixed implementation often uses both.

Gauss–Lobatto quadrature

Gauss–Lobatto nodes include endpoints, facilitating element coupling. With N+1 nodes the rule is exact through degree 2N-1.

Collocated quadrature may diagonalize the mass matrix but under-integrate nonlinear and geometric products. This creates aliasing debt.

Conforming and discontinuous formulations

Conforming methods identify interface traces strongly and form a global continuous space. DG methods retain independent traces and impose a numerical flux.

DG stability depends on upwinding or penalty parameters. Conforming methods use fewer interface states but are less flexible for nonmatching meshes. The two formulations use different interface quotients.

Static condensation

Interior and boundary unknowns are partitioned. Eliminating interior variables produces a Schur complement on element interfaces. Interior information is encoded in the Schur operator rather than discarded.

This reduces the global system and permits parallel local factorization. Recovery of interior states remains part of liftback.

Sum factorization

Tensor-product basis evaluation is performed through successive one-dimensional contractions. This reduces cost from naive multidimensional scaling to near O(d N^(d+1)).

Matrix-free sum-factorized kernels are central to high-order performance. Data motion and cache layout often dominate arithmetic.

Curved elements and geometric aliasing

Curved maps produce variable metric and Jacobian factors. Their products with solution and flux polynomials may exceed quadrature exactness. Aliasing can violate conservation and free-stream preservation.

Over-integration, metric projection, and compatible geometric identities control this residue. Geometry is part of the numerical operator.

h-, p- and hp-adaptivity

h-refinement localizes nonsmoothness; p-refinement exploits smoothness. hp methods choose between them using error indicators and modal decay.

Exponential modal decay supports p-refinement. Persistent algebraic decay or localized residuals support h-refinement. The decision is an xSCD clock transition.

Mortar and interface coupling

Mortar spaces connect nonmatching interfaces through weak projection. The interface operator should preserve constants, flux balance, and orientation.

Poor projection can inject or remove energy. Interface error propagates globally and must be included in the ledger.

Parallel and GPU-oriented implementation

High-order methods offer high arithmetic intensity through tensor contractions. Efficient execution uses matrix-free element kernels, batched transforms, contiguous data, and minimized synchronization.

Performance should be reported as time to certified error. Raw degrees of freedom per second can favor an inaccurate method.

Week 16 — Stability on curved and nonlinear systems

Metric identities and discrete geometric conservation laws

Continuous coordinate transformations satisfy metric identities ensuring that constant physical states remain constant. The discrete metric representation must satisfy analogous relations.

A free-stream test is an executable counterkernel: a constant state should produce zero discrete divergence. Moving meshes additionally require consistency between Jacobian evolution and grid velocity.

Summation-by-parts structure

A derivative matrix and mass matrix satisfy H D + D^T H = B. This is the discrete integration-by-parts identity. It transfers continuum energy analysis to the grid.

Boundary and interface terms remain visible through B. SBP structure is an algebraic certificate, not merely a differentiation approximation.

Split forms and entropy stability

Nonlinear product rules fail discretely under finite quadrature. Split forms symmetrize products to recover kinetic-energy or entropy balances. Flux differencing uses two-point entropy-conservative fluxes and adds controlled dissipation.

The semidiscrete entropy inequality is a structural certificate. It depends on symmetry, consistency, and discrete telescoping.

De-aliasing on mapped elements

Nonlinear fluxes, material coefficients, and geometric factors each increase polynomial degree. Over-integration or polynomial projection must address all components.

Dealiasing only the physical flux while leaving metric products under-resolved does not restore the full energy identity. Each residue is audited separately.

Shock detection and spectral viscosity

Modal decay, interelement jumps, or entropy residuals identify loss of smoothness. Stabilization is activated locally through spectral viscosity, filtering, or artificial diffusion.

The added dissipation is a controlled model modification. Its spatial profile, modal range, and scaling determine whether accuracy survives in smooth regions.

Positivity and realizability

High-order polynomials can violate positivity between nodes. Limiters rescale the high-order correction toward a positive average. For systems, admissibility may require positive density and pressure or a realizable moment vector.

Positivity is part of the state identity relation. A small norm error does not make an inadmissible state physically equivalent to an admissible one.

Incompressible Navier–Stokes projection methods

Projection methods compute an intermediate velocity and solve a pressure Poisson problem to remove divergence. Boundary conditions for the correction influence accuracy and may create splitting layers.

The certificate includes post-projection divergence, boundary consistency, pressure gauge, and energy balance. Nonlinear de-aliasing remains essential.

Transition, turbulence and under-resolved computation

Under-resolved turbulence cannot be certified by pointwise convergence. Relevant observers include spectra, transfer rates, dissipation, coherent structures, and statistics.

Stabilization, subgrid modeling, and numerical dissipation must be distinguished. A visually plausible field can still have incorrect energy transport.

Part VI — Modern extensions and verification

Week 17 — Scientific machine learning and spectral operators

Solution operators between function spaces

A solution operator maps input functions to output functions, for example G:a(.) -> u(.). Training samples are finite representations of this map. The operator itself is defined only after input and output spaces, norms, boundary conditions, and parameter ranges are declared.

Discretization independence requires consistency under refinement. A learned array-to-array map on one grid is not automatically a continuum operator.

Fourier and spectral neural operators

A Fourier neural layer applies learned modewise matrices between FFT transforms and combines them with local mixing and nonlinear activation. Polynomial or spherical variants replace Fourier modes with other spectral bases.

Nonlinear activations generate unresolved frequencies and can alias. Boundary conditions and geometry require explicit architectural or residual treatment. Learned multipliers should not be confused with exact PDE diagonalization.

Learning in coefficient space

Inputs and outputs are expanded in a declared basis and the network learns coefficient maps. This separates representation error from regression error and permits mode-dependent weighting.

Coefficient normalization, truncation, parity, reality conditions, and basis metadata must travel with the model. An untyped coefficient vector has no stable meaning.

Spectral residual losses

A predicted solution gives residual R=L u_theta-f. A coefficient-space loss is sum_k w_k|R_hat_k|². The weights define the norm; derivative-weighted choices approximate Sobolev residuals.

A small residual implies small error only when the inverse operator is stable in the corresponding spaces. Boundary and initial residuals remain separate.

Parseval-based training objectives

In an orthonormal basis, ||u-u_ref||²=sum_k|u_k-u_ref,k|². Sobolev losses weight modes by powers of the operator eigenvalues.

Relative mode weighting can overemphasize tiny noisy coefficients. The loss should reflect the physical observable and certification norm.

Aliasing and resolution-transfer failure

Spectral neural networks perform nonlinear operations that generate unresolved modes. Without padding or projection, aliasing becomes part of the learned map.

A model may exploit a grid-specific aliasing pattern and fail at another resolution. A commuting refinement test compares projection of the continuum prediction with prediction by the refined model.

Hybrid solver–learner architectures

Learned components can supply initial guesses, preconditioners, closures, or coarse corrections inside a residual-controlled classical solver. The outer method supplies convergence tests and rejection.

This is safer than replacing the entire solver because the learned module remains subordinate to a certified correction loop.

Training-data and discretization dependence

Training data inherit the errors of the generating solver, mesh, quadrature, and physical model. A network can learn those artifacts.

Cross-discretization testing and independently generated references reveal whether the model approximates the continuum operator or one numerical implementation.

Comparison with classical spectral solvers at equal accuracy

A fair comparison fixes a target residual and error. Total cost includes data generation, training, inference, correction, setup, and repeated solves. The relevant quantity is time to certified accuracy.

Robustness, conservation, resolution transfer, and out-of-distribution behavior must be compared alongside runtime.

Stability, conservation and posterior certification

Architectural constraints can preserve symmetry, mass, positivity, or contractivity. They do not by themselves certify the complete model.

After inference, residual, boundary defect, invariant error, and uncertainty determine whether the prediction is accepted, corrected by a classical solve, or rejected.

Week 18 — Error certification and computational audit

A priori versus a posteriori error

A priori estimates predict approximation from regularity and resolution. A posteriori estimates use the computed residual, jumps, boundary defect, or coefficient tail.

A priori results guide carrier selection. A posteriori results guide adaptation and acceptance. Both require meaningful constants and norms.

Residual, truncation and quadrature decomposition

The total defect can be written as data approximation plus operator/quadrature error plus algebraic residual. Additional components include aliasing, geometry, time integration, and modeling.

Terms expressed in different norms must be propagated through stability estimates rather than added blindly. The ledger records both magnitude and carrier.

Roundoff and conditioning

Backward error measures the perturbation for which the computed answer is exact. Forward error is backward error amplified by conditioning. Spectral differentiation and high-degree endpoint evaluation can magnify roundoff.

Increasing N beyond the approximation–roundoff crossover can worsen accuracy. Scaling, orthonormal bases, preconditioning, and higher precision move the floor.

Manufactured solutions

Choose u_exact, derive f=L u_exact, and generate compatible boundary data. This tests assembly and convergence.

Manufactured tests should include variable coefficients, mixed conditions, nonlinear terms, curved geometry, analytic solutions, and finite-regularity cases. A polynomial exactly representable by the basis is useful for debugging but not for convergence validation.

Grid and polynomial-order refinement

h-refinement rates are estimated from successive mesh errors. Spectral convergence is diagnosed by semilog or log-log plots versus N, depending on geometric or algebraic decay.

Refinement continues until the asymptotic regime and error floor are visible. All coupled parameters, including time step and quadrature order, must be controlled.

Cross-formulation verification

Solve the same problem using independent formulations, such as Galerkin and collocation or direct and mixed methods. Agreement within separate error estimates is stronger than refinement of one code path.

Shared quadrature or boundary routines can create common-mode errors, so genuine independence matters.

Conservation and stability audits

For a continuous balance dE/dt=-D+B, compute the discrete defect dE_N/dt+D_N-B_N. Mass, momentum, energy, entropy, and fluxes are tracked according to the model.

Transient growth and long-time drift are included. Final boundedness alone is not a stability audit.

Reproducibility and benchmark design

A benchmark fixes equations, parameters, geometry, reference values, norms, tolerances, precision, and software environment. Reference solutions require independent validation.

Benchmarks should include smooth, singular, stiff, non-normal, nonlinear, and geometry-sensitive regimes rather than only favorable analytic cases.

Cost versus accuracy

Performance is measured as error versus wall time, memory, and energy. High-order methods have larger per-degree cost but may dominate at stringent tolerances.

Setup, transforms, factorization reuse, communication, and number of right-hand sides affect the break-even point. Complexity without certified error is incomplete.

When spectral convergence fails

Failure sources include nonsmoothness, endpoint singularities, layers, shocks, geometry irregularity, aliasing, inconsistent boundaries, conditioning, and roundoff. The response is diagnosis followed by carrier mutation.

Possible mutations are domain decomposition, singularity factoring, coordinate mapping, weighted bases, over-integration, filtering, viscosity, preconditioning, or increased precision. Blindly increasing N is prohibited.

Week 19 — Final projects and research presentations

Navier–Stokes or magnetohydrodynamics on periodic or curved domains

The project must specify divergence enforcement, pressure gauge, nonlinear form, dealiasing, metric identities, time integration, and energy or magnetic-energy diagnostics. Curved geometry requires a free-stream test. MHD additionally requires control of div B=0 and magnetic flux.

High-order Schrödinger eigenvalue computations

The project formulates the operator domain, potential singularity, boundary conditions, and eigenpair certificate. Shift-invert or contour methods may target interior eigenvalues. Coulomb singularities require basis adaptation or domain decomposition.

Residual, backward error, spectral convergence, and comparison with independent references are mandatory.

Fractional or nonlocal PDEs

The project declares the fractional derivative or kernel, exterior condition, singularity structure, quadrature, and fast application. Weighted Jacobi or Fourier carriers should be justified from the operator.

Validation includes asymptotic behavior, refinement, and independent quadrature or transform comparison.

Spectral-element wave propagation

The project examines dispersion, numerical fluxes, CFL restrictions, curved elements, PML or absorbing boundaries, and matrix-free performance. Phase error over many wavelengths is more important than a single short-time norm.

Non-normal operator pseudospectra

The project computes eigenvalues, resolvent norms, pseudospectral contours, and transient growth for a spectrally discretized operator. Resolution continuation distinguishes continuum features from matrix artifacts.

Adaptive ultraspherical solvers

The project implements sparse differentiation, conversion, boundary rows, adaptive QR, residual evaluation, and automatic resolution selection. Smooth, layered, and singular examples test mutation logic.

Fast Jacobi connection transforms

The project compares direct and fast transforms, measures scaling and round-trip error, and demonstrates an operator whose sparse implementation depends on parameter transport. Endpoint and high-degree stability are central.

Spectral methods for kinetic equations

The project addresses phase-space dimension, velocity truncation, positivity, conservation, filamentation, and tensor compression. Vlasov–Poisson, Fokker–Planck, or reduced Boltzmann models are suitable carriers.

Structure-preserving Cahn–Hilliard or Hamiltonian systems

For Cahn–Hilliard, mass conservation and energy decay are audited. For Hamiltonian systems, symplectic or Poisson structure and long-time invariants are measured. Equal short-time accuracy comparisons expose structural differences.

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

The project builds both methods for one parameterized PDE, generates independently certified references, and compares total cost, break-even query count, residual, boundary error, conservation, resolution transfer, and failure rate. The neural output must be subjected to posterior correction or rejection.

Every final project ends with an xSCD manifest containing the continuum carrier, primitive failure, numerical carrier, transport maps, debt ledger, counterkernel tests, certificates, liftback statement, mutation triggers, and terminal status. A project is locally certified only when the claimed continuum result follows from measured residuals and a declared stability mechanism.