Search

This paper presents a novel algorithm, based on use of rational approximants of randomly scalarized boundary integral resolvents, for the evaluation of acoustic and electromagnetic resonances in open and closed cavities; for simplicity we restrict treatment to cavities in two-dimensional space. The desired open cavity resonances (also known as ``eigenvalues'' for interior problems, and ``scattering poles'' for exterior and open problems) are obtained as the poles of associated rational approximants; both the approximants and their poles are obtained by means of the recently introduced AAA rational-approximation algorithm. In fact, the proposed resonance-search method applies to any nonlinear eigenvalue problem (NEP) associated with a given function $F: U \to \mathbb{C}^{n\times n}$, wherein a complex value $k$ is sought for which $F_kw = 0$ for some nonzero $w\in \mathbb{C}^n$. For the cavity problems considered in this paper, $F_k$ is taken as a spectrally discretized version of a Green function-based boundary integral operator at spatial frequency $k$. In all cases, the scalarized resolvent is given by an expression of the form $u^* F_k^{-1} v$, where $u,v \in \mathbb{C}^n$ are fixed random vectors. A variety of numerical results are presented for both scattering resonances and other NEPs, demonstrating the accuracy of the method even for high frequency states.

The AAA algorithm, introduced in 2018, computes best or near-best rational approximations to functions or data on subsets of the real line or the complex plane. It is much faster and more robust than previous algorithms for such problems and has been used in many applications since its appearance, including the numerical solution of Laplace, Poisson, and biharmonic PDE problems in irregular domains. AAA has also been extended in new directions and seems likely to be a tool of lasting importance in the future.

In this dissertation we have applied Trefethen and Gopal's Lightning Method to solve the Helmholtz equation in the exterior of two dimensional piecewise smooth domains. The background theory motivating the method is presented, and we explore the optimal method implementation for the unit square, which is subsequently used to give a guide on parameter selection for a general region. The behaviour of the computed solutions is verified to act in accordance with our intuition and current understanding of wave propagation, and we show that the wave decays to approximately 0 in the shadow region., Comment: 51 pages, 42 figures

Low Reynolds number fluid flows are governed by the Stokes equations. In two dimensions, Stokes flows can be described by two analytic functions, known as Goursat functions. Brubeck and Trefethen (2022) recently introduced a lightning Stokes solver that uses rational functions to approximate the Goursat functions in polygonal domains. In this paper, we present the "LARS" algorithm (Lightning-AAA Rational Stokes) for computing 2D Stokes flows in domains with smooth boundaries and multiply-connected domains using lightning and AAA rational approximation (Nakatsukasa et al., 2018). After validating our solver against known analytical solutions, we solve a variety of 2D Stokes flow problems with physical and engineering applications. Using these examples, we show rational approximation can now be used to compute 2D Stokes flows in general domains. The computations take less than a second and give solutions with at least 6-digit accuracy.

AAA rational approximation has normally been carried out on a discrete set, typically hundreds or thousands of points in a real interval or complex domain. Here we introduce a continuum AAA algorithm that discretizes a domain adaptively as it goes. This enables fast computation of high-accuracy rational approximations on domains such as the unit interval, the unit circle, and the imaginary axis, even in some cases where resolution of singularities requires exponentially clustered sample points, support points, and poles. Prototype MATLAB (or Octave) and Julia codes aaax, aaaz, and aaai are provided for these three special domains; the latter two are equivalent by a Moebius transformation. Execution is very fast since the matrices whose SVDs are computed have only three times as many rows as columns. The codes include a AAA-Lawson option for improvement of a AAA approximant to minimax, so long as the accuracy is well above machine precision. The result returned is pole-free in the approximation domain.

In this short, conceptual paper we observe that essentially the same mathematics applies in three contexts with disparate literatures: (1) sigmoidal and RBF approximation of smooth functions, (2) rational approximation of analytic functions near singularities, and (3) $hp$ mesh refinement for solution of PDEs. The relationship of (1) and (2) is as simple as the change of variables $s = \log(x)$, and our informal mnemonic for this relationship is ``sigmoid = log(ratapprox).''

Results on the rational approximation of functions containing singularities are presented. We build further on the ''lightning method'', recently proposed by Trefethen and collaborators, based on exponentially clustering poles close to the singularities. Our results are obtained by augmenting the lightning approximation set with either a low-degree polynomial basis or poles clustering towards infinity, in order to obtain a robust approximation of the smooth behaviour of the function. This leads to a significant increase in the achievable accuracy as well as the convergence rate of the numerical scheme. For the approximation of $x^\alpha$ on $[0,1]$, the optimal convergence rate as shown by Stahl in 1993 is now achieved simply by least-squares fitting.

M\"untz's theorem asserts, for example, that the even powers $1, x^2, x^4,\dots$ are dense in $C([0,1])$. We show that the associated expansions are so inefficient as to have no conceivable relevance to any actual computation. For example, approximating $f(x)=x$ to accuracy $\varepsilon = 10^{-6}$ in this basis requires powers larger than $x^{280{,}000}$ and coefficients larger than $10^{107{,}000}$. We present a theorem establishing exponential growth of coefficients with respect to $1/\varepsilon$.

We propose AAA rational approximation as a method for interpolating or approximating smooth functions from equispaced data samples. Although it is always better to approximate from large numbers of samples if they are available, whether equispaced or not, this method often performs impressively even when the sampling grid is fairly coarse. In most cases it gives more accurate approximations than other methods.

Often the easiest way to discretize an ordinary or partial differential equation is by a rectangular numerical method, in which n basis functions are sampled at m>>n collocation points. We show how eigenvalue problems can be solved in this setting by QR reduction to square matrix generalized eigenvalue problems. The method applies equally in the limit "m=infinity" of eigenvalue problems for quasimatrices. Numerical examples are presented as well as pointers to some related literature.

Gopal and Trefethen recently introduced "lightning solvers" for the 2D Laplace and Helmholtz equations, based on rational functions with poles exponentially clustered near singular corners. Making use of the Goursat representation in terms of analytic functions, we extend these methods to the biharmonic equation, specifically to 2D Stokes flow. Solutions to model problems are computed to 10-digit accuracy in less than a second of laptop time. As an illustration of the high accuracy, we resolve two or more counter-rotating Moffatt eddies near a singular corner.

A two-step method for solving planar Laplace problems via rational approximation is introduced. First complex rational approximations to the boundary data are determined by AAA approximation, either globally or locally near each corner or other singularity. The poles of these approximations outside the problem domain are then collected and used for a global least-squares fit to the solution. Typical problems are solved in a second of laptop time to 8-digit accuracy, all the way up to the corners, and the conjugate harmonic function is also provided. The AAA-least squares combination also offers a new method for avoiding spurious poles in other rational approximation problems, and for greatly speeding them up in cases with many singularities. As a special case, AAA-LS approximation leads to a powerful method for computing the Hilbert transform or Dirichlet-to-Neumann map.

The standard design principle for quadrature formulas is that they should be exact for integrands of a given class, such as polynomials of a fixed degree. We show how this principle fails to predict the actual behavior in four cases: Newton-Cotes, Clenshaw-Curtis, Gauss-Legendre, and Gauss-Hermite quadrature. Three further examples are mentioned more briefly.

This article is about both approximation theory and the numerical solution of partial differential equations (PDEs). First we introduce the notion of {\em reciprocal-log} or {\em log-lightning approximation} of analytic functions with branch point singularities at points $\{z_k\}$ by functions of the form $g(z) = \sum_k c_k /(\log(z-z_k) - s_k)$, which have $N$ poles potentially distributed along a Riemann surface. We prove that the errors of best reciprocal-log approximations decrease exponentially with respect to $N$ and that exponential or near-exponential convergence (i.e., at a rate $O(\exp(-C N / \log N))$) also holds for near-best approximations with preassigned singularities constructed by linear least-squares fitting on the boundary. We then apply these results to derive a "log-lightning method" for numerical solution of Laplace and related PDEs in two-dimensional domains with corner singularities. The convergence is near-exponential, in contrast to the root-exponential convergence for the original lightning methods based on rational functions., Comment: 20 pages, 11 figures

Rational approximations of functions with singularities can converge at a root-exponential rate if the poles are exponentially clustered. We begin by reviewing this effect in minimax, least-squares, and AAA approximations on intervals and complex domains, conformal mapping, and the numerical solution of Laplace, Helmholtz, and biharmonic equations by the "lightning" method. Extensive and wide-ranging numerical experiments are involved. We then present further experiments showing that in all of these applications, it is advantageous to use exponential clustering whose density on a logarithmic scale is not uniform but tapers off linearly to zero near the singularity. We give a theoretical explanation of the tapering effect based on the Hermite contour integral and potential theory, showing that tapering doubles the rate of convergence. Finally we show that related mathematics applies to the relationship between exponential (not tapered) and doubly exponential (tapered) quadrature formulas. Here it is the Gauss--Takahasi--Mori contour integral that comes into play.

Vandermonde matrices are exponentially ill-conditioned, rendering the familiar "polyval(polyfit)" algorithm for polynomial interpolation and least-squares fitting ineffective at higher degrees. We show that Arnoldi orthogonalization fixes the problem.

New algorithms are presented for numerical conformal mapping based on rational approximations and the solution of Dirichlet problems by least-squares fitting on the boundary. The methods are targeted at regions with corners, where the Dirichlet problem is solved by the "lightning Laplace solver" with poles exponentially clustered near each singularity. For polygons and circular polygons, further simplifications are possible.

Analytic continuation is ill-posed, but becomes merely ill-conditioned (although with an infinite condition number) if it is known that the function in question is bounded in a given region of the complex plane. In an annulus, the Hadamard three-circles theorem implies that the ill-conditioning is not too severe, and we show how this explains the effectiveness of Chebfun and related numerical methods in evaluating analytic functions off the interval of definition. By contrast, we show that analytic continuation is far more ill-conditioned in a strip or a channel, with exponential loss of digits of accuracy at the rate $\exp(-\pi x/2)$ as one moves along. The classical Weierstrass chain-of-disks method loses digits at the faster rate $\exp(-e\kern .3pt x)$.

Rational minimax approximation of real functions on real intervals is an established topic, but when it comes to complex functions or domains, there appear to be no algorithms currently in use. Such a method is introduced here, the {\em AAA-Lawson algorithm,} available in Chebfun. The new algorithm solves a wide range of problems on arbitrary domains in a fraction of a second of laptop time by a procedure consisting of two steps. First, the standard AAA algorithm is run to obtain a near-best approximation and a set of support points for a barycentric representation of the rational approximant. Then a "Lawson phase" of iteratively reweighted least-squares adjustment of the barycentric coefficients is carried out to improve the approximation to minimax.

A new method is introduced for solving Laplace problems on 2D regions with corners by approximation of boundary data by the real part of a rational function with fixed poles exponentially clustered near each corner. Greatly extending a result of D. J. Newman in 1964 in approximation theory, we first prove that such approximations can achieve root-exponential convergence for a wide range of problems, all the way up to the corner singularities. We then develop a numerical method to compute approximations via linear least-squares fitting on the boundary. Typical problems are solved in < 1s on a laptop to 8-digit accuracy, with the accuracy guaranteed in the interior by the maximum principle. The computed solution is represented globally by a single formula, which can be evaluated in tens of microseconds at each point.

New numerical algorithms based on rational functions are introduced that can solve certain Laplace and Helmholtz problems on two-dimensional domains with corners faster and more accurately than the standard methods of finite elements and integral equations. The new algorithms point to a reconsideration of the assumptions underlying existing numerical analysis for partial differential equations.

The traditional view in numerical conformal mapping is that once the boundary correspondence function has been found, the map and its inverse can be evaluated by contour integrals. We propose that it is much simpler, and 10-1000 times faster, to represent the maps by rational functions computed by the AAA algorithm. To justify this claim, first we prove a theorem establishing root-exponential convergence of rational approximations near corners in a conformal map, generalizing a result of D. J. Newman in 1964. This leads to the new algorithm for approximating conformal maps of polygons. Then we turn to smooth domains and prove a sequence of four theorems establishing that in any conformal map of the unit circle onto a region with a long and slender part, there must be a singularity or loss of univalence exponentially close to the boundary, and polynomial approximations cannot be accurate unless of exponentially high degree. This motivates the application of the new algorithm to smooth domains, where it is again found to be highly effective.

At the ANZIAM conference in Hobart in February, 2018, there were several talks on the solution of Laplace problems in multiply connected domains by means of conformal mapping. It appears to be not widely known that such problems can also be solved by the elementary method of series expansions with coefficients determined by least-squares fitting on the boundary. (These are not convergent series; the coefficients depend on the degree of the approximation.) Here we give a tutorial introduction to this method, which converges at an exponential rate if the boundary data are sufficiently well-behaved. The mathematical foundations go back to Runge in 1885 and Walsh in 1929. One of our examples involves an approximate Cantor set with up to 2048 components.

Let $E_{kk}^{(n)}$ denote the minimax (i.e., best supremum norm) error in approximation of $x^n$ on $[\kern .3pt 0,1]$ by rational functions of type $(k,k)$ with $k

Computing rational minimax approximations can be very challenging when there are singularities on or near the interval of approximation - precisely the case where rational functions outperform polynomials by a landslide. We show that far more robust algorithms than previously available can be developed by making use of rational barycentric representations whose support points are chosen in an adaptive fashion as the approximant is computed. Three variants of this barycentric strategy are all shown to be powerful: (1) a classical Remez algorithm, (2) a "AAA-Lawson" method of iteratively reweighted least-squares, and (3) a differential correction algorithm. Our preferred combination, implemented in the Chebfun MINIMAX code, is to use (2) in an initial phase and then switch to (1) for generically quadratic convergence. By such methods we can calculate approximations up to type (80, 80) of $|x|$ on $[-1, 1]$ in standard 16-digit floating point arithmetic, a problem for which Varga, Ruttan, and Carpenter required 200-digit extended precision., Comment: 29 pages, 11 figures

The trigonometric interpolants to a periodic function $f$ in equispaced points converge if $f$ is Dini-continuous, and the associated quadrature formula, the trapezoidal rule, converges if $f$ is continuous. What if the points are perturbed? With equispaced grid spacing $h$, let each point be perturbed by an arbitrary amount $\le \alpha h$, where $\alpha\in [\kern .5pt 0,1/2)$ is a fixed constant. The Kadec 1/4 theorem of sampling theory suggests there may be be trouble for $\alpha\ge 1/4$. We show that convergence of both the interpolants and the quadrature estimates is guaranteed for all $\alpha<1/2$ if $f$ is twice continuously differentiable, with the convergence rate depending on the smoothness of $f$. More precisely it is enough for $f$ to have $4\alpha$ derivatives in a certain sense, and we conjecture that $2\alpha$ derivatives is enough. Connections with the Fej\'er--Kalm\'ar theorem are discussed.

We introduce a new algorithm for approximation by rational functions on a real or complex set of points, implementable in 40 lines of Matlab and requiring no user input parameters. Even on a disk or interval the algorithm may outperform existing methods, and on more complicated domains it is especially competitive. The core ideas are (1) representation of the rational approximant in barycentric form with interpolation at certain support points and (2) greedy selection of the support points to avoid exponential instabilities. The name AAA stands for "adaptive Antoulas--Anderson" in honor of the authors who introduced a scheme based on (1). We present the core algorithm with a Matlab code and nine applications and describe variants targeted at problems of different kinds. Comparisons are made with vector fitting, RKFIT, and other existing methods for rational approximation.

A theorem is proved concerning approximation of analytic functions by multivariate polynomials in the $s$-dimensional hypercube. The geometric convergence rate is determined not by the usual notion of degree of a multivariate polynomial, but by the {\it Euclidean degree,} defined in terms of the 2-norm rather than the 1-norm of the exponent vector $\bf k$ of a monomial $x_1^{k_1}\cdots \kern .8pt x_s^{k_s}$.

Chebfun and related software projects for numerical computing with functions are based on the idea that at each step of a computation, a function $f(x)$ defined on an interval $[a,b]$ is "rounded" to a prescribed precision by constructing a Chebyshev series and chopping it at an appropriate point. Designing a chopping algorithm with the right properties proves to be a surprisingly complex and interesting problem. We describe the chopping algorithm introduced in Chebfun Version 5.3 in 2015 after many years of discussion and the considerations that led to this design.

Algorithms and underlying mathematics are presented for numerical computation with periodic functions via approximations to machine precision by trigonometric polynomials, including the solution of linear and nonlinear periodic ordinary differential equations. Differences from the nonperiodic Chebyshev case are highlighted., Comment: 20 pages

Although flow in a circular pipe is stable to infinitesimal perturbations, it can be excited to turbulence by finite perturbations whose minimal amplitude shrinks as R goes to infinity (R = Reynolds number). Laboratory experiments have appeared to disagree with one another and with theoretical predictions about the dependence of this minimal amplitude on $R$, with published results ranging approximately from $R^{-1/4}$ to $R^{-3/2}$. Here it is shown that these discrepancies can be explained by the use of different definitions of amplitude by different authors. An attempt is made to convert the existing results to a uniform definition of amplitude, the nondimensionalized $L^2$ definition common in the theoretical literature. Although subtleties in the physics raise some questions, agreement appears to be reached on a minimal amplitude that scales as $R^{-3/2 \pm 0.3}$.

There has been much recent interest, initiated by work of the physicists Hatano and Nelson, in the eigenvalues of certain random non-Hermitian periodic tridiagonal matrices and their bidiagonal limits. These eigenvalues cluster along a "bubble with wings" in the complex plane, and the corresponding eigenvectors are localized in the wings, delocalized in the bubble. Here, in addition to eigenvalues, pseudospectra are analyzed, making it possible to treat the non-periodic analogues of these random matrix problems. Inside the bubble, the resolvent norm grows exponentially with the dimension. Outside, it grows subexponentially in a bounded region that is the spectrum of the infinite-dimensional operator. Localization and delocalization correspond to resolvent matrices whose entries exponentially decrease or increase, respectively, with distance from the diagonal. This article presents theorems that characterize the spectra, pseudospectra, and numerical range for the four cases of finite bidiagonal matrices, infinite bidiagonal matrices ("stochastic Toeplitz operators"), finite periodic matrices, and doubly infinite bidiagonal matrices ("stochastic Laurent operators")., Comment: withdrawn