Rational Functions and Beyond

By Lloyd N. Trefethen

Real analysis, complex analysis, and numerical analysis all start from polynomials. Since Newton at least, polynomials and their limits as Taylor series have captured the local behavior of functions, and the approximations become global with numerical methods such as data fitting, quadrature formulas, spectral methods, “proxy” rootfinding, and Chebfun. For smooth functions on bounded domains, this is often all you need — with the footnote that if a function is periodic, it is advantageous to switch from polynomials to Fourier series.

A rational function is a quotient \(r(z) = p(z)/q(z)\) of polynomials. We say that \(r\) is of degree \(n\) if it can be written in this way, where \(p\) and \(q\) are of degree at most \(n\). Clearly \(r\) may have \(n\) zeros and \(n\) poles, real or complex, and a polynomial is nothing more than a rational function whose poles are constrained to all lie at \(z=\infty\). But \(p/q\) is often not the best representation conceptually or numerically. A better start for many purposes is to write \(r\) in partial fractions:

\[r(z) = c_0 + \sum_{k=1}^n {c_k\over z-z_k}.\tag1\]

This generic representation cannot treat poles of order higher than \(1\) exactly—or poles at infinity—but it can approximate them arbitrarily closely.

Formula \((1)\) highlights the fact that a rational function can have its poles \(\{z_k\}\) anywhere. In particular, if \(r\) approximates a function \(f\) with a singularity at \(z=z_0\), it may achieve great accuracy by clustering its poles near \(z_0\). This effect hit the headlines with Donald Newman’s 1964 paper about the approximation of \(f(x) = |x|\) on \([-1,1]\), or equivalently, the approximation (with half the degree) of \(f(x) = \sqrt x\) on \([0,1]\). Newman showed that by clustering their poles and zeros exponentially near the singularity, rational functions can achieve root-exponential convergence: errors of order \(O(\exp(-C\sqrt n))\) for some \(C > 0\). For six-digit accuracy in this example, a rational function of degree \(n=26\) is enough, whereas a polynomial would need \(n \approx 280,000\). 

All of this is well established and extensively studied by rational approximation theorists. But can one apply these ideas to solve computational problems? Tradition regards rational functions as a computational minefield because of complications like “Froissart doublets” (pole-zero pairs), both theoretical and numerical. For example, in their work during the 1980s and 1990s, Richard Varga, Arden Ruttan, and Amos Carpenter used up to 200-digit precision. But in the past four years and in collaboration with Yuji Nakatsukasa, Abinand Gopal, and others, I have been part of some developments that appear to be changing this picture — all in ordinary floating-point arithmetic.

First came the AAA algorithm (adaptive Antoulas-Anderson), a method of unprecedented speed and robustness that finds near-best rational approximations [4]. The magic of AAA comes from its barycentric representation of rational functions in a third fashion, as a quotient of two partial fractions:

\[r(z) = \sum_{k=0}^n {a_k\over z-t_k} \bigg/\sum_{k=0}^n {b_k\over z-t_k}.\tag2\]

The numbers \(\{t_k\}\) are not the poles of \(r\), but rather a set of arbitrary support points that are chosen to enable numerical stability, even when the poles and zeros are exponentially clustered. For example, suppose that we execute the following commands in Chebfun:

    Z = rand(1000,1) + 1i*rand(1000,1); 
    F = sqrt(Z-Z.^2)./(Z-2);
    [r,poles] = aaa(F,Z);
    plot(Z,’.k’), hold on

    phaseplot(r,[-0.5 2.5 -1.5 1.5])

In \(1/50\) seconds (s) on my laptop, AAA computes a degree-20 rational approximation that matches \(F\) to 14 digits at the 1,000 random data points. Figure 1 shows how the poles and zeros of \(r\) simulate a circular branch cut that connects the branch points \(0\) and \(1\). There is also a pole at \(z = 2.00000077\), whose six-digit agreement with the pole of \(f\) demonstrates the power of rational functions for extrapolation and analytic continuation. In fact, most existing methods for estimation of poles and acceleration of convergence (Padé, Aitken, eta, epsilon, etc.) are based on rational functions.

Next came lightning approximation [2]. AAA has free poles that are adaptively determined, and there is no guarantee that they will not fall in a region where one wants analyticity. Moreover, we do not know how to use AAA to approximate harmonic functions, i.e., real parts of analytic functions, which is what one needs to solve Laplace problems via function approximation. The new idea is to leverage our understanding of exponential clustering at singularities by prescribing the poles \(\{z_k\}\) a priori — simply placing them in a configuration with exponential clustering at the corners, which are likely to be branch points of the solution. A polynomial term to handle “the smooth part of the problem” is included as well. Now one has the linear problem of finding good coefficients \(c_k\) for an approximation \((1)\) with known points \(\{z_k\}\), which is readily solved by least-squares fitting. If you want to fit the real part to solve a Laplace boundary value problem, this makes no significant difference to the calculation. Figure 2 depicts the solution to a Laplace problem on an octagon computed by this method. In \(0.13\) s on a laptop, the code laplace.m has computed a solution with \(138\) poles that is accurate to six digits, all the way up to the corner singularities. Being just the real part of a rational function, the solution can be evaluated in \(11\) \(\mu\textrm{s}\) per point. For 10-digit accuracy, these figures change to \(0.9\) s, \(286\) poles, and \(26\) \(\mu\textrm{s}\). Similar computations are possible for biharmonic problems, which are reducible to harmonic functions, and Helmholtz problems, where the poles of a rational function become center points of shifted Hankel functions [1].

I want to finish with a third development that turned up unexpectedly just a few months ago [3]. What about the use of functions other than rationals, with singularities other than the poles of \((1)\)? In a new development, we have found that reciprocal-log or log-lightning approximations of the form

\[r(z) = c_0 + \sum_{k=1}^n {c_k\over \log(z-z_k)-s_k}\tag3\]

can speed up convergence from root-exponential to exponential or exponential-minus-log, i.e., \(O(\exp(-C n/\log n))\). The approximations take advantage of analyticity on a Riemann surface and can be used for analytic continuation to other Riemann sheets beyond the branch cuts. Figure 3 shows an approximation to the function \(z^{1/3}\) of this kind with \(n=30\), \(z_k = 0\), and \(s_k = \frac{1}{2} n (1+it_k)^2\),  \(t_k = -\pi+ 2\pi (k-\frac{1}{2})/n\), \(1\le k\le n\). 

It seems that a new era of numerical computation with rational functions and other functions with singularities is arriving. This short essay is confined to scalar problems, but there are also exciting ongoing developments that involve rational functions in large-scale linear algebra. Some key names are Athanasios Antoulas, Christopher Beattie, Bernhard Beckermann, Peter Benner, Vladimir Druskin, Serkan Güğercin, Stefan Güttel, Leonid Knizhnerman, Eric Polizzi, Valeria Simoncini, Alex Townsend, Heather Wilber, and Karen Willcox.

---

This article is based on Nick Trefethen’s John von Neumann Prize Lecture at the 2020 SIAM Annual Meeting, which took place virtually this July. Trefethen’s presentation is available on SIAM’s YouTube Channel. 

The figures in this article were provided by the author.

References
[1] Gopal, A., & Trefethen, L.N. (2019). New Laplace and Helmholtz solvers. Proc. Nat. Acad. Sci., 116, 10223.
[2] Gopal, A., & Trefethen, L.N. (2019). Solving Laplace problems with corner singularities via rational functions. SIAM J. Numer. Anal., 57, 2074-2094.
[3] Nakatsukasa, Y., & Trefethen, L.N. (2021). Reciprocal-log approximation and planar PDE solvers. SIAM J. Numer. Anal. Submitted.
[4] Nakatsukasa, Y., Sète, O., & Trefethen, L.N. (2018). The AAA algorithm for rational approximation. SIAM J. Sci. Comput., 40, A1494-A1522.
[5] Wegert, E. (2012). Visual Complex Functions: An Introduction with Phase Portraits, New York, NY: Springer Birkhäuser.

Nick Trefethen is Professor of Numerical Analysis at the University of Oxford. He served as president of SIAM during 2011-2012.