Publications

Many data-driven algorithms in dynamical systems rely on ergodic averages that converge painfully slowly. One simple idea changes this: taper the ends. Weighted Birkhoff averages can converge much faster (sometimes superpolynomially, even exponentially) and can be incorporated seamlessly into existing methods. We demonstrate this with five weighted algorithms: weighted Dynamic Mode Decomposition (wtDMD), weighted Extended DMD (wtEDMD), weighted Sparse Identification of Nonlinear Dynamics (wtSINDy), weighted spectral measure estimation, and weighted diffusion forecasting. Across examples ranging from fluid flows to El Niño data, the message is clear: weighting costs nothing, is easy to implement, and often delivers markedly better results from the same data.

Nonlinear spectral problems arise across a range of fields, including mechanical vibrations, fluid-solid interactions, and photonic crystals. Discretizing infinite-dimensional nonlinear spectral problems often introduces significant computational challenges, particularly spectral pollution and invisibility, which can distort or obscure the true underlying spectrum. We present the first general, convergent computational method for computing the spectra and pseudospectra of nonlinear spectral problems. Our approach uses new results on nonlinear injection moduli and requires only minimal continuity assumptions: specifically, continuity with respect to the gap metric on operator graphs, making it applicable to a broad class of problems. We use the Solvability Complexity Index (SCI) hierarchy, which has recently been used to resolve the classical linear problem, to systematically classify the computational complexity of nonlinear spectral problems. Our results establish the optimality of the method and reveal that Hermiticity does not necessarily simplify the computational complexity of these nonlinear problems. Comprehensive examples -- including nonlinear shifts, Klein--Gordon equations, wave equations with acoustic boundary conditions, time-fractional beam equations, and biologically inspired delay differential equations -- demonstrate the robustness, accuracy, and broad applicability of our methodology.

Koopman operators provide a linear framework for data-driven analyses of nonlinear dynamical systems, but their infinite-dimensional nature presents major computational challenges. In this article, we offer an introductory guide to Koopman learning, emphasizing rigorously convergent data-driven methods for forecasting and spectral analysis. We provide a unified account of error control via residuals in both finite- and infinite-dimensional settings, an elementary proof of convergence for generalized Laplace analysis—a variant of filtered power iteration that works for operators with continuous spectra and no spectral gaps—and review state-of-the-art approaches for computing continuous spectra and spectral measures. The goal is to provide both newcomers and experts with a clear, structured overview of reliable data-driven techniques for Koopman spectral analysis.

The fluctuation-dissipation theorem is a cornerstone result in statistical mechanics that can be used to translate the statistics of the free natural variability of a system into information on its forced response to perturbations. By combining this viewpoint on response theory with the key ingredients of Koopmanism, it is possible to deconstruct virtually any response operator into a sum of terms, each associated with a specific mode of natural variability of the system. This dramatically improves the interpretability of the resulting response formulas. We show here on three simple yet mathematically meaningful examples how to use the extended dynamic mode decomposition algorithm on an individual trajectory of the system to compute with high accuracy correlation functions as well as Green functions associated with acting forcings. This demonstrates the great potential of using Koopman analysis for the key problem of evaluating and testing the sensitivity of a complex system.

We establish a connection between quantum mechanics and computation, revealing fundamental limitations for algorithms computing spectra, especially in non-Hermitian settings. Introducing the concept of locally trivial pseudospectra, we show such assumptions are necessary for spectral computation. Locally trivial pseudospectra adapt dynamically to system energies, enabling spectral analysis across a broad class of challenging non-Hermitian problems. Exploiting this framework, we overcome a longstanding obstacle by computing the eigenvalues and eigenfunctions of the imaginary cubic oscillator \(H_{\mathrm{B}} = p^2 + i x^3\) with error bounds and no spurious modes—yielding, to our knowledge, the first such error-controlled result. We confirm, for instance, the 100th eigenvalue as \(627.6947122484365113526737029011536\ldots\). Here, truncation-induced \(\mathcal{PT}\)-symmetry breaking causes spurious eigenvalues—a pitfall our method avoids, highlighting the link between truncation and physics. Finally, we illustrate the approach's generality via spectral computations for a range of physically relevant operators. This letter provides a rigorous framework linking computational theory to quantum mechanics and offers a precise tool for spectral calculations with error bounds.

We introduce the Rigged Dynamic Mode Decomposition (Rigged DMD) algorithm, which computes generalized eigenfunction decompositions of Koopman operators. By considering the evolution of observables, Koopman operators transform complex nonlinear dynamics into a linear framework suitable for spectral analysis. While powerful, traditional Dynamic Mode Decomposition (DMD) techniques often struggle with continuous spectra. Rigged DMD addresses these challenges with a data-driven methodology that approximates the Koopman operator's resolvent and its generalized eigenfunctions using snapshot data from the system's evolution. At its core, Rigged DMD builds wave-packet approximations for generalized Koopman eigenfunctions and modes by integrating Measure-Preserving Extended Dynamic Mode Decomposition with high-order kernels for smoothing. This provides a robust decomposition encompassing both discrete and continuous spectral elements. We derive explicit high-order convergence theorems for generalized eigenfunctions and spectral measures. Additionally, we propose a novel framework for constructing rigged Hilbert spaces using time-delay embedding, significantly extending the algorithm's applicability (Rigged DMD can be used with any rigging). We provide examples, including systems with a Lebesgue spectrum, integrable Hamiltonian systems, the Lorenz system, and a high-Reynolds number lid-driven flow in a two-dimensional square cavity, demonstrating Rigged DMD's convergence, efficiency, and versatility. This work paves the way for future research and applications of decompositions with continuous spectra.

We introduce a simple, general, and convergent scheme to compute generalized eigenfunctions of self-adjoint operators with continuous spectra on rigged Hilbert spaces. Our approach does not require prior knowledge about the eigenfunctions, such as asymptotics or other analytic properties. Instead, we carefully sample the range of the resolvent operator to construct smooth and accurate wave packet approximations to generalized eigenfunctions. We prove high-order convergence in key topologies, including weak-* convergence for distributional eigenfunctions, uniform convergence on compact sets for locally smooth generalized eigenfunctions, and convergence in seminorms for separable Fréchet spaces, covering the majority of physical scenarios. The method's performance is illustrated with applications to both differential and integral operators, culminating in the computation of spectral measures and generalized eigenfunctions for an operator associated with Poincaré’s internal waves problem. These computations corroborate experimental results and highlight the method's utility for a broad range of spectral problems in physics.

Koopman operators are infinite-dimensional operators that linearize nonlinear dynamical systems, facilitating the study of their spectral properties and enabling the prediction of the time evolution of observable quantities. Recent methods have aimed to approximate Koopman operators while preserving key structures. However, approximating Koopman operators typically requires a dictionary of observables to capture the system's behavior in a finite-dimensional subspace. The selection of these functions is often heuristic, may result in the loss of spectral information, and can severely complicate structure preservation. This paper introduces Multiplicative Dynamic Mode Decomposition (MultDMD), which enforces the multiplicative structure inherent in the Koopman operator within its finite-dimensional approximation. Leveraging this multiplicative property, we guide the selection of observables and define a constrained optimization problem for the matrix approximation, which can be efficiently solved. MultDMD presents a structured approach to finite-dimensional approximations and can more accurately reflect the spectral properties of the Koopman operator. We elaborate on the theoretical framework of MultDMD, detailing its formulation, optimization strategy, and convergence properties. The efficacy of MultDMD is demonstrated through several examples, including the nonlinear pendulum, the Lorenz system, and fluid dynamics data, where we demonstrate its remarkable robustness to noise.

This paper examines algorithmic nonlinear dictionary enhancements for Dynamic Mode Decomposition (DMD) based analysis of an industrially relevant flow problem. Several nonlinear variants of DMD are used to analyze flow interacting with two tandem cylinders in a bi-stable configuration. The study investigates the coexistence of two vortex-shedding regimes, distinguished by their characteristic frequencies. Using DMD, we extract coherent structures associated with each regime and consider their relation to surface pressure and far-field noise at these frequencies. Simultaneous experimental measurements of velocity fields, unsteady surface pressure, and far-field acoustic pressure are performed using high-speed Particle Image Velocimetry (PIV), a remote microphone technique, and a far-field microphone. Two DMD algorithms are applied to investigate the flow. First, we apply kernelized extended DMD (kEDMD) to demonstrate how a nonlinear dictionary can improve modal reliability compared with the usual linear alternative. The extracted modes are validated by analyzing their modal residuals (ResDMD), which are further used to explore phenomena at other frequencies and harmonics while examining the impact of dictionary depth and characteristics. Second, we study the flow physics of each shedding frequency via the Rigged DMD algorithm, which provides detailed spectral insights for dynamical systems with continuous spectra. A new kernelized Rigged DMD (kRigged DMD) is introduced, showing improved accuracy in resolving generalized eigenfunctions and spectral measures compared to linear dictionaries.

Sharpness is an almost generic assumption in continuous optimization that bounds the distance from minima by objective function suboptimality. It leads to the acceleration of first-order methods via restarts. However, sharpness involves problem-specific constants that are typically unknown, and restart schemes typically reduce convergence rates. Moreover, such schemes are challenging to apply in the presence of noise or approximate model classes (e.g., in compressive imaging or learning problems), and typically assume that the first-order method used produces feasible iterates. We consider the assumption of approximate sharpness, a generalization of sharpness that incorporates an unknown constant perturbation to the objective function error. This constant offers greater robustness (e.g., with respect to noise or relaxation of model classes) for finding approximate minimizers. By employing a new type of search over the unknown constants, we design a restart scheme that applies to general first-order methods and does not require the first-order method to produce feasible iterates. Our scheme maintains the same convergence rate as when assuming knowledge of the constants. The rates of convergence we obtain for various first-order methods either match the optimal rates or improve on previously established rates for a wide range of problems. We showcase our restart scheme on several examples and point to future applications and developments of our framework and theory.

The first step when solving an infinite-dimensional eigenvalue problem is often to discretize it. We show that one must be extremely careful when discretizing nonlinear eigenvalue problems. Using examples from the NLEVP collection, we demonstrate that discretization can lead to several issues, including: (1) introduction of spurious eigenvalues, (2) omission of spectra, (3) severe ill-conditioning, and (4) emergence of ghost essential spectra. While many eigensolvers are available for solving finite matrix nonlinear eigenvalue problems, we propose InfBeyn, a solver for general holomorphic infinite-dimensional nonlinear eigenvalue problems that circumvents these discretization issues. We prove that InfBeyn is stable and converges. Furthermore, we provide an algorithm that computes the problem's pseudospectra with explicit error control, enabling verification of computed spectra. Both algorithms and numerical examples are publicly available in the \texttt{infNEP} software package, which is written in MATLAB.

We study the convergence of Hermitian Dynamic Mode Decomposition (DMD) to the spectral properties of self-adjoint Koopman operators. Hermitian DMD is a data-driven method that approximates the Koopman operator associated with an unknown nonlinear dynamical system, using discrete-time snapshots. This approach preserves the self-adjointness of the operator in its finite-dimensional approximations. We prove that, under suitably broad conditions, the spectral measures corresponding to the eigenvalues and eigenfunctions computed by Hermitian DMD converge to those of the underlying Koopman operator. This result also applies to skew-Hermitian systems (after multiplication by \(i\)), applicable to generators of continuous-time measure-preserving systems. Along the way, we establish a general theorem on the convergence of spectral measures for finite sections of self-adjoint operators, including those that are unbounded, which is of independent interest to the wider spectral community. We numerically demonstrate our results by applying them to two-dimensional Schrödinger equations.

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 a posteriori 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.

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.

Residual Dynamic Mode Decomposition (ResDMD) offers a method for accurately computing the spectral properties of Koopman operators. It achieves this by calculating an infinite-dimensional residual from snapshot data, thus overcoming issues associated with finite truncations of Koopman operators (e.g., Extended Dynamic Mode Decomposition), such as spurious eigenvalues. Spectral properties computed by ResDMD include spectra, pseudospectra, spectral measures, Koopman mode decompositions, and dictionary verification. In scenarios where the number of snapshots is fewer than the dictionary size, particularly for exact DMD and kernelized Extended DMD, ResDMD has traditionally been applied by dividing snapshot data into a training set and a quadrature set. We demonstrate how to eliminate the need for two datasets through a novel computational approach of solving a dual least-squares problem. We analyze these new residuals for exact DMD and kernelized Extended DMD, demonstrating ResDMD’s versatility and broad applicability across various dynamical systems, including those modeled by high-dimensional and nonlinear observables. The utility of these new residuals is showcased through three diverse examples: the analysis of a cylinder wake, the study of airfoil cascades, and the compression of transient shockwave experimental data. This approach not only simplifies the application of ResDMD but also extends its potential for deeper insights into the dynamics of complex systems.

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 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\times10^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.

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.

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 well-conditioned 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 integers \(K > 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.

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 beam equations.

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.

Sharpness conditions directly control the recovery performance of restart schemes for first-order optimization methods without the need for restrictive assumptions such as strong convexity. However, they are challenging to apply in the presence of noise or approximate model classes (e.g., approximate sparsity). We provide a first-order method: Weighted, Accelerated and Restarted Primal-dual (WARPd), based on primal-dual iterations and a novel restart-reweight scheme. Under a generic approximate sharpness condition, WARPd achieves stable linear convergence to the desired vector. Many problems of interest fit into this framework. For example, we analyze sparse recovery in compressed sensing, low-rank matrix recovery, matrix completion, TV regularization, minimization of \(\|Bx\|_{l^1}\) under constraints (\(l^1\)-analysis problems for general \(B\)), and mixed regularization problems. We show how several quantities controlling recovery performance also provide explicit approximate sharpness constants. Numerical experiments show that WARPd compares favorably with specialized state-of-the-art methods and is ideally suited for solving large-scale problems. We also present a noise-blind variant based on a Square-Root LASSO decoder. Finally, we show how to unroll WARPd as neural networks. This approximation theory result provides lower bounds for stable and accurate neural networks for inverse problems and sheds light on architecture choices. Code and a gallery of examples are available online as a MATLAB package.

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 length scales 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 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 parallelisable, 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 \(\mathcal{O}(N)\). 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.

Modal decomposition techniques are widely employed to analyze turbulent flows by identifying coherent structures, which are often linked to key flow features, such as unsteady loading and pressure fluctuations. Modal decomposition offers a powerful framework for uncovering these structures and interpreting their role in flow dynamics. The most popular methods include POD (proper orthogonal decomposition), SPOD (spectral proper orthogonal decomposition), and DMD (dynamic mode decomposition). In this paper, we explore the flow interaction with tandem cylinders arranged in a bi-stable configuration. The velocity field, unsteady surface pressure, and far-field acoustic pressure of a tandem cylinder configuration are measured simultaneously using a high-speed PIV system, a remote microphone technique, and a far-field microphone, respectively. We utilize POD, SPOD and kEDMD (kernelized extended DMD) to identify coherent structures in the vorticity field. These three methods are compared to relate the vertical structures associated with the key flow regimes to surface and far-field pressure measurements. Results show that POD and SPOD offer insight into the important frequencies in an efficient manner, whilst the versatility of DMD-based methods offers further clarity into the dynamics of the system via the Koopman modes.

In this paper, we investigate large-scale superstructures over a k-type homogeneous rough wall (k_s+ = 209) at a speed of 1.3 × 10^6 Rem−1, under varying pressure gradient conditions characterized by Clauser parameters ranging from beta = −0.6 (Favorable Pressure Gradient) to beta = 0.98 (Adverse Pressure Gradient). Time-resolved stereoscopic wall-parallel PIV measurements, combined with two-point spatial and temporal correlation analyses, reveal the presence of predominantly streamwise-aligned large-scale structures. While these correlations qualitatively capture the spatial and temporal coherence of the flow, only subtle variations are observed across different pressure gradient conditions. To gain deeper insights, kernelized Residual Dynamic Mode Decomposition (kResDMD) was employed to extract coherent modes that are not superstructures themselves but served as the skeletal representation of these structures. This modal characterization enabled us to quantify pressure gradient effects on structure distribution, organization, orientation, and frequency content, providing a more comprehensive description of superstructure dynamics in rough-wall turbulent boundary layers.

The Hypothalamic-Pituitary-Adrenal (HPA) axis is a major neuroendocrine system, and its dysregulation is implicated in various diseases. This system also presents interesting mathematical challenges for modeling. We consider a nonlinear delay differential equation model and calculate pseudospectra of three different linearizations: a time-dependent Jacobian, linearization around the limit cycle, and dynamic mode decomposition (DMD) analysis of Koopman operators (global linearization). The time-dependent Jacobian provided insight into experimental phenomena, explaining why rats respond differently to perturbations during corticosterone secretion's upward versus downward slopes. We developed new mathematical techniques for the other two linearizations to calculate pseudospectra on Banach spaces and apply DMD to delay differential equations, respectively. These methods helped establish local and global limit cycle stability and study transients. Additionally, we discuss using pseudospectra to substantiate the model in experimental contexts and establish bio-variability via data-driven methods. This work is the first to utilize pseudospectra to explore the HPA axis.

In recent years, algorithm unrolling has emerged as a promising methodology in various signal-processing applications, intending to combine the strengths of iterative algorithms and deep learning. Despite its success, unrolling greedy sparse recovery algorithms, particularly Orthogonal Matching Pursuit (OMP), has yet to receive much attention. The primary challenge is the non-differentiable nature of the argsort operator, a key component in greedy algorithms, which hinders gradient backpropagation during training. To address this, we introduce ‘OMPNet’ by utilizing softsorting to approximate the argsort operator in a differentiable manner. Our numerical and theoretical analysis shows that under certain conditions, the approximation error is minimal, and the performance of the approximated OMP, which we call ‘Soft-OMP’, closely matches the original. We then incorporate Soft-OMP into the feedforward neural network’s layers, integrating learnable weight parameters to connect our approach to weighted sparse recovery. Our numerical results demonstrate that this network is trainable and can surpass the performance of the original OMP in certain scenarios.

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.

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.

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.