Computational approximation theory on higher-genus domains

Abstract

We develop computational approximation theory on higher-genus domains, finite unions of disjoint subintervals of the real axis that arise naturally in orthogonal polynomial theory, iterative numerical linear algebra, and the study of integrable nonlinear wave equations. In each case, exploiting this underlying structure to choose appropriate basis functions (orthogonal polynomials) is essential to obtain a fast and efficient numerical method. A unifying theme is the numerical solution of Riemann--Hilbert problems on these domains. First, we develop a numerical framework and apply it to compute orthogonal polynomials. For Chebyshev-like weight functions supported on a union of disjoint intervals, we use the Fokas--Its--Kitaev Riemann--Hilbert problem~\cite{Fokas1992} to produce an algorithm that computes recurrence coefficients and allows for the pointwise evaluation of Cauchy transforms in optimal ($\mathrm{O}(1)$) complexity. Using these orthogonal polynomials, we develop the Akhiezer iteration, a new iterative method for solving indefinite linear systems and computing matrix functions. The iteration uses orthogonal polynomial recurrence coefficients to efficiently compute the action of a matrix polynomial on a vector. It is inner-product-free, features an a priori computable convergence rate, and is often faster in practice than standard Krylov subspace methods. A particular extension yields a fast inverse-free iteration for solving large-scale Sylvester matrix equations. Finally, we present a class of numerical methods for computing and evolving soliton gasses, ensembles of a very large or infinite number of interacting solitons, in the context of the Korteweg--de Vries and modified Korteweg--de Vries equations. While traditionally difficult to simulate due to ill-conditioning, the methods are based on the Riemann--Hilbert formulation of~\cite{primitive-potentials} and the analysis of~\cite{Girotti2021,Girotti2023}. By using basis functions that encode the correct endpoint asymptotics, this approach enables the first accurate simulation of such solutions in standard (double) precision arithmetic. This allows for fast, pointwise evaluation before asymptotic estimates are valid.