Reciprocal-log approximation and planar PDE solvers

In this paper we introduce a new problem in approximation theory. On the interval $[-\infty,0\kern 0.5pt]$, it is known that functions $F(s)$ like $e^{as}$ ($\hbox{\rm\kern 0.3ptRe\kern 1.0pt}a>0$) with essential singularities at $s=-\infty$ can be approximated by rational functions $r(s)$ with exponential convergence as a function of $n$, the degree of the rational function. Our starting point in section 2 is the change of variables $z=e^{s}$ from $s\in[-\infty,0\kern 0.5pt]$ to $z\in[\kern 0.5pt0,1]$. This gives us exponentially convergent approximations on $[\kern 0.5pt0,1]$ of functions $f(z)$ like $z^{a}$ ($a>0$) with branch point singularities at $z=0$ by functions of the form $g(z)=r(\log(z))$. In the generic case where the poles $\{s_{k}\}$ of $r$ are distinct, the approximation can be written in the partial fractions form

We show that approximations of this kind can be computed by linear least-squares fitting, and we note the equioscillatory characterization of the error in the case of minimax (best $\infty$-norm) approximations. With both best and near-best approximations, one encounters the startling property that the oscillations typically cluster double-exponentially near the singularity, readily coming as close (in theory) as a distance of, say, $10^{-1000}$. This is in contrast to the well-known case of rational approximation, where the clustering is just exponential [22].

Section 3 generalizes the discussion to domains in the complex $z$-plane with $m\geq 1$ singularities, typically at corners. Here we consider approximations of the generic form

where $p_{0}$ is a polynomial of degree $n_{0}$ (possibly 0), with total number of degrees of freedom $N=1+\sum n_{j}$. The least-squares method extends to this case too and again gives exponentially or near-exponentially convergent approximations. (By near-exponential, we mean at a rate $O(\exp(-CN/\log N))$ with $C>0$.) This and other claims are first explained heuristically and illustrated numerically. Then they are established rigorously by a succession of theorems in section 4, deriving accuracy estimates from Hermite integrals, and section 5, showing how a problem with several singularities can be decomposed into single-singularity problems by Cauchy integrals over open arcs.

As is typical with computations based on redundant expansions, numerical instability resulting from the use of ill-conditioned matrices may be a problem; we stabilize our computations by the Arnoldi orthogonalization technique presented in [1]. Section 6 considers additional fine points associated with the “$10^{-1000}\kern 1.0pt$” effect mentioned above, proposing in particular that it may be a good idea to constrain the singular part of a reciprocal-log approximation to approach zero more rapidly at each singularity.

Up to this point, the paper is all approximation theory and algorithms. The motivation, however, is numerical solution of PDEs by a “log-lightning” analogue of the lightning solvers for Laplace, Helmholtz, and biharmonic equations presented in [2], [6], and [7]. A proof of concept is briefly presented in section 7 to establish that, as intended, these approximations can be used to derive exponentially convergent new methods for solving PDEs in planar domains with corner singularities.

In a word, this is a paper about approximation by functions with logarithmic branch points rather than just pole singularities. What starts out as a seemingly arbitrary rearrangement of familiar rational approximation theory turns out to have potentially significant consequences for scientific computing.

A publication of Driscoll and Fornberg in 2001 proposed the introduction of a logarithmic term to approximate functions with jumps, with coefficients computed in the manner of quadratic Hermite–Padé approximation [4]. Their technique is different from what is being proposed here, but the use of logarithms to resolve singularities gives a family resemblance.

A classic problem in approximation theory is the rational minimax approximation of $F(s)=e^{s}$ on the interval $[-\infty,0\kern 0.5pt]$. To be precise, for each $n$, let us consider rational functions $r(s)$ of degree $n$, meaning that $r(s)$ can be written as $p(s)/q(s)$ for degree $n$ polynomials $p$ and $q$. In the generic case where the poles are distinct and finite, $r$ can be written in partial fractions form as

with coefficients $\{c_{k}\}$ and poles $\{s_{k}\}$. Of course, some rational functions have confluent poles (i.e., poles of order $>1$) and cannot be written in this form, and we shall sometimes make use of such approximations.

Cody, Meinardus, and Varga showed in 1969 [3] that approximations $r(s)$ to $F(s)=e^{s}$ exist with exponential convergence, meaning $\|F-r\|=O(C^{-n})$ as $n\to\infty$ for some $C>1$, where $\|\cdot\|$ is the $\infty$-norm on $[-\infty,0\kern 0.5pt]$. The optimal constant is $C=1/H=9.28903\dots,$ where $H$ is a number related to elliptic integrals known as Halphen’s constant; for a discussion with references, see chapter 25 of [19]. Exponential convergence of best rational approximations at the same asymptotic rate applies to $F(s)=e^{as}$ for any $a$ with $\hbox{\rm\kern 0.3ptRe\kern 1.0pt}a>0$ and to certain other functions $F(s)$ related to $e^{as}$ by rational transformations [16]. Moreover, convergence at this rate occurs not just on $[-\infty,0\kern 0.5pt]$ but throughout the complex plane [16].

By the change of variables $z=e^{s}$ and the definitions $f(z)=F(s)$ and $g(z)=r(s)$, we transplant this theory to a new setting. Now the domain is $z\in[\kern 0.5pt0,1]$, and $f$ is a function such as $z^{a}$ that may have a branch point singularity at $z=0$. What is new is that (3) now takes the form

where $\log$ denotes the standard branch of the logarithm. Thus, by the transplantation of existing theory, we know that a wide class of functions on $[\kern 0.5pt0,1]$ with branch point singularities at $z=0$ can be approximated by functions of the form (4) with exponential convergence at the rate $O((9.28903\dots)^{-n})$.

More important for scientific computing are approximations that are not best but readily computable, still with exponential convergence though not at the optimal rate. One way to derive such approximations is by expressing $F(s)$ as an integral along a Hankel contour winding around $(-\infty,0\kern 0.5pt]$ in the $s$-plane. If the integral is approximated by a truncation of the transplanted equispaced trapezoidal rule on the real axis, one obtains rational approximations (3) to $F(s)$, hence reciprocal-log approximations (4) to $f(z)$, which converge at rates $O(C^{-n})$ with values of $C$ on the order of 3 or 4 depending on details of contours and parameters. The leader in developing these methods has been Weideman [24, 28]; for a review, see section 15 of [23].

The approach we shall make use of is a variation on this theme motivated by the lightning PDE solvers of [2, 6, 7] and section 7.  Rather than deriving coefficients from the trapezoidal rule, we will obtain them by solving a linear least-squares problem. Typically this will be posed on a discrete subset of $[\kern 0.5pt0,1]$, with exponential clustering of sample points at $z=0$. Experiments show that an effective choice of singularities $\{s_{k}\}$ for such approximations lies on a parabolic contour in the style of [24, 28],

For simplicity, however, we will also sometimes use a configuration in which all the points $s_{k}$ are taken to be confluent at a single value $s_{0}$ of scale $O(n)$; we normally take $s_{0}=n/2$. In this case (4) changes to the confluent form

where $p$ is a polynomial of degree $n$. As we shall prove in section 4, these methods give exponential convergence as $n\to\infty$.

We illustrate these ideas with a MATLAB code segment whose fifth line realizes the distribution (5):

          
          
                M = 1000;
                z = logspace(-24,0,M)’;
                f = sqrt(z);
                n = 10;
                kk = 1:n; tk = -pi + 2*pi*(kk-.5)/n; sk = (n/4)*(1+1i*tk).^2;
                A = [z.^0 1./(log(z)-sk)];
                c = A\f;
                g = A*c;
          

The second row of Figure 1 shows corresponding results for the confluent-pole approximation (6) with $s_{0}=n/2$. These poles are further from optimal, and the exponential convergence is about 15% slower. Certain fine points associated with confluent computations are discussed in section 6. To minimize the effects very near the singularity discussed there and make the comparison of Figure 1 a fair one, the least-squares fitting for this figure was done on a grid of $2000$ points logarithmically spaced in $[10^{-100},1]$.

Figure 2 shows error curves $g(z)-\sqrt{z}$, $z\in[10^{-20},1]$, for approximations (4). On the left these are the approximations computed by least-squares fitting with the singularities $\{s_{k}\}$ of (5), as in the first row of Figure (1). Note that the maximal error decreases exponentially with $n$ and that the error is much smaller for $z\approx 0$ than for $z\approx 1$, because the clustered grid puts more weight there. The right column shows error curves for approximations with the same $\{s_{k}\}$ whose coefficients have been modified to achieve $L^{\infty}$ optimality by 20 steps of a linear Lawson iteration, also known as iteratively reweighted least-squares [9]. These equioscillatory error curves are elegant, but we doubt whether the “optimal” approximations are superior in practice. Their maximal error is lower by only about a factor of 2, at the cost of greatly worsening the accuracy for $z\approx 0$. See chapter 16 and Myth 3 of Appendix A in [19].

The poles of our approximations $g(z)$ lie exponentially close to the origin in the $z$-plane, at least when the formula (5) is used. Exponential clustering of poles is familiar in rational approximation theory [22]. Another phenomenon appears in Figure 2 that is unfamiliar, however, and that is doubly-exponential clustering of oscillating extrema. To see this, consider the leftmost minimum of each error curve in the right column of the figure. With $n=2,$ $4$, and $8$, these minima lie at approximately $z=10^{-1.5}$, $10^{-4.5}$, and (as shown by computations on a longer interval) $10^{-25}$. A rough model of such behavior would be $z_{\min{}}=\exp(-1.7^{\kern 0.3ptn})$, which suggests values for $n=16$ and 32 on the order of $10^{-2000}$ and $10^{-10{,}000{,}000}$. In the figure, we see nothing like this, since the computations are carried out on $[10^{-20},1]$, and it would be impossible to resolve much beyond $[10^{-300},1]$ in standard IEEE floating point arithmetic. However, this probably does not matter, as the approximants on $[\kern 0.5pt0,1]$ will not differ much from those associated with a truncated interval such as $[10^{-20},1]$. We discuss these matters further in section 6.

All these approximants are linear, based on poles $s_{k}$ fixed a priori. What about true minimax approximations, with $\{s_{k}\}$ chosen optimally to minimize the error? For the current setting of $[\kern 0.5pt0,1]$, everything is known about such approximations since they are transplants of minimax rational approximants on $[-\infty,0\kern 0.5pt]$ [19, chap. 24]. As illustrated in Figure 3, which was computed by the method outlined on pp. 217–219 of [19], the errors will be about the squares of what we have seen before, and now there are $2n+2$ equioscillatory extreme points. We shall not discuss such approximations further, because they become impractical once one moves to planar domains.

We now turn to domains in the complex plane. In the basic case one has a closed Jordan region $E$ with piecewise analytic boundary $\Gamma$ and a function $f$ continuous in $E$ and analytic in the interior. We suppose that $f$ is analytic on $\Gamma$ except for branch point singularities at a finite set of points $\{z_{j}\}$, $1\leq j\leq m$, which may for example be corners. As mentioned in the introduction, we consider approximations to $f$ of the generic form

where $p_{0}$ is a degree $n_{0}$ polynomial, with $N=1+\sum n_{j}$. Each function $\log(z-z_{j})$ denotes a fixed choice of a branch of the logarithm that is continuous in $E$. Since standard hardware puts the branch of the logarithm on the negative real axis, it is often convenient in a computer code to work with rotated terms of the form $\log(e^{i\theta_{j}}(z-z_{j}))$ for certain real numbers $\theta_{j}$. We also employ the confluent variant

where each $p_{j}$ is a polynomial of degree $n_{j}$. In both (7) and (8), we often take $n_{0}=0$, so that $p_{0}$ reduces to a constant.

The essential point, to be proved in the next two sections, is that exponential or near-exponential convergence is guaranteed so long as the numbers $n_{j}$ (for $j\geq 1$) increase in proportion to $N$ as $N\to\infty$, in the absence of rounding errors of course. In this section we will explore just one example defined on the domain shown in Figure 5, a two-thirds bite of the unit disk. The function is

with $\omega=e^{\pi i/3}$, with branch point singularities at each of the corners. The boundary $\Gamma$ is discretized by 500 points along each segment exponentially clustered down to distances from the corners of about $10^{-14}$.

As indicated in Figure 4, we computed three approximations to $f$ for values of $N$ up to $150$. One curve shows the accuracy of AAA rational approximations of degree $N-1$. (AAA approximations are rational approximations that are computed very quickly by an algorithm described in [12], implemented in the aaa command of Chebfun; we calculated them with the “cleanup” option turned off.) The AAA approximations are close to optimal when the degree is not too large, and this curve shows root-exponential convergence down to an error on the order of $10^{-8}$. Another curve shows lightning rational approximations [6] of the form

with $\sum n_{j}=N-1$, where $p_{0}$ is a polynomial of degree $n_{0}$; we take $n_{0}=N/4-1$ and $n_{j}=N/4$ for $j=1,2,3$. (The choice $n_{0}=0$, which is often a good one for (7) and (8), tends to slow down convergence significantly for (10).) The poles were fixed at distances $d_{jk}=\exp(4(\sqrt{\vphantom{n_{j}}k}-\sqrt{\vphantom{k}n_{j}}\kern 1.0pt))$, $1\leq k\leq n_{j}$ from the corners $z_{j}$ [6, eq. (3.2)]. Note that the convergence is again root-exponential, at a rate five or six times slower than with AAA with its adaptively determined poles. Clearly AAA has some advantages, but it has drawbacks of sometimes greater sensitivity to rounding errors and occasional placement of poles inside a domain where they are not wanted. More fundamentally, as we shall discuss in section 7, lightning methods generalize immediately to approximation of the real part for solving Laplace problems, whereas AAA does not.

The final curve in Figure 4 shows log-lightning approximations (8) with $s_{j}=n_{j}/3$, now with $n_{0}=0$ and $n_{j}=(N-1)/3$ for $j=1,2,3$. The convergence appears exponential to around $10^{-9}$, after which it slows for reasons at least partly related to rounding errors.

Figure 5 displays phase portraits [27] for $f$ and its three approximations with $N=100$. The rational approximations simulate branch cuts by strings of poles, whereas the reciprocal log approximation incorporates true branch cuts. For $f$ itself, the shapes of the branch cuts are arbitrary. If we had coded the function differently, their directions could have been altered, for example, but this would have had no effect on the values of $f$ on $\Gamma$ and therefore on the three approximations and the other three images. For the lightning and log-lightning approximations, the forms of the branch cuts are fixed by our computation: they have been set to be straight lines bisecting the exterior angles, the natural default choice in the absence of other information. For the AAA approximation, the branch cuts are presumably more nearly optimal, curved in a fashion analyzed by Stahl for the case of Padé approximation [15].

In this section we consider the reciprocal-log approximation of a function $f(z)$ with a branch point singularity at $z=0$. By arguments adapted from those of section 2 of [6], we show that if $f$ can be analytically continued along any contour in the complex plane avoiding $z=0$, then exponentially convergent approximations exist with the singularities $\{s_{k}\}$ extending a distance $O(n)$ into the right half-plane (Thm. 2). If the analyticity assumption holds just in a bounded neighborhood of $0$, then near-exponentially convergent approximations exist with $\{s_{k}\}$ of bounded positive real part (Thm. 6).

To make these ideas precise the theorems of this and the next section use the following definition, which covers the familiar singularities of type $z^{a}(\log z)^{b}$ [10, 26]. The first condition, (11), essentially asserts that $f$ can be analytically continued along arbitrary contours in a punctured neighborhood of each singularity $z_{k}$. But all we really need is continuation along curves that spiral in logarithmically to $z_{k}$, so that after the $s=\log(z-z_{k})$ change of variables they lie in a $V$-shaped wedge in the $s$-plane, and that is why (11) takes the form it does. The second condition, (12), asserts that $f$ is Hölder continuous at $z_{k}$, so that with $s=\log(z-z_{k})$ we get exponential decay as $\hbox{\rm\kern 0.3ptRe\kern 1.0pt}s\to-\infty$. A function like $(z-z_{k})^{a}\log(z-z_{k})$, however, will grow slowly as $z$ winds around $z_{k}$, and that is why (12) includes the factor involving $\hbox{\rm\kern 0.3ptarg\kern 1.0pt}(z-z_{k})$.

For our first theorem, it is enough to use confluent singularities and consider approximations of the form (6) with $s_{0}=\sigma n$,

where $p$ is a polynomial of degree $n$ and $\sigma>0$ is a constant. This result is sharp in the sense that faster than exponential convergence is in general not possible, as we know from the equivalence of the problem of reciprocal-log approximation of $z^{a}$ on $[\kern 0.5pt0,1]$ to that of rational approximation of $e^{as}$ on $[-\infty,0\kern 0.5pt]$. The condition (14) looks almost like a repeat of (12), but (12) is a local condition, only needing to apply for $z$ near $z_{k}$, whereas (14) applies even as $z\to\infty$. With $s=\log z$, it is designed to ensure that $F(s)$ is under control in the right half of the $s$-plane as well as the left.

By the change of variables $s=\log z$ and subtraction of $f(0)$, Theorem 1 is equivalent to the following theorem on rational approximation. As just mentioned, (16) constrains $F(s)$ in the left half $s$-plane, and (17) in the right. The function $e^{s}+s\kern 0.3pte^{2s}$, for example, satisfies the conditions with $a<1$, $a^{\prime}>2$, and $b,b^{\prime}>1$.

for some $C>0$, where $\|\cdot\|$ is the supremum norm on $E$.

Thus our job is to prove Theorem 3. The proof will be based on the Hermite integral formula for rational interpolation [25, Thm. 2, Chap. 8]. Neither the poles $\{s_{k}\}$ nor the interpolation points $\{\alpha_{k}\}$ in the statement below need to be distinct, and an interpolation point of multiplicity $\nu>1$ is interpreted, as usual, to mean interpolation at that point of $f,f^{\prime},\dots,f^{(\nu-1)}$. By a Hankel contour, we mean a continuous curve winding counterclockwise around $(-\infty,0\kern 0.5pt]$ from $\infty$ to $\infty$.

In proving Theorem 3, we will suppose $E=(-\infty,0\kern 0.5pt]$ without loss of generality; the estimates change only by constant factors for other domains since the arguments involve a geometry that scales with $n$ as $n\to\infty$. We will use this choice of interpolation points in $(-\infty,0\kern 0.5pt]$:

These are derived from potential theory [11, 25], where poles and interpolation points are interpreted as point charges of opposite signs, and a minimal-energy charge configuration gives asymptotically optimal approximations as $n\to\infty$. The function

maps the unit disk in the $t$-plane to the $s$-plane slit along $(-\infty,0\kern 0.5pt]$, with $t=0$ mapping to $s=n\kern 0.5pt\sigma$. Thus the equilibrium distribution of interpolation points on the unit circle corresponding to poles at $t=0$, namely the uniform distribution, maps to the distribution determined by (23) in the $s$-plane slit along $(-\infty,0\kern 0.5pt]$ with poles at $s=n\kern 0.5pt\sigma$, as in the function (18). (The $\{\alpha_{j}\}$ are called Fejér–Walsh points [17].) The following lemma, illustrated in Figure 6, was provided to us by Peter Baddoo.

Proof of Theorem 3. As mentioned above, we take $E=(-\infty,0\kern 0.5pt]$ without loss of generality, and we choose $\Gamma$ for Lemma 4 as a V-shape passing through $n\kern 0.5pt\sigma$ at a fixed acute angle, as sketched in Figure 7. (Any acute angle will suffice.) We divide $\Gamma$ into a “head” in the right half-plane, a “tail” with $\hbox{\rm\kern 0.3ptRe\kern 1.0pt}t\leq-n$, and a “middle” with $-n\leq\hbox{\rm\kern 0.3ptRe\kern 1.0pt}t\leq 0$, and show that their contributions $I_{\scriptsize\rm head}$, $I_{\scriptsize\rm tail}$, and $I_{\scriptsize\rm middle}$ to (20) are each exponentially small. (A quantity depending on $n$ is exponentially small if it is $O(\exp(-Cn))$ as $n\to\infty$ for some $C>0$ and exponentially large if its reciprocal is exponentially small.) Note that the denominator $t-s$ of (20) is bounded below in absolute value, $|\phi(s)|$ is bounded by $1$ by (25), and $2\pi$ is just a constant. Thus it is enough to show that each of the three integrals

is exponentially small. In the tail, $F(t)$ is exponentially small and decreases exponentially as $\hbox{\rm\kern 0.3ptRe\kern 1.0pt}t\to-\infty$ by (16), whereas $1/|\phi(t)|$ is bounded by Lemma 5; thus $I_{\scriptsize\rm tail}$ is exponentially small. In the middle, $F(t)$ grows at most algebraically with $n$ by (16), whereas $1/|\phi(t)|$ is exponentially small by Lemma 5; thus $I_{\scriptsize\rm middle}$ is exponentially small. In the head, $F(t)$ is exponentially large, but by (17), at most $O(\exp(A\kern 0.7pt\sigma n))$. On the other hand, by Lemma 5 again, $1/|\phi(t)|$ is exponentially small. In particular, $|\phi(t)|>\exp(d\kern 0.3ptn)$ uniformly in this region for some $d>0$ that is independent of $\sigma$. It follows that if $\sigma<d/A$, then $I_{\scriptsize\rm head}$ too is exponentially small.

We now turn to the second theorem about approximation near $z=0$, assuming analyticity only in a bounded domain. As before, our prototypical example for discussion without loss of generality is $E=[\kern 0.5pt0,1]$. For the argument we again use a V-shaped contour $\Gamma$ for the Hermite integral, as shown in Figure 8, but now, $\Gamma$ must be fixed rather than growing with $n$. This necessitates a change in the choice of poles. If the poles are all at the apex, then the same argument as before ensures exponentially small contributions to the error from $I_{\scriptsize\rm head}$ and $I_{\scriptsize\rm middle}$, but no longer so small (just root-exponential) from $I_{\scriptsize\rm tail}$. Instead it is necessary to distribute poles along $\Gamma$, and a successful strategy is to distribute them along segments extending a distance $\sigma L=\sigma\rho\kern 1.0ptn$ from the apex on the two sides of $\Gamma$, for any constant $\rho>0$.

Specifically, one way to distribute poles is by solving a potential theory problem in the infinite $V$ domain. As suggested by the contour lines $u=0.1,0.2,\dots,0.9$ in Figure 8, let $u$ be the unique harmonic function in the infinite slit wedge taking values $u=0$ on $(-\infty,0\kern 0.5pt]$ and $u=1$ on the segments of length $\sigma L$ just mentioned, with a homogeneous Neumann condition on the remainder of the boundary. This function $u$ can be calculated by means of a conformal map of the lower half of the slit wedge onto a rectangle $-K<\hbox{\rm\kern 0.3ptRe\kern 1.0pt}w<K$, $0<\hbox{\rm\kern 0.3ptIm\kern 1.0pt}w<K^{\prime}$ in the complex $w$-plane, with the Dirichlet boundary components mapping to the vertical sides and the Neumann components to the top and bottom; the ratio $K^{\prime}/2K$ is known as the conformal module of the rectangle. If the interior angle of the wedge is $\pi/\mu$ for some $\mu>1$, the map has the explicit representation

where sn is the Jacobi elliptic sine function and

The numbers $K$ and $K^{\prime}$ are the complete elliptic integrals of the first kind with parameters $M^{-2}$ and $1-M^{-2}$, respectively [14]. Our Fejér–Walsh choice of poles and interpolation points for the theorem will be

and

Theorems 2 and 6 are local assertions, establishing exponential and near-exponential resolution of isolated singularities. It remains to show that global approximations can be constructed by adding together these local pieces. We do this in the style of Theorem 6, requiring analyticity just in a neighborhood of $E$. The argument is adapted from the discussion around Theorem 2.3 in [6] and illustrated schematically in Figure 9.

As mentioned in the penultimate paragraph of section 2 and illustrated in Figure 3, reciprocal-log approximations have surprising properties near the singularities. Whereas error curves for rational approximations vary over a scale exponentially close to the singularity, for reciprocal-log approximations this becomes doubly-exponential. Fortunately, so far as we are aware, these effects need not cause difficulties in using these approximations, and in particular, it appears that they do not necessitate the use of extended-precision arithmetic. We shall now explain our understanding of these matters with reference to Figure 10, which presents four approximations of $\sqrt{z}$ for $z\in[\kern 0.5pt0,1]$. In each case the approximation is computed over $[10^{-20},1]$ and then the absolute value of the error is plotted over $[10^{-80},1]$. For comparison, fine dots show the corresponding results for computations over $[10^{-80},1]$.

Image (a) shows the best (minimax) reciprocal-log approximation with $n=4$. We see here that the best approximation over $[10^{-20},1]$ is by no means optimal over $[10^{-80},1]$. The error is 10 times larger there (6.2e-4 vs. 5.0e-5), and this ratio will grow rapidly with increasing $n$. One might think this implies that reciprocal-log approximations must be ineffective at approximating all the way up to the singularity, but the small dots in the figure, corresponding to the best approximation over $[10^{-80},1]$, contradict this expectation, showing an error of just 7.6e-5. As can be seen in Figure 3, the error in best approximation over all of $[\kern 0.5pt0,1]$ is also not much larger, just 8.7e-5.

Image (b) shows the error for reciprocal-log approximation (4) with poles positioned on a Hankel contour according to (5), now with $n=8$. Again the error in the left-hand part of $[10^{-80},1]$ is larger than it needs to be, but this may be less important for these approximations since the error is dominated by larger values of $z$.

Image (c) shows a very different situation for reciprocal-log approximation (6) with confluent poles at $s_{0}=n/2$. Now the approximation over $[10^{-20},1]$ gives an error over $[10^{-80},1]$ that is orders of magnitude too large.

It is at this point that we wish to pause and ask, what might be the practical implications of such behavior? It certainly seems disturbing if an approximation is much less accurate on $[10^{-80},1]$ than on $[10^{-20},1]$. There are three rather disparate observations to be made, which combine to an encouraging picture.

First we note that in IEEE floating point arithmetic, whereas one can represent numbers near $0$ smaller than $10^{-300}$, the resolution near other points $z_{k}$ is only on the order of 16 digits. Thus even if we wanted to compute approximations by least-squares fitting on a grid in $z_{k}+(0,10^{-20}\kern 0.8pt]$, we could not do so. This might seem worrying, if one takes Figure 10 to suggest that such domains will be required.

The second observation is that as a purely practical matter, they should not be required. On a planar domain $E$ such as that of Figure 5, for example, if a function $f$ is approximated to a satisfactory precision everywhere except at distances $<10^{-20}$ from the vertices $\{z_{k}\}$, what difference does it make? What is wrong with an approximation that is mathematically inaccurate in theory, but only at points where it will never be evaluated?

And yet the use of such approximations would be troubling. This brings us to the third observation, which is the basis of image (d) of Figure 10: it may be possible to make an approximation accurate on $z_{k}+(0,10^{-20}\kern 0.8pt]$ without doing any arithmetic there. The trick is to forget least-squares gridding in $z_{k}+(0,10^{-20}\kern 0.8pt]$ and instead impose additional conditions at $z_{k}$, forcing the singular part of the approximation to decay there at a rate $O(1/(\log|z-z_{k}|)^{J})$ for some power $J>1$. This is natural since in floating point arithmetic, all of $z_{k}+(0,10^{-20}\kern 0.8pt]$ rounds to $z_{k}$ anyway. Image (d) corresponds to an approximation of this kind with $J=2$: we approximate in the $n$-dimensional space spanned by $1/(\log(z)-n/2)^{j}$ with $2\leq j\leq n+2$. The work and the number of degrees of freedom are unchanged, and now least-squares fitting on $[10^{-20},1]$ gives an approximation that is accurate on all of $[\kern 0.5pt0,1]$. We call this a pinned approximation, since certain terms of the approximation are constrained to be zero at the singularity. In this example $J=2$ is effective, but for larger $n$ a choice closer to $J=n/2$ will be called for. To estimate this number, one can use analysis such as that of (22) or (32) to monitor how many interpolation points might be expected to fall closer to $z_{k}$ than about machine precision. For the problem of image (d) with $n=8$, (22) gives distances $\exp(\alpha_{j})\approx 0.8,0.6,0.3,0.06,0.003,10^{-5},10^{-13},10^{-56}$, and our pinned approximation might be thought of as effectively replacing the last two of these numbers by zero. There are many mathematical questions here and it would be interesting to investigate them.

The motivation for this work has been the development of numerical methods for the solution of partial differential equations on planar domains with corners, beginning with the Laplace equation with Dirichlet boundary conditions. For these problems we do not know an analytic function on $E$, just its real part on the boundary. This is no difficulty for linear least-squares fitting, however, where all that matters is having an efficient approximation space. The idea of solving Laplace problems in this fashion was presented in [6], where root-exponential convergence of approximations based on rational functions is established theoretically and experimentally. Here we move to reciprocal-log approximations, and the theorems of sections 4 and 5 assert that the convergence should improve to exponential or near-exponential.

Figure 11 presents a pair of curves confirming this expectation. The solutions computed here are for the problem posed on NA Digest in December 2019, defined by the L-shaped region with vertices $0,2,2+i,1+i,1+2i,2i$ and boundary data $u(z)=\hbox{\rm\kern 0.3ptRe\kern 1.0pt}(z^{2})$ [18]. As usual, $N$ denotes the total number of degrees of freedom, which in this plot is a number of the form $14n+1$ with $n=2,4,6,\dots.$ For each $n$, the approximating function consists of the real part of a polynomial of degree $n$ together with $n$ singularities $\{s_{k}\}$ at each of the six corners as defined by (5), though with a factor $n/3$ instead of $n/4$ in front. This gives $7(n+1)$ complex degrees of freedom, hence $14n+2$ real ones, and the number reduces to $14n+1$ since the imaginary part of the constant term of the polynomial does not affect the fit.

This experiment is a long way from the software that has been developed for lightning approximation [20]. The boundary has been resolved simply by 500 exponentially clustered points on each side, whereas adaptive software can make the gridding both more efficient and more careful. More importantly, each corner has been allocated the same number $n$ of singularities, whereas adaptive software will put more singularities at vertices with stronger singularities. So although the log-lightning convergence curve of Figure 11 is promising, the reliable fast way to solve Laplace problems at this date is still by means of the laplace code at [20]. This code also addresses a number of variations such as Neumann boundary conditions, discontinuous data, and piecewise smooth curved boundaries.

The number of degrees of freedom $N$ is not a direct measure of computer time. For reasons of linear algebra, the work actually increases in proportion to $N^{2}$ or $N^{3}$, depending on whether the grid is fixed or refined with $N$, which tends to make the difference between the lightning and log-lightning computations greater than it may appear in Figure 11. A consideration pushing in the other direction is that computer evaluation of a complex logarithm $\log z$ is typically slower than that of a reciprocal $1/z$ by a factor such as $3$ or $4$.

It was shown in [7] that the original lightning Laplace solver can be generalized to a lightning Helmholtz solver by replacing poles by Hankel functions. Presumably there is a Helmholtz generalization of the log-lightning method too, but we have not investigated this.

Exponential or near-exponential convergence of approximations to branch point singularities is a new phenomenon. Reciprocal-log approximations are interesting in theory and may have consequences in practice if good software can be developed.

Philosophical questions are attached to this kind of approximation. It is an old idea that polynomials and rational functions have a special status since “the only operations that can really be carried out numerically are the four elementary operations of addition, subtraction, multiplication and division.”¹¹1This quote comes from the PhD thesis of Hilbert’s student Kirchberger in 1902 [19, pp. 196–197]. This point of view has contributed to the dominance of polynomials and rational functions in approximation theory, with other classes of approximating functions seeming less fundamental. We are not ourselves immune to the impression that there is something contrived about approximations of the forms (1) and (2). But the philosophical distinction that may have seemed sharp before the arrival of computers seems less sharp today.

A fascinating aspect of reciprocal-log approximations is their exploitation of the behavior of a function on a Riemann surface. As one of the smallest consequences of this feature, these approximations have no difficulty at all in treating domains with slits, and the prospects for further adventures with multivalued functions seem enticing.

The research for this article was carried out in regular communication with André Weideman, whose many suggestions we are happy to acknowledge. Another key contribution was Lemma 5, from Peter Baddoo. In addition we are grateful for helpful comments from Marco Fasondini, Abi Gopal, Alex Townsend, and Elias Wegert.