Source-linked AI summary
Dedalus: A Flexible Framework for Numerical Simulations with Spectral Methods
Keaton J. Burns, Geoffrey M. Vasil, Jeffrey S. Oishi, Daniel Lecoanet, Benjamin P. Brown
TL;DR
Many PDEs require numerical discretization, while conventional methods can be difficult to implement for complex equations and converge relatively slowly. Dedalus addresses this gap with a symbolic, parallel spectral framework that produces banded systems and demonstrates efficient scaling and adaptability across physical models.
Problem
Many PDEs require numerical computation, but conventional discretizations can be difficult to implement for complex equations and relatively slow to converge, limiting straightforward simulation of non-standard equation sets.
Method
Dedalus lets users symbolically specify custom PDEs on products of spectral series, then uses a first-order tau formulation, basis choices, and preconditioning to generate banded systems that are automatically parallelized.
Results
Dedalus scales efficiently to thousands of cores, typically maintains over 50% parallel efficiency down to 8 pencils per core, and demonstrates adaptability across diverse physical models.
Takeaways & Limitations
Dedalus provides an open-source Python framework for constructing, solving, and analyzing custom PDE problems from individual laptops to supercomputers with tens-of-thousands of cores.
Abstract
from arXiv · showhide
Numerical solutions of partial differential equations enable a broad range of scientific research. The Dedalus Project is a flexible, open-source, parallelized computational framework for solving general partial differential equations using spectral methods. Dedalus translates plain-text strings describing partial differential equations into efficient solvers. This paper details the numerical method that enables this translation, describes the design and implementation of the codebase, and illustrates its capabilities with a variety of example problems. The numerical method is a first-order generalized tau formulation that discretizes equations into banded matrices. This method is implemented with an object-oriented design. Classes for spectral bases and domains manage the discretization and automatic parallel distribution of variables. Discretized fields and mathematical operators are symbolically manipulated with a basic computer algebra system. Initial value, boundary value, and eigenvalue problems are efficiently solved using high-performance linear algebra, transform, and parallel communication libraries. Custom analysis outputs can also be specified in plain text and stored in self-describing portable formats. The performance of the code is evaluated with a parallel scaling benchmark and a comparison to a finite-volume code. The features and flexibility of the codebase are illustrated by solving several examples: the nonlinear Schrodinger equation on a graph, a supersonic magnetohydrodynamic vortex, quasigeostrophic flow, Stokes flow in a cylindrical annulus, normal modes of a radiative atmosphere, and diamagnetic levitation. The Dedalus code and the example problems are available online at http://dedalus-project.org/.
I. INTRODUCTION
PDE simulation needs numerical methods that balance accuracy, geometric flexibility, and implementation cost. Dedalus addresses this by combining sparse spectral discretization with symbolic equation specification and automatic parallelization.
- Motivation: Most PDEs require discretization into finite algebraic systems, while finite-element, finite-volume, and finite-difference methods can be difficult to implement for complex equations and converge relatively slowly.These methods accommodate complex geometries but derive relations over cells or points.
- Motivation: Spectral methods expand variables in basis functions, yielding rapidly converging solutions for smooth functions but typically restricting applications to simple geometries.Chebyshev and Fourier bases can provide exponential convergence when the represented functions are smooth.
- Motivation: Many optimized simulation packages hardcode equations, variables, controls, integration schemes, and analysis outputs, limiting their flexibility for new models.This specialization can leave novel equations underserved when development costs are difficult to justify.
- Dedalus Project: Dedalus bridges flexible model specification and parallel spectral simulation by letting users enter custom PDE systems symbolically and automatically producing sparse, parallelized discretizations.The project is open-source, modular, and designed for a broad range of PDEs.
- Spectral methods: Spectral representations truncate basis expansions after N modes, with smooth-function coefficients typically decaying exponentially and determining highly accurate approximations.Quadrature rules provide transforms between spectral coefficients and sampled function values.
- Dedalus Project: Dedalus’s general sparse tau approach uses different test and trial bases to represent differential operators with sparse and banded matrices for general systems and boundary conditions.The approach is intended to accommodate general equations and boundary conditions automatically.
2. Banded boundary conditions
Dedalus converts general PDE systems into sparse, fully banded matrices by combining a first-order tau formulation, basis recombination, and variable-mode interleaving. This structure supports efficient solution of coupled systems while accommodating boundary conditions and smooth non-constant coefficients.
- 2. Banded boundary conditions: Dirichlet preconditioning makes boundary conditions depend only on the first two expansion coefficients.The method is also called basis recombination.
- 2. Banded boundary conditions: Choosing φ_n = D_n, ψ_n = U_n, and P = U_N−1 produces fully banded matrices for first-order systems with Dirichlet conditions.Sparse/banded algorithms can then solve the resulting systems efficiently.
- 2. Banded boundary conditions: First-order reformulation expresses Neumann and global integral conditions as Dirichlet conditions on first-order variables.This allows the same banded formulation to handle broader boundary-condition types and other orthogonal polynomial bases.
- 2. Banded boundary conditions: Smooth non-constant coefficients retain low-bandwidth multiplication matrices when their spectral expansions are truncated to retained terms.For smooth coefficients, the retained-term count N′ is much smaller than N.
- 2. Banded boundary conditions: Interleaving variables and coefficients changes coupled-system bandwidth from O(SN) to O(S), enabling practical solves at large mode counts.The system matrix acts on interleaved variable coefficients after reversing the grouping order.
- 2. Banded boundary conditions: The Poisson construction becomes highly sparse and completely banded after T-to-U conversion, boundary treatment, Dirichlet recombination, and mode grouping.The stages progress from dense upper-triangular differentiation matrices to sparse boundary rows and a final banded matrix.
B. Dependencies
Dedalus is implemented as an open-source Python package built around extensible spectral-basis objects and optimized numerical libraries. Its generalized design keeps the codebase compact while supporting transforms, parallel computation, and reproducible use.
- B. Dependencies: Dedalus is an open-source Python3 package that relies on optimized C libraries for much of the computation.Primary dependencies include NumPy, SciPy, Cython, FFTW, MPI/mpi4py, and HDF5/h5py.
- B. Dependencies: Generalized algorithms let Dedalus provide broad functionality in roughly 10,000 lines of Python.The paper reports that accommodating many equations and domains produced a compact, maintainable codebase.
- B. Dependencies: Dedalus supports public development, tutorials, example problems, package releases, and installation channels for single machines and clusters.The project maintains documentation, source repositories, and community support infrastructure.
- B. Dependencies: Basis objects define one-dimensional spectral series, coordinate mappings, resolutions, dealiasing, and forward and backward transforms.Transform plans or matrices can be cached by basis objects, and multiple libraries can share a unified interface.
- B. Dependencies: Fourier bases use evenly spaced grids on periodic intervals, while Chebyshev and Legendre bases use finite-interval polynomial grids.Fourier coefficients are computed with FFTs; real-valued Fourier fields store only nonnegative wavenumbers using conjugate symmetry.
B. Sine/Cosine basis
Dedalus supplies sine/cosine and polynomial spectral bases for periodic, finite, and unbounded domains. Their grids and transforms are matched to symmetry, endpoint behavior, and the geometry of the computational interval.
- B. Sine/Cosine basis: SinCos bases represent endpoint-symmetric periodic functions with either cosine or sine series on [0, π].Even-parity functions use cosine series, while odd-parity functions use sine series.
- B. Sine/Cosine basis: SinCos transforms use type-II forward and type-III backward DCTs or DSTs on evenly spaced interior grids.The transforms are implemented through both SciPy and FFTW and rescaled so coefficients represent mode amplitudes.
- B. Sine/Cosine basis: Chebyshev bases use Gauss-Chebyshev interior nodes that cluster quadratically near finite-interval endpoints.Endpoint clustering permits resolution of small structures.
- B. Sine/Cosine basis: Chebyshev rational functions can cover half-lines and the whole real line, but Dedalus requires users to include the coordinate changes manually.These functions are not implemented explicitly in the package.
- B. Sine/Cosine basis: Legendre bases use Legendre polynomials and Gauss-Legendre quadrature nodes, with quadrature-based matrix-multiply transforms.The nodes and weights are computed with scipy.special.roots_legendre.
- B. Sine/Cosine basis: Hermite bases represent functions on the whole real line using Hermite polynomials or enveloped Hermite functions.Their affine mapping is specified by center and stretch parameters, and grids use Gauss-Hermite nodes.
F. Laguerre basis
Dedalus combines Laguerre bases for half-line problems with compound polynomial bases, scalable transforms, dealiasing, and distributed-memory domain layouts. These components support localized resolution and parallel spectral computations.
- F. Laguerre basis: Laguerre bases represent problems on (0, ∞) with standard or enveloped Laguerre functions.Their mapping uses edge and stretch parameters, and negative stretch can span the negative half-line.
- F. Laguerre basis: Compound bases concatenate adjacent polynomial segments with independent resolutions and lengths.Interior grids avoid overlapping interface points, and continuity need not be imposed a priori.
- F. Laguerre basis: Compound expansions can reduce the modes needed for nonsmooth functions when nondifferentiable locations are known.They also place higher resolution at fixed interior locations through endpoint clustering on subsegments.
- F. Laguerre basis: Transforms map N_c coefficients to grids of size N_g = sN_c by truncating coefficients for s < 1 or zero-padding for s > 1.Padding supports spectral interpolation and dealiasing in grid-space nonlinear operations.
- F. Laguerre basis: The 3/2 rule requires s ≥ 3/2 to properly dealias quadratic nonlinearities evaluated on the grid.Dedalus sets each basis’s dealias scale when the basis is instantiated.
- F. Laguerre basis: Domains distribute multidimensional fields across MPI processes, while Distributor and Layout objects coordinate transforms and global data transpositions.The coefficient-to-grid sequence uses D local transforms and P distributed array transpositions.
B. Transpose routines
Dedalus distributes multidimensional spectral data through ordered layouts, using local transforms and distributed transpositions to move between coefficient and grid representations. Its field and operator abstractions hide these layout details while supporting efficient parallel manipulation and linear algebra.
- Transpose routines: Distributed transpositions switch which array dimensions are distributed over each process-mesh axis through subgroup communication and data reordering.MPI implementations use subgroup all-to-all calls, while FFTW provides an optimized distributed-transpose interface.
- Transpose routines: Multiple arrays can be concatenated into one planned transposition, reducing latency from initiating separate communications.The grouped operation expands the reduced subgroup shape along its first dimension by the number of arrays.
- Transpose routines: The layout system constructs ordered transform/distribution states for arbitrary domain and process-mesh dimensions, while accounting for empty-process edge cases.Layout objects expose global and local shapes, block sizes, coordinates, and slices for distributed fields.
- Transpose routines: Field methods move data sequentially toward grid or coefficient space, or directly to a requested layout, without requiring users to manage intermediate transforms and transpositions.Local grid and coefficient arrays can be initialized robustly across changing parallelization scenarios.
- Transpose routines: Operator trees represent composed field expressions, supporting deferred repeated evaluation and matrix construction for linear operations.These capabilities provide the basis for implicitly or explicitly evaluating arbitrary expressions in general PDE solvers.
D. Linear spectral operators
Dedalus represents linear spectral operations as basis-specific operators that can be evaluated on fields or converted into diagonal or sparse matrix forms. Its operator system covers differentiation, integration, interpolation, Hilbert transforms, and user-defined functions.
- D. Linear spectral operators: Linear operators provide matrix forms for spectral coefficients: separable operators are diagonal, whereas coupled operators use sparse matrices.Operators generally act on coefficient-space data, with coupled operators requiring the relevant axis to be local.
- 1. Differentiation: Differentiation preserves sparsity by using diagonal Fourier/SinCos operators, naturally banded Hermite/Laguerre operators, or test-basis conversions for polynomial bases.Chebyshev, Legendre, and Laguerre forward differentiation uses O(Nc) recurrence relations rather than dense matrices.
- D. Linear spectral operators: Integration returns constants and marks integrated axes as constant, while indefinite integration solves a simple boundary value problem to fix the integration constant.Basis-specific integration vectors are rescaled for affine coordinate maps.
- D. Linear spectral operators: Interpolation constructs basis-function evaluation vectors at specified problem-coordinate positions, including left, center, and right endpoints.Compound bases use the sub-basis containing the interpolation position, choosing the first at an interface.
- D. Linear spectral operators: Hilbert transforms are implemented as separable operators for Fourier and SinCos bases, with factory methods for mixed and higher-order transforms.Constant operands along the transformed axis return zero without creating an operator.
- D. Linear spectral operators: General Python functions can be wrapped as deferred operators, and operator expressions support recursive algebraic transformations such as expansion, splitting, and replacement.These operations enable symbolic manipulation of expressions before evaluation or matrix construction.
VIII. PROBLEMS
Dedalus problem classes turn plain-text equations and boundary conditions into structured operator systems for initial-value, boundary-value, and eigenvalue problems. The parser validates linearity, derivative order, metadata compatibility, and problem-specific constraints before constructing solver representations.
- A. Problem creation: Problem classes support LBVP, NLBVP, EVP, and IVP formulations, parsing plain-text equations into sparse linear operators and explicit nonlinear operator trees.Variables, domains, parameters, substitutions, and boundary conditions define each problem instance.
- A. Problem creation: Problem creation requires first-order systems, with auxiliary variables added when needed; polynomial bases are restricted to at most one and must occupy the last axis.IVPs use an optional temporal-variable name that defaults to 't'.
- D. Substitutions: Substitutions provide composable string-replacement rules and lambda-like mathematical functions that simplify complex equation entry.They can define aliases for computed quantities or reusable expressions such as kinetic energy and Laplacians.
- A. Problem creation: All problem types require compatible metadata, a nonzero linear LHS, first-order coupled derivatives, and boundary conditions constant along the last axis.These checks are applied after equations are split into LHS and RHS operator expressions.
- VIII. PROBLEMS: Linear boundary-value problems require RHS expressions independent of problem variables, allowing inhomogeneous forcing while the LHS remains linear.The linear expressions are expanded and converted into canonical linear form.
- 2. Nonlinear boundary value problems: Nonlinear boundary-value problems construct the Frechet differential of the nonlinear RHS, but Dedalus supports them only in one dimension because variable-dependent coefficients generally couple horizontal modes.The differential is constructed symbolically from perturbations of the problem variables.
4. Initial value problems
Dedalus discretizes initial and boundary-value problems into preconditioned, banded pencil systems and solves them independently where possible. Nonlinear problems use Newton iterations whose convergence depends on the starting state.
- Problem formulation: Initial-value problems require time-independent, first-order-in-time LHS terms and a RHS without time derivatives.They are represented as M · ∂tX + L · X = F(X, t).
- Matrix construction: Dedalus constructs operator matrices per pencil, applying basis and Dirichlet preconditioning to render derivatives and boundary conditions banded.Boundary and interface rows enforce external conditions, continuity, and constant-direction constraints before matrices are joined.
- Matrix construction: Processed pencil matrices are interleaved by mode, with bandwidth set by the number of coupled variables and the largest sub-block bandwidth.The resulting sparse matrices are stored using SciPy sparse representations.
- Linear solves: Linear boundary-value solves reduce to independent sparse, banded matrix problems for each pencil and recover the state through right preconditioning.Parallel processes solve local pencils, with cost scaling as O(NcD) in the stated formulation.
- Nonlinear solves: Nonlinear boundary-value solves use Newton updates obtained by solving a Frechet-derivative system around the current state.The matrix changes each iteration, and convergence can depend sensitively on the initial state; sufficiently good starting points yield quadratic convergence.
D. Eigenvalue solver
The solver framework supports generalized eigenvalue problems alongside time-dependent integration and configurable stopping, timestep, and output controls. Dense solves compute all modes but become impractical at large resolutions, whereas sparse shift-and-invert targets selected modes.
- Eigenvalue solver: Eigenvalue solvers construct preconditioned M and L matrices and solve one pencil at a time, storing eigenvalues and eigenvectors.Computed eigenmodes can be loaded into the solver state for visualization or further computation.
- Eigenvalue solver: The dense eigenvalue method solves the full generalized problem, but its dense memory and computational cost become prohibitive at large resolutions.It uses SciPy’s dense eigensolver after converting the matrices to dense arrays.
- Eigenvalue solver: The sparse method uses shift-and-invert with ARPACK to compute selected eigenmodes near a target eigenvalue, assuming the shifted operator is full rank.A sparse linear operator applies M and solves against L + σT M.
- Initial-value integration: Initial-value integration uses mixed implicit-explicit schemes, treating LHS terms implicitly and RHS terms explicitly while allowing user-controlled timestepping.Different pencils are coupled through the RHS terms, and the mass matrix may be rank deficient because of constraints and boundary conditions.
- Timestep control: Adaptive CFL timesteps are based on the maximum summed grid frequency, then modified by safety, absolute, relative-change, and threshold limits.These controls can avoid excessive overstepping, instability from abrupt timestep changes, and unnecessary matrix factorizations.
2. Runge-Kutta IMEX integrators
Dedalus implements reusable IMEX Runge-Kutta infrastructure that advances stiff linear terms implicitly and nonlinear terms explicitly. Its integrator hierarchy supports adaptive timesteps, factorization reuse, and straightforward switching or extension.
- Runge-Kutta IMEX integrators: Runge-Kutta IMEX schemes construct stages from preconditioned mass, linear, and RHS operators, then solve one system per stage.Each stage evaluates the RHS, forms the stage system, applies preconditioning, and updates the state.
- Runge-Kutta IMEX integrators: Globally stiffly accurate methods identify the final stage with the advanced solution and do not require a full-rank mass matrix.This accommodates algebraic constraints and boundary conditions common in Dedalus problems.
- Runge-Kutta IMEX integrators: Stage updates solve preconditioned systems for each pencil and recover the state with right preconditioning.If the timestep is unchanged, stored matrix factorizations can reduce solve time.
- Integrator interface: Dedalus provides first-, second-, and third-order IMEX methods and lets users switch integrators or implement new schemes through shared base classes.Runge-Kutta methods do not depend on previous state iterations, supporting full-order initial steps and adaptive timestepping.
- Linear algebra: A generic matrix-solver interface wraps direct sparse, banded, and block-diagonal routines and supports reusable factorizations.The interface is extensible, and solver choices can be specified per Solver object or through configuration.
C. Post-processing
Dedalus couples symbolic analysis tasks with parallel, self-describing output and post-processing tools. Its benchmarks and examples demonstrate scalable execution, competitive accuracy, and applicability across diverse PDE settings.
- C. Post-processing: Distributed process files can be merged into global spatial files or combined space-time datasets, although merging sets can produce very large files.High-level plotting and parallel visualization helpers simplify interaction with simulation products.
- Demonstration problems: Dedalus examples span graph-based nonlinear Schrödinger dynamics, compressible MHD, quasigeostrophic flow, cylindrical Stokes flow, atmospheric waves, and diamagnetic levitation.The examples exercise networks, shocks, polar coordinates, eigenmodes, immersed boundaries, and nonlocal conditions.
- Performance: Parallel efficiency typically remains above 50% down to 8 pencils per core, while weak scaling follows the expected 1/log N behavior for FFT-based computations.The benchmark covers 32–2048 cores and reports linear execution-time scaling with problem variables.
- Accuracy benchmark: 163842 finite-volume cells were required to match Dedalus accuracy at 20482 modes in the Kelvin–Helmholtz benchmark.At low-to-moderate resolution, finite-volume errors triggered spurious secondary instabilities that disappeared at increased resolution.
- Accuracy benchmark: For smooth, low-Mach-number flows in simple geometries, rapid spectral convergence can outweigh the higher per-iteration cost of spectral methods.Dedalus also follows shock interactions in high-Mach-number MHD when diffusion regularizes the shocks.
E. Quasigeostrophic Flow
The QG example uses Dedalus to simulate a general finite three-dimensional layer with boundaries and non-constant coefficients, while illustrating a physically motivated first-order formulation. The resulting flow becomes statistically stationary, dominated by compact eddies with nontrivial depth-dependent potential vorticity.
- Model scope: QG is used as a demanding example because the general problem requires a three-dimensional layer with boundaries and non-constant coefficients.The formulation can also be simplified to two-dimensional or triply periodic cases.
- Model formulation: Using upwelling simplifies the boundary conditions without increasing computational cost relative to the traditional potential-vorticity formulation.The first-order formulation avoids solving a nonlinear buoyancy-advection equation on the boundary.
- Model formulation: The alternative formulation solves for upwelling w and pressure p, with buoyancy and velocity defined diagnostically from them.Dedalus substitutions encode the geostrophic and hydrostatic diagnostic balances.
- Results: Figure 13 shows surface and lateral slices of potential vorticity and buoyancy perturbation at t = 200 in the statistically saturated state.Potential vorticity appears in the top images and buoyancy perturbation in the bottom images.
- Simulation setup: The simulation uses a linearly unstable thermal-wind profile with β = 0.1, ν4 = κ4 = 10^-6, γ− = 0.16, and 256 × 128 × 32 spectral modes.The horizontal domain is 40 × 20 Rossby radii, with 3/2 dealiasing and no top friction.
- Results: After roughly 75 dynamical time units, the system reaches statistical stationarity and develops compact eddies that form, merge, and break up.The flow is mostly two-dimensional but retains nontrivial depth-dependent potential-vorticity variations.
F. Stokes Flow
The Stokes-flow example models reversible tracer evolution in a cylindrical annulus using sparse spectral discretization and reveals how shear dispersion and diffusion shape the final dye patterns.
- Stokes Flow: At low Reynolds number, Taylor-Couette flow reverses after the inner cylinder rotates forward and backward, causing dye tracers to mix and then unmix.The Peclet number controls tracer mass diffusion.
- Stokes Flow: A Fourier–Chebyshev direct product accurately discretizes the annulus, while multiplying equations by r or r^2 keeps non-constant coefficient matrices banded.This avoids the slowly converging Chebyshev expansion of r^-2.
- Stokes Flow: The simulation uses 512×512 modes and runs approximately four rotations forward and backward, with nrot = 8.It uses the Runge-Kutta 443 timestepper with fixed dt = 0.005.
- Stokes Flow: For diffusion times longer than the simulation, the required Peclet number is Pe ≳ 4 × 10^5; with Pe = ∞, the sheets return to the original Gaussian apart from slight timestep errors.This establishes a diffusion-controlled boundary between reversible-looking and visibly altered tracer evolution.
- Additional examples: The example demonstrates Dedalus’s ability to connect complex atmospheric structure with detailed eigenvalue and eigenfunction solves, including gravity-mode trapping and acoustic-mode penetration.The same framework is illustrated across radiative atmospheres and diamagnetic levitation.
XII. CONCLUSION & OUTLOOK
Dedalus is presented as a public, open-source Python framework for accurate, reproducible PDE simulations using spectral methods. Its symbolic specification, custom domains, automated parallelization, and broad examples support diverse physical models and future community-driven extensions.
- Dedalus is introduced as a public Python framework for solving PDEs with spectral methods and an interdisciplinary user and developer community.
- Users can construct product spectral domains, symbolically specify equations and boundary conditions, and perform custom data analysis.
- Dedalus supports initial-value, eigenvalue, and linear and nonlinear boundary-value problems with MPI-based automatic parallelization.
- A diverse set of examples demonstrates the code’s adaptability to many different physical models.
- Planned extensions include coupled dimensions, coordinate-free equation entry, enhanced analysis, and non-direct-product bases for curvilinear tensorial quantities.
Appendix A: Spectral operator matrices
The appendix describes sparse spectral operator matrices for Fourier, trigonometric, Chebyshev, Legendre, Hermite, and Laguerre bases. Recurrences and basis conversions produce sparse or banded differentiation operators suitable for preconditioning and efficient computation.
- Fourier and trigonometric bases: Fourier differentiation is represented as a separable matrix indexed by each mode’s signed wavenumber.
- Fourier and trigonometric bases: Sine/cosine differentiation is separable and flips the parity of its operand.
- Polynomial bases: Chebyshev and Legendre differentiation uses recurrence relations and sparse conversions to construct operators for left-preconditioning differential equations.
- Hermite bases: Hermite polynomial and enveloped-Hermite differentiation matrices are banded, so conversion operators or left-preconditioners are unnecessary.
- Laguerre bases: Laguerre differentiation matrices become sparse through conversion to first generalized Laguerre polynomials, enabling left-preconditioning of Laguerre differential equations.
2. Dirichlet recombination matrices
Dedalus uses basis recombinations and sparse conversion matrices to localize Dirichlet boundary support while preserving sparse equation matrices. The appendix also describes band-limited multiplication operators and current limitations for enveloped-function products.
- Dirichlet recombination: Chebyshev and Legendre endpoint conditions allow recombinations that isolate Dirichlet boundary support to the first two modes.
- Dirichlet recombination: For Chebyshev and Legendre bases, sparse conversion matrices right-precondition boundary rows while maintaining equation-matrix sparsity.
- Dirichlet recombination: Laguerre endpoint conditions allow recombination that isolates Dirichlet boundary support to the first mode.
- Dirichlet recombination: A sparse Laguerre conversion matrix right-preconditions polynomial and enveloped-function problems while compressing Dirichlet boundary rows.
- Multiplication operators: Chebyshev, Legendre, Hermite, and Laguerre multiplication matrices are constructed from polynomial identities or recurrence relations for band-limited spectral operators.
- Multiplication operators: Multiplication between enveloped Hermite or Laguerre functions is not implemented because it is not band-limited in the corresponding enveloped bases.