Numerical analytic continuation

Let f be an analytic function on a simply-connected compact continuum E of the complex z-plane. This might be an interval of the real line, where f might be real analytic. How can we calculate good estimates of the analytic continuation of f to other points z∈C\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z\in {\mathbb C}$$\end{document}? How can we estimate the locations of real or complex singularities of f? We review both the theory and the practice of some existing methods for these problems and propose that excellent results can be obtained from the computation of rational approximations of f by the AAA algorithm. In the case of analytic functions of two or more variables, the rational approximations are applied along line segments or other analytic arcs.


Introduction
Let f be an analytic function defined on a simply-connected compact continuum E in the complex z-plane, which might be an interval of the real axis.Analyticity means that at each point z ∈ E, the Taylor series for f exists and converges to f in a neighborhood of that point.We consider the problem, how can we estimate values of the analytic continuation of f at other real or complex points z ∈ C, or locate nearby real or complex singularities?
Mathematically, the values of an analytic continuation of f are determined throughout some open set ⊆ C; a thorough presentation can be found in [68, chap. 8, vol. 3].As explained in [58,103], the analytic continuation problem is ill-posed a priori, since the result depends discontinuously on the data, but becomes well-posed, though with infinite condition numbers, if we confine attention to a simply-connected subset of such as a disk or a strip where it is known that f is bounded.
In the literature of analytic continuation going back many years, the dominant numerical tool has been Padé approximation, that is, approximation by rational functions constructed to match initial terms of the Taylor series of f at a point [8].This can be an excellent method for analytic continuation of functions whose Taylor series are known analytically, but it quickly runs into trouble with functions that are just known by their values.A more robust method in many situations should be to approximate f by a rational function over a domain E bigger than a single point, but until recently, the tools available for such approximations were limited.This has changed with the appearance of the AAA algorithm [33,70], which works remarkably well as a black box near-best rational approximator on real or complex sets.In this paper, we recommend AAA approximation as a method for univariate and multivariate analytic continuation.Sections 2 and 3 present some examples, and Sect. 4 proposes the "one-wavelength principle": in 16-digit arithmetic, rational approximations can analytically continue an oscillatory function about one wavelength beyond its sampling domain.As a converse, Section 5 points out that unrelated function branches on disjoint domains can be interpolated by a single smooth analytic function, to 16-digit accuracy, if the domains are well separated.Sections 6-8 present aspects of the theory of rational approximation in the context of analytic continuation; these results are mainly not new but have not been collected together like this before.Section 9 applies the AAA method to various problems from the literature involving ordinary and partial differential equations (ODEs and PDEs).In particular several problems are examined of the numerical determination of singularities, related to fluid and solid dynamics, where there is an established tradition going back many years; see [9,15,19,27,60,67,89,90,98,107,108,112,115,116].
Section 10 turns to the analytic continuation of a real analytic bivariate function f (x, y) [62].There is not much literature on this topic, but we are aware of two relevant themes.One is the general problem of "function extension," in which, for any of a variety of purposes, it is desired to extend a function smoothly beyond its initial domain of definition.This literature uses many tools including partitions of unity and radial basis functions (RBFs), usually without an emphasis on analyticity; see [41,55,61,81,117].The other is the more specific topic of analytic continuation of Helmholtz fields, that is, solutions of the Helmholtz equation u + k 2 u = 0; see [1,10,11,25,63].The method we propose for bivariate problems is a new one.Given data f (x, y) on a region E ⊆ , one could attempt to compute a bivariate rational approximation r (x, y) and then use this for analytic continuation.An algorithm for such computations is presented in [2].However, at present there are no methods for bivariate rational approximation that have anything like the speed and flexibility of AAA approximation in the univariate case, and accordingly, our new method is based on AAA.Specifically, we propose the construction of rational approximations r to f defined on one-dimensional curves-typically line segments.This leads to a method that is extremely fast, and trivial to implement using the existing Chebfun aaa command [32] or the MATLAB or Julia codes available with [33].
Most of the existing literature on analytic continuation is related either to polynomials, whose effectiveness is limited, or to Padé approximants when it comes to rational functions.The view of the present paper is that Padé approximation is a very special case, often presenting accuracy problems requiring the use of extended precision arithmetic or other special tools; and that for many applications, a simpler starting point is to approximate on a continuum E rather than at a single point.
There are also a number of methods for numerical analytic continuation that are not based (at least explicitly) on polynomial or rational approximation.One fruitful idea is to impose some kind of regularization to combat ill-conditioning [28,37,84].It is also important to emphasize that the best method of analytic continuation may be analytic rather than numerical, and a review of the subject from this point of view can be found in Chapter 3 of [38]. 1 A famous example is the Riemann zeta function, whose representation ζ(z) = ∞ k=1 k −z converges only for Rez > 1 but which has other representations valid for other values z ∈ C (see Fig. 4

below).
We will not give details of the AAA algorithm here, because for the purposes of this paper, it is a black box [33,70], and if another algorithm with similar speed and approximation power came along, we would expect it to be equally useful for analytic continuation.We will only mention the basic feature that AAA involves manipulations (of so-called barycentric support points and barycentric coefficients) only on the approximation domain, nowhere else in the complex plane.In particular, unlike some algorithms, it does not work with poles of rational approximations; the poles are only computed after the rational approximation is obtained.Thus, for example, the approximate symmetries of poles in Figures 1, 2, 5, and 6 are not explicitly imposed on the approximations, but are just a result of their near-optimality.This is in contrast to most of the theory of rational approximation, to be discussed in some detail in Sects.6-8, which starts from an explicit consideration of poles by means of the Hermite integral formula, as developed originally by Walsh.Concerning alternatives to AAA, we note that for approximation on real intervals, the methods of Chebyshev-Padé (nonperiodic) and Fourier-Padé (periodic) approximation, or simply rational interpolation in Chebyshev or equispaced points, respectively, can be effective, as can the gold standard of minimax (best supremumnorm) approximation; see [51,115,116] and the Chebfun commands chebpade, ratinterp, and minimax.AAA approximation is exceptionally fast, robust, and easy to work with, however, and most of the figures presented in this paper are generated by codes just a few lines long that run in a fraction of a second.Most of our examples invoke AAA in its default mode as implemented in the aaa command of Chebfun, where the approximation degree is determined adaptively to achieve a relative supremum norm error tolerance of 10 −13 .In one case (Fig. 7) we apply AAA in 128-digit arithmetic in Julia.
Following standard terminology, by a rational function of type (m, n), we mean a function that can be written as a quotient p/q of polynomials of degrees (at most) m and n, respectively.A rational function of type (n, n) is also said to be of degree n.Representing rational functions by quotients of polynomials is numerically unstable, however, and one of the reasons AAA approximation is so robust and accurate is that it employs barycentric representations (and sometimes simpler partial fractions representations) instead.

Comparison of polynomials and rational functions
In this section we focus on the most common case for applications, where f is a real or complex analytic function defined on a real interval E, which without loss of generality we take to be [−1, 1].Such a function has a Chebyshev series that converges absolutely and uniformly on [−1, 1], where T k is the degree k Chebyshev polynomial, and analyticity implies that the coefficients a k decrease exponentially as k → ∞.Specifically, f must be continuable to a bounded analytic function in the Bernstein ellipse E ρ in the complex z-plane for some number ρ > 1, defined as the open set bounded by the ellipse with foci z = ±1 and semimajor and semiminor axis lengths summing to ρ.
which implies that the partial sums of the Chebyshev series converge at the same exponential rate, where The results of this paragraph are due to Bernstein in 1912 [14] (except the interpolation) and are presented in Chapter 8 of [102].
Bernstein also proved a converse to these results [102,Thm. 8.3].We display this as a theorem since it is can be interpreted as an assertion about numerical analytic continuation.
Theorem 1 (Analytic continuation by polynomials) Suppose f is any function defined on [−1, 1] for which there exist degree n polynomials p n , n = 0, 1, 2, . . ., for some ρ > 1 and C > 0. Then f can be analytically continued to an analytic function in the open Bernstein ellipse E ρ .Moreover, the polynomials p n converge exponentially to f throughout E ρ , uniformly on compact subsets.
In this theorem, and throughout the paper, we use the following terminology: The "big O" notation has its usual precise meaning as an upper bound as n → ∞.
Theorem 1 together with the bounds (3)-( 4) tells us that Chebyshev series, interpolants, or minimax approximants defined on [−1, 1] can be used for analytic continuation to every point within the smallest Bernstein ellipse on which f has a singularity.(A point z 0 ∈ C is a singularity of f if f can be analytically continued to points z arbitrarily close to z 0 but not to any neighborhood of z 0 .)For a discussion of the practical, numerical side of this effect in the context of Chebfun and "Chebfun ellipses," see [103, sec. 6].The phenomenon of convergence beyond the domain of approximation was called overconvergence by Walsh [111].It is a generalization of the familiar phenomenon of convergence of Taylor series throughout the disk bounded by the singularity of minimal modulus, which is the starting point of Weierstrass's "chain of disks" procedure used to define multi-valued analytic functions in complex analysis.(Chains of disks are mainly a theoretical idea, not of much use for computation [38,103]; see the beautiful illustrations of the idea in [113, sec. 3.6].) Outside the Bernstein ellipse, on the other hand, polynomial approximations determined from data on [−1, 1] will be useless.For Taylor series and disks, this observation is brought into focus by Jentzsch's theorem [57], which asserts that the zeros of the Taylor partial sums accumulate on the circle of convergence.For Chebyshev series on [−1, 1], a corresponding theorem involving zeros on the Bernstein ellipse is due to Walsh [110]; see [102, p. 140].
Rational approximations are an entirely different story, as illustrated in Fig. 1.Here the function f (z) = tanh(z) is approximated on [−1, 1] by a Chebyshev interpolant of degree 30 on the left and a AAA rational function of degree 7 on the right (or one can use minimax approximations; it makes little difference).The function f has poles at ±πik/2, where k ranges over the odd integers, and in the image on the left, we see that polynomial approximation gives nothing useful outside the Bernstein ellipse that passes through the first pair of these poles, at ±πi/2.The image on the right, by contrast, shows good approximation much further out, and in particular, the AAA poles approximate those of f .Six of the seven AAA poles lie at positions that would correspond to the values k = ±1, ±3, ±5 being adjusted to ±1.00000000057, ±3.0066, and ±5.99.Thus the first pair of singularities is located by the rational approximant to an accuracy of better than 10 −9 and the second to about 0.2%.(The seventh pole is at about 977.1i; it would be at ∞ for the minimax approximant, since f is an odd function.) Figure 2 presents another pair of images in the same format, but now for the function which in addition to poles at z = ±3i has an essential singularity at z = −2 and a branch point at z = 2. Despite this non-meromorphic structure, the rational approximation on [−1, 1] again does a good job of approximating f further out in the complex z-plane.Theorems 5 and 6 of Sect.7 will explain this success.The approximation r (z) is of degree 9, with two poles at −0.00014 ± 3.00012i, three poles clustering near the branch point at 3.99, 2.57, and 2.12, three poles near the essential singularity at −2.31 and −2.13 ± 0.20i, and one pole at 14.61.
One should not conclude from Figs. 1 and 2 that rational functions will always be more effective than polynomials for analytic continuation.As a counterexample, consider the amber function A(z) defined in [56] by the Chebyshev series   8), which we believe has a natural boundary along the Bernstein 2-ellipse.The polynomial approximation on the left has degree 46, and the rational approximation on the right has degree 23, so they have equal numbers of free parameters, 47.Inside the ellipse, rational approximations have no advantage over polynomials.Outside, the function f and thus the error f − r are not defined, so no color is plotted.The reason rational functions are so important for analytic continuation is that functions like this one are uncommon in applications In Sects.6-8 we will consider what can be said theoretically about rational approximants.Here in this section, we make a pair of observations about the more straightforward polynomial case.These considerations are fundamental if one wants to understand how numerical analytic continuation is possible.The first observation is that, in stating in this section that polynomial approximants only converge in a Bernstein ellipse, we do not mean to say that no polynomials exist that would be accurate further out.Indeed, Runge's theorem of 1885 guarantees that polynomials accurate on arbitrary compact subsets of do exist if is simply connected [42,85].The point is that one will not find them by computing approximations on [−1, 1].
The second observation pertains to the question, if p n is chosen to approximate f on [−1, 1], why is it also accurate further out in the Bernstein ellipse?In fact, there exist plenty of polynomials p n that approximate f on [−1, 1] without having any accuracy elsewhere.However, such approximations can not converge at the optimal rate O(ρ −n ).The crux of Bernstein's result is that the only way that polynomial approximations can be optimal or near-optimal on [−1, 1] is if they approximate f also on a larger domain.And thus approximation automatically leads to analytic continuation. 2

Complex examples
Analytic continuation also makes sense with complex functions and complex domains, and in this section we look at two examples.The first concerns the Riemann zeta function.As was mentioned in the introduction, ζ(z) is defined a priori by the Dirichlet series This series converges just for Rez > 1, but in fact, ζ(z) is analytic throughout the complex z-plane except for a simple pole at z = 1.An enormous amount is known about its numerical evaluation, including algorithms so fast as to enable computation around the 10 36 th zero on the critical line Rez = 1 2 [17].Generic numerical analytic continuation of the kind discussed in this paper can hardly compete with these specialized algorithms, but it is still interesting to get an idea of how it behaves.Suppose we evaluate ζ(z) using (9) at 300 points in the z-plane linearly spaced from 4 − 100i to 4 + 100i.Taking 40,000 terms in the series suffices for 14-digit accuracy of the evaluations, and the code segment zeta = @(z) sum((4e4:-1:1).ˆ(-z),2);Z = linspace(4-100i,4+100i,300).';r = aaa(zeta,Z); phaseplot(r,50*[-1 1 -1 1]) produces a rational approximation r (z) ≈ ζ(z) in a fraction of a second, whose phase portrait is shown in Fig. 4. (The color at each point of a phase portrait indicates the complex argument of the function there, with red corresponding to the positive real axis, yellow to argument π/3, and so on [113].)The striped region of the figure closely matches the phase portrait of ζ(z) itself [95], revealing the pole at z = 1, the zeros at z = −2, −4, −6, −8, and the first 10 pairs of zeros along the critical line.The accuracy of the zeros is 8-10 digits throughout the region shown in the figure .For example, the first two zeros of r in the upper half-plane are Fig. 5 AAA analytic continuation of a conformal map from boundary data.The black curve is the superellipse defined by x 4 + y 4 = 1, whose boundary correspondence function f (z) for a conformal map to the unit disk has been computed numerically via an integral equation.The inverse map f −1 from the unit circle to is then fitted by a rational function r , and the red and blue curves are images under r of a polar grid inside and outside the unit circle, respectively.The red dots are poles of a rational approximation to the forward map f 0.499999999946+14.134725141712i,0.499999999940+21.022039638788i. ( It is obvious in these numbers how closely the real parts of the zeros of r and ζ agree, and the accuracy of the imaginary parts is similar. Our second example comes from a 1986 paper by Reichel on numerical analytic continuation [84].Let be the region in C bounded by a Jordan curve , let f be a conformal map of onto the unit disk, and suppose we know the boundary correspondence function, that is, the values of f on .In practice, this information would traditionally be found by numerical solution of an integral equation [104,114].Based on these data, how can we evaluate f (z) for z in the interior , or, assuming is smooth enough for this to make sense, how can we analytically continue it to the exterior?As proposed in [52], an excellent method is AAA approximation, as illustrated in Fig. 5 for Reichel's example of a superellipse defined by the curve x 4 + y 4 = 1.The black curve is , the red curves are images in the interior of of a polar grid in the unit disk, and the blue curves are images of a polar grid outside the unit circle.The core of the computation is carried out by the commands Z = exp(2i*pi*(0:100)/100); g = chebfun2(@(x,y) x.ˆ4+y.ˆ4,1.1*[-1 1 -1 1]); superellipse = conj(roots(g-1)); f = conformal(superellipse); F = f(Z); r = aaa(F,Z); The superellipse is defined in lines 2-3 by a bit of Chebfun trickery having no connection to this paper, the boundary correspondence function f ( ) is computed via the Kerzman-Stein integral equation in line 4, and line 5 computes the rational approximation that serves to evaluate r in the exterior of and continue it analytically outside.A further computation [˜,pol] = aaa(Z,F) finds poles of a rational approxi-mation of the inverse map from the unit disk to , and it is these poles that are plotted as red dots in the figure.

The one-wavelength principle
This section proposes a principle that has both experimental and (a little) theoretical support: One-wavelength principle.In 16-digit arithmetic, rational approximations can analytically continue an oscillatory function about one wavelength beyond its sampling domain.
Two wavelengths would be possible in 32-digit arithmetic, three wavelengths in 48digit arithmetic, and so on.Related results for Fourier-based extrapolation methods are presented in [13].
Let us explore this idea experimentally.Figure 6 shows a pair of computations in the pattern of the last section for the analytic continuation of f (z) = cos(2π z) to z > 0 from its values for z ∈ [−2.6, 0] (the number 2.6 is arbitrary).The wavelength of this function is 1, and with the grid lines as a guide, we focus on the question of how far along the positive x-axis a minimax polynomial (left) or rational approximation (right) retains some accuracy.With degrees n = 28 and 14, respectively, both approximations have accuracy about 7 × 10 −14 on [−2.6, 0].The upper pair of plots are in the same format as before, and the lower pair show the data (black) and the approximations (red) along the x-axis.It is apparent that the polynomial approximation extrapolates to z ≈ 0.7 and the rational approximation to z ≈ 1.
Here is what we mean by the one-wavelength principle.When an oscillatory function is approximated by a minimax, AAA, or other rational function to an accuracy of, say, 10-16 digits, one often observes that the approximation retains some accuracy at a distance of about one wavelength from the approximation domain.Of course, this principle would not be of much use if it only applied to sines and cosines.One sees similar behavior more generally, however, for functions that have something of a smooth oscillatory character near the region of interest.Several of our later examples will illustrate this, and in fact, one can see the one-wavelength effect in Fig. 1 in the approximation of tanh(z).Like the cosine, this is an exactly periodic function, with period π in the Im z direction.The colored contours on the right of Fig. 1 show that the rational approximation retains some accuracy up to the second pole in the upper half-plane.
A rough empirical model of the relationship between accuracy of a rational approximation and the number of wavelengths of extrapolation is this: no. of wavelengths = C 10 • (no. of digits of accuracy), C 10 ≈ 0.075, (11) or equivalently, no. of wavelengths = C 2 • (no. of bits of accuracy), C 2 ≈ 0.023.( 12) The function f is entire, whereas g is analytic in the strip −1 < Im z < 1. (Another way to write g to display its analytic structure explicitly is The data show approximately linear increase of the number of wavelengths of analytic continuation with the number of digits of approximation accuracy.In this figure, the number of wavelengths of analytic continuation of f is measured by the first point z > 0 at which | f (z) − r (z)| ≥ 0.1, and for g, which is smaller in amplitude on the real line, it is the first point at which |g(z)−r (z)| ≥ 0.001.
I do not have a formulation of the one-wavelength principle precise enough to be cast as a theorem or a conjecture.An approximate indication of where the phenomenon comes from will be presented in Sect.8.3.Meanwhile, here is the closest I have found to a similar statement in the literature, from a 1986 sampling theory paper by Henry Landau [64].The one-wavelength aspect of the matter comes in via his assumption that the data are band-limited.
Theorem 2: When sample measurements are accurate only to within ε > 0 in amplitude or in total energy, good extrapolation is possible for only a bounded distance (having an order of magnitude − log ε) beyond the interval of observation, regardless of the amount of data used.From τ it is straightforward to construct "blending-to-zero" [23] or partition-ofunity functions with similar properties, (with L and R standing for left and right), as shown in Fig. 8.These analytic functions add up exactly to 1 while having numerical support (−∞, 1] and [−1, ∞), respectively.With these tools one can piece disjoint analytic functions f L and f R together with the formula as illustrated for two pairs of functions in Fig. 9.
The examples of Fig. 9 are extension problems of what Boyd calls the "first kind," where f L and f R are known and analytic throughout the desired blending region [20].Sometimes, however, f L or f R may have singularities ("second kind") or be unknown ("third kind") in the blending region.These cases can also be handled by the tools we have introduced.We reduce the third kind situation to the first or second kind by AAA extrapolation, which will normally provide extensions of f L and f R at least The stretching process of ( 18)-( 17) applied to two functions analytic for x ≤ 0.5 but not x ≤ 1, shown in grey.In each case the portion of the function in [0, 0.5] is analytically stretched to [0, 2].On (−∞, 0], the stretched function matches the original to 16 digits some distance beyond their initial intervals.If the extensions reach across the blending region without singularities, then we are in the first kind setting and the problem is done. For problems where f L and f R , after AAA extension, remain unknown or have singularities in the blending region, we can use τ to "stretch" these functions analytically so that they have the necessary property.To formulate this with some generality, let us shift the intervals and suppose that a function f (x) is given that is analytic for x ∈ [−L, a] for some a, L > 0 and we want to find another function g(t) that is analytic for . This is mathematically impossible in general, but with the help of τ , it is achievable to machine precision.We define g(t) = f (x(t)), where x(t) is an analytic homeomorphism of [−L, b] to [−L, a] that matches the identity to 16 digits, i.e. x(t) ≈ t, with If x(t) is the inverse function (readily computed with Chebfun), this gives us Two examples are shown in Fig. 10.
Once functions have been analytically continued by AAA and/or analytically stretched as just described, we are back in the "first kind" situation and can combine them as illustrated in Fig. 9. Figure 11 shows as an example an analytic blending of the upper function of Fig. 10 (shifted left 1 unit) and the lower one (flipped and then shifted right one unit).
Normally, a branch cut is a curve across which an analytic function discontinuously changes from one branch to another.As suggested in Figs.11 and 12, our subject in this section might be whimsically called fat branch cuts, that is, domains of positive width across which an analytic function smoothly changes from one branch to another.The change is so smooth that we may regard the result as a single globally analytic function, with the understanding that it matches the original functions not exactly but to a prescribed accuracy such as 16 digits.
The blending method we have described for computing fat branch cuts is not the only possibility of this kind.More familiar choices would make use of entire functions such as Gaussians and the error function, which are even smoother than ( 14) but approach their limiting values more slowly.I do not claim that ( 14) is necessarily better than these methods, merely that it is a choice worth considering.In the context of numerical quadrature, a great deal is known about the comparison of double expenoential, erf, and other formulas based on variable transformations.See for example Theorem 4.1 of [97].
As a very different approach to blending unrelated analytic functions, one might simply calculate a least-squares fit to the data on the given disjoint domains by a global polynomial, which will then give some kind of analytic interpolant in the inbetween region.In general it is necessary to use Stieltjes orthogonalization, also known as "Vandermonde with Arnoldi," to generate a well-conditioned polynomial basis for this process [22].An example involving two real intervals is given in Figure 3.1 of [22], and Fig. 12 shows an example involving two disks.Compare Figure 8 of [34].
Fig. 12 Phase portraits of approximations of sign(Re z) on a pair of disks in the left and right halves of the complex z-plane.The AAA approximant on the left generates an approximation to a branch cut close to the imaginary axis, an effect discussed further in Section 7.3.The degree 80 polynomial least-squares fit on the right (which matches the data with the large maximal error of 0.056) illustrates another approach to fat branch cuts Before moving on we should mention that sometimes one is in the fortunate situation of trying to blend two analytic functions that are in fact the same-that is, the problem of missing data or data with gaps.In such a case a global AAA interpolant can be very effective, as is shown in Fig. 2 of [119] and is also essentially the theme of [56].This is related to the easy problem of analytic continuation inwards from a closed boundary, as in the contours inside the superellipse in Fig. 5.

Hermite integrals and potential theory
To say anything precise about rational approximations, one usually makes use of a set of tools built on the Hermite integral formula [66,106,111].The origins go back to Cauchy in the 1820 s and Hermite in the 1870 s, as well as other works including Runge's 1901 paper on equispaced interpolation.It was Joseph Walsh at Harvard in the early and mid-20th century who used the Hermite integral extensively for analyzing rational approximations [111], and then Gonchar, Bagby, Rakhmanov, and others later in the 20th century who connected these tools with notions of potential theory [5, 6, 46-48, 65, 77, 83, 86].
Walsh's idea was to analyze a type and n interpolation points where r which are distinct from the poles.If the poles are distinct, the rational function can be written for some coefficients {a k }, and if the interpolation points are distinct, the interpolation conditions are We will use ( 21)-( 22) for simplicity, but in fact the theory applies equally if some of the poles and/or interpolation points are confluent, in which case one has poles of higher order and/or interpolation conditions involving derivatives as well as function values.
There is no reason why r has to interpolate f at exactly n points.In fact, r and f could agree at any number of points from 0 to ∞ (countable or uncountable).However, the number n is special since it is the number of coefficients {a k } in (21), and it is readily proved that if n poles (19) and n interpolation points (20) are fixed arbitrarily, then there is a unique interpolant ( 21)- (22) for any data values f (z 1 ), . . ., f (z n ) [111,Thm. 8.1].The formulas of the new few pages are valid with any choice of n interpolation points, regardless of what additional interpolation points may also be present.
The Hermite integral formula starts from the function which implies The numerator of ( 24) is the geometric mean distance of z to the interpolation points {z k }, and the denominator is the geometric mean distance to the poles {π k }.Thus |φ(z)| and this is a harmonic function (i.e., it satisfies the Laplace equation) away from {z k } and {π k }.We can interpret u as the potential generated by n negative point charges of strength n −1 at {z k } and n positive point charges of strength −n −1 at {π k }, with And now we come to the contour integral.Suppose is a continuous closed contour enclosing {z k } but not {π k }, and suppose f is analytic in the region bounded by and continuous up to the boundary.More generally, could have several components, so that has some holes.The Hermite integral takes this form [111, Thm.8.2]: Theorem 2 (Hermite integral formula) Under the circumstances just described, the rational interpolant r satisfies Here is the crucial implication of ( 27): Much of rational approximation theory is built on this principle.For example, consider Fig. 13.On the innermost black level curve we have |φ(z)| 1/n = 0.3, hence |φ(z)| < 10 −5 .On the outermost, we have |φ(t)| 1/n = 0.9, hence |φ(t)| > 10 −1 .Suppose we take the contour in (27) to be this outermost curve.Then for any z ∈ [−1, 1], and any f analytic within , we find In words, the degree 10 rational interpolant to f can have error no greater than order 10 −4 .If n is doubled, this will improve to about 10 −8 , and so on: exponential convergence.
The Hermite integral tells us something about analytic continuation, too.Between the blue dots and the red ones in Fig. 13, we are still guaranteed good accuracy, though with progressively weaker constants as one moves outward.This is an analogue for rational approximations of Bernstein's result for polynomials, explaining in a general way the success of rational approximation as observed, for example, in Figs. 1, 2 and 3.There is a logical gap, however, in that whereas Theorem 1 applies to any polynomial approximations of a function f , the Hermite formula assumes that the poles lie at prescribed locations.So although it is suggestive, it gives us no immediate information about other best or near-best approximations, whose poles are free.We will take up that subject in Sect.8.
Given an approximation domain E such as [−1, 1] and a domain of analyticity such as the disk bounded by = {z : |z| ≤ 2}, what is the best choice of {z k } and {π k }?A natural idea is to aim to minimize the ratio φ(z)/φ(t) of ( 27) by minimizing the quotient The infimum of (30) will be approximately achieved by a minimal energy configuration in which {π k } and {z k } are distributed on and E, respectively, in such a way is to minimize the energy functional Determining an optimal configuration (for the given choice of ) is challenging for finite n, and would not be truly optimal in any case since f (t) and t − z also affect the value in (27).In the limit n → ∞, however, the optimality problem becomes very clean, just a planar Laplace problem.To treat this limit, we imagine continua of interpolation points and poles defined by a signed measure μ supported on E, where it is nonpositive with total mass −1, and on , where it is nonnegative with total mass 1.The energy functional becomes Fig. 14 Level curves |φ(z)| 1/n = 0.4, 0.5, . . ., 0.9 in the limit n → ∞ of ( 23) (from inside out), or equivalently, u = log(0.4), . . ., log(0.9) for the optimal charge configuration determined by solving a Laplace problem for a potential u(z) taking constant values on both boundaries.The contours hug the circle a little more closely than in Fig. 13 because that configuration was not quite optimal.The light grey orthogonal curves, equispaced with respect to the normal derivative u n on the boundaries, indicate the optimal charge density distributions (interpolation points, poles) on the inner and outer boundaries with associated potential function It can be shown that there is a unique measure μ that minimizes I (μ), and that this corresponds to a configuration where u is a harmonic function taking constant values u E < 0 on E and u = 0 on .The minimum I min = inf μ I (μ), which is positive, is the logarithm of the modulus of the pair (E, ), which we denote by mod(E, ).
Let us return to the problem of Fig. 13, taking E = [−1, 1] and as the circle of radius 2. The charge configuration we chose there, with {π k } equispaced and {z k } equal to Chebyshev points, is a good one as n → ∞, giving exponential convergence, but it is not quite optimal.The optimal potential for n → ∞ concentrates charge a little more near z = ±2 on the outer boundary, and is shown in Fig. 14.
We have just described the Laplace problem that was solved to produce Fig. 14, but we did not prescribe it fully since we did not specify u E .This constant is determined by the condition that there is −1 unit of charge on E or, equivalently, 1 unit on .(The charges must be equal and opposite because u = 0 on the outer boundary .)The charge density along is 1/2π times the normal derivative u n defined with respect to the outward-pointing normal.Along ∂ E the essence of the matter is the same, but if E is an interval rather than a Jordan domain, one must think of it as having two sides, each with half the charge density.Thus the additional condition to fully specify u is or equivalently, subject to the qualification just mentioned, In both cases the contours or ∂ E are traversed in the usual positive sense, with the domain always lying to the left.All together, the Laplace problem in the region between E and that determines asymptotic accuracy of Walsh rational interpolants is with the constant u E being determined implicitly by (36) or (37).
Once the equilibrium potential u(z) for a condenser (E, ) has been found, the associated convergence rates for rational interpolants follow essentially from ( 26) and (27).Let us say that {r n }, n = 0, 1, 2, . . ., is a family of Walsh approximants of f in (E, ) if it is a family of rational interpolants ( 21)-( 22) defined by poles {π n } and interpolation points {z k } asymptotically distributed according to the potential u(z).(The notion of asymptotic distribution can be made precise; we omit details.)The next theorem is due to Walsh [111, chap. 9].

Theorem 3 (Asymptotic accuracy of rational interpolants) Any family of Walsh approximants {r n } for an analytic function f in a condenser (E, ) satisfies
for each z enclosed by .For z ∈ E this becomes Note that ( 39) and ( 40) are inequalities, not equalities.The reason for this is that they are derived from bounding the Hermite integral (27) via the absolute value of its integrand, and f (z) − r (z) might be smaller if there is cancellation.This effect can be arbitrarily strong, but in many cases of best and near-best approximations it is a factor of exactly 2 for reasons of orthogonality, as we will discuss in the next section: exp(−1/cap(E, )) becomes exp(−2/cap(E, )).

Rational approximation on E
How fast can rational approximations on a simply-connected compact continuum E converge to a function f that is analytic there?Certainly exponential convergence is always possible, since polynomials can achieve this, as established in the case E = [−1, 1] by ( 3) and ( 4); for more general sets E ⊆ C this result is due to Runge [85].The convergence rate for rational approximations, however, is generally much faster than polynomial approximations can achieve, and it depends on what portion of C it is possible to analytically continue f to.In important cases the rate improves to super-exponential.
We now apply the theory of the last section to study these convergence rates.As always we assume that E is a simply-connected compact continuum and is a Jordan curve or collection of Jordan curves containing E inside which f is analytic.
Before considering the various cases, we make a comment about best vs. near-best rational approximation.It is known that f has a best rational approximation of each degree n, i.e., a function r * n that minimizes the supremum norm f − r n E [111].(In general r * n need not be unique, though it is unique for real approximation of real functions on an interval.)In the theoretical literature, many discussions are framed in terms of best approximations, whether on a set E or in the limiting case of Padé approximation at a point.In numerical practice, however, and certainly when using the AAA algorithm, we are normally dealing with near-best rather than exactly best approximations.Therefore in this paper, where possible, we avoid speaking of r * n .

Natural boundaries
The first case we consider is a simple one conceptually, but rarely turns up in applications.A natural boundary or natural barrier for an analytic function is a curve across which no analytic continuation is possible.For example, a function f analytic in the open unit disk has a natural boundary on the unit circle if there exists no function g that takes the same values as f within the disk and is also analytic in some larger connected domain.One source of examples of functions with natural boundaries is Taylor or Chebyshev series with random coefficients, as we have seen with the amber function A(z) of Fig. 3 and is discussed in the monograph by Kahane [59].
Taking the most straightforward case of this kind, let us suppose that f extends from E to an analytic function throughout the region bounded by E and a Jordan curve or set of Jordan curves , and that is a natural boundary of f .Then convergence at least as fast as the rate determined by the condenser capacity cap(E, ) is guaranteed.

Theorem 4 For an analytic function f with a natural boundary , there exist rational approximations {r n } with exponential convergence,
The bound (41) is sharp in the sense that, given E and , there exist functions f of the prescribed class for which the inequality is an equality.However, it is far from sharp in that there also exist other functions of this class with rational approximants that converge much faster, indeed arbitrarily fast.We can illustrate this by an example.Consider the function f (z) defined for |z| < 1 by what looks like an infinite series of simple poles on the unit circle, where {c k } is an absolutely convergent series of nonzero coefficients, i.e., |c k | < ∞.
It is easily seen that f is analytic in the unit disk.But since the numbers e ik are distinct, it can also be shown that f cannot be analytically continued to any point e ik , where k is a positive integer. 3Since these points are dense on the unit circle, the unit circle must be a natural boundary for f .(Thus the points e ik are not actually poles of f .Curiously, the same formula (42) also defines an analytic function f outside the unit circle, but f and f are not analytically related.)Now pick coefficients in ( 42) that decay as fast as you like, such as By truncating the series, we see that on any disk E = {z : |z| ≤ τ }, 0 < τ < 1, rational approximations of degree n = 1, 2, 3, . . .can achieve errors bounded by 1.1/(1 − τ ) times 10 −10 , 10 −100 , 10 −1000 , . . . .Thus the convergence can be arbitrarily fast, even though f has a natural boundary.

Entire functions
At the other extreme, suppose there exists an analytic continuation of f to all of C: f is entire.Here, (3) and ( 4) imply that polynomials, let along more general rational functions, can achieve super-exponential convergence.Another way to say it is that if a function is entire, it has super-exponentially convergent rational approximants even if you constrain all their poles to lie at ∞.We do not state this as a theorem as it is so basic and also a special case of the next theorem.E Γ Fig. 15 Sketch of the proof of Theorem 5 as a consequence of Theorem 3.For a function f with isolated singularities, we may choose defined by an arbitrarily large outer contour and arbitrarily small inner loops around each singularity.The associated convergence rate 1/cap(E, ) can accordingly be made as fast as we like

Isolated singularities
If f has singularities in the complex plane, Theorem 1 implies that polynomial approximations will be limited to exponential convergence.Rational functions, however, can still achieve super-exponential convergence if the singularities of f are isolated.For a start, suppose f can be continued to all of C except for poles: f is meromorphic in C. In such a case each pole is necessarily isolated from the others, and their number is finite or at most countably infinite.Here, it is obvious that rational approximations can converge super-exponentially, since they can "peel off poles" one after another and thus effectively approximate over bigger and bigger regions of analyticity.
In fact, a kind of "peeling off singularities" works for essential singularities too, and the theory of the last section explains why; the idea is sketched in Fig. 15.If f can be continued to all of C apart from isolated singularities, then it is analytic within curves and unions of curves whose reciprocal condenser capacity 1/cap(E, ) can be made as large as we want by making hug the singularities more and more closely.We enclose E in a large circuit of , and around each essential singularity within this outer circuit, we put a component of in the form of a small loop.This gives a convergence estimate based on the capacity cap(E, ).Taking larger and larger outer circuits and smaller and smaller loops around the singularities, we can make cap(E, ) as small as we like and thus get any exponential rate of convergence.Theorem 5 For a function f analytic in C except for isolated singularities, there exist rational approximations {r n } with super-exponential convergence,

Branch points
There is another kind of structure we need to consider that appears all the time in applications: branch points.A branch point is a singularity in a neighborhood of which f cannot be analytic and single-valued, like the point z = 2 for the function As mathematicians, we are used to introducing branch cuts to handle such functions, or to thinking of f as a multivalued function in the complex plane or a single-valued function a Riemann surface.But a rational approximation is defined everywhere in C and single-valued.What will happen when the analytic continuation involves branch points?
The theory of the last section points the way, and the idea is sketched in Fig. 16.The theory applies to any condenser domain (E, ) in which f is analytic (in the standard single-valued sense).This means that we are free to make use of any choice that divides up C in such a way as to make the analytic continuation of f single-valued.For any such choice, Theorem 3 will hold, with its bound (40) involving the modulus exp(1/cap(E, )).So we are free to adjust to optimize this constant.The supremum of all such moduli is called the extremal modulus for f and E: where the supremum runs over all contours inside which f is analytic.The extremal modulus is ∞ for the cases considered in the last two subsections, but it is always finite if f has branch points.
Theorem 6 For a function f analytic in C except for isolated singularities and branch points, there exist rational approximations {r n } with exponential convergence at the rate determined by the extremal modulus As a converse of Theorem 6, it is known that if rational functions r n exist such that lim n→∞ f − r n 1/n E = 0, or even lim inf n→∞ f − r n 1/n E = 0, then the analytic extension of f -even though it may live in just a small portion of C-cannot have any branch points [48].
The extremal modulus corresponds to a boundary (no longer strictly composed of Jordan curves) constituting an optimal branch cut or set of branch cuts for rational approximation of f on E. Thus although the analytic continuation of a function f intrinsically only has branch points, not branch cuts, its optimal rational approximations do indeed have a well-defined set of branch cuts.I have heard them called "God's branch cuts," and a more sedate term might be "Stahl's branch cuts" since they were powerfully analyzed by Stahl for the case of Padé approximations [92,94].Stahl's results were later extended by Buslaev [26] to multipoint Padé approximation.
When rational approximations are computed for functions with branch points, the theory we have just described emerges beautifully.One sees that they feature approximate branch cuts delineated by interlacing zeros and poles.To illustrate this effect, Fig. 17   Figure 18 shows analogous images for f (z) = √ z 2 − 2z + 2, which has a pair of branch points at z = 1 ± i. Now the rational approximation has an approximate branch cut along a circular arc connecting these two points.For other functions f , more complicated curves will appear.
The poles and zeros of approximate branch cuts constructed by rational approximations tend to be exponentially clustered near the branch points.(In Fig. 17, the distances of the poles from z = 2 are about 123.3, 12.5, 3.72, 1.37, 0.472, 0.103, and the distances of the zeros are 29.8, 6.50, 2.24, 0.822, 0.245, 0.025.)This clustering effect goes back at least to Newman in 1964 [72] and is investigated in detail in [106].

The factor of 2
All our discussion so far has concerned convergence at the rate determined by mod = exp(1/cap(E, )), which we can express loosely as In actuality, however, best and near-best rational functions often converge at twice this rate: For example, Fig. 3 shows that this is the case with the amber function.The polynomial convergence behavior is that of (48), where is the Bernstein 2-ellipse.This is a natural boundary, so Theorem 4 would lead us to expect convergence of rational approximants at the same rate.In fact, however, the rate is more like (49), as is reflected in the figure showing comparable accuracies for polynomial degree 46 and rational degree 23.The factor of 2 can be explained intuitively by rational functions of a given degree n having twice as many parameters as polynomials of the same degree.At a deeper level it is related to certain orthogonality properties [82], analogous to the factor of 2 boost that Gauss quadrature obtains from orthogonality in comparison with other quadrature formulas [102].Yet (49) does not always hold, and here is an example.Consider approximation on the disk E = {z : |z| ≤ τ }, 0 < τ < 1, of the function which has a natural boundary on the unit circle.Polynomial truncations obviously converge at the rate τ n , and it can be shown that in the lim-sup sense, rational approximations can do no better. 4xamples like (50) are quite extreme: the approximation difficulty only arises for a certain subset of degrees n.In fact there is a counter-result of Parfenov [74], confirming a conjecture of Gonchar, that shows that the factor of 2 always applies if we look at the lim-inf instead of the lim-sup [66, Eq. ( 21)]: More satisfyingly, the doubled rate actually does hold in the lim-sup sense for a class of problems that includes most applications: functions whose only non-isolated singularities are algebraic branch points.More precisely, we assume that f is a (multivalued) analytic function that can be analytically continued along arbitrary curves in C so long as they avoid a fixed set S of capacity zero.We say that f is continuable along all curves excluding S. The notion of a set of capacity zero (also known as a polar set) is a standard one [83].
Theorem 7 (Gonchar-Stahl ρ 2 theorem) Let f be an analytic function on a compact set E which is continuable along all curves in C excluding a set S of capacity 0. Then the best degree n rational approximations r * n to f on E satisfy

Analytic continuation beyond E
The last section explored the exponential or super-exponential convergence that is possible for rational approximants to analytic functions on a simply-connected compact continuum E. We considered Walsh approximants, constructed on the basis of the asymptotically optimal distribution of poles and interpolation points as determined by potential theory, and for these rational functions, convergence outside E is also assured, with exponential rates diminishing as one moves outward from E toward .Numerical analytic continuation, however, is not based on prescribed poles, for we would not know where to put them.On the contrary, an approximation r ≈ f on E is computed somehow based on just information on E. In this situation, can we expect its behavior for z ∈ C\E to be comparable to that of Walsh approximants?
The rough answer is yes.In practice, best and near-best approximations, and AAA approximants in particular, usually behave at most points z ∈ C\E about as the theory of the last section would lead us to expect.Exponential or super-exponential convergence is typically observed as expected.Most remarkably, as we have already demonstrated by examples in Figs. 17 and 18, rational approximants to functions with branch points typically converge in the region defined by the optimal branch cuts associated with the extremal modulus.For high accuracy very close to the branch points, however, or for continuation across branch cuts onto different Riemann sheets, different methods are called for that are based on approximants which themselves have branch cuts.A starting point is quadratic Padé approximation [36,91], which generalizes to Hermite-Padé approximation at higher orders [8,21,96], and AAA analogues of such methods can also be devised.Another approach is reciprocal-log approximation [4,71].
All this encouraging news, however, as we said, is "rough."Rational approximants are not guaranteed always to converge cleanly, and anomalies often arise.In the next two subsections we will describe the difficulty, which goes by the name of Froissart doublets, followed by the theoretical notion that has been devised to cope with it, namely convergence in capacity.
As a further observation of surprising effects that may appear with rational approximation, we note that no matter how fast approximants r n converge to an analytic function f on a set E, this will never guarantee that f is analytic at any particular point of C\E.To see this, consider again the function (42) with a set of coefficients decaying very fast as in (43).By truncating the series, we see that on any disk E = {z : |z| ≤ τ }, 0 < τ < 1, rational approximations of degree n = 1, 2, 3, . . .can achieve errors bounded by 1.1/(1 − τ ) times 10 −10 , 10 −100 , 10 −1000 , . . . .Thus no rational approximation convergence rate, no matter how fast, implies that f can be analytically continued to any domain beyond the unit circle.This is in marked contrast to the situation with polynomials as asserted in Theorem 1. Natural boundaries may be impassable, yet they can be surprisingly flimsy!

Froissart doublets
Rational approximations sometimes have poles where one doesn't expect them, far from any singularities of f , which may destroy the quality of an approximation at least as measured in the ∞-norm.Figure 19 illustrates this effect.The left image repeats the calculation of Fig. 1, showing in a closeup two of the seven poles of this AAA rational approximation.The image on the right is the same, except that now, the approximation has been computed from data corrupted by random noise of amplitude on the order of 10 −13 , which is the default tolerance of the Chebfun AAA code.The striking result is the appearance of 41 new poles, 27 of them on the real axis and 7 conjugate pairs.To judge by the figure, it would seem that we now have a useless approximation.
The first thing to say about bad poles is that they are often not as bad as they look.The new poles in this experiment have extremely small residues, none bigger than 1.18 × 10 −11 . 5Thus each of these new poles corresponds to a contribution r k to r that Fig. 19 On the left, computed poles as in Fig. 1 for the AAA approximation of f (z) = tanh(z) in 100 points of [−1, 1].On the right, the same but with the approximation data perturbed by noise on the scale of 10 −13 .Forty-one poles with negligible residues have appeared on or near [−1, 1], which contribute nothing useful to the approximation behaves locally for z ≈ z k like with |d k | < 1.2 × 10 −11 .Such a pole with an extremely small residue, far from singularities of f , is called a spurious pole.When the residue is so small, this means that r k is of negligible size, for most purposes, except when z is extremely close to z k .Thus, although spurious poles make the supremum-norm error of an approximation over a region of interest infinite, they may hardly show up at all when the approximation is applied in practice.
If r has a pole with residue ε 1, it probably also has a zero nearby, at a distance O(ε).For example, in the right image of Fig. 19, each of the 41 new poles has a zero at a distance no greater than 10 −10 , and in all but two cases the distance is less than 10 −14 .A pole-zero pair like this is called a Froissart doublet, and this term highlights another way to interpret the limited effect of spurious poles on approximations.Except very close to z k , the pole and zero nearly cancel, making r behave for many purposes almost as if they were not there.In this paper we speak interchangeably of spurious poles and Froissart doublets; for consistency, the term we usually choose is Froissart doublets.
In numerical rational approximation, Froissart doublets appear almost invariably when one tries to approximate to an accuracy close to the level of rounding errors or other noise.For discussions of this issue, many with fascinating numerical plots, see [40, 43-45, 50, 51, 76, 96].Sometimes robust algorithms can be devised to limit the effect, for example by applying the singular value decomposition [50,51].Broadly speaking, however, it is usually necessary to loosen the tolerance so that one is approxeigenvalue problems involving arrowhead matrices [70].The actual computation of poles, residues, and zeros is carried out by the Chebfun function prz, which has just 12 executable lines.
imating to an accuracy well above the noise level.The Froissart doublets of Fig. 19 go away, for example, if aaa is called with a tolerance of 10 −12 instead of the default 10 −13 .With AAA approximation on a discretization of a continuum E like [−1, 1], Froissart doublets also often appear if the discretization is too coarse, in which case it may be enough simply to take more sample points.This issue becomes particularly important in problems with poles exponentially clustered near E, in which case it is important that the sample points be exponentially clustered too [52,106].
In a purely mathematical analysis, when rounding errors or other noise are not an issue, Froissart doublets arise less commonly.Nevertheless, they may still appear.It is a fact of rational approximation that optimality of approximations on a set E, or at a point in the case of Padé approximation, may imply the appearance of pole-zero pairs that seem to contribute nothing useful to an approximation and yet are mathematically correct.This effect has been known since Perron early in the 20th century, who showed that functions exist whose Padé approximations have poles appearing infinitely often on a dense set of points in the complex plane [75, §78].Later Wallin constructed a function whose Padé approximations are unbounded for every z = 0 [109].A general analysis of the phenomenon of spurious poles can be found in [93].

Convergence in capacity
In the face of this mathematical reality, it has of course been impossible to prove that optimal rational approximations must converge beyond the domain of approximation E as one would like.However, a slightly modified form of convergence can still be established.As n → ∞, the residues of unwanted poles of certain degree n Padé approximants are guaranteed to decrease fast enough that the portion of the complex plane contaminated by them gets smaller and smaller.The initial idea, pursued by Nuttall in 1970 [73], is that for any ε, the measure of points in the complex plane where | f (z) − r (z)| > ε must shrink to 0 as n → ∞.Shortly afterwards it was shown by Pommerenke [78] that the same conclusion holds with measure replaced by capacity, and this stronger and sharper result is now called the Nuttall-Pommerenke theorem [8].Although the theorem applies a priori just to Padé approximation of functions in a certain class, it provides a model of what we may hope for more generally: that nearoptimal rational approximants will converge in capacity as n → ∞ on compact subsets of C disjoint from singularities of f .We say that r n converges in capacity to f on if for any ε, η > 0, there are sets Theorem 8 (Nuttall-Pommerenke theorem) Let f be a single-valued analytic function in C\S, where S has capacity 0, which is analytic at z = ∞.Then the Padé approximants to f at z = ∞ converge super-exponentially to f in capacity.
If f has branch points, then a similar result holds with exponential convergence, as proved by Stahl [92,94].For the following theorem, we assume that f is a (multivalued) analytic function that can be analytically continued along arbitrary curves in C so long as they avoid a fixed set S of capacity zero.As before, we say that such a function is continuable along all curves excluding S. Theorem 9 (Stahl's theorem for convergence with branch points) Let f be a multivalued analytic function in C\S, where S has capacity 0, which is analytic at z = ∞ and continuable along all curves excluding S. Then there is a certain branch-cut set S ⊆ C such that the Padé approximants to f at z = ∞ converge exponentially to f in capacity in C\S.
As commented above, mathematical Froissart doublets, in the absence of rounding errors or other noise, are not as common as these worst-case results might suggest, and they are certainly unlikely to be found at much of a distance from the approximation set E. The reason is that a pole with a very small residue can contribute very little to the accuracy of an approximation on E. Therefore a pole of this kind can only appear in a marginal situation where a different pole with a bigger residue could not be useful.In particular, if a degree n rational approximation exhibits a Froissart doublet, then an almost equally good approximation must exist of degree n −1.Thus Froissart doublets are associated with stagnation of convergence curves. 6

Return to the one-wavelength principle
I do not know a formulation of the one-wavelength principle precise enough to be called a conjecture, so of course I do not have a proof.However, an indication can be given of the kind of balance that leads in this direction.
Theorem 2 presented the Hermite integral for the accuracy of rational approximations, As emphasized in the boxed formula (28) and the subsequent discussion of sections 6-8, good rational approximations are associated with poles lying on contours far enough from the approximation set E that |φ(t)| is very big for t ∈ .If f is entire, say, then could be taken as large as desired, but there is a tradeoff because f will then be very large.For the model problem of degree n rational approximation on [−1, 1] of the function f (z) = sin(z), taking as the circle |z| = R gives the rough estimate size of f on : e R , size of φ on : R n .
Dividing these quantities gives size of f size of φ on : e R−n log R , and differentiating the exponent with respect to R shows that this ratio takes its minimum value with R = n.So we expect R, and hence the number of wavelengths of f that are approximated, to grow roughly linearly with Since the convergence will also be roughly exponential with n, this suggests a roughly linear relationship between digits of approximation accuracy and wavelengths of analytic continuation.We hasten to add that this argument is only a very crude first approximation to a phenomenon that is certainly more complicated.
The asymptotically linear growth of the scale of poles and zeros of rational approximations as functions of degree was investigated comprehensively for Padé approximation of e z by Saff and Varga [87], extending an earlier result for zeros of Taylor polynomial approximations due to Szegő [99].

ODE and PDE applications
Throughout the mathematical sciences one finds problems involving ordinary and partial differential equations (ODEs and PDEs) posed, for the most part, in terms of real variables.Typically the solutions are real analytic, and it may be of interest to know, what happens in the complex plane?For example, many problems feature rapid transitions, and these tend to be associated with complex singularities close to the real axis.Some authors have argued, for problems ranging from the Navier-Stokes equations to the nonlinear Schrödinger equation, that a full understanding of the dynamics is only possible through an analysis of complex singularities.For example, see [27,108].
One-dimensional analytic continuation problems may involve a time variable t and/or a space variable x.In the present section we will focus on cases where a solution has been computed numerically, and it is desired to explore how the solution changes if a variable becomes complex.Sometimes it may be possible to compute the solution in the complex plane directly by a suitable extension of the numerical discretization.When this is not feasible, or just for speed or simplicity, an alternative is to use numerical analytic continuation.

Lorenz equations
For an example involving just t-dependence, consider the Lorenz equations, a system of three coupled nonlinear ODEs, The upper plot of Fig. 20 shows the trajectory u(t) for 0 ≤ t ≤ 5, numerically computed from the initial data u(0) = v(0) = −15, w(0) = 20.One sees a few cycles of the familiar chaotic oscillation, and gray lines have been drawn to mark local maxima of |u(t)|.
For the lower half of the figure, a AAA approximation r (t) to this trajectory has been computed based on 500 equispaced samples (degree 60, 0.2 s of desktop time), using Fig. 20 Above, a short segment of a chaotic trajectory u(t) of the Lorentz equations (55).Below, poles in the complex plane of a AAA approximation of u(t), giving an indication of where singularities of the analytic continuation lie.[100,108,112] a AAA tolerance of 10 −8 to roughly match the accuracy of the numerical trajectory.The poles of r have then been plotted as red dots in the complex plane. 7We see from this computation that u(t) appears to be analytic about the real axis in a complex strip of half-width about 0.2, and within this strip, r (t) will presumably be an accurate approximation of the true analytic continuation of u(t).Note that the singularities come close to the real axis at the times where |u(t)| reaches its maxima.This is typical of all kinds of dynamical processes: rapid transitions for real t reflect nearby singularities for complex t.
From an image like Fig. 20 alone, one does not get much information about whether the complex singularities of the problem are poles, branch points, or more complicated.Analysis of the Lorenz equations by Tabor and Weiss and later Viswanath and Sahatoglu has indicated that in this case, the singularities are branch points formed from powers of logarithms ("psi series") [100,108].
For an earlier experiment involving rational approximation as in Fig. 20, though based on least-squares fitting in Chebyshev points (Chebfun ratinterp) instead of AAA, see the paper by Webb [112], who also applies rational extrapolation techniques to the Lotka-Volterra equations and choreography solutions of the three-body problem.

Blasius equation
The Blasius equation is a 3rd-order nonlinear ODE originating in the analysis of boundary layers of fluid flows [16].The standard formulation poses the equation on [0, ∞) with two boundary conditions at x = 0 and an inhomogeneous Neumann condition at x = ∞: (An alternative solution method, about twice as fast in Chebfun, is to fix u (0) to an arbitrary nonzero value and solve an initial-value problem, which can then be rescaled to satisfy u (L) = 1.) Figure 21 presents a contour plot of the absolute value of the rational approximation r computed by the code segment just listed.In the complex plane, we see approximate three-fold symmetry; it is known that the true solution has this property exactly.The pole on the negative real axis is at about −5.505, and the conjugate pair in the right half-plane have moduli about 5.679; the exact value is about 5.690.These numbers are not known analytically, but have been computed to high accuracy along with many other quantities in a paper by Boyd [19], which includes this charming comment: The smoothness and monotonicity for real x belie a complex-plane structure which is rather complicated.No nonlinear differential equation relevant to engineering, it seems, is too simple to be uncomplicated off the real axis.

Burgers equation
A problem whose behavior in the complex plane has had a good deal of attention is the viscous Burgers equation, where the subscripts denote partial derivatives and ν > 0 is a viscosity parameter.For computational simplicity, periodic boundary conditions are often chosen, though there is no fundamental need for this restriction.
The idea of numerical analytic continuation of numerical solutions of (57) into the complex x-plane goes back to Sulem, Sulem, and Frisch in 1983 [98] and Kida shortly thereafter [60].Their method involved Fourier series, which, as a periodic analogue of polynomials, are mostly useful for locating the closest complex singularity to the solution of interest.(A multivariate generalization can be found in [67].)In other words, Fourier series used in this fashion give an approximation to u(x, t) in a strip of analyticity around the real x-axis, but not beyond.Later authors including Caflisch, et al. [27] and Weideman [116] took the step to rational approximations, using Fourier-Padé approximation to fit the periodic data [8].(Some authors just call the technique Padé approximation, since it amounts to Padé approximation with respect to the variable z = e i x .)Most recently, AAA approximation has been used for this problem by VandenHeuvel et al. [107]; see also [30] for another application.We do not claim that AAA is superior to Fourier-Padé for this application, merely that it is a fast and universal tool applicable in almost any context.Figure 22 shows an example of AAA approximation for ( 57)-( 58) with ν = 0.075 based on the initial condition u(x, 0) = − sin(x).Although exact solutions for this equation are possible via the Cole-Hopf transformation [107], the PDE was solved by a 128-point Fourier spectral discretization in space coupled with the MATLAB solver ode113 with tolerance 10 −8 in time.(The computational part of the code is just a dozen lines long.)At each time t = 0, 0.3, 1.5, 4, and 10, an AAA fit with tolerance 10 −6 has been computed to the 128-point data.The poles are plotted together with a label indicating the total number of poles of the AAA approximation, about half of them off-screen.(Note that in this experiment, ordinary non-periodic AAA approximation is used, even though the problem is periodic and there exists a trigonometric variant of AAA approximation for periodic approximations, with a corresponding Chebfun code AAAtrig [3].)The computation is very fast, and works well for a real time movie with a new AAA computation at each time step.
There is much to be said about the poles of Burgers solutions, and of their AAA approximations as shown in Fig. 22. Theoretically, for this initial data, it is known that for each t > 0, the analytic continuation of u(x, t) is a meromorphic function in the complex x-plane with an infinite sequence of conjugate pole pairs extending along the imaginary axis to ±i ∞.(The Painlevé property that all singularities are poles holds for the Burgers equation generally for t > 0, for any initial data [107].)As Fig. 22 In the upper row, numerically computed solutions u(x, t) of the periodic Burgers equation ( 57)-( 58) with ν = 0.075 at five times t, starting from initial data u(x, 0) = − sin(x).In the lower row, poles in the complex x-plane of AAA approximations to these five solutions.As t increases, the shock sharpens up and the poles move in from ∞ to approach the real axis.As it increases further, the shock weakens and the poles move slowly back to ∞.The singularities of the exact analytic continuation of solutions of the Burgers equation are known to consist of an infinite number of poles the approximate shock steepens, these poles move toward the real axis, and then they diverge off to ∞ again as t → ∞ and the shock loses amplitude.This is true for any positive viscosity ν in (57), but as ν decreases, the poles come closer together, and in the inviscid case ν = 0, we have a pair of branch points, and the analytic continuation of u(x, t) is no longer single-valued in C [15].At the critical time T c = 1, a true shock forms and the branch points coalesce on the real axis.There has been a great deal of study of these effects [15,27,36,60,89,90,98,107,116].For example, in the recent paper [107], a case is considered where the initial function already has singularities in the complex x-plane, namely u(x, 0) = (1 + x 2 ) −1 .It is found that as soon as t becomes positive, there are an infinity of nearby poles.
Figure 23 shows a phase portrait of the degree 12 AAA rational approximation to u(x, t) at t = 1.5, now on larger axes −2π ≤ Re x, Im x ≤ 2π , with poles marked by white dots and zeros by black dots.The figure emphasizes that of course, the rational approximation captures only a few of the infinity of poles of the exact solution.From a theorist's point of view this may seem unfortunate, but as always in analytic continuation, the poles far from the real x-axis are of negligible importance to the approximation, and would change completely if the initial condition were perturbed even very slightly.We see this sensitivity to perturbations in the fact that the approximate branch cut of the AAA approximation is curved rather than lying exactly on the imaginary axis.The curvature has only a small effect on the accuracy of the approximation on the real axis.It would be hard to argue that the exact poles of u(x, t) far from the real x-axis are of much physical or mathematical significance.

Analytic continuation in two variables
In this final major section we consider the numerical continuation of a function f (x, y) that is real analytic on a domain E in the real x-y plane.Analyticity at a point (x 0 , y 0 ) ∈ E means that the bivariate Taylor series exists at that point and converges to f for (x, y) in a real neighborhood of (x 0 , y 0 ).This will imply that the series also converges in a complex neighborhood of (x 0 , y 0 ).We speak of two variables for concreteness, but everything carries over immediately to higher dimensions.
The theme of this article has been the power of rational approximation, and for a bivariate function f , one might expect a bivariate rational approximation to be effective.This may be true in principle [31], and I have been sent some encouraging experiments by Anthony Austin (private communication) based on the recently developed methods presented in [2].However, it appears that no algorithm for multivariate rational approximation is yet known that approaches the speed of AAA approximation in the univariate case.For this reason, our discussion will be restricted to taking advantage of standard, univariate AAA.
The idea is simply to apply AAA along line segments or other analytic arcs.If S is a line segment in E, we can parametrize S as (x(t), y(t)) with x(t) = at + b and y(t) = ct + d for t ∈ [0, 1] and some real constants a, b, c, d.This gives us a univariate function F(t) = f (x(t), y(t)) that is analytic in t and can be analytically continued by AAA approximation in the usual way.The same applies if S is a segment of a circle or any other analytic arc.Univariate extensions for bivariate functions have certainly been applied by authors previously; see e.g.[12, p. 137] and [23].

Oscillatory function in 2D
For a starting example, consider the function  The resulting plot nicely illustrates the one-wavelength principle.

Absolute value function
Next we look at an example with a single singularity, that is, the 2-norm of the vector (x, y) T .This function is analytic throughout the real (x, y)-plane except at the point (0, 0).(If x and y are viewed as complex variables, f is analytic everywhere in C 2 except on the planes x = iy and x = −iy.)Following the same pattern as in Fig. 24, Fig. 25 illustrates the method of AAA approximation for f .We see that the accuracy is good along lines that are well separated from the singularity at (0, 0).This example illustrates a point mentioned in section 2 for polynomials.There exist bivariate rational functions r (x, y) that approximate f accurately all around the singularity; it is just that one can not expect to find them by approximation of data at is radially symmetric and matches f with maximal error about 0.026 over the domain −3 ≤ x, y ≤ 3, producing a figure almost indistinguishable from that on the left of Fig. 25 (not shown).This approximation is derived from the type (2, 2) minimax approximation of √ s on [0, 9], exploiting the fact that f (x, y) depends only on s = x 2 + y 2 .

Singular values and pseudospectra
The absolute value function (60) can be interpreted as the minimal singular value of a matrix, in the trivial case where A is the zero matrix of dimension 1 × 1 (or indeed of any dimension).Singular values of parametrized matrices are real analytic functions [18,24].In particular, the level curves of σ min (A − z I ) are the boundaries of the pseudospectra of A in the complex z-plane [105].Thus the pseudospectra of a matrix A may be mathematically determined throughout the complex plane by the behavior in a small neighborhood, though as always, computing the analytic continuation in practice will be challenging.Figure 26 applies numerical analytic continuation in the case where A is the random 4×4 matrix generated by the MATLAB commands rng(1), A = (randn(n)+1i* Singular value surfaces offer thought-provoking examples of real analyticity.Wherever the singular value of σ k (A − (x + iy)I ) is simple, the dependence on x and y is analytic, but where singular values coalesce, one must be more careful.Any normal matrix, starting with a 2 × 2 case like A = [0 0; 0 1], will have singular values that change their order at points equidistant between eigenvalues.Analytic continuation still makes sense, but ceases to give the information one might hope for.In particular, σ 1 now typically becomes σ 2 ; see [18] and [24].If A is a normal matrix, one cannot learn anything from analytic continuation of pseudospectra near one eigenvalue about any of the other eigenvalues.
Figure 27 considers the same function as in Fig. 26 from two more angles.In the first image, AAA approximants are computed based on data from both sides, both −3 ≤ y ≤ −2 and 2 ≤ y ≤ 3. It is not surprising that one now gets a good approximation to much of the surface, apart from near the eigenvalues of A. As mentioned in the final paragraph of section 5, if this function were globally analytic, this would be a problem of missing data, and we would expect the two-sided approximation to be accurate everywhere.In the second image, based on the same data values with −3 ≤ y ≤ −2 and 2 ≤ y ≤ 3, the procedure illustrated in Fig. 11 64) of ( 63) the x-axis of the x-y-plane.The AAA approximation, based on 150 equispaced data points with −1 ≤ x ≤ 1, gives good accuracy for a little more than one wavelength at both ends of the interval −2 ≤ y ≤ 2 as described in Sect. 5 and patched together by the formulas ( 14)-( 16) (with 3.2x adjusted to 1.6 y).This gives us a function of (x, y) that can be expected to be smooth with respect to both x and y, though analyticity with respect to the x variable is not actually guaranteed since AAA approximations are involved that are not guaranteed to vary analytically with x.

Helmholtz wave field
Our remaining three examples illustrate different aspects of analytic continuation of solutions to the Helmholtz equation where k is a constant.(The function of Fig. 24 was also a Helmholtz solution.)This generalization of the Laplace equation governs wave propagation in acoustics, elasticity, and electromagnetics.Figure 28 shows a plot of a solution to (63) with k = 40, where J 1 is the Bessel function of the first kind of order 1.This corresponds to a field of concentric waves around a point source at x = 2, y = 1.In the figure, AAA approximation has been applied on 150 equispaced points with y = 0 and x ∈ [−1, 1] to extrapolate the solution to other values of x with y = 0. Good accuracy appears to be maintained for a little more than one wavelength at both ends, near x = −1 and x = 1.15.19725, 19.73921, 29.52148, 31.91264, 41.47451. (65) Figure 29 shows the corresponding eigenmodes, and in addition, each contour plot has been analytically continued into the upper-right quadrant, which is not part of the L shape.This continuation has been computed by AAA approximation along vertical lines from the lower-right leg of the L, as suggested by the white arrows in the first image.For this problem, the exact analytic continuation is just a reflection because of the Schwarz reflection principle.Near the right edge of each image, the numerical continuation is evidently accurate, as we would expect in view of the one-wavelength principle.Toward the center, however, the continuation loses accuracy because of the singularity at the center point, much as in Figs. 25, 26, and 27.Only the third eigenmode shows no inaccuracy, because in this case the domain decouples into three squares and there is no singularity.
Fig. 30 Numerical extension of a 2D function across a boundary, after [41].The function is given inside the curve, where it is sampled on a polar grid.These samples are then used as data for AAA approximations, one for each angle of the grid.About one wavelength of successful extension is achieved.This computation of 80 AAA approximants takes about 0.1 s on a laptop

2D function extension
Our final example comes from a paper of Fryklund, Lehto, and Tornberg [41], representative of a wider literature.For a variety of reasons, especially related to the numerical solution of PDEs by integral equations and other methods, it is often desirable to extend a smooth function in the interior of a domain locally into the exterior [23,41,61].Various methods have been proposed, and [41] presents a local 2D radial basis function (RBF) approximation, which can then be evaluated on the other side of the boundary and also made global via a partition of unity decomposition.The computation of Fig. 30, modeled on the example shown in Figure 5 of [41], is based on a polar grid with 80 radial lines each sampled at 90 equispaced points.

Discussion
Analytic continuation is a large and old subject, and there are many other aspects of it that might have been discussed here.A particularly interesting one is the Schwarz function, which is the analytic function S defined in a neighborhood of an analytic arc that satisfies S(z) = z for z ∈ [29].As discussed in [11], S is a natural starting point for understanding the analytic continuation of analytic and real-analytic functions across curved boundaries, and it can be computed numerically by AAA approximation.I hope to consider the Schwarz function in a future paper.
The speed and power of AAA approximation for analytic continuation were illustrated by an example that came up in my work recently.In Matlab, I needed to evaluate the reciprocal gamma function f (z) = 1/ (z) for complex values of z in the range 0 ≤ Rez ≤ 3, −1.5 ≤ Im z ≤ 1.5.The Matlab gamma command, however, only works for real arguments.At first, I assumed I would have to track down some gamma function software for complex arguments.Then I realized that, for mid-range accuracy at any rate, this was a routine problem of analytic continuation.(Indeed, rational approximations are an established method for computing the gamma function [88].)The code segment X = linspace(-3,6,50); r = aaa(1./gamma(X),X);runs on my laptop in a couple of milliseconds and produces a function handle r for a rational function degree 14 whose maximum error |r (z) − f (z)| over the region of interest is less than 10 −8 .Evaluations of r take a fraction of a microsecond per point z, so I was able to evaluate f (z) to 8-digit accuracy by this method at a million points in one second.

Fig. 1 Fig. 2
Fig. 1 Error contours in the complex z-plane (log base 10) for approximations of f (z) = tanh(z) for z ∈ [−1, 1].On the left, degree 30 polynomial minimax approximation or Chebyshev interpolation, with poles of f marked by small black dots.(The figure is computed with Chebyshev interpolation, but minimax looks almost the same.)On the right, degree 7 rational minimax or AAA approximation based on 100 equispaced sample points in [−1, 1], with poles of r marked by thicker red dots.(The figure is computed with AAA, but again minimax looks almost the same.)Polynomial approximations are only useful for analytic continuation within a Bernstein ellipse, whereas rational approximation may be effective further out

Fig. 3
Fig.3Repetition of Fig.1again but now for the "amber function" A(z) of (8), which we believe has a natural boundary along the Bernstein 2-ellipse.The polynomial approximation on the left has degree 46, and the rational approximation on the right has degree 23, so they have equal numbers of free parameters, 47.Inside the ellipse, rational approximations have no advantage over polynomials.Outside, the function f and thus the error f − r are not defined, so no color is plotted.The reason rational functions are so important for analytic continuation is that functions like this one are uncommon in applications

Fig. 4
Fig. 4 Phase portrait of the AAA approximation r (z) to the Riemann zeta function ζ(z) based on sample values at 300 points linearly spaced in the complex z-plane from 4 − 100i to 4 + 100i.The striped region of the figure closely matches the phase portrait for ζ(z) itself.Its left boundary lies roughly along the curve defined by |ζ(z)| = 10 7

Fig. 6
Fig. 6 In the upper row, a repetition of Fig. 1 for f (z) = cos(2π z), a cosine function with wavelength 1.The values at 150 points of [−2.6, 0] represent several wavelengths of sample data, and the approximations have accuracy 7 × 10 −14 on this interval.For z > 0, successful analytic continuation by somewhat less than one wavelength is evident in the left image, and by about one wavelength on the right.The lower curves confirm these estimates, showing the data in black and the approximations p(z) and r (z) in red.With sample data on a larger interval [−L, 0], the polynomial continuation would get worse for reasons explicable by Bernstein ellipses, whereas the rational one would not change much

Fig. 7
Fig.7Extended precision Julia computations of AAA approximants of the functions f and g of (13) in 500 Chebyshev points scaled to [−4.7, 0].As the degree n of the approximations increases from about 11 to 47, the number of digits of accuracy on [−4.7, 0] increases steadily.The number of digits of successful analytic continuation along the positive real axis also increases steadily in a fashion consistent with the approximate model(11)

Fig. 9
Fig. 9 Two examples of blending by (16) of one analytic function f L on (−∞, −1] with another unrelated one f R on [1, ∞).The global function f is analytic in a strip around R and matches f L and f R on their subintervals to 16-digit accuracy

Fig. 11
Fig.11Analytic blending of the two functions of Fig.10, flipped and shifted as described in the text.The blended function is analytic and matches the given functions to 16 digits.We may think of the red part of the curve as a "fat branch cut"

Fig. 13
Fig. 13 Level curves |φ(z)| 1/n = 0.3, 0.4, . . ., 0.9 of (23) (from inside out) for a rational approximation in which the poles {π k } are n = 10 equally spaced points on the circle |z| = 2, marked in red (the nth roots of 1024), and the interpolation points {z k } are n Chebyshev points on [−1, 1], marked in blue ) Physically, we think of {π k } as a set of point charges of amplitude n −1 that repel each other pairwise with forces of magnitude n −2 /|π k − π j |, and similarly {z k } are point charges of amplitude −n −1 repelling each other pairwise with forces of magnitude n −2 /|z k − z j |; meanwhile each pair π k and z j are attracted by a force of magnitude n −2 /|z j − π k |.If all these charges are distributed in an equilibrium fashion on and E, this will minimize I .
shows a AAA approximant r to (45) on [−1, 1].On the left are plotted the

Fig. 16 Fig. 17
Fig.16 If f has branch points, then Walsh rational approximants can only be constructed in a domain bounded by curves that force f to be single-valued.Optimizing over all such curves implicitly defines one or more optimal branch cuts for rational approximation of f over E

Fig. 18
Fig.18 Same as Fig.17but for the function f (z) = z 2 − 2z + 2, which has branch points at z = 1 ± i. Now the approximate branch cut lies on a circle connecting the two branch points, as suggested schematically in Fig.16and explained rigorously for Padé and multipoint Padé approximation by Stahl[92] and Buslaev[26]

Fig. 21
Fig. 21 AAA extrapolation into C of the solution to the Blasius problem (56) computed numerically on the interval [0, 8] (marked in blue).The extrapolant catches the 3-fold symmetry of the solution, with the triad of closest to zero having moduli within a few percent of the true value 5.69003805 . . .[19].The absolute value |u(z)| is plotted

Fig. 23
Fig.23 Phase portrait in the complex x-plane of the AAA rational approximation of the numerically computed solution u(x, t) of Fig.22at t = 1.5.The poles near the real axis are good approximations to those of the exact analytic continuation

Fig. 24
Fig. 24 Bivariate AAA analytic continuation of the smooth oscillatory function f (x, y) = sin(2x + y).On the left, level curves | f (x, y)| = −0.75,−0.5, . . ., 0.75; the negative levels are marked in red.On the right, the same plot for the approximation g(x, y) computed by AAA approximation along the horizontal lines marked by black dots

Figure 24
shows numerical analytic continuation from data with −4 ≤ x ≤ −2 into the region −2 ≤ x ≤ 4. The AAA algorithm is applied along lines y = const., with the approximation in each case based on data at 21 points (x, y) with y fixed and x equispaced from −4 to −2.All together, this amounts to the construction of 81 different AAA approximations (taking about 0.1 s all together).

Fig. 25
Fig. 25 Bivariate AAA analytic continuation of the absolute value function f (x, y) = x 2 + y 2 .On the left, level curves | f (x, y)| = 0.25, 0.5, . . ., 3.0.On the right, the same levels for the approximation g(x, y) computed by AAA approximation along the horizontal lines marked by black dots

Fig. 26
Fig. 26 Bivariate AAA analytic continuation of the function σ min (A − (x + iy)I ) for a 4 × 4 random matrix A, showing in each case levels 0.1, 0.3, 0.5, . . ., 1.9 (from inside out) on the domain −3 ≤ x, y ≤ 3. The first image corresponds to A and the second to analytic continuation along vertical lines from data points with −3 ≤ y ≤ −2.The blue dots are the eigenvalues of A Fig. 27 On the left, an image as in the right side of Fig. 26 but now with AAA approximations based on two-sided data at 2 ≤ y ≤ 3 as well as −3 ≤ y ≤ 2. On the right, the procedure of Section 5 is applied to produce a global smooth blend between the same two stripes of data values has been applied.First, AAA approximation is used to extrapolate values from −3 ≤ y ≤ −2 to −2 < y ≤ −1 and from 2 ≤ y ≤ 3 to 1 ≤ y < 2. These values are then analytically stretched to

Fig. 28
Fig.28 Analytic continuation of the solution (64) of (63) the x-axis of the x-y-plane.The AAA approximation, based on 150 equispaced data points with −1 ≤ x ≤ 1, gives good accuracy for a little more than one wavelength at both ends of the interval

Fig. 29
Fig. 29 Eigenmodes of an L-shaped drum.In each case the lower-right portion of the eigenmode has been continued by AAA approximation upward into the upper-right quadrant, as suggested by the white arrows.The third image is special in that the center point is nonsingular • [−1,1] denotes the supremum norm over [−1, 1].It follows that minimax (i.e., best supremum-norm) polynomial approximants also have errors satisfying the same bound.The Chebyshev interpolants to f , that is, the degree n polynomials defined by interpolation of f in the n + 1 Chebyshev points z j = cos( jπ/n), 0 ≤ j ≤ n, also converge at the rate O(ρ −n ), with the upper bound differing from (3) by just a factor of 2: