Exponential node clustering at singularities for rational approximation, quadrature, and PDEs
aa r X i v : . [ m a t h . NA ] J u l Numer. Math. manuscript No. (will be inserted by the editor)
Exponential node clustering at singularities for rationalapproximation, quadrature, and PDEs
Lloyd N. Trefethen · Yuji Nakatsukasa · J. A. C. Weideman
Received: date / Accepted: date
Abstract
Rational approximations of functions with singularities can con-verge at a root-exponential rate if the poles are exponentially clustered. Webegin by reviewing this effect in minimax, least-squares, and AAA approxima-tions on intervals and complex domains, conformal mapping, and the numericalsolution of Laplace, Helmholtz, and biharmonic equations by the “lightning”method. Extensive and wide-ranging numerical experiments are involved. Wethen present further experiments showing that in all of these applications, itis advantageous to use exponential clustering whose density on a logarithmicscale is not uniform but tapers off linearly to zero near the singularity. We givea theoretical explanation of the tapering effect based on the Hermite contourintegral and potential theory, showing that tapering doubles the rate of con-vergence. Finally we show that related mathematics applies to the relationshipbetween exponential (not tapered) and doubly exponential (tapered) quadra-
L. N. TrefethenMathematical InstituteUniversity of OxfordOxford OX2 6GG, UKE-mail: [email protected]. NakatsukasaMathematical InstituteUniversity of OxfordOxford OX2 6GG, UKE-mail: [email protected]. A. C. WeidemanApplied MathematicsStellenbosch UniversityPrivate Bag X1Matieland 7602, South AfricaE-mail: [email protected] Lloyd N. Trefethen et al. ture formulas. Here it is the Gauss–Takahasi–Mori contour integral that comesinto play.
Keywords rational approximation · lightning PDE solvers · potentialtheory · tanh and tanh-sinh quadrature Mathematics Subject Classification (2010) · · Analytic functions can be approximated by polynomials with exponential con-vergence, i.e., k f − p n k = O (exp( − Cn )) for some C > n → ∞ . Here n isthe polynomial degree and k · k is the ∞ -norm on an approximation domain E ,which may be a closed interval of the real axis or more generally a simply con-nected compact set in the complex plane. This result is due to Runge [32,60]and explains the exponential convergence of many numerical methods whenapplied to analytic functions, including Gauss and Clenshaw–Curtis quadra-ture [53,55] and spectral methods for ordinary and partial differential equa-tions [52,54]. It is also the mathematical basis of Chebfun [7].If f is not analytic in a neighborhood of E , then Bernstein showed in 1912that exponential convergence of polynomial approximations is impossible [3,55]. Bernstein also showed that in approximation of functions with derivativediscontinuities such as f ( x ) = | x | on [ − , O ( n − ) [4]. Now from the beginning, going back to Chebyshev inthe mid-19th century, approximation theorists had investigated approximationby rational functions as well as polynomials. Yet it was not until fifty yearsafter these works by Bernstein that it was realized that for this problem of ap-proximating | x | on [ − , root-exponential convergence, that is, k f − r n k = O (exp( − C √ n )) for some C >
0. This result was published by Newman in 1964 [26], who also showedthat faster convergence is not possible. With hindsight, it can be seen thatthe root-exponential effect was implicit in the results of Chebyshev’s studentZolotarev nearly a century earlier [10,23,38,65], but this was not noticed.Newman’s theorem has been a great stimulus to further research in ratio-nal approximation theory [10,11,12,20,33,36,38,59]. It has not, however, hadmuch impact on scientific computing until very recently with the discoverythat it can be the basis of root-exponentially converging numerical methodsfor the solution of partial differential equations (PDEs) in domains with cornersingularities [6,13,14,15,56]. The aim of this paper is to contribute to buildingthe bridge between approximation theory and numerical computation.In particular, we shall focus on the key feature that gives rational approxi-mations their power: the exponential clustering of poles near singularities. (Thezeros are also exponentially clustered, typically interlacing the poles, with thealternating pole-zero configuration serving as proxy for a branch cut.) Thishas been a feature of the theory since Newman’s explicit construction. Ouraim is, first, to show how widespread this effect is, not only with minimax xponential node clustering at singularities 3 approximations (i.e., optimal in the ∞ -norm), the focus of most theoreticalstudies, but also for other kinds of approximations that may be more usefulin computation. Section 2 explores this effect in a wide range of applications.In section 3 we turn to a new contribution of this paper, the observationthat good approximations tend to make use of poles which, although exponen-tially clustered, have a density on a logarithmic scale that tapers to zero atthe endpoint. Specifically, the distances of the clustered poles to the singular-ity appear equally spaced when the log of the distance is plotted against thesquare root of the index. We show experimentally that this scaling appearsnot just with minimax approximations but more generally.To explain this effect, we begin with a review in section 4 of the Hermitecontour integral, which is the basis of the application of potential theory inapproximation. We show how this leads to the idea of condenser capacity forthe analysis of rational approximation of analytic functions. Section 5 thenturns to functions with singularities and explains the tapering effect. In thiscase the condenser is short-circuited, and it is not possible to estimate theHermite integral by considering the ∞ -norm of the factors of its integrand, butthe 1-norm gives the required results. Analysis of a model problem shows howthe tapered exponential clustering of poles enables better overall resolution,potentially doubling the rate of convergence. These arguments are related tothose developed in the theoretical approximation theory literature by Stahland others [36,38,40], but we believe that section 5 of this paper is the firstto connect this theory with numerical analysis.Finally in section 6 we turn to a different problem, the quadrature of func-tions with endpoint singularities on [ − , R n denotes the set of rational functions of degree n ,that is, functions that can be written as r ( x ) = p ( x ) /q ( x ) where p and q arepolynomials of degree n . The norm k·k is the ∞ -norm on E , but, as mentionedabove, other measures will come into play in sections 5 and 6, and indeed, atheme of our discussion is that certain aspects of rational approximation areoften concealed by too much focus on the ∞ -norm.The numerical experiments in this paper are a major part of the contribu-tion; we are not aware of comparably detailed studies elsewhere in the liter-ature. Our emphasis is on the results, not the algorithms, but our numericalmethods are briefly summarized in the discussion section at the end. In this section we explore the convergence of a variety of rational approxi-mations to analytic functions with boundary branch point singularities. Ourstarting point is Fig. 1, which presents results for six kinds of approximations
Lloyd N. Trefethen et al. of f ( x ) = √ x on [0 ,
1] by rational functions of degrees 1 ≤ n ≤
20. (By thesubstitution x = t , this is equivalent to Newman’s problem of approximationof | t | on [ − , f is not special; as we shall illustrate in Figs. 2and 4, other functions with endpoint singularities give similar results.First, the big picture. The upper-left image of the figure shows ∞ -normerrors k f − r n k plotted on a log scale as functions not of n but of √ n . Withthe exception of the erratic case labeled AAA, all the curves plainly approachstraight lines as n → ∞ : root-exponential convergence. (The shapes would beparabolas if we plotted against n .) The upper-right image shows the absolutevalues of the 20 poles for the approximations with n = 20, that is, theirdistances from the singularity at x = 0. On this logarithmic scale the polesare smoothly distributed: exponential clustering. This clustering is furthershown in the lower images, for the approximation labeled minimax, by a phaseportrait [62] of the square root function (the standard branch) and its degree20 rational approximation after an exponential change of variables.The top four approximations have preassigned poles, making the approxi-mation problems linear; indeed the Stenger, trapezoidal, and Newman approxi-mations are given by explicit formulas. The AAA and minimax approximationsare nonlinear, with poles determined during the computation. Although it istempting to rank these candidates from worst at the top to best at the bot-tom (the minimax approximation is best by definition), this is not the point.All these approximations converge root-exponentially, and the differences inefficiency among them amount to constant factors of order 10, which can infact be improved in most cases by introducing a scaling parameter or two. Inparticular, minimax and other nonlinear approximations can approximatelydouble the rate of convergence of the linear approximations [28]. All these ap-proximations can achieve accuracy 10 − with degrees n ≈ n = 140,085.We comment now on the individual approximations of Fig. 1. The New-man approximation comes from the explicit formula presented in his four-pagepaper [26]. The approximation is r ( x ) = √ x ( p ( √ x ) − p ( −√ x )) / ( p ( √ x ) + p ( −√ x )), where p ( t ) = Q n − k =0 ( t + ξ k ) and ξ = exp( − / √ n ); this canbe shown to be a rational function in x of degree n . The asymptotic con-vergence rate is exp( −√ n ) [64]. This can be improved to approximatelyexp( − ( π/ √ n ) by defining ξ = exp( − ( π/ / √ n ), an example of the scal-ing parameters mentioned in the last paragraph (these values are conjecturedto be optimal based on numerical experiments).The trapezoidal approximation originates with Stenger’s investigationsof sinc functions and associated approximations [43,44,45]. Following p. 211of [55], we approximate √ x by starting from the identity √ x = 2 x/π R ∞ ( t + x ) − dt , which with the change of variables t = e s becomes √ x = 2 xπ Z ∞−∞ e s dse s + x . (1) xponential node clustering at singularities 5 -8 -6 -4 -2 -12 -8 -4 -16 -12 -8 -4 0-0.500.5 -16 -12 -8 -4 0-0.500.5 Fig. 1
Root-exponential convergence of six kinds of degree n rational approximations of f ( x ) = √ x on [0 ,
1] as n → ∞ . On the upper-left, the asymptotically straight lines on thislog scale with √ n on the horizontal axis (except for AAA) show the root-exponential effect.On the upper-right, the distances of the poles in ( −∞ ,
0) from the singularity at x = 0 showthe exponential clustering. Below, phase portraits in the complex plane of the square rootfunction (the standard branch) and its degree 20 minimax approximation on [0 , z instead of e z to enable comparisonwith the axis labels in the images above. For n ≥
1, we approximate this integral by an equispaced n -point trapezoidalrule with step size h > r ( x ) = 2 hxπ ( n − / X k = − ( n − / e kh e kh + x . (2)(If n is even, the values of k are half-integers.) There are n terms in the sum,so r is a rational function of degree n with simple poles at the points p k = − exp(2 kh ). Two sources of error make r ( x ) differ from √ x . The terminationof the sum at n < ∞ introduces an error of the order of exp( − nh/ Lloyd N. Trefethen et al. the finite step size introduces an error on the order of exp( − π /h ), since theintegrand is analytic in the strip around the real s -axis of half-width π/ h ≈ π p /n andapproximation error k r − √ x k ≈ exp( − π p n/ ∞ as well as at 0, and indeed, it converges root-exponentially not just on [0 ,
1] but on any interval [0 , L ] with
L > / √ n canbe expected to lead to a root-exponential result. Such effects are familiar inthe analysis of hp discretizations of partial differential equations when thestep sizes h and orders p of multiscale discretizations are balanced to achieveoptimal rates of convergence near corners [34].A drawback of the trapezoidal approximation is that its derivation dependson the precise spacing of the poles, since it relies on the property that thetrapezoidal rule is exponentially accurate in this special case [58]. The curveslabeled Stenger in Fig. 1 come from a more flexible alternative approach, alsoproposed by Stenger [44], where we fix n distinct poles p k ∈ ( −∞ , ≤ k ≤ n and n + 1 interpolation points x k ∈ [0 , ≤ k ≤ n , and then take r to be theunique rational function of degree n with these poles that interpolates f ( x ) inthese points. The theory of rational interpolation with preassigned poles wasdeveloped by Walsh [60] and will be discussed in section 4. For our problemof approximation on [0 ,
1] with a singularity at x = 0, a good choice is to take x = 0 and x k = − p k for k ≥
1. In particular, our
Stenger approximant isthe rational function r resulting from the choices − p k = x k = exp( − ( k − h ) , ≤ k ≤ n, (3)with h = O (1 / √ n ). Figure 1 takes h = π/ √ n .Interpolation is important for theoretical analysis, but for practical com-putation, least-squares fitting is more robust and more accurate, since it doesnot require knowledge of good interpolation points. The least-squares dataof Fig. 1 come from fixing the same exponentially clustered poles as in (3),but now choosing approximation coefficients by minimizing the least-squareserror f − r on a discretization of [0 ,
1] by standard methods (MATLAB back-slash). As always when discretizing near singularities, we use an exponentiallygraded mesh ( logspace(-12,0,2000 )), and a weight function w ( x ) = √ x isintroduced in the discrete least-squares problem so that it approximates auniformly weighted problem on the continuum. The error curve r ( x ) − √ x for x ∈ [0 ,
1] for this approximation (not shown) approximately equioscillates Stenger considered rational approximations of this kind, though not in this precise set-ting of a finite interval with just one endpoint singularity.xponential node clustering at singularities 7 -12 -8 -4 -16 -12 -8 -4 Fig. 2
Four more minimax approximations, showing the same root-exponential convergenceand exponential clustering of poles as in Fig. 1. Two involve the functions x /π and x log x on [0 , x on [0 ,
1] but with the ∞ -norm weighted by x , and one involves √ z on the disk about of radius . In the right image, n takes its final value from the leftimage for each problem, 14 for the weighted approximation and 20 for the other cases. between n + 2 extrema, indicating that it is a reasonable approximation to thebest L ∞ approximation with these fixed poles.The minimax data in Fig. 1 correspond to the true optimal (real) ap-proximations, rational approximations with free poles. Here the error curveequioscillates between 2 n + 2 extrema [55], and the error is approximatelysquared; the asymptotic convergence rate is exp( − π √ n ) [40,59].Computing minimax approximations, however, can be challenging [8], andon a complex domain they need not even be unique [16]. This brings us tothe data in the figure for AAA (adaptive Antoulas–Anderson) approximation,a fast method of near-best rational approximation introduced in [24]. AAAapproximation is at its least robust on real intervals, as reflected in the erraticdata of the figure, but for more complicated problems and in the complexplane, it is often the most practical method for rational approximation.This concludes our discussion of Fig. 1. The next figure, Fig. 2, illustratesthat these effects are not confined to approximation on a real interval or to thefunction √ x . The figure presents data for four further examples of minimaxapproximations. One set of curves shows approximation of x /π on [0 , /π chosen to dispel any thought that rational exponents might bespecial. This problem requires poles particularly close to the singularity sincethe exponent is so small. Another shows approximation of x log x on [0 , √ x again, but now it is weightedminimax approximation, with a weight function x (and the error measured isnow the weighted error, notably smaller than before). Finally the fourth set ofdata shows minimax approximation of √ z on the complex disk { z : | z − | < } . Lloyd N. Trefethen et al. -0.5 0 0.5 1-1-0.500.51
A BC -6 -4 -2 ABC
Fig. 3
The conformal map of a circular pentagon onto the unit disk has been computedand then approximated numerically by a rational function of degree 70 [13,57] by the AAAalgorithm. The poles cluster exponentially at the corners, where the map is singular.
Figure 3 turns to our first problem of scientific computing. Following meth-ods presented in [13] and [57], a region E in the complex plane bounded bythree line segments and two circular arcs has been conformally mapped ontothe unit disk, and the map has then been approximated to about eight dig-its of accuracy by AAA approximation, which finds a rational function with n = 70. This process is entirely adaptive, based on no a priori informationabout corners or singularities, yet it clusters the poles near the corners justas in Figs. 1 and 2. Many poles cluster at the strong singularity A and onlya few at the weak singularity B . Note that the poles lie asymptotically onthe bisectors of the external angles. This effect is well known especially fromthe theory of Pad´e approximation as worked out initially by Stahl [41,46].Optimal approximations line up their poles along curves which balance thenormal derivatives of a potential gradient on either side, and evidently theAAA method comes close enough to optimal for the same effect to appear.We finish this section with a look at lightning solvers for PDEs in two-dimensional domains, introduced in 2019 and applied to date to Laplace [14,15,56], Helmholtz [15], and biharmonic equations (Stokes flow) [6]. In the basiccase of a Laplace problem ∆u = 0, the idea is to represent the solution on adomain E as u ( z ) ≈ Re r ( z ), the real part of a rational function with no polesin E that approximates the boundary data to an accuracy typically of 6–10digits. The rational functions have preassigned poles that cluster exponentiallyat the corners, where the solution will normally have singularities [19,61], andthe name “lightning” alludes to this exploitation of the same mathematics thatmakes lightning strike objects at sharp corners. Coefficients for the solutionare found by least-squares fitting, making this an approximation process ofthe same structure as in the least-squares example of Fig. 1. The difference isthat the approximations are now applied to give values of u ( z ) in the interior xponential node clustering at singularities 9 sqrt(DoF) -16 -12 -8 -4 e rr o r convergence solve time = 2.31 secseval time = 22.4 microsecs per pt dim(A) = 3246 x 1073 -1 0 1-101 -0.4-0.3-0.2-0.10 Fig. 4
Example of the lightning Laplace solver [14,15] as implemented in the code laplace.m [56]. For each number of degrees of freedom (DoF), poles are clustered exponen-tially near the 12 corners of the domain E , and the numbers are increased until a solutionto 10-digit accuracy is obtained in the form of a rational function with 480 poles. This takes2 . µ s per point, with the accuracy ofeach evaluation guaranteed by the maximum principle. of the domain E , where it is not known a priori. See Fig. 4 for an example ona “snowflake” with boundary data log | z | .Lightning solvers have been generalized to the Helmholtz equation ∆u + k u = 0 [15] and the biharmonic equation ∆ u = 0 [6], as illustrated in Fig. 5.In the Helmholtz case, poles ( z − z k ) − of rational functions become singular-ities of complex Hankel functions H ( k | z − z k | ) exp( ± i arg( z − z k )), and thebiharmonic case is handled by the Goursat reduction u ( z ) = Im( zf ( z ) + g ( z ))to a coupled pair of analytic functions f and g , each of which is approxi-mated by its own rational function. The mathematics of lightning methods forHelmholtz and biharmonic problems has not yet been worked out fully, andthe analysis given in section 5 applies just to the Laplace case.Although it is not the purpose of this article to give details about lightningPDE solvers, they are at the heart of our motivation. Usually in approximationtheory, minimax approximations are investigated as an end in themselves, andthe locations of their poles may be examined as an outgrowth of this process;a magnificent example is [39]. Here, the order is reversed. Our aim is to exploitan understanding of how poles cluster to construct approximations on the flyto solve problems of scientific computing. In the last section, 13 plots were presented of the distances of poles to singu-larities on a log scale, the right-hand images of Figs. 1, 2, and 3. All showedexponential clustering, and all but three showed a further effect which wecall tapered exponential clustering , the main subject of the rest of this paper:
Fig. 5
Lightning solvers have been generalized to the two-dimensional Helmholtz (left) [15]and biharmonic equations (right) [6]. The Helmholtz image shows a plane wave incidentfrom the left scattered from a sound-soft equilateral triangle. The biharmonic image showscontours of the stream function for Stokes flow in a cavity driven by a quarter-circular bound-ary segment rotating at speed 1 and with zero velocity on the remainder of the boundary.The black contours in the corners, representing the stream function value ψ = 0, delimitcounter-rotating Moffatt vortices. Tapered exponentially clustered singularities are used inboth computations. on the log scale, the spacing of the poles grows sparser near the singularity.This was also colorfully evident in the phase portrait at the bottom of Fig. 1.The three exceptions were the Stenger, least-squares, and trapezoidal approx-imations of Fig. 1, all of which are based on poles preassigned with strictlyuniform exponential clustering. These examples illustrate that tapering of thepole distribution is not necessary for root-exponential convergence. A fourthset of data in Fig. 1 also involves preassigned poles, the Newman data, andsome tapering is apparent in this case.Figure 6 shows the nine remaining examples of exponential clustering ofpoles from Figs. 1–3, the ones with free poles, presenting the distances { d k } ofthe poles from their nearest singularities on a log scale. What is immediatelyapparent is that all the curves look straight for smaller values of k . Note thatfive of them stop at n = 20, one at n = 14, and the remaining three, from theapproximation of a conformal map of Fig. 3, at different values determinedadaptively by the AAA algorithm.Yet the horizontal axis in Fig. 6 is not k but √ k . Plotted against k (notshown), the data would look completely different. Evidently in a wide range ofrational approximations, both best and near-best, the distances { d k } of polesto singularities is well approximated by the formulalog d k ≈ α + σ √ k (4)for some constants α and σ , that is, d k ≈ β exp( σ √ k ) (5)for some β and σ . xponential node clustering at singularities 11 -15 -10 -5 ABC
Fig. 6
Tapered exponential clustering of poles near singularities for the nine examples withfree poles from Figs. 1–3 of the last section. The crucial feature is that the curves appearstraight with this horizontal axis marking √ k rather than k , where { d k } are the sorteddistances of the poles from the singularities. The data for the poles at vertex A of Fig. 3have been deemphasized to diminish clutter (black dots), since they lie at such a differentslope from the others. To make sense of the √ k scaling, let us remove the exponential from theproblem by defining a distance variable s = log d , thereby transplanting aninterval such as d ∈ [0 ,
1] to s ∈ ( −∞ , ρ ( s ) of poles with respect to s ? If ρ ( s ) were constant, this wouldcorrespond to a uniform exponential distribution of poles, requiring an infinitenumber of poles since s goes to −∞ . So some kind of cutoff of ρ ( s ) to 0 mustoccur as s → −∞ . An abrupt cutoff, as with the Stenger, trapezoidal, andleast-squares distributions of Fig. 1, leads to a linear cumulative distribution,as shown in the left column of Fig. 7. By contrast, a linear cutoff gives aquadratic cumulative distribution, as shown in the right column, and whenthis is inverted, the result is the √ k distribution we have observed.Thus the straight lines of Fig. 6 can be explained if pole density functions ρ ( s ) for good rational approximations tend to take the form sketched in theupper-right of Fig. 7. (Aficionados of deep learning may call this the “ReLU”shape.) In section 5 we will explain why this is the case and continue the storyof Fig. 7 in Fig. 11.We have not presented data in this section for lightning PDE solutions,but it was in this context that we first became aware of the importance of UNIFORM TAPERED
Fig. 7
The algebra of exponential clustering. With respect to the variable s = log d , where d is the distance to the singularity, the simplest exponential clustering of poles would haveuniform density ρ ( s ) down to a certain value and then cut off abruptly (left column). Atapered distribution cuts off linearly instead (right column), resulting in poles exponentiallyclustered in the √ k fashion seen in Fig. 6. tapered exponential clustering. In the course of the work leading to [14], thefirst author noticed that although straight exponential spacing of preassignedpoles gave root-exponential convergence, better efficiency could be achieved ifthe resulting approximations were re-approximated a second time by the AAAalgorithm. On examination it was found that the AAA approximations hadpoles in a tapered distribution, just like cases A – C of Fig. 6. The model (4)–(5)was developed empirically in this context, with σ ≈ The basic tool for estimating accuracy of rational approximations is the Her-mite integral formula [20,60]. In this section we review how this formula leadsto the use of potential theory [29], and in particular the quantity known asthe condenser capacity, for approximations of analytic functions. Building onthe work of Walsh [60], these ideas began to be developed by Gonchar andRakhmanov in the Soviet Union not long after the appearance of Newman’spaper [11,12].The following statement is adapted from Thm. 8.2 of [60].
Theorem 1
Let Ω be a simply connected domain in C bounded by a closedcurve Γ , and let f be analytic in Ω and extend continuously to the boundary.Let distinct interpolation points x , . . . , x n ∈ Ω and poles p , . . . , p n anywhere xponential node clustering at singularities 13 in the complex plane be given. Let r be the unique degree n rational functionwith simple poles at { p k } that interpolates f at { x k } . Then for any x ∈ Ω , f ( x ) − r ( x ) = 12 πi Z Γ φ ( x ) φ ( t ) f ( t ) t − x dt, (6) where φ ( z ) = n Y k =0 ( z − x k ) (cid:30) n Y k =1 ( z − p k ) . (7)To see how this theorem is applied, let Ω be a simply connected domainbounded by a closed curve Γ , as indicated in Fig. 8 (see also Fig. 9 in the nextsection), and let f be analytic in Ω and extend continuously to Γ . Suppose f is to be approximated on a compact set E ⊂ Ω , which in this section we taketo be disjoint from Γ . Theorem 1 implies that for any x ∈ E , | f ( x ) − r ( x ) | ≤ Cτ, (8)where C is a constant independent of n and τ is the ratio τ = max z ∈ E | φ ( z ) | min z ∈ Γ | φ ( z ) | . (9)If φ is much smaller on E than on Γ , then τ and hence f − r must be small.Figure 8 gives an idea of how this can happen. In each image, the red dotson Γ represent a good choice of poles { p k } and the blue dots on the boundaryof E a corresponding good choice of interpolation points { x k } . Consider firstthe upper-left image, where E and Γ define a circular annulus. The equispacedconfigurations of { p k } and { x k } ensure that τ will decrease exponentially as n → ∞ . To see this, in view of (7), we define u ( z ) = n − n X k =0 log | z − x k | − n − n X k =1 log | z − p k | . (10)This is the potential function generated by n + 1 negative point chargesof strength n − at the interpolation points and n positive point charges ofstrength − n − at the poles. Then exp( nu ( z )) = | φ ( z ) | , and therefore τ = exp( − n [min z ∈ Γ u ( z ) − max z ∈ E u ( z )]) . (11)For τ to be small, we want u to be uniformly bigger on Γ than on E . Findingthe best such configuration is an extremal problem that will be approximatelysolved if the points are placed in an energy-minimizing equilibrium position.In each of the images of Fig. 8, the points are close to such an equilibrium.Each charge is attracted to the charges of the other color, but repelled bycharges of its own color.Finding an optimal configuration (for the given choice of Γ ) is complicatedfor finite n , but the problem becomes cleaner in the limit n → ∞ , and this is Fig. 8
Potential theory and rational approximation. In each image, the shaded region is anapproximation domain E for a function f analytic in the region Ω bounded by Γ . If we thinkof the poles of an approximation r ≈ f as positive point charges and the interpolation pointsas negative point charges, then a minimal-energy equilibrium distribution of the chargesgives a favorable configuration for approximation. This is a discrete problem of potentialtheory that becomes continuous in the limit n → ∞ , enabling one to take advantage ofinvariance under conformal maps. In these images E and Γ are disjoint and the convergenceis exponential, but the third domain and its close-up illustrate the clustering effect, whichwill become more pronounced as the gap shrinks to zero. The pairs of interpolation pointsand poles marked by hollow dots delimit one half of the total, highlighting how both sets ofpoints accumulate close to the singularity. where the power of potential theory is fully revealed. We now imagine continuaof 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 nonnegativewith total mass 1. It can be shown that there is a unique measure of this kindthat minimizes the energy I ( µ ) = − Z Z log | z − t | dµ ( z ) dµ ( t ) , (12)with associated potential function u ( z ) = − Z log | z − t | dµ ( t ) , (13)and u takes constant values u E < E and u Γ = 0 on Γ . The minimum I min = inf µ I ( µ ) is known to be positive, and for minimax degree n rational xponential node clustering at singularities 15 approximations r ∗ n one has exponential convergence as n → ∞ at a corre-sponding rate: lim sup k f − r ∗ n k /n ≤ exp( − I min ) . (14)(The actual rate is in fact twice as fast as this, exp( − I min ), for functionswhose singularities in the complex plane are just isolated algebraic branchpoints [20, p. 93], [36].)The reciprocal of I min is known as the condenser capacity for the ( E, Γ )pair, a term that reflects an electrostatic interpretation of the approximationproblem. In electronics, capacitance is the ratio of charge to voltage difference.A capacitor has high capacitance if its positive and negative plates are closeto one another, so that the attraction of charges of opposite sign enables agreat deal of charge to be accumulated on them without the need for much ofa voltage difference. For fast-converging rational approximation, on the otherhand, we want Γ and E to be far apart, corresponding to a small amountof charge relative to the voltage difference, hence small numbers of poles andinterpolation points needed to achieve a given ratio τ .We can now see how the second and third images of Fig. 8 were drawn.They were obtained by conformal transplantation, exploiting the invarianceof problems of potential theory under conformal maps. The eccentric domainof the second image comes from a M¨obius transformation, and the pincheddomain of the third image comes from a further squaring. The blue and redpoints obtained as conformal images of equispaced points in the symmetricannulus are known as Fej´er–Walsh points [42].One might wonder, for arguments of this kind, is it necessary to place thepoles of r on the boundary of the region of analyticity of f ? In fact, Γ doesnot have to lie as far out as that boundary, nor do the poles have to be on Γ , for as stated in Theorem 1, the integral representation (6) is valid for anyplacement of the poles. Asymptotically as n → ∞ , however, it is known thatthe convergence rate cannot be improved by placing poles beyond the regionof analyticity of f [20]. A special choice is to put all the poles at x = ∞ , inwhich case rational approximation reduces to polynomial approximation, stillwith exponential convergence though at a lower rate than in (14). Now we examine how the analysis of the last section must change for approxi-mations with singularities. There is a considerable specialist literature here byauthors including Aptekarev, Saff, Stahl, Suetin, and Totik [2,20,31,33,40,41,46], which investigates certain best approximations in detail. Our emphasis ison the broad ideas applicable to near-best approximations too.From Fig. 8 it is clear that potential theory should give some insight when f has a singularity on the boundary of E . The lower pair of images shows clus-tering of poles where Γ has a cusp close to the boundary of E , and as the cuspis brought closer to E , the clustering will grow more pronounced. However,the argument we have presented breaks down when Γ actually touches E . The Fig. 9
Two kinds of problems of rational approximation of a function f on a domain E .On the left (section 4), f is analytic on E and poles can be placed on a contour Γ enclosing E in the region of analyticity: convergence is exponential with accuracy on the order ofexp( − nδ ) for a constant δ >
0. On the right (section 5), f has a singularity at a point z c on the boundary of E , and Γ must touch E at z c : convergence is root-exponential withaccuracy of order exp( − nδ ) again, but now with δ diminishing at the rate 1 / √ n as n → ∞ .In the circled region, the potential makes the transition from u Γ = 0 to u E = − δ . situation is sketched in Fig. 9. Physically, this would be a capacitor of infi-nite capacitance, implying that an equipotential distribution u with a nonzerovoltage difference would require an infinite quantity of charge. Mathematically,the estimate (8) fails because τ cannot be smaller than 1.To see what happens in such cases, we can examine the function φ com-puted numerically for an example problem. The left column of Fig. 10 showserror curves in type (9 ,
10) minimax approximation of √ x on [0 . ,
1] (above)and [0 ,
1] (below). (Type ( m, n ) means numerator degree at most m and de-nominator degree at most n ; we choose these parameters rather than ( n, n )to make the plots slightly cleaner.) The curves each equioscillate between m + n + 2 = 21 extrema, and in the lower curve, on the semilogx scale, wesee the wavelength increasing as x →
0. As a minimax approximation withfree poles, this rational function has m + n + 1 = 20 points of interpolationrather than the standard number n = 10 for an approximation with preas-signed poles, so for the cleanest display of the potential function φ in the rightcolumn we have picked out just half of these, marked by the red dots.The right column of Fig. 10 shows the function | φ | plotted on the approx-imation interval (the lower blue curve) and on the important portion [ − , Γ (the upper red curve). (To be precise, for theseplots the numerator of (7) ranges over just the interpolation points x , . . . , x n marked by red dots.) In the upper image, for [0 . , τ of (9) is far below 1, and the estimate (8)serves to bound the approximation error. (The actual error is about the squareof this bound since we have omitted half the interpolation points.)The lower image, which is a centerpiece of this paper, tells a strikinglydifferent story. Here again the blue curve is flat, showing the even dependence xponential node clustering at singularities 17 -10 -5 -505 10 -12 -10 -5 -8 -4 -10 -5 -1-0.500.51 10 -5 -10 -5 -4 Fig. 10
On the left, error curves in type (9 ,
10) minimax approximation of √ x on [0 . , , φ ( z ) as defined by (7) on these approximation intervalsand on [ − , − ,
0] (upper curve) and on [0 . ,
1] (lower curve), but in the lower-right image, nothing like constant behavior of | φ ( z ) | on [ − ,
0] is evident. We explain thisby noting that what matters to the accuracy of an approximation is the integral (15) of | φ ( x ) /φ ( t ) | with respect to t ∈ Γ , not its maximum. Taking advantage of this property, polesand interpolation points distribute themselves more sparsely near the singularity, freeingmore of them to contribute to the approximation further away—the phenomenon of taperedexponential clustering. on x we expect in a minimax approximation. The red curve for | φ ( z ) | on[ − , t moves from − − ,
0] is not at all a curve of constant | φ | .To understand the linearly closing gap in Fig. 10, we note that what fails inthe analysis of the last section for an approximation problem with a singularityis not the Hermite integral, f ( x ) − r ( x ) = 12 πi Z Γ φ ( x ) φ ( t ) f ( t ) t − x dt, (15)but the estimate (8) we derived from it. Implicitly (8) came from bounding(15) by H¨older’s inequality, | f ( x ) − r ( x ) | ≤ π (cid:13)(cid:13)(cid:13)(cid:13) φ ( x ) φ ( t ) (cid:13)(cid:13)(cid:13)(cid:13) ∞ (cid:13)(cid:13)(cid:13)(cid:13) f ( t ) t − x (cid:13)(cid:13)(cid:13)(cid:13) ∞ k k , (16) where the ∞ - and 1- norms are defined over t ∈ Γ . (The norm k k is equalto the length of Γ .) When Γ and E are disjoint, the first ∞ -norm in (16)is exponentially small as n → ∞ and the second is bounded. However, theseproperties fail as Γ and E touch. We can rescue the argument by notingthat | φ ( x ) /φ ( t ) | does not have to be small for all t so long as its integral issmall. More precisely, the quantity f ( t ) / ( t − x ) of (16) may not be bounded as t, x → z c but f ( t ) | t − z c | − α / ( t − x ) will be bounded if we assume f ( t − z c ) = O ( | t − z c | α ) for some constant α . So what actually matters is that the integralof | t − z c | α − | φ ( x ) /φ ( t ) | should be small, and we accordingly replace (16) bythe alternative H¨older estimate | f ( x ) − r ( x ) | ≤ π (cid:13)(cid:13)(cid:13)(cid:13) φ ( x ) φ ( t ) | t − z c | α − (cid:13)(cid:13)(cid:13)(cid:13) (cid:13)(cid:13)(cid:13)(cid:13) f ( t ) t − x | t − z c | − α (cid:13)(cid:13)(cid:13)(cid:13) ∞ . (17)For simplicity let us assume that the singularity lies at z c = 0 and the partof the contour Γ that matters is [0 , | φ ( x ) | ≈ x ∈ E . We want | φ ( t ) | to be at most 1 for t < t/ε ) α for t > ε , where ε is the distance of the closest pole to the singularity.See the upper-right image of Fig. 11. Defining u ( t ) = n − log φ ( t ) leads usto the model problem sketched in the image below this in the figure: find aharmonic function u ( t ) in the upper half t -plane such that u ( t ) = (cid:26) , t ≤ ε , αn − log( t/ε ) , t > ε . (18)We now make the change of variables s = log t , which transplants the Laplaceproblem to the infinite strip S = { s ∈ C : 0 < Im s < π } , as sketched in the(3 ,
2) position of the figure: find a harmonic function u ( s ) in S satisfying u ( s ) = , Im s = π ,0 , Im s = 0 and Re s ≤ log ε , αn − ( s − log ε ) , Im s = 0 and Re s > log ε . (19)This change of variables is convenient mathematically, and it is also importantconceptually, since it is well known that influences on harmonic functionsdecay exponentially with distance along a strip. Consequently, if ε is small,the solution to a Laplace problem for log ε ≪ Re s ≪ t variable, where behavior for | t | of order ε or less is unimportant because itcontributes negligibly to the integral (15) and behavior for | t | of order 1 ormore is unimportant because it is far from the singularity under investigation.So we address our attention to (19). An exact solution can be obtained viathe Poisson integral formula for an infinite strip [63], u ( x + iy ) = αn − π Z ∞ ξ sin( y )cosh( ξ − ( x − log ε )) − cos( y ) dξ, (20) xponential node clustering at singularities 19 UNIFORM TAPERED
Fig. 11
The potential theory of exponential clustering (continuation of Fig. 7). The firsttwo rows (right column) show the function | φ ( t ) | of (7) and the associated potential u ( t ) = n − log | φ ( t ) | for the model problem (18)–(19). The third row shows the behavior alongthe real s -axis after the change of variables to s = log t ; the domain is now the infinitestrip 0 < Im s < π , with u = 0 for Im s = π . The final row shows the charge density ρ ( s ) = nu n ( s ) /π , where u n is the normal derivative of u on the boundary of the strip. Theintervals that matter (emphasized by solid rather than dashed lines) are ε < | t | < t variable and log ε < Re s < s variable. Smaller values of | t | and s contributenegligibly to the integral (15), and larger values are far from the singularity. where we have set s = x + iy . However, we do not need exactly this since theregion where our model applies is log ε ≪ Re s ≪
0. In this region, the bilinearharmonic function u ( x + iy ) = αn − (1 − yπ )( x − log ε ) (21) satisfies the boundary conditions and is accordingly a good approximation tothe solution to (19). The corresponding pole density distribution on the real x axis is ( n/π ) times the normal derivative, ρ ( s ) = − nπ ∂∂y u ( x + iy ) = απ ( x − log ε ) . (22)This linear growth, sketched in the bottom-right image of Fig. 11, is just whatwe set out to explain in Fig. 7.Let us now look at the quantitative implications of this argument, com-paring uniform exponential clustering (left column of Fig. 11) with taperedexponential clustering (right column). According to our model, the integral ofthe solid portions of the ρ ( s ) curves in the bottom should be equal to n , thetotal number of poles. For uniform clustering the integral is ( − α log ε ) /π ,leading to the estimatesClosest pole: ε ≈ exp( − π p n/α ) , Accuracy: ε α ≈ exp( − π √ αn ) . (23)For tapered clustering the integral is ( − α log ε ) /π , leading to the estimatesClosest pole: ε ≈ exp( − π p n/α ) , Accuracy: ε α ≈ exp( − π √ αn ) . (24)Thus, as mentioned in the abstract, our model leads to the prediction of afactor of 2 speedup with tapered clustering. It would be interesting to inves-tigate whether, for certain problems, exactly this ratio could be establishedtheoretically in the limit n → ∞ .As an example of a problem in which we may make such a comparisonnumerically, consider Fig. 12. These data show ∞ -norm errors for rationallinear-minimax approximations of even degrees n from 2 to 50 with preassignedexponentially clustered poles. That is, the approximations are optimal in the ∞ -norm among rational functions in R n with simple poles at the prescribedpoints; they are characterized by error curves equioscillating between n +2 extrema. The upper curves correspond to uniformly clustered poles p k = − exp( − πk/ √ n ), 0 ≤ k ≤ n −
1, and the lower curves to tapered poles p k = − exp( √ π ( √ k − √ n )), 1 ≤ k ≤ n . The asymptotic errors appear to be aboutexp( −√ . n ) for uniform clustering and exp( −√ . n ) for tapered clustering.With α = 1 / f ( x ) = √ x , the corresponding estimates (23) and (24) areexp( −√ . n ) and exp( −√ . n ).Analyses related to the argument we have presented were published byStahl for rational minimax approximation of | x | on [ − ,
1] and x α on [0 ,
1] [37,38,39,40]. For x α Stahl gives the resultAccuracy: ε α ≈ exp( − π √ αn ) , (25)which is not just an estimate but a theorem concerning the limit n → ∞ (assuming α is not an integer), with precise constants. This is exactly whatone would expect based on (24), since, as mentioned earlier, the effective valueof n is doubled in the case of true minimax approximants [28]. xponential node clustering at singularities 21 -6 -4 -2 -12 -8 -4 Fig. 12
Linear-minimax approximation of f ( x ) = √ x on [0 ,
1] with preassigned expo-nentially clustered poles in [ − , n = 2 , , . . . ,
50. Tapering the distribution makes theconvergence rate approximately double, as predicted by the model of Section 5.
Stahl worked essentially in the variable t rather than s , so his boundaryconditions involved logarithms, as in the second image of the right column ofFig. 11. Whenever one has a Laplace problem with Dirichlet boundary data,one can interpret it as the problem of finding an equipotential distributionin the presence of an external field defined by that boundary data, and thisinterpretation has been carried far in approximation theory [33]. From thispoint of view one can say that tapered exponential clustering results from polesand zeros being slightly pushed away from a singular point by a logarithmicpotential field. In this final section we turn to another problem where exponential clusteringappears. Let f be a continuous function on [ − , f by a linear combination I n = n X k =1 w k f ( x k ) , (26)where { x k } are distinct nodes in [ − ,
1] and { w k } are corresponding weights,in such a way that the accuracy is good even if f has branch point singularitiesat the endpoints. To this end, we introduce a change of variables g ( s ) fromthe real line to [ − , I = Z − f ( x ) dx = Z ∞−∞ f ( g ( s )) g ′ ( s ) ds, (27) and we apply the equispaced trapezoidal rule. This involves an infinity ofsample points in principle, but if g ′ ( s ) decays rapidly, we may truncate theseto an n -point rule like (2): I n = h ( n − / X k = − ( n − / f ( g ( kh )) g ′ ( kh ) . (28)Quadrature formulas of this kind were introduced around 1970 by Mori, Taka-hasi, and other Japanese researchers and also in the analysis of sinc methodsby Stenger. See [17,18,43,45,49,58], as well as [21] for the history as told byMori himself. The standard “exponential” choice of g is g ( s ) = tanh( s ) , g ′ ( s ) = sech ( s ) , (29)with which (28) becomes the tanh formula . As in section 2, we estimate thetruncation error as of order exp( − nh ) and the discretization error of orderexp( − π /h ). (The latter could be worse if f has additional singularities near( − , h ≈ π/ √ n , with convergence rate of orderexp( − π √ n ). An estimate of this form is valid for any H¨older continuous branchpoint singularity; see [43, Thm. 3.4], [51, Thm. 2.1], and [58, Thm. 14.1].Root-exponential convergence! This is much better than any algebraic or-der, but for practical applications on one-dimensional domains, methods ofthis kind often seem very wasteful, with almost all the points being used upin resolving the singularity (100% of them, in the limit n → ∞ ) [30]. A yearor two after the first exponential formulas appeared, it was realized that onecan do better with “double exponential” formulas. We focus on the tanh-sinhformula proposed by Takahasi and Mori in [50] and subsequently used andanalyzed by many others including Okayama, Sugihara, and Tanaka as well asBailey and Borwein [1,22,27,47,51]. Here (29) is replaced by g ( s ) = tanh( π sinh( s )) , g ′ ( s ) = π cosh( s ) sech ( π sinh( s )) . (30)Under suitable assumptions we can now estimate the truncation and discretiza-tion errors as of orders exp( − ( π/
2) exp( nh/ − π /h ). The firstof these estimates is the big improvement, for this quantity can be almost-exponentially small with a much smaller value of h than before, of orderlog( n ) /n rather than 1 / √ n . By almost-exponential, we mean of orderexp( − Cn/ log n ) for some C >
0. With this reduced value of h , the secondestimate becomes almost-exponentially small too.Figure 13 shows data for the tanh and tanh-sinh formulas. (We usedthe empirical choices h = π/ √ n and h = 1 . πn ) /n , respectively.) Theleft image plots | I n − I | against √ n for n from 1 to 40 for the integrationof f ( x ) = √ x . The tanh curve appears straight, confirming the root-exponential convergence, and the tanh-sinh curve bends downward, confirmingthat its rate is faster. The unexpected image is on the right, a plot of distancesof the nodes from the endpoint x = −
1. For tanh quadrature, these distancesare uniformly exponentially spaced, appearing as a parabola on these axes. xponential node clustering at singularities 23 -16 -12 -8 -4 -16 -12 -8 -4 Fig. 13
On the left, root-exponential convergence of the tanh quadrature formula appliedto integration of √ x (note the √ n axis as usual); the tanh-sinh formula converges muchfaster down to machine precision. On the right, the distances of nodes from poles (with a √ k axis) show uniform exponential clustering for the tanh formula with n = 40 and taperedexponential clustering for tanh-sinh. The curve for tanh-sinh quadrature, however, is almost perfectly straight. Itwould seem that tanh-sinh quadrature exploits tapered exponential clustering!It surprised us when we first saw curves like this. Why is there a resemblancebetween the tanh and tanh-sinh quadrature formulas and the phenomena ofrational approximation discussed in the earlier sections of this article?Some steps toward an answer come from a beautiful connection introducedby Gauss and exploited by Takahasi and Mori [9,48,53]: every quadratureformula can be associated with a rational approximation. Suppose first that f can be analytically continued to a neighborhood Ω of [ − ,
1] bounded by acontour Γ . Then the integral can be written I = Z − f ( x ) dx = 12 πi Z Γ f ( t ) ϕ ( t ) dt, (31)where the characteristic function ϕ is defined by ϕ ( t ) = Z − dxt − x = log t + 1 t − . (32)On the other hand the quadrature sum (26) can be written I n = 12 πi Z Γ f ( t ) r ( t ) dt, (33)where r is the degree n rational function defined by r ( t ) = n X k =1 w k t − x k . (34) Subtracting (33) from (31) gives what we call the
Gauss–Takahasi–Mori (GTM)contour integral, I − I n = 12 πi Z Γ f ( t )( ϕ ( t ) − r ( t )) dt (35)and the corresponding error bound | I − I n | ≤ π k f k ∞ k ϕ − r k ∞ k k , (36)which we have written in the style of (16), with the norms defined over Γ .Equations (35) and (36) relate accuracy of a quadrature formula to anapproximation problem: if the nodes and weights are such that ϕ − r is smallon the boundary Γ of a region where f is analytic, then | I − I n | must be small.This reasoning was applied by Takahasi and Mori to a range of quadratureformulas [48]. Now ϕ is an analytic function in the extended complex planeminus the segment [ − , Γ is disjoint from [ − , ϕ that converge exponentially on Γ as n → ∞ . In particular, this holds for the rational functions associated withGauss and Clenshaw–Curtis quadrature [53], where it is convenient to take Γ in the form of an ellipse about [ − ,
1] with foci ±
1. It follows that both thesequadrature formulas converge exponentially as n → ∞ for analytic integrands(cf. [55, Thm. 19.3]).But what if f has endpoint singularities? Now Γ must touch [ − ,
1] at theendpoints, and (36) fails just as (16) did in such a case. In fact, this failureis more severe, since k ϕ − r k ∞ = ∞ for any r because of the logarithmicsingularities of ϕ . The last section, however, suggests a solution. Instead of(36), we can derive from (35) the bound | I − I n | ≤ π k f k ∞ k ϕ − r k . (37)The switch from the ∞ - to the 1-norm changes the problem of rational approxi-mation of ϕ profoundly. Since the dominant effects just concern approximationof a logarithmic singularity near the singular point, the essential question be-comes, how fast can log t be approximated by rational functions over bothsides of the interval [ − ,
0] in the 1-norm?As we did with Fig. 10, let us get some insight by looking at the details ofthe approximation problem. The rational function (34) for the tanh rule is r ( t ) = h ( n − / X k = − ( n − / sech ( kh ) t − tanh( kh ) , (38)and for the tanh-sinh rule, it is r ( t ) = h ( n − / X k = − ( n − / π cosh( kh ) sech ( π sinh( kh )) t − tanh( π sinh( kh )) . (39) xponential node clustering at singularities 25 -15 -10 -5 -15 -10 -5 Fig. 14
Error | ϕ ( t ) − r ( t ) | as a function of distance to the left from t = − n = 40 with h = π/ √ n and h = 1 . πn ) /n , respectively.By symmetry, the same behavior would appear to the right from t = 1. Compare Fig. 10,where the ratio of the blue and red values in the lower-right image is closely analogous tothe blue curve here. The 1-norms of the approximation errors over [ − , −
1] are indicated.The slight irregularities at the left are the result of rounding error.
Figure 14 plots | ϕ ( t ) − r ( t ) | for these two approximations.For tanh quadrature, we know that | ϕ ( t ) − r ( t ) | must diverge to ∞ as t →− ϕ at t = −
1. Yet the singularity is so weakthat the divergence only shows up as a gentle upward drift in the blue curve atthe left. Over the main part of the plot, | ϕ ( t ) − r ( t ) | decreases steadily downto around 10 − . The 1-norm, measured here over t ∈ [ − , − n → ∞ , this 1-norm decays root-exponentially.For tanh-sinh quadrature, again no approximation of ϕ is possible in the ∞ -norm. In the 1-norm, however, one might expect that the convergence willnow be almost-exponential. Indeed, | ϕ − r | decays almost-exponentially as n → ∞ over any domain bounded away from the singularity. But the 1-normdecay over the whole interval is in fact just root-exponential, as is suggestedby the number listed being barely smaller than before. The following reasoningsuggests why this must be. Consider approximation of f ( x ) = log x on [0 , g ( x ) = x log x − x with faster than root-exponential convergencein the ∞ -norm, which would contradict the evidence of Fig. 2.If k ϕ − r k decreases only root-exponentially as n → ∞ , how does thequadrature formula converge almost-exponentially? It appears that this de-pends on additional properties that go beyond rational approximation, involv-ing analytic continuation of the integrand onto an infinitely-sheeted Riemannsurface in exponentially small neighborhoods of the endpoints [47,51].There remains the phenomenon of tapered exponential clustering, so vividlyevident in Fig. 13. We do not yet have an explanation for this, nor a view of whether an approximate √ k dependence is genuine or just an artifact. Thisis a topic for ongoing research, where it would be good to investigate also thedistributions of exponentially clustered nodes, also apparently tapered, thatarise with the “universal quadrature” formulas of Bremer, et al. [5,35]. Exponential clustering of poles at singularities has been part of the landscapeof rational approximation for half a century, but we believe this is the firststudy to focus on this effect. Our motivation is that this clustering is whatmakes rational approximations so powerful, and understanding it enables oneto improve existing numerical algorithms and develop new ones. We find thesephenomena fascinating, especially the tapered clustering effect, and discover-ing that tapering also appears in double exponential quadrature was a bonus.The elucidation of these matters with the help of a sometimes seemingly end-less program of numerical experiments will forever be associated in our mindswith the Covid-19 shutdowns of 2020.Here are some details of our computations. Figures 1, 2, 6 and 10 made useof the Chebfun minimax command [8], principally due to Silviu Filip, and Filipalso kindly provided us with a modified code for the weighted minimax ap-proximations of Figs. 2 and 6. For successful results in some of these problems,we applied a M¨obius transformation of [0 ,
1] to itself to weaken the singularitywhile preserving the space R n . For the approximations of Figs. 2 and 6 on acomplex disk, the AAA-Lawson algorithm was used as implemented in Cheb-fun [7,25], again with a M¨obius transformation. Figure 3 was produced withthe confmap code available at [56], which in turn calls aaa from Chebfun [24]and laplace from [56]. The aaa code was also used directly in Figs. 1 and 6,and laplace in Fig. 5. The Stokes and Helmholtz results of Fig. 5 were pro-duced by experimental codes that are not yet publicly available developed withAbi Gopal and Pablo Brubeck, respectively. In Fig. 12, a least-squares problemwas extended by a Lawson iteration (iteratively reweighted least-squares) tocompute minimax approximations with preassigned poles. All the remainingresults are based on straightforward computations in MATLAB and Chebfun. Acknowledgements
We have benefited from helpful advice from Bernd Beckermann,Pablo Brubeck, Silviu Filip, Abi Gopal, Stefan G¨uttel, Arno Kuijlaars, Andrei Mart´ınez-Finkelshtein, Ed Saff, Kirill Serkh, Alex Townsend, and Heather Wilber.
References
1. Bailey, D. H., Borwein, J.: Hand-to-hand combat with thousand-digit integrals, J. Com-put. Sci. | x | par des polynˆomes de degr´esdonn´es, Acta Math.
8. Filip, S.-I., Nakatsukasa, Y., Trefethen, L. N., Beckermann, B.: Rational minimax approx-imation via adaptive barycentric representations, SIAM J. Sci. Comput. A2427–A2455(2018)9. Gauss, C. F.: Methodus nova integralium valores per approximationem inveniendi, Com-ment. Soc. Reg. Scient. Gotting. Recent. 39–76 (1814)10. Gonˇcar, A. A.: Estimates of the growth of rational functions and some of their applica-tions, Math. USSR-Sbornik A1494–A1522 (2018)25. Nakatsukasa, Y., Trefethen, L. N.: An algorithm for real and complex rational minimaxapproximation, SIAM J. Sci. Comput., to appear (2020)26. Newman, D. J.: Rational approximation to | x | , Mich. Math. J. ρ -theorem and associated directions in thetheory of rational approximations of analytic functions, Math. Sbornik p - and hp -Finite Element Methods, Clarendon Press (1999)35. Serkh, K.: personal communication, September 201836. Stahl, H.: General convergence results for rational approximants, in C. K. Chui et al.,eds., Approximation Theory VI, Academic Press, pp. 605–634 (1989)37. Stahl, H.: Best uniform rational approximation of | x | on [ − , x α on [0 , | x | , Constr. Approx. x α on [0 , https://people.maths.ox.ac.uk/trefethen/lightning , 202057. Trefethen, L. N.: Numerical conformal mapping with rational functions, Comp. Meth.Funct. Thy. to appear (2020)58. Trefethen, L. N., Weideman, J. A. C.: The exponentially convergent trapezoidal rule,SIAM Rev. | x | by rational functions, Sov.Math. Dokl. | x | by Newman’srational operators, Acta Math. Hungar.30