Dynamic Mode Decomposition (DMD) is a popular data-driven analysis technique used to decompose complex, nonlinear systems into a set of modes, revealing underlying patterns and dynamics through spectral analysis. This review presents a comprehensive and pedagogical examination of DMD, emphasizing the role of Koopman operators in transforming complex nonlinear dynamics into a linear framework. A distinctive feature of this review is its focus on the relationship between DMD and the spectral properties of Koopman operators, with particular emphasis on the theory and practice of DMD algorithms for spectral computations. We explore the diverse “multiverse” of DMD methods, categorized into three main areas: linear regression-based methods, Galerkin approximations, and structure-preserving techniques. Each category is studied for its unique contributions and challenges, providing a detailed overview of significant algorithms and their applications as outlined in Table 1. We include a MATLAB® package with examples and applications to enhance the practical understanding of these methods. This review serves as both a practical guide and a theoretical reference for various DMD methods, accessible to both experts and newcomers, and enabling readers to explore their areas of interest in the expansive field of DMD.

Koopman operators linearize nonlinear dynamical systems, making their spectral information of crucial interest. Numerous algorithms have been developed to approximate these spectral properties, and Dynamic Mode Decomposition (DMD) stands out as the poster child of projection-based methods. Although the Koopman operator itself is linear, the fact that it acts in an infinite-dimensional space of observables poses challenges. These include spurious modes, essential spectra, and the verification of Koopman mode decompositions. While recent work has addressed these challenges for deterministic systems, there remains a notable gap in verified DMD methods for stochastic systems, where the Koopman operator measures the expectation of observables. We show that it is necessary to go beyond expectations to address these issues. By incorporating variance into the Koopman framework, we address these challenges. Through an additional DMD-type matrix, we approximate the sum of a squared residual and a variance term, each of which can be approximated individually using batched snapshot data. This allows verified computation of the spectral properties of stochastic Koopman operators, controlling the projection error. We also introduce the concept of variance-pseudospectra to gauge statistical coherency. Finally, we present a suite of convergence results for the spectral information of stochastic Koopman operators. Our study concludes with practical applications using both simulated and experimental data. In neural recordings from awake mice, we demonstrate how variance-pseudospectra can reveal physiologically significant information unavailable to standard expectation-based dynamical models.

Koopman operators are infinite-dimensional operators that globally linearize nonlinear dynamical systems, making their spectral information valuable for understanding dynamics. However, Koopman operators can have continuous spectra and infinite-dimensional invariant subspaces, making computing their spectral information a considerable challenge. This paper describes data-driven algorithms with rigorous convergence guarantees for computing spectral information of Koopman operators from trajectory data. We introduce residual dynamic mode decomposition (ResDMD), which provides the first scheme for computing the spectra and pseudospectra of general Koopman operators from snapshot data without spectral pollution. Using the resolvent operator and ResDMD, we compute smoothed approximations of spectral measures associated with general measure-preserving dynamical systems. We prove explicit convergence theorems for our algorithms (including for general systems that are not measure-preserving), which can achieve high-order convergence even for chaotic systems when computing the density of the continuous spectrum and the discrete spectrum. Since our algorithms have error control, ResDMD allows aposteri verification of spectral quantities, Koopman mode decompositions, and learned dictionaries. We demonstrate our algorithms on the tent map, circle rotations, Gauss iterated map, nonlinear pendulum, double pendulum, and Lorenz system. Finally, we provide kernelized variants of our algorithms for dynamical systems with a high-dimensional state space. This allows us to compute the spectral measure associated with the dynamics of a protein molecule with a 20,046-dimensional state space and compute nonlinear Koopman modes with error bounds for turbulent flow past aerofoils with Reynolds number \(>10^5\) that has a 295,122-dimensional state space.

Koopman operators globally linearize nonlinear dynamical systems and their spectral information is a powerful tool for the analysis and decomposition of nonlinear dynamical systems. However, Koopman operators are infinite-dimensional, and computing their spectral information is a considerable challenge. We introduce measure-preserving extended dynamic mode decomposition (mpEDMD), the first Galerkin method whose eigendecomposition converges to the spectral quantities of Koopman operators for general measure-preserving dynamical systems. mpEDMD is a data-driven algorithm based on an orthogonal Procrustes problem that enforces measure-preserving truncations of Koopman operators using a general dictionary of observables. It is flexible and easy to use with any pre-existing DMD-type method, and with different types of data. We prove convergence of mpEDMD for projection-valued and scalar-valued spectral measures, spectra, and Koopman mode decompositions. For the case of delay embedding (Krylov subspaces), our results include the first convergence rates of the approximation of spectral measures as the size of the dictionary increases. We demonstrate mpEDMD on a range of challenging examples, its increased robustness to noise compared with other DMD-type methods, and its ability to capture the energy conservation and cascade of a turbulent boundary layer flow with Reynolds number \(> 6\times 10^4\) and state-space dimension \(>10^5\).

Dynamic mode decomposition (DMD) describes complex dynamic processes through a hierarchy of simpler coherent features. DMD is regularly used to understand the fundamental characteristics of turbulence and is closely related to Koopman operators. However, verifying the decomposition, equivalently the computed spectral features of Koopman operators, remains a significant challenge due to the infinite-dimensional nature of Koopman operators. Challenges include spurious (unphysical) modes and dealing with continuous spectra, which both occur regularly in turbulent flows. Residual dynamic mode decomposition (ResDMD), introduced by Colbrook & Townsend (Rigorous data-driven computation of spectral properties of Koopman operators for dynamical systems. 2021. arXiv:2111.14889), overcomes such challenges through the data-driven computation of residuals associated with the full infinite-dimensional Koopman operator. ResDMD computes spectra and pseudospectra of general Koopman operators with error control and computes smoothed approximations of spectral measures (including continuous spectra) with explicit high-order convergence theorems. ResDMD thus provides robust and verified Koopmanism. We implement ResDMD and demonstrate its application in various fluid dynamic situations at varying Reynolds numbers from both numerical and experimental data. Examples include vortex shedding behind a cylinder, hot-wire data acquired in a turbulent boundary layer, particle image velocimetry data focusing on a wall-jet flow and laser-induced plasma acoustic pressure signals. We present some advantages of ResDMD: the ability to resolve nonlinear and transient modes verifiably; the verification of learnt dictionaries; the verification of Koopman mode decompositions; and spectral calculations with reduced broadening effects. We also discuss how a new ordering of modes via residuals enables greater accuracy than the traditional modulus ordering (e.g. when forecasting) with a smaller dictionary. This result paves the way for more significant dynamic compression of large datasets without sacrificing accuracy.

Topological insulators (TIs) are renowned for their remarkable electronic properties: quantised bulk Hall and edge conductivities, and robust edge wave-packet propagation, even in the presence of material defects and disorder. Computations of these physical properties generally rely on artificial periodicity (the supercell approximation, which struggles in the presence of edges), or unphysical boundary conditions (artificial truncation). In this work, we build on recently developed methods for computing spectral properties of infinite-dimensional operators. We apply these techniques to develop efficient and accurate computational tools for computing the physical properties of TIs. These tools completely avoid such artificial restrictions and allow one to probe the spectral properties of the infinite-dimensional operator directly, even in the presence of material defects, edges, and disorder. Our methods permit computation of spectra, approximate eigenstates, spectral measures, spectral projections, transport properties, and conductances. Numerical examples are given for the Haldane model, and the techniques can be extended similarly to other TIs in two and three dimensions.

Deep learning (DL) has had unprecedented success and is now entering scientific computing with full force. However, current DL methods typically suffer from instability, even when universal approximation properties guarantee the existence of stable neural networks (NNs). We address this paradox by demonstrating basic wellconditioned problems in scientific computing where one can prove the existence of NNs with great approximation qualities, however, there does not exist any algorithm, even randomised, that can train (or compute) such a NN. For any positive integersK > 2 and L, there are cases where simultaneously: (i) no randomised training algorithm can compute a NN correct to K digits with probability greater than 1/2, (ii) there exists a deterministic training algorithm that computes a NN with K−1 correct digits, but any such (even randomised) algorithm needs arbitrarily many training data, (iii) there exists a deterministic training algorithm that computes a NN with K − 2 correct digits using no more than L training samples. These results imply a classification theory describing conditions under which (stable) NNs with a given accuracy can be computed by an algorithm. We begin this theory by establishing sufficient conditions for the existence of algorithms that compute stable NNs in inverse problems. We introduce Fast Iterative REstarted NETworks (FIRENETs), which we both prove and numerically verify are stable. Moreover, we prove that only \(\mathcal{O}(\log(\epsilon^{-1}))\) layers are needed for an \(\epsilon\)-accurate solution to the inverse problem.

Instability is AI's Achilles’ heel. We show the following paradox: there are cases where stable and accurate AI exists, but it can never be trained by any algorithm. We initiate a foundations theory for when AI can be trained - such a programme will shape political and legal decision-making in the coming decades, and have a significant impact on markets for AI technologies.

Computing spectra is a central problem in computational mathematics with an abundance of applications throughout the sciences. However, in many applications gaining an approximation of the spectrum is not enough. Often it is vital to determine geometric features of spectra such as Lebesgue measure, capacity or fractal dimensions, different types of spectral radii and numerical ranges, or to detect essential spectral gaps and the corresponding failure of the finite section method. Despite new results on computing spectra and the substantial interest in these geometric problems, there remain no general methods able to compute such geometric features of spectra of infinite-dimensional operators. We provide the first algorithms for the computation of many of these longstanding problems (including the above). As demonstrated with computational examples, the new algorithms yield a library of new methods. Recent progress in computational spectral problems in infinite dimensions has led to the Solvability Complexity Index (SCI) hierarchy, which classifies the difficulty of computational problems. These results reveal that infinite-dimensional spectral problems yield an intricate infinite classification theory determining which spectral problems can be solved and with which type of algorithm. This is very much related to S. Smale's comprehensive program on the foundations of computational mathematics initiated in the 1980s. We classify the computation of geometric features of spectra in the SCI hierarchy, allowing us to precisely determine the boundaries of what computers can achieve (in any model of computation) and prove that our algorithms are optimal. We also provide a new universal technique for establishing lower bounds in the SCI hierarchy, which both greatly simplifies previous SCI arguments and allows new, formerly unattainable, classifications.

The problem of computing spectra of operators is arguably one of the most investigated areas of computational mathematics. However, the problem of computing spectra of general bounded infinite matrices has only recently been solved. We establish some of the foundations of computational spectral theory through the Solvability Complexity Index (SCI) hierarchy, an approach closely related to Smale's program on the foundations of computational mathematics and McMullen's results on polynomial root finding with rational maps. Infinite-dimensional problems yield an intricate infinite classification theory, determining which spectral problems can be solved and with what types of algorithms. We provide answers to many longstanding open questions on the existence of algorithms. For example, we show that spectra can be computed, with error control, from point sampling operator coefficients for large classes of partial differential operators on unbounded domains. Further results include: computing spectra of (possibly unbounded) operators on graphs and separable Hilbert spaces with error control; determining if the spectrum intersects a compact set; the computational spectral gap problem and computing spectral classifications at the bottom of the spectrum; and computing discrete spectra, multiplicities, eigenspaces and determining if the discrete spectrum is non-empty. Moreover, the positive results with error control can be used in computer-assisted proofs. In contrast, the negative results preclude computer-assisted proofs for classes of operators as a whole. Our proofs are constructive, yielding a library of new algorithms and techniques that handle problems that before were out of reach. We demonstrate these algorithms on challenging problems, giving concrete examples of the failure of traditional approaches (e.g., ``spectral pollution'') compared to the introduced techniques.

Robust edge transport can occur when charged particles in crystalline lattices interact with an applied external magnetic field. Such systems have a spectrum composed of bands of bulk states and in-gap edge states. For quasicrystalline systems, we still expect to observe the basic characteristics of bulk states and current-carrying edge states. We show that, for quasicrystals in magnetic fields, there is an additional third option—bulk localized transport (BLT) states. BLT states share the in-gap nature of the well-known edge states and can support transport, but they are fully contained within the bulk of the system, with no support along the edge. Thus, transport is possible along the edge and within distinct regions of the bulk. We consider both finite-size and infinite-size systems, using rigorous error controlled computational techniques that are not prone to finite-size effects. BLT states are preserved for infinite-size systems, in stark contrast to edge states. This allows us to observe transport in infinite-size systems, without any perturbations, defects, or boundaries being introduced. We confirm the in-gap topological nature of BLT states for finite- and infinite-size systems by computing the Bott index and local Chern marker (common topological measures). BLT states form due to magnetic aperiodicity, arising from the interplay of lengthscales between the magnetic field and the quasiperiodic lattice. BLT could have interesting applications similar to those of edge states, but now taking advantage of the larger bulk of the lattice. The infinite-size techniques introduced here, especially the calculation of topological measures, could also be widely applied to other crystalline, quasicrystalline, and disordered models.

We develop an algorithm that computes strongly continuous semigroups on infinite-dimensional Hilbert spaces with explicit error control. Given a generator \(A\), a time \(t>0\), an arbitrary initial vector \(u_0\) and an error tolerance \(\epsilon>0\), the algorithm computes \(\exp(tA)u_0\) with error bounded by \(\epsilon\). The algorithm is based on a combination of a regularized functional calculus, suitable contour quadrature rules, and the adaptive computation of resolvents in infinite dimensions. As a particular case, we show that it is possible, even when only allowing pointwise evaluation of coefficients, to compute, with error control, semigroups on the unbounded domain \(L^2(\mathbb{R}^d)\) that are generated by partial differential operators with polynomially bounded coefficients of locally bounded total variation. For analytic semigroups (and more general Laplace transform inversion), we provide a quadrature rule whose error decreases like \(\exp(−cN/\log(N))\) for \(N\) quadrature points, that remains stable as \(N\rightarrow\infty\), and which is also suitable for infinite-dimensional operators. Numerical examples are given, including: Schrödinger and wave equations on the aperiodic Ammann-Beenker tiling, complex perturbed fractional diffusion equations on \(L^2(\mathbb{R})\), and damped Euler–Bernoulli beamequations.

We develop a rapid and accurate contour method for the solution of time-fractional PDEs. The method inverts the Laplace transform via an optimised stable quadrature rule, suitable for infinite-dimensional operators, whose error decreases like \(\exp(-cN/\log(N))\) for \(N\) quadrature points. The method is parallisable, avoids having to resolve singularities of the solution as \(t\downarrow 0\), and avoids the large memory consumption that can be a challenge for time-stepping methods applied to time-fractional PDEs. The ODEs resulting from quadrature are solved using adaptive sparse spectral methods that converge exponentially with optimal linear complexity. These solutions of ODEs are reused for different times. We provide a complete analysis of our approach for fractional beam equations used to model small-amplitude vibration of viscoelastic materials with a fractional Kelvin-Voigt stress-strain relationship. We calculate the system's energy evolution over time and the surface deformation in cases of both constant and non-constant viscoelastic parameters. An infinite-dimensional "solve-then-discretise" approach considerably simplifies the analysis, which studies the generalisation of the numerical range of a quasi-linearisation of a suitable operator pencil. This allows us to build an efficient algorithm with explicit error control. The approach can be readily adapted to other time-fractional PDEs and is not constrained to fractional parameters in the range \(0<\nu<1\).

Structured compressed sensing takes signal structure into account, thereby outperforming earlier compressed sensing methods. However, results are usually based on sampling in the Fourier domain, such as in Magnetic Resonance Imaging. In the time domain, the benefits of structured compressed sensing are still unknown. This paper introduces concepts that incorporate the signal structure into both the acquisition and reconstruction of compressed sensing in time and image domain applications. First, a stratified-random sampling pattern is proposed to improve the recovery of the dominant low-frequency range of natural signals. A heuristic decay of primes criterion is developed to evaluate the properties of the sensing matrix and is used to optimize the sampling pattern. Second, the sparsity of the Fourier transform as the representation domain is improved by estimating the signal structure in a preprocessing step, and then adapting the grid of the Fourier transform. In contrast to existing methods, grid stretching is integrated into the fast Fourier transform to reduce computational complexity. Both structured acquisition and reconstruction are evaluated using simulations, as well as two real-world applications: wireless sensor networks in structural health monitoring and electron microscopy. Results show that both reconstruction errors and robustness can significantly be improved by incorporating structure into the acquisition and reconstruction.

Using the resolvent operator, we develop an algorithm for computing smoothed approximations of spectral measures associated with self-adjoint operators. The algorithm can achieve arbitrarily high-orders of convergence in terms of a smoothing parameter for computing spectral measures of general differential, integral, and lattice operators. Explicit pointwise and \(L^p\)-error bounds are derived in terms of the local regularity of the measure. We provide numerical examples, including a partial differential operator, a magnetic tight-binding model of graphene, and compute one thousand eigenvalues of a Dirac operator to near machine precision without spectral pollution. The algorithm is publicly available in \(\texttt{SpecSolve}\), which is a software package written in MATLAB.

Spectral measures arise in numerous applications such as quantum mechanics, signal processing, resonance phenomena, and fluid stability analysis. Similarly, spectral decompositions (into pure point, absolutely continuous and singular continuous parts) often characterise relevant physical properties such as the long-time dynamics of quantum systems. Despite new results on computing spectra, there remains no general method able to compute spectral measures or spectral decompositions of infinite-dimensional normal operators. Previous efforts have focused on specific examples where analytical formulae are available (or perturbations thereof) or on classes of operators that carry a lot of structure. Hence the general computational problem is predominantly open. We solve this problem by providing the first set of general algorithms that compute spectral measures and decompositions of a wide class of operators. Given a matrix representation of a self-adjoint or unitary operator, such that each column decays at infinity at a known asymptotic rate, we show how to compute spectral measures and decompositions. We discuss how these methods allow the computation of objects such as the functional calculus, and how they generalise to a large class of partial differential operators, allowing, for example, solutions to evolution PDEs such as the linear Schrödinger equation on \(L^2(\mathbb{R}^d)\). Computational spectral problems in infinite dimensions have led to the Solvability Complexity Index (SCI) hierarchy, which classifies the difficulty of computational problems. We classify the computation of measures, measure decompositions, types of spectra, functional calculus, and Radon-Nikodym derivatives in the SCI hierarchy. The new algorithms are demonstrated to be efficient on examples taken from orthogonal polynomials on the real line and the unit circle (giving, for example, computational realisations of Favard's theorem and Verblunsky's theorem, respectively), and are applied to evolution equations on a two-dimensional quasicrystal.

This paper presents a rapid numerical method for predicting the aerodynamic noise generated by foam-like porous aerofoils. In such situations, particularly for high-frequency noise sources, Darcy's law may be unsuitable for describing the pressure jump across the aerofoil. Therefore, an inertial Forchheimer correction is introduced. This results in a non-linear boundary condition relating the pressure jump across the material to the fluid displacement. We aim to provide a quick, semi-analytical model that incorporates such non-linear effects without requiring a full turbulent simulation. The numerical scheme implemented is based on local Mathieu function expansions, leading to a semi-analytical boundary spectral method that is well-suited to both linear and non-linear boundary conditions (including boundary conditions more general than the Forchheimer correction). In the latter case, Newton's method is employed to solve the resulting non-linear system of equations for the unknown coefficients. Whilst the physical model is simplified to consider just the scattering by a thin porous aerofoil with no background flow, when the non-linear inertial correction is included good agreement is seen between the model predictions and both experimental results and large eddy simulations. It is found that for sufficiently low-permeability materials, the effects of inertia can outweigh the noise attenuation effects of viscosity. This helps explain the discrepancy between experimental results and previous (linear) low-fidelity numerical simulations or analytical predictions, which typically overestimate the noise reduction capabilities of porous aerofoils.

This paper considers the effects of smoothly varying chordwise porosity of a finite perforated plate on turbulence-aerofoil interaction noise. The aeroacoustic model is made possible through the use of a novel Mathieu function collocation method, rather than a traditional Wiener--Hopf approach which would be unable to deal with chordwise varying quantities. The main focus is on two bio-inspired porosity distributions, modelled from air flow resistance data obtained from the wings of barn owls (tyto alba) and common buzzards (buteo buteo). Trailing-edge noise is much reduced for the owl-like distribution, but, perhaps surprisingly, so too is leading-edge noise, despite both wings having similar porosity values at the leading edge. A general monotonic variation is then considered indicating that there may indeed be a significant acoustic impact of how the porosity is distributed along the whole chord of the plate, not just its values at the scattering edges. Through this investigation, it is found that a plate whose porosity continuously decreases from the trailing edge to a zero-porosity leading edge can, in fact, generate lower levels of trailing-edge noise than a plate whose porosity remains constant at the trailing-edge value.

Many problems in fluid mechanics and acoustics can be modelled by Helmholtz scattering off poro-elastic plates. We develop a boundary spectral method, based on collocation of local Mathieu function expansions, for Helmholtz scattering off multiple variable poro-elastic plates in two dimensions. Such boundary conditions, namely the varying physical parameters and coupled thin-plate equation, present a considerable challenge to current methods. The new method is fast, accurate, and flexible, with the ability to compute expansions in thousands (and even tens of thousands) of Mathieu functions, thus making it a favourable method for the considered geometries. Comparisons are made with elastic boundary element methods, where the new method is found to be faster and more accurate. Our solution representation directly provides a sine series approximation of the far-field directivity and can be evaluated near or on the scatterers, meaning that the near field can be computed stably and efficiently. The new method also allows us to examine the effects of varying stiffness along a plate, which is poorly studied due to limitations of other available techniques. We show that a power-law decrease to zero in stiffness parameters gives rise to unexpected scattering and aeroacoustic effects similar to an acoustic black hole metamaterial.

Kernel density estimation on a finite interval poses an outstanding challenge because of the well-recognized bias at the boundaries of the interval. Motivated by an application in cancer research, we consider a boundary constraint linking the values of the unknown target density function at the boundaries. We provide a kernel density estimator (KDE) that successfully incorporates this linked boundary condition, leading to a non-self-adjoint diffusion process and expansions in nonseparable generalized eigenfunctions. The solution is rigorously analyzed through an integral representation. The new KDE possesses many desirable properties, such as consistency, asymptotically negligible bias at the boundaries, and an increased rate of approximation, as measured by the AMISE. We apply our method to the motivating example in biology and provide numerical experiments with synthetic data, including comparisons with state-of-the-art KDEs (which currently cannot handle linked boundary constraints). Results suggest that the new method is fast and accurate. Furthermore, we demonstrate how to build statistical estimators of the boundary conditions satisfied by the target function without a priori knowledge. Our analysis can also be extended to more general boundary conditions that may be encountered in applications.

We provide the first significant extension of the unified transform (also known as the Fokas method) applied to elliptic boundary value problems, namely, we extend the method to curvilinear polygons and partial differential equations (PDEs) with variable coefficients. This is used to solve the generalized Dirichlet-to-Neumann map. The central component of the unified transform is the coupling of certain integral transforms of the given boundary data and of the unknown boundary values. This has become known as the global relation and, in the case of constant coefficient PDEs, simply links the Fourier transforms of the Dirichlet and Neumann boundary values. We extend the global relation to PDEs with variable coefficients and to domains with curved boundaries. Furthermore, we provide a natural choice of global relations for separable PDEs. These generalizations are numerically implemented using a method based on Chebyshev interpolation for efficient and accurate computation of the integral transforms that appear in the global relation. Extensive numerical examples are provided, demonstrating that the method presented in this paper is both accurate and fast, yielding exponential convergence for sufficiently smooth solutions. Furthermore, the method is straightforward to use, involving just the construction of a simple linear system from the integral transforms, and does not require knowledge of Green’s functions of the PDE. Details on the implementation are discussed at length.

There is an increasing literature on random non-self-adjoint infinite matrices with motivations ranging from condensed matter physics to neural networks. Many of these operators fall into the class of "pseudoergodic" operators, which allows the elimination of probabilistic arguments when studying spectral properties. Parallel to this is the increased awareness that spectral properties of non-self-adjoint operators, in particular stability, may be better captured via the notion of pseudospectra as opposed to spectra. Although it is well known that the finite section method applied to these matrices does not converge to the spectrum, it is often found in practice that the pseudospectrum behaves better with appropriate boundary conditions. We make this precise by giving a simple proof that the finite section method with periodic boundary conditions converges to the pseudospectrum of the full operator. Our results hold in any dimension (not just for banded bi-infinite matrices) and can be considered as a generalisation of the well-known classical result for banded Laurent operators and their circulant approximations. Furthermore, we numerically demonstrate a convergent algorithm for the pseudospectrum, including cases where periodic boundary conditions converge faster than the method of uneven sections. Finally, we show that the result carries over to pseudoergodic operators acting on \(l^p\) spaces for \(p\in [1,\infty ]\).

This paper considers the use of compliant boundary conditions to provide a homogenized model of a finite array of collinear plates, modelling a perforated screen or grating. While the perforated screen formally has a mix of Dirichlet and Neumann boundary conditions, the homogenized model has Robin boundary conditions. Perforated screens form a canonical model in scattering theory, with applications ranging from electromagnetism to aeroacoustics. Interest in perforated media incorporated within larger structures motivates interrogating the appropriateness of homogenized boundary conditions in this case, especially as the homogenized model changes the junction behaviour considered at the extreme edges of the screen. To facilitate effective investigation we consider three numerical methods solving the Helmholtz equation: the unified transform and an iterative Wiener–Hopf approach for the exact problem of a set of collinear rigid plates (the difficult geometry of the problem means that such methods, which converge exponentially, are crucial) and a novel Mathieu function collocation approach to consider a variable compliance applied along the length of a single plate. We detail the relative performance and practical considerations for each method. By comparing solutions obtained using homogenized boundary conditions to the problem of collinear plates, we verify that the constant compliance given in previous theoretical research is appropriate to gain a good estimate of the solution even for a modest number of plates, provided we are sufficiently far into the asymptotic regime. We further investigate tapering the compliance near the extreme endpoints of the screen and find that tapering with tanh functions reduces the error in the approximation of the far field (if we are sufficiently far into the asymptotic regime). We also find that the number of plates and wavenumber has significant effects, even far into the asymptotic regime. These last two points indicate the importance of modelling end effects to achieve highly accurate results.

Computing the spectra of operators is a fundamental problem in the sciences, with wide-ranging applications in condensed-matter physics, quantum mechanics and chemistry, statistical mechanics, etc. While there are algorithms that in certain cases converge to the spectrum, no general procedure is known that (a) always converges, (b) provides bounds on the errors of approximation, and (c) provides approximate eigenvectors. This may lead to incorrect simulations. It has been an open problem since the 1950s to decide whether such reliable methods exist at all. We affirmatively resolve this question, and the algorithms provided are optimal, realizing the boundary of what digital computers can achieve. Moreover, they are easy to implement and parallelize, offer fundamental speed-ups, and allow problems that before, regardless of computing power, were out of reach. Results are demonstrated on difficult problems such as the spectra of quasicrystals and non-Hermitian phase transitions in optics.

Spectral computations of infinite-dimensional operators are notoriously difficult, yet ubiquitous in the sciences. Indeed, despite more than half a century of research, it is still unknown which classes of operators allow for the computation of spectra and eigenvectors with convergence rates and error control. Recent progress in classifying the difficulty of spectral problems into complexity hierarchies has revealed that the most difficult spectral problems are so hard that one needs three limits in the computation, and no convergence rates nor error control is possible. This begs the question: which classes of operators allow for computations with convergence rates and error control? In this paper, we address this basic question, and the algorithm used is an infinite-dimensional version of the QR algorithm. Indeed, we generalise the QR algorithm to infinite-dimensional operators. We prove that not only is the algorithm executable on a finite machine, but one can also recover the extremal parts of the spectrum and corresponding eigenvectors, with convergence rates and error control. This allows for new classification results in the hierarchy of computational problems that existing algorithms have not been able to capture. The algorithm and convergence theorems are demonstrated on a wealth of examples with comparisons to standard approaches (that are notorious for providing false solutions). We also find that in some cases the IQR algorithm performs better than predicted by theory and make conjectures for future study.

This paper implements the unified transform to problems in unbounded domains with solutions having corner singularities. Consequently, a wide variety of mixed boundary condition problems can be solved without the need for the Wiener–Hopf technique. Such problems arise frequently in acoustic scattering or in the calculation of electric fields in geometries involving finite and/or multiple plates. The new approach constructs a global relation that relates known boundary data, such as the scattered normal velocity on a rigid plate, to unknown boundary values, such as the jump in pressure upstream of the plate. By approximating the known data and the unknown boundary values by suitable functions and evaluating the global relation at collocation points, one can accurately obtain the expansion coefficients of the unknown boundary values. The method is illustrated for the modified Helmholtz and Helmholtz equations. In each case, comparisons between the traditional Wiener–Hopf approach, other spectral or boundary methods and the unified transform approach are discussed.

Recent work has given rise to a novel and simple numerical technique for solving elliptic boundary value problems formulated in convex polygons in two dimensions. The method, based on the unified transform, involves expanding the unknown boundary values in a Legendre basis and determining the expansion coefficients by evaluating the so-called global relation at appropriate points in the complex Fourier plane (spectral collocation). In this paper we provide a significant advancement of this numerical technique by providing a fast and efficient method to evaluate the solution in the domain interior. The use of a Legendre basis allows the relevant integrals to be computed efficiently and accurately using Chebyshev interpolation, even for large degree. For the particular case of the Laplace equation this allows an explicit expansion in the domain interior in terms of hypergeometric functions. Evaluation in the interior is found to converge more rapidly than the approximation of the unknown boundary values, allowing accurate approximation of solutions with weak corner singularities. For stronger singularities, the method can be combined with global singular functions for rapid convergence. Numerical examples are provided, showing that the method compares well against standard spectral methods and opens up the possibility of applying the method to general curvilinear domains.

This paper presents a new approach to solving acoustic scattering problems: the Unified Transform method. This spectral, boundary-based collocation method can be readily applied to acoustic scattering by disjoint two-dimensional structures, and, for the purposes of this paper, is illustrated in the case of multiple flat plates, which also addresses the additional difficulty of mathematical singularities in the scattered field due to diffraction at sharp edges. Fluid-structure interaction may also be incorporated into the method, such as plate elasticity, which when applied to aerofoil trailing edges, is known to reduce aerodynamic noise. While a range of examples are illustrated to show the versatility of the method, attention is in particular given to the scattering of quadrupole sources by rigid plates with finite elastic extensions. It is seen that whilst a fully elastic plate is most beneficial acoustically, plates with only small extensions can considerably reduce the far-field sound power versus a fully rigid plate.

This paper employs the unified transform, also known as the Fokas method, to solve the advection-dispersion equation on the half-line. This method combines complex analysis with numerics. Compared to classical approaches used to solve linear partial differential equations (PDEs), the unified transform avoids the solution of ordinary differential equations and, more importantly, constructs an integral representation of the solution in the complex plane which is uniformly convergent at the boundaries. As a consequence, such solutions are well suited for numerical computations. Indeed, the numerical evaluation of the solution requires only the computation of a single contour integral involving an integrand which decays exponentially fast for large values of the integration variable. A novel contribution of this paper, with respect to the solution of linear evolution PDEs in general, and the implementation of the unified transform in particular, is the following: using the advection-dispersion equation as a generic example, it is shown that if the transforms of the given data can be computed analytically, then the unified transform yields a fast and accurate method that converges exponentially with the number of evaluations N yet only has complexity. Furthermore, if the transforms are computed numerically using \(\mathcal{O}(M)\) evaluations, the unified transform gives rise to a method with complexity \(\mathcal{O}(MN)\). Results are successfully compared to other existing solutions.

There exists a growing literature on using the Fokas method (unified transform method) to solve Laplace and Helmholtz problems on convex polygonal domains. We show here that the convexity requirement can be eliminated by the use of a ‘virtual side’ concept, thereby significantly increasing the flexibility and utility of the approach. We also show that the inclusion of singular functions in the basis to treat corner singularities can greatly increase the rate of convergence of the method. The method also compares well with other standard methods used to cope with corner singularities. An example is given where this inclusion leads to exponential convergence. As well as this, we give new results on several additional issues, including the choice of collocation points and calculation of solutions throughout domain interiors. An appendix illustrates the algebraic simplicity of the methodology by showing how the core part of the present approach can be implemented in only about a dozen lines of MATLAB code.

Passive-scalar mixing (metals, molecules, etc.) in the turbulent interstellar medium (ISM) is critical for abundance patterns of stars and clusters, galaxy and star formation, and cooling from the circumgalactic medium. However, the fundamental scaling laws remain poorly understood in the highly supersonic, magnetized, shearing regime relevant for the ISM. We therefore study the full scaling laws governing passive-scalar transport in idealized simulations of supersonic turbulence. Using simple phenomenological arguments for the variation of diffusivity with scale based on Richardson diffusion, we propose a simple fractional diffusion equation to describe the turbulent advection of an initial passive scalar distribution. These predictions agree well with the measurements from simulations, and vary with turbulent Mach number in the expected manner, remaining valid even in the presence of a large-scale shear flow (e.g. rotation in a galactic disc). The evolution of the scalar distribution is not the same as obtained using simple, constant ‘effective diffusivity’ as in Smagorinsky models, because the scale dependence of turbulent transport means an initially Gaussian distribution quickly develops highly non-Gaussian tails. We also emphasize that these are mean scalings that apply only to ensemble behaviours (assuming many different, random scalar injection sites): individual Lagrangian ‘patches’ remain coherent (poorly mixed) and simply advect for a large number of turbulent flow-crossing times.

Spanwise homogeneous turbulent boundary layers convect groups of flow structures that move with similar momentum, often called coherent structures. These coherent structures remain statistically similar as they get convected downstream. The length scales of such coherent structures vary with the wall-normal distance of the boundary layer. Large-scale motions on the order of ~20\(\delta\), often termed superstructures, are conjectured to be present in the outer layer. However, Particle Image Velocimetry (PIV) data acquired at the bottom of the logarithmic layer of the boundary layer provides evidence of the existence of superstructures in the overlap region. Their presence in the close vicinity of the wall makes them a strong candidate to be the prime contributors to the low-wavenumber pressure fluctuations in a smooth wall. A discussion on the fundamental nature of these superstructures in the smooth-wall turbulent boundary layer is presented, followed by the flow field's decomposition into its streamwise wavenumber and frequency components. The experimental data is also analyzed using a novel dynamic mode decomposition technique, ResDMD. This decomposition approach differs from other known techniques due to error control, verification (e.g., of dictionaries) and convergence theorems. This technique was used to uncover the transient behavior within the system and identify turbulence events based on their length scales and convection velocities. The extent of the streamwise flow homogeneity is also discussed. Evidence of high spectral levels at low-wavenumbers confirms the role of superstructures in containing a significant fraction of the turbulence energy. The streamwise wavenumber-frequency spectra of the pressure fluctuations at sub-convective wavenumbers support this.

Self-adjoint operators on infinite-dimensional spaces with continuous spectra are abundant but do not possess a basis of eigenfunctions. Rather, diagonalization is achieved through spectral measures. The SpecSolve package [SIAM Rev., 63(3) (2021), pp. 489–524] computes spectral measures of general (self-adjoint) differential and integral operators by combining state-of-the-art adaptive spectral methods with an efficient resolvent-based strategy. The algorithm achieves arbitrarily high orders of convergence in terms of a smoothing parameter, allowing computation of both discrete and continuous spectral components. This article extends SpecSolve to two important classes of operators: singular integro-differential operators and general operator pencils. Essential computational steps are performed with off-the-shelf spectral methods, including spectral methods on the real line, the ultraspherical spectral method, Chebyshev and Fourier spectral methods, and the (hp-adaptive and sparse) ultraspherical spectral element method. This collection illustrates the power and flexibility of SpecSolve's "discretization-oblivious" paradigm.

We consider a finite perforated plate and study the effects of smoothly varying chordwise porosity on turbulence-aerofoil interaction noise. To study this problem, we use a novel Mathieu function collocation method, rather than a traditional Wiener-Hopf approach which would be unable to deal with chordwise varying quantities. Our main focus is on two bio-inspired porosity distributions, which are modelled based on air flow resistance data from the wings of barn owls and common buzzards. As expected, trailing-edge noise is much reduced for the owl-like distribution. However, and perhaps surprisingly, so too is leading-edge noise, despite both wings having similar porosity values at the leading edge. We then consider a general monotonic variation. Our study indicates that there may be a significant acoustic impact of how the porosity is distributed along the whole chord of the plate (i.e. not just its values at the scattering edges). Indeed, a plate whose porosity continuously decreases from the trailing edge to a zero-porosity leading edge can, in fact, generate lower levels of trailing-edge noise than a plate whose porosity remains constant at the trailing-edge value.

We present a fast numerical method for predicting the aerodynamic noise generated by foam-like porous aerofoils. Darcy's law, describing the pressure jump across the aerofoil, may be unsuitable in such situations, particularly for high-frequency noise sources where the unsteady velocity fluctuations are large relative to the porous structure size. We therefore introduce an inertial Forchheimer correction that results in a non-linear boundary condition relating the pressure jump across the material to the fluid displacement. Based on local Mathieu function expansions, we provide a semi-analytical boundary spectral method that is well-suited to both linear and non-linear boundary conditions. In the latter case, Newton's method is employed to solve the resulting non-linear system of equations. The outcome is a fast semi-analytical model that incorporates such non-linear effects without requiring a full turbulent simulation. Whilst we consider only the simplified case of scattering by a thin porous aerofoil with no background flow, when the non-linear inertial correction is included good agreement is seen between the model predictions and experimental results. A key conclusion is that for sufficiently low-permeability materials, the effects of inertia can outweigh the noise attenuation effects of viscosity. This helps explain the discrepancy between experimental results and previous (linear) low-fidelity numerical simulations or analytical predictions, which typically overestimate the noise reduction capabilities of porous aerofoils.

Deep learning has emerged as a competitive new tool in image reconstruction. However, recent results demonstrate such methods are typically highly unstable - tiny, almost undetectable perturbations cause severe artefacts in the reconstruction, a major concern in practice. This is paradoxical given the existence of stable state-of-the-art methods for these problems. Thus, approximation theoretical results nonconstructively imply the existence of stable and accurate neural networks. Hence the fundamental question: Can we explicitly construct/train stable and accurate neural networks for image reconstruction? We prove the answer is yes and construct such networks. Numerical examples of the competitive performance are also provided.

This paper employs the unified transform to present a boundary-based spectral collocation method suitable for solving acoustic scattering problems. The method is suitable for both interior and exterior scattering problems, and may be extended to three dimensions. A number of simple two-dimensional examples are presented to illustrate the versatility of this method, and, upon comparison with other spectral methods or boundary-based methods, the approach presented in this paper is shown to be very competitive.

Spectral computations in infinite dimensions are ubiquitous in the sciences. However, their many applications and theoretical studies depend on computations which are infamously difficult. This thesis, therefore, addresses the broad question, "What is computationally possible within the field of spectral theory of separable Hilbert spaces?" The boundaries of what computers can achieve in computational spectral theory and mathematical physics are unknown, leaving many open questions that have been unsolved for decades. This thesis provides solutions to several such long-standing problems.

To determine these boundaries, we use the Solvability Complexity Index (SCI) hierarchy, an idea which has its roots in Smale's comprehensive programme on the foundations of computational mathematics. The Smale programme led to a real-number counterpart of the Turing machine, yet left a substantial gap between theory and practice. The SCI hierarchy encompasses both these models and provides universal bounds on what is computationally possible. What makes spectral problems particularly delicate is that many of the problems can only be computed by using several limits, a phenomenon also shared in the foundations of polynomial root-finding as shown by McMullen. We develop and extend the SCI hierarchy to prove optimality of algorithms and construct a myriad of different methods for infinite-dimensional spectral problems, solving many computational spectral problems for the first time.

For arguably almost any operator of applicable interest, we solve the long-standing computational spectral problem and construct algorithms that compute spectra with error control. This is done for partial differential operators with coefficients of locally bounded total variation and also for discrete infinite matrix operators. We also show how to compute spectral measures of normal operators (when the spectrum is a subset of a regular enough Jordan curve), including spectral measures of classes of self-adjoint operators with error control and the construction of high-order rational kernel methods. We classify the problems of computing measures, measure decompositions, types of spectra (pure point, absolutely continuous, singular continuous), functional calculus, and Radon–Nikodym derivatives in the SCI hierarchy. We construct algorithms for and classify; fractal dimensions of spectra, Lebesgue measures of spectra, spectral gaps, discrete spectra, eigenvalue multiplicities, capacity, different spectral radii and the problem of detecting algorithmic failure of previous methods (finite section method). The infinite-dimensional QR algorithm is also analysed, recovering extremal parts of spectra, corresponding eigenvectors, and invariant subspaces, with convergence rates and error control. Finally, we analyse pseudospectra of pseudoergodic operators (a generalisation of random operators) on vector-valued lp spaces.

All of the algorithms developed in this thesis are sharp in the sense of the SCI hierarchy. In other words, we prove that they are optimal, realising the boundaries of what digital computers can achieve. They are also implementable and practical, and the majority are parallelisable. Extensive numerical examples are given throughout, demonstrating efficiency and tackling difficult problems taken from mathematics and also physical applications.

In summary, this thesis allows scientists to rigorously and efficiently compute many spectral properties for the first time. The framework provided by this thesis also encompasses a vast number of areas in computational mathematics, including the classical problem of polynomial root-finding, as well as optimisation, neural networks, PDEs and computer-assisted proofs. This framework will be explored in the future work of the author within these settings.