Schedule

Oscillatory integrals arise in a broad range of models for wave-based phenomena. The method of Numerical Steepest Descent (NSD) combines complex contour deformation with numerical quadrature to provide an efficient and accurate approach for evaluating such integrals. Unless the phase function governing the oscillation is particularly simple, the application of Numerical Steepest Descent requires a significant amount of a priori analysis and expert user input. In particular, naive implementations of NSD are known to break down in the presence of coalescing saddle points, which occur in various applications — for example, approximation of special functions, such as the Airy Functions and the Catastrophe Integrals. PathFinder is an algorithm implemented in Matlab/Octave on top of C code, which automatically deforms the integration range/contour and applies a modified version of NSD to evaluate integrals of the form $$ I = \int_a^b f(x) \exp(\mathrm{i}\omega g(x))~\mathrm{d}x, $$ where the endpoints $a$ and $b$ can be complex-valued, even infinite; $\omega>0$ is the frequency parameter, $f$ is an entire function and $g$ is a polynomial. PathFinder remains accurate with a bounded computational cost for all $\omega>0$ and is robust in the presence of coalescing saddle points. It can be easily used by non-experts, without understanding the underlying complex analysis. In this talk, I will begin by introducing NSD, explaining when it works well and when it breaks down. I will then explain why the PathFinder algorithm works when NSD fails.

Interference is one of the most universal phenomena in nature. As ubiquitous as interference and interference patterns are, they are generally hard to compute. The oscillatory integrals involved are often only conditionally and not absolutely convergent, meaning they converge slowly and artefacts such as dependence on unphysical cutoffs may be hard to avoid. Likewise, if the integrals are performed iteratively, as is often the only practicable method, conditional convergence is in general insufficient to guarantee uniqueness, since the order in which partial integrals are taken can affect the result. Using Picard-Lefschetz theory, a general, exact approach to multidimensional oscillatory integrals based upon saddle point and steepest descent techniques, I will present ways to define and evaluate multi-dimensional oscillatory integrals efficiently in the complex plane. I will demonstrate its use in several applications.

Periodic operators are a fundamental object of study in many application fields, notably condensed matter physics. Their spectral properties are conveniently turned into integrals of periodic functions over the Brillouin zone using the Bloch-Floquet transform. Of interest are ground state properties of insulators (smooth integrands), metals (discontinuous integrands) and spectral properties of materials (singular integrands). I will review the broader context, and present in particular recent results on an integration method based on a multidimensional generalization of contour deformation.

We consider the application of the generalized Convolution Quadrature (gCQ) to approximate the solution of an important class of sectorial problems. The gCQ is a generalization of Lubich's Convolution Quadrature (CQ) that allows for variable steps. In many applications of interest the data is not smooth enough and the original CQ, with uniform steps, is known to present an order reduction close to the singularity. We present recent results about the stability and convergence of the gCQ in the low regularity setting and show how to choose optimal time meshes according to the behaviour of the data. An important advantage of the gCQ method is that it allows for a fast and memory-reduced (oblivious) implementation. The fast and oblivious gCQ relies on efficient special quadratures for appropriate contour integral representations of the numerical solution. We illustrate our results with several numerical experiments.

Despite the development of increasingly successful methods for the approximation of highly oscillatory integrals in recent years, such methodologies are arguably not yet as versatile as methods for classical quadrature. In this talk we will focus on improving the Filon and the Levin method which are amongst the most successful methods for approximating such integrals. Filon methods generally are limited by the moment-problem - the necessity to efficiently compute frequency-dependent quadrature moments. In particular this restricts their applicability in the case of singular, oscillatory integrals appearing for example in high-frequency wave scattering problems. We will address this problem by exploiting the banded nature of certain differential operators when acting on relevant interpolation bases to compute Filon moments rapidly for a broader class of oscillators. Moment-free methods such as the Levin method are a suitable alternative to this problem, but in turn require the solution of a dense collocation system and therefore are often computationally fairly expensive. Repeatedly exploiting the banded structure of a modified Levin operator allows us to significantly accelerate the solution of this collocation system in Levin–Clenshaw–Curtis methods. This approach can be combined with existing approaches improving the versatility of Levin methods thus taking a step towards a flexible and efficient technique for highly oscillatory quadrature. This is joint work with Arieh Iserles and Nigel Peake.

One of the key ingredients for computing special functions is integral representations (when available). In many cases, these integrals are highly oscillatory, and methods inspired by asymptotic analysis can be used to obtain stable representations. This concept is illustrated with some examples.

We will illustrate how numerical integration can be used to: (1) compute the complicated coefficients in uniform asymptotic expansions in a stable manner; (2) determine the location of zeros, poles or other singularities that appear on the boundary of the sections of validity of asymptotic approximations.

Catastrophe theory leads to a well-established hierarchy of integrals that represent the local behaviour of many physical phenomena where contributions are dominated by contributions from coalescences of features of their integrands. Normally such integrals involve exponents of poly/multi-nomial functions. In this talk we review work over the past few years that has found that the presence of logarithmic terms in these exponents is both common and essential in understanding certain wave problems near event horizons. Such inclusions gives rise to subtle, but significant changes to changes to asymptotic evaluation of the catastrophe integrals, with notable physical consequences. We provide examples from bores, acoustics and Bose-Einstein condensates.

In this talk, I will give an overview of recent developments linking wave theory and multidimensional complex analysis. I will explain how a procedure of complex deformation of the integration surface of Fourier-like highly oscillatory double integrals can lead to closed-form far-field asymptotics results in wave diffraction theory. Each far-field component will be shown to be connected to a special point on the singularity set of the integrand. The procedure will be illustrated through the three-dimensional problem of plane wave diffraction by a quarter-plane and the two-dimensional problem of plane wave diffraction by a penetrable wedge. I will also show how it can be used to shed some light on wave propagation in periodic structures. Moreover, I’ll discuss how the surface deformation procedure can be used explicitly for numerical integration of double oscillatory integrals and I’ll mention a link between this procedure and the Hybrid-Numerical-Asymptotic (HNA) boundary element method.

A 2D Fourier integral is studied, whose transformant is a function having algebraic growth at infinity and holomorphic everywhere except some polar and branching sets of complex codimension 1. The integration surface is the real plane, possibly slightly shifted to avoid the singularities of the transformant. The aim of the talk is to build the deformation of the integration surface into some other surface to make the exponential factor of the Fourier integral decaying as fast as possible on it. According to the multidimensional Cauchy's theorem, a deformation (homotopy) not hitting the singular sets does not change the value of the integral. At the same time, a proper deformation makes the integral more suitable for numerical evaluation or for asymptotical investigation. A general procedure of deformation of the integration surface is proposed in the talk. The resulting surface is a sum of components stemming from the special points of the singularities of the transformant, such as the saddle on singularities or the crossings (see [1]). The components have topological structure dictated by the special point type. Thus, the field becomes splitted explicitly into terms corresponding to the special points of the singularities. The work is being done in collaboration with Raphael C. Assier and Andrey I. Korolkov from the University of Manchester. References 1. R.C.Assier, A.V.Shanin, A.I.Korolkov, A contribution to the mathematical theory of diffraction. Part I: A note on double Fourier integrals // QJMAM, 76(2):211-241, 2022.

I will discuss a class of solvers for linear scalar ordinary differential equations which run in time bounded independent of frequency. They operate by producing exponential representations of a basis in the space of solutions of the equation. These exponential representations can be used to rapidly evaluate any desired solution of the differential equation at any point in the solution domain with accuracy on the order of the condition number of the problem. I will also present a theorem which bounds the complexity of the exponential representations as a function of a measure of the complexity of the equation's coefficients. If time permits, I will discuss how these methods can be used to solve certain systems of linear ordinary differential equations in time independent of frequency, as well as various applications of this work. Applications I might discuss include the rapid computation of barycentric interpolation nodes and weights, a method for the computation of the first $N$ eigenvalues and eigenfunctions of a singular Sturm-Liouvile problem in $O(N)$ time and methods for the rapid application of special function transforms. The theorem I will present was developed in collaboration with Vladimir Rokhlin (Yale), the work on the rapid application of special function transforms was is joint with Haizhao Yang (U. of Maryland) and the work on Sturm-Liouville problems is joint with my student, Richard Chow (U. of Toronto).

We investigate, by a mixture of analysis and rigorous computation, a conjecture of Carlos Kenig, that the essential spectral radius of the classical double-layer potential integral operator on the boundary of a general Lipschitz domain is $< 1/2$. Our computations are for 2D domains that are piecewise analytic and locally dilation invariant in the sense of Chandler-Wilde, Hagger, Perfekt, and Virtanen, Numer. Math. (2023). This domain class permits infinitely many oscillations of arbitrarily large amplitude adjacent to singular points on the boundary. Such domains are promising candidates as counter-examples since it is already known (Chandler-Wilde and Spence, Numer. Math. 2022) that, for these domains, the essential numerical radius can be arbitrarily large. Via localisation results and a Floquet-Bloch transform we reduce the computation of the essential spectral radius to a problem of estimating the spectral radii of a family of integral operators with periodic analytic kernels, for which the trapezium rule and associated error bounds are perfectly suited.

Gaussian quadrature is typically associated with orthogonal polynomials. In this talk we explore different perspectives. We start by showing that orthogonality is not an essential aspect of Gaussian quadrature and that in fact Gaussian quadrature rules (with positive weights, points inside the integration interval and only half as many points as expected) can be devised for a variety of non-polynomial spaces. This could include functions with singularities. Next, we find that the non-polynomial perspective carries over to higher-dimensional integrals more easily than that of orthogonal polynomials, suggesting novel ways of constructing cubature rules for multivariate domains.

It is well known that matrices with low Hessenberg-structured displacement rank enjoy fast algorithms for certain matrix factorizations. We will show how n x n principal finite sections of the Gram matrix for the orthogonal polynomial measure modification problem has such a displacement structure, unlocking a collection of fast algorithms for computing connection coefficients (as a Cholesky factor) between a known orthogonal polynomial family and the modified family. In general, the $O(n^3)$ complexity is reduced to $O(n^2)$, and if the symmetric Gram matrix has upper and lower bandwidth b, then the $O(b^2n)$ complexity for a banded Cholesky factorization is reduced to $O(b n)$. In the case of modified Chebyshev polynomials, it is known that the Gram matrix has Toeplitz-plus-Hankel structure, and if the modified Chebyshev moments satisfy a decay condition, then an off-diagonal low-rank structure is observed in the Gram matrix, enabling a further reduction in the complexity of an approximate Cholesky factorization. Finally, for Koornwinder's class of bivariate orthogonal polynomials of total degree n, we will explore how the original displacement-based algorithms can be adapted to reduce an $O(n^6)$ complexity down to $O(n^4)$, and in case the Gram matrix has bandwidth $O(b n)$, from $O(b^2 n^4)$ down to $O(b n^3)$. These fast modification algorithms provide yet another computationally feasible avenue to investigate oscillatory and singular quadrature rules.

Hermite approximations are widely used in solving PDEs in unbounded domains. However, it is difficult to find some sharp convergence results for analytic functions. In this talk, I will introduce recent progress in rigorous error analysis of Hermite approximations for analytic functions.

We present and analyse quadrature rules for the evaluation of weakly singular integrals over self-similar fractal measures. A special case is where integration is with respect to an appropriate Hausdorff measure restricted to the fractal attractor of an iterated function system of contracting similarities satisfying the open set condition, such as a Cantor set or a Koch curve. For standard power law and logarithmic singularities we show that using the homogeneity properties of the integrand and the self-similarity of the measure, one can express singular integrals in terms of regular integrals with smooth integrands, which can then be evaluated using any appropriate quadrature rule, such as a composite barycentre rule or a higher order alternative. We apply our quadrature rules in the context of Galerkin integral equation methods for simulating acoustic scattering by sound-soft fractal scatterers. We present numerical results for a range of examples including Cantor sets and dusts, Koch curves and snowflakes, and Sierpinski triangles and tetrahedra. We make our software available as a Julia code. This is joint work with António Caetano (Aveiro), Simon Chandler-Wilde (Reading), Xavier Claeys (Sorbonne), Andrew Gibbs (UCL), Botond Major (UCL) and Andrea Moiola (Pavia).

In this talk we consider the problem of computing an integral of a restriction of a regular function to a self-similar set (e.g. Sierpinski sieve), with respect to the Hausdorff measure. This mathematical problem has applications notably in the boundary element methods for fractal antenna engineering, see Caetano et el. 2024, Hewett at al. 2023. Our method is based on two key ideas: 1) use of the values of the function outside of the self-similar set and 2) exploiting a self-similar structure of the set in order to evaluate the quadrature weights. The first idea enables us to use tensor product quadratures, which inherit their convergence properties from the underlying 1D quadratures. The second idea yields the characterization of the weights as an eigenvector of a certain matrix that encodes information about the quadrature points and self-similarity properties of the geometry. Depending on a structure of a self-similar set, the computed weights coincide with exact weights (integrals of Lagrange polynomials) or approximate them. Surprisingly, even the inexact weights integrate the $Pk$-polynomials exactly. We discuss further properties of the quadrature weights (e.g. questions related to the positivity) and illustrate our findings with numerical experiments. This is a joint work with Patrick Joly and Zois Moitier (INRIA, France).

A quadrature is an approximation of the definite integral of a function by a weighted sum of function values at specified points, or nodes, within the domain of integration. Gaussian quadratures are constructed to yield exact results for any polynomials of degree $2r-1$ or less by a suitable choice of r nodes and weights. Cubature is a generalization of quadrature in higher dimension. Constructing a cubature amounts to find a linear form $ \Lambda : R[x] \to R$, $p \mapsto \sum_{j=1}^r a_j \, p(\xi_j)$ from the knowledge of its restriction to $R[x]_{\leq d}$. The unknowns to be determined are the weights $a_j$ and the nodes $\xi_j$. An approach based on moment matrices was proposed in [1]. We give a basis-free version in terms of the Hankel operator $\mathcal{H}$ associated to $\Lambda$. The existence of a cubature of degree $d$ with $r$ nodes boils down to conditions of ranks and positive semidefiniteness on $\mathcal{H}$. We then recognize the nodes as the solutions of a generalized eigenvalue problem. Standard domains of integration are symmetric under the action of a finite group. It is natural to look for cubatures that respect this symmetry [2]. They are exact for all anti-symmetric functions beyond the degree of the cubature. Introducing adapted bases obtained from representation theory, the symmetry constraint allows to block diagonalize the Hankel operator $\mathcal{H}$. The size of the blocks is explicitly related to the orbit types of the nodes. From the computational point of view, we then deal with smaller-sized matrices both for securing the existence of the cubature and computing the nodes. [1] L. Fialkow and S. Petrovic. A moment matrix approach to multivariable cubature. Integral Equations Operator Theory, 52(1):85–124, 2005. [2] R. Cools. Constructing cubature formulae: the science behind the art. Acta numerica, 6:1–54, 1997. Joint work with Mathieu Collowald (Université Côte d'Azur & Inria) partly publised as Studies in Applied Mathematics, 141:501-546 (2018). doi:10.1111/sapm.12240 See also the preprint https://hal.science/hal-01188290

An effective method for computing singular integrals with logarithmic and Cauchy kernels is to expand the integrand in a basis of weighted orthogonal polynomials and to use simple three-term recurrences that the integral operators applied to the basis satisfy. Using the classical theory of minimal solutions we can compute these extremely fast for any point of evaluation, including arbitrary close to the integration domain. This utilises a simple algorithm built by carefully choosing between forward recurrence and (F. W. J.) Olver’s algorithm. We review this approach for logarithmic kernels, extend it to power law kernels, discuss application to more general kernels like Bessel kernel, and finally discuss extension to logarithmic kernels on rectangles.

We devise a solver for a class of equations, posed on $\mathbb{R}^d, d∈{1,2}$, involving the fractional Laplacian – a singular integral operator that has generated exceptional interest due to its nonlocal properties, particularly in equations capturing anomalous diffusion. By leveraging exceptionally elegant formulae for the fractional Laplacian applied to weighted classical orthogonal polynomials, one can diagonalize the equation. The problem is then reduced from a solve of a nonlocal equation to the expansion of the right-hand side in a non-standard family of functions. However, by appealing to frame theory, we show that a robust technique for finding the expansion coefficients is via a truncated SVD projection. We apply our solver to the fractional heat equation (utilizing up to a 6th-order Runge–Kutta time discretization), a fractional heat equation with a time-dependent exponent $s(t)$ in the fractional Laplacian, and a two-dimensional problem. We observe spectral convergence in the spatial dimension for sufficiently smooth data.

A large toolbox of numerical schemes for dispersive equations has been established, based on different discretization techniques such as discretizing the variation-of-constants formula (e.g., exponential integrators) or splitting the full equation into a series of simpler subproblems (e.g., splitting methods). In many situations these classical schemes allow a precise and efficient approximation. This, however, drastically changes whenever non-smooth phenomena enter the scene such as for problems at low regularity and high oscillations. Classical schemes fail to capture the oscillatory nature of the solution, and this may lead to severe instabilities and loss of convergence. In this talk I present a new class of resonance based schemes. The key idea in the construction of the new schemes is to tackle and deeply embed the underlying nonlinear structure of resonances into the numerical discretization. As in the continuous case, these terms are central to structure preservation and offer the new schemes strong geometric properties at low regularity.

Layer potentials represent solutions to partial differential equations in an integral equation formulation. When numerically evaluating layer potentials at target points close to the domain boundary, specialized quadrature techniques are required for accuracy because of the rapid variations in the integrand. With the goal of efficiently achieving results within a specified error tolerance, we introduce an adaptive quadrature method, with automatic parameter adjustment, facilitated through error estimation. This approach is tailored for axisymmetric surfaces, employing a trapezoidal rule in the azimuthal angle and a Gauss-Legendre quadrature rule in the polar angle. Notably, while each surface must be axisymmetric, the layer density need not be, thus enabling application to complex geometries featuring multiple axisymmetric shapes. The adaptive quadrature method utilizes a so-called interpolatory semi-analytical quadrature in conjunction with a singularity swap technique in the azimuthal angle. In the polar angle, such a technique is utilized as needed, depending on the integral kernel, in combination with an adaptive subdivision of the integration interval. Estimates for both quadrature and interpolation errors are derived by complex analysis to complement existing error estimates. For each evaluation point close to the surface, an error estimate is used to to determine if the special quadrature method is needed to meet the error tolerance, and if so, to set the adaptive subdivision of panels, and at some instances, to decide on which approach should be utilized to evaluate a certain quantity. Numerical examples are presented to elucidate the method's efficacy.

We present a method for accurately solving boundary integral equations on implicitly defined surfaces in $\mathbb{R}^d$. The method relies on combining a dimension-independent technique for generating a high-order surface quadrature on level-set surfaces, with the general-purpose density interpolation method for handling the singular and nearly singular integrals ubiquitous in boundary integral formulations. The proposed methodology, based on a Nystrom discretization scheme, bypasses the need for generating a body-conforming mesh for the implicit surface, allowing in principle for an efficient coupling between a robust level-set representation of deforming surfaces, and boundary integral equation solvers. Particular attention will be paid to the computation of singular integrals when only a surface quadrature is available (i.e. in the absence of an actual mesh). We believe such techniques could prove useful in applications involving microscopic flows governed by the Stokes equations; in particular, the simulation of micro-swimmers and droplet microfluidics.

In studying higher order finite-size corrections to various limit laws at the soft edge, we discovered a particularly simple, yet perplexing structure: the correction terms are linear combinations of higher order derivatives of the limit law with rational polynomial coefficients. Given the highly nonlinear nature of perturbations of operator determinants, such a result is unexpected and points to a further layer of integrability. Here, it is the unique solvability of certain rectangular linear systems over the ring of rational polynomials with the solution of Painlevé II and its derivative added as indeterminates. We will present the ideas and evidence in the simplest example, the large matrix limit of Gaussian ensembles.

A standard step in electronic-structure calculations for crystals is the computation of Wannier functions, which describe localized molecular orbitals. From a mathematical viewpoint, this amounts to finding a common phase to smoothly "glue" together eigenfunctions of a parameter-dependent Schrodinger operator. In this talk, I will describe a new procedure for doing this for 1D crystalline material. This is joint work with Hanwen Zhang.

Enhancing the predictive accuracy of climate and weather models necessitates a more accurate representation of ice crystal scattering. Current models of ice crystal scattering exhibit inconsistencies across the electromagnetic spectrum, particularly in the far-infrared, microwave, and sub-millimetre ranges. To address this issue, the development of new ice crystal scattering models that adhere to observed mass-size and area-size power laws is proposed. These models aim to enhance the consistency of ice crystal optical properties across the spectrum, thereby rectifying the inconsistencies present in current models. We have employed the boundary element method in the development of such a model. This paper will provide examples of the application of this new model and discuss the need for more realistic representations of scattering by surface-roughened ice crystals and multi-phase particles. Here, we underscore the importance of improving ice crystal scattering models to enhance the reliability and accuracy of climate and weather predictions.

In this talk I will review an emerging approach to spectral methods for time-dependent PDEs, using T-systems and W-systems. This leads to a raft of quadrature problems, often involving singularities and high oscillation.

Analytic continuation is impossible or at least highly ill-conditioned in theory, depending how you frame the problem, yet sometimes surprisingly tractable in practice. The outline of the talk will be as follows:

1. Theorems (pretty gloomy)
2. Blending functions and "fat branch cuts"
3. Experiments (pretty encouraging)
4. The one-wavelength principle
5. The Schwarz function