### Page Navigation

## Journal Articles

#### (*denotes first authors where applicable, others alphabetical)

(paper prizes shown in orange )

### Spectral Computations

Computing spectral measures of self-adjoint operatorsMatthew Colbrook, Andrew Horning, Alex Townsend SIAM Review, to appear. pdfUsing 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. |

Computing spectral measures and spectral typesMatthew Colbrook* Communications in Mathematical Physics, 2021. pdf journal linkSpectral 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. |

How to compute spectra with error controlMatthew Colbrook*, Bogdan Roman, Anders Hansen Physical Review Letters, 2019, 122(25), 250201. (chosen as cover article) pdf journal linkComputing 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. |

On the infinite-dimensional QR algorithmMatthew Colbrook*, Anders Hansen Numerische Mathematik, 2019, 143(1), 17-83. pdf journal linkSpectral 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. |

Pseudoergodic operators and periodic boundary conditionsMatthew Colbrook* Mathematics of Computation, 2020, 89(322), 737-766. pdf journal linkThere 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 ]\). |

### Spectral Methods

A Mathieu function boundary spectral method for diffraction by multiple variable poro-elastic plates, with applications to metamaterials and acousticsMatthew Colbrook*, Anastasia Kisil Proceedings of the Royal Society A, 2020, 476(2241), 20200184.
pdf journal linkMany 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. |

The unified transform for mixed boundary condition problems in unbounded domainsMatthew Colbrook*, Lorna Ayton, Athanassios Fokas Proceedings of the Royal Society A, 2019, 475(2222), 20180605. pdf journal linkThis 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. |

Extending the unified transform: curvilinear polygons and variable coefficient PDEsMatthew Colbrook* IMA Journal of Numerical Analysis, 2020, 40(2), 976-1004. pdf journal linkWe 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. |

A hybrid analytical-numerical technique for elliptic PDEsMatthew Colbrook*, Thanasis Fokas, Parham Hashemzadeh SIAM Journal on Scientific Computing, 2019, 41(2), A1066-A1090. pdf journal linkRecent 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. |

On the Fokas method for the solution of elliptic problems in both convex and non-convex polygonal domainsMatthew Colbrook*, Natasha Flyer, Bengt Fornberg Journal of Computational Physics, 2018, 374, 996-1016. pdf journal linkThere 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. |

### Numerical Solution of Scattering Problems

Do we need non-linear corrections? On the boundary Forchheimer equation in acoustic scatteringMatthew Colbrook*, Lorna Ayton Journal of Sound and Vibration, 2021, 495, 115905.
pdf journal linkThis 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. |

Reducing aerofoil-turbulence interaction noise through chordwise-varying porosityLorna Ayton*, Matthew Colbrook, Thomas Geyer, Paruchuri Chaitanya, Ennes Sarradj Journal of Fluid Mechanics, 2020, 906, p. A1. DOI: 10.1017/jfm.2020.746. pdf journal linkThis 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 ( |

Fast and spectrally accurate numerical methods for perforated screens(with applications to Robin boundary conditions) Matthew Colbrook, Matthew Priddin IMA Journal of Applied Mathematics, 2020, 85(5), 790-821. pdf journal linkThis 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. |

A spectral collocation method for acoustic scattering by multiple elastic platesMatthew Colbrook*, Lorna Ayton Journal of Sound and Vibration, 2019, 461, 114904. pdf journal linkThis 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. |

### Computational Fluid Mechanics

A hybrid analytical-numerical method for solving advection-dispersion problems on a half-lineFelipe de Barros, Matthew Colbrook, Athanassios Fokas International Journal of Heat and Mass Transfer, 2019, 139, 482-491. pdf journal linkThis 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. |

Scaling laws of passive-scalar diffusion in the interstellar mediumMatthew Colbrook*, Xiangcheng Ma, Philip Hopkins, Jonathan Squire Monthly Notices of the Royal Astronomical Society, 2017, 467(2), 2421-2429. pdf journal linkPassive-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. |

### Computation and Analysis of PDEs

Kernel density estimation with linked boundary conditionsMatthew Colbrook*, Zdravko Botev, Karsten Kuritz, Shev MacNamara Studies in Applied Mathematics, 2020, 145(3), 357-396. pdf journal linkKernel 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 given by the unified transform (or Fokas method). 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. |

## Refereed Conference Articles

On the existence of stable and accurate neural networks for image reconstructionMatthew Colbrook*, Vegard Antun, Anders Hansen Signal Processing with Adaptive Sparse Structured Representations (SPARS), 2019. pdf journal linkDeep 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. |

The unified transform: a spectral collocation method for acoustic scatteringLorna Ayton, Matthew Colbrook, Athanassios Fokas AIAA/CEAS Aeroacoustics, 2019. pdf journal linkThis 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. |

## Selected Submitted Articles/Preprints

The foundations of spectral computations via the solvability complexity index hierarchy: Part IMatthew Colbrook*, Anders Hansen arXiv pdf |

On the computation of geometric features of spectra of linear operators on Hilbert spacesMatthew Colbrook* |

Computing Spectra - On the solvability complexity index hierarchy and towers of algorithmsJonathan Ben-Artzi, Matthew Colbrook, Anders Hansen, Olavi Nevanlinna, Markus Seidel arXiv pdf |

Can stable and accurate neural networks be computed? - On the barriers of deep learning and Smale's 18th problemVegard Antun, Matthew Colbrook, Anders Hansen arXiv pdf |

Computing semigroups with error controlMatthew Colbrook* |

Stratified sampling based compressed sensing for structured signalsTheresa Loss*, Matthew Colbrook, Anders Hansen |

## Publication Profiles and Code

Researchgate

Some of my code can be found on my GitHub profile.

Code for \(\texttt{SpecSolve}\) can be found here.

Code for \(\texttt{FIRENETs}\) can be found here.