%% 23. Nonlinear approximation: why rational functions? ATAPformats %% % Up to now, this book has been about polynomials, or in the last chapter, % their transplants. The final six chapters of the book are about rational % functions, which have been a mainstay of approximation theory from the % beginning. Why do rational approximations occupy such a large place in % the literature? Polynomials are familiar and comfortable, but rational % functions seem complex and specialized. Is their position in % approximation theory justified, or is it an artifact of history, perhaps % a holdover from the pre-computer era? In this chapter we attempt to % answer these questions, and in doing so we shall find ourselves % considering the broader question of what the uses are of the whole % subject of approximation theory. %% % I think the answer is this. Although rational functions indeed became an % established part of approximation theory long before computers and many % of today's practical applications, their place in the subject is % deserved. Their importance stems from a conjunction of two facts. On the % one hand, rational functions are more powerful than polynomials at % approximating functions _near singularities_ and on _unbounded domains._ % On the other hand, for various reasons related for example to partial % fraction decompositions, they are easier to work with than their % nonlinearity might suggest---indeed, sometimes no more complicated than % polynomials. %% % A rational function is the ratio of two polynomials, and in particular, % given $m\ge 0$ and $n\ge 0$, we say that $r$ is a rational function of % _type_ $(m,n)$ if it can be written as a quotient % $p_m/q_n$ with $p_m\in {\cal P}_m$ and $q_n\in {\cal P}_n$. The set of % all rational functions of type $(m,n)$ is denoted by ${\cal R}_{mn}$, and % any $r\in {\cal R}_{mn}$ can be written in the form % $$r(x) = \sum_{k=0}^m a_k x^k\kern -4pt \left/ % \sum_{k=0}^n b_k x^k\right. \eqno (23.1)$$ % for some real or complex coefficients $\{a_k\}$ and $\{b_k\}$. The % degrees need not be exact, i.e., there is no requirement that $a_m^{}$ or % $b_n^{}$ must be nonzero. Nor do we require that the numerator and % denominator are relatively prime, that is, that they have no common % zeros. %% % Suppose, however, that for some nonzero $r\in {\cal R}_{mn}$, we choose a % representation with relatively prime numerator and denominator. Define % $\mu \le m$ to be the index of the highest degree nonzero numerator % coefficient and similarly $\nu\le n$ for the denominator, and further % normalize the coefficients by requiring $b_\nu = 1$. Then we can write % $$r(x) = \sum_{k=0}^\mu a_k x^k\kern -4pt \left/ % \sum_{k=0}^\nu b_k x^k \right. , \quad a_\mu\ne 0, ~b_\nu=1. \eqno (23.2)$$ % In this case $r$ has exactly $\mu$ finite zeros and $\nu$ finite poles, % counted with multiplicity: we say that $r$ is of _exact type_ $(\kern 1pt % \mu,\nu)$. (The case in which $r$ is identically zero is a special one, % with no nonzero coefficients in the numerator, and we say it has exact % type $(-\infty,0)$.) If $\mu>\nu$, then $r$ has a pole at $x=\infty$ of % order $\mu-\nu$, and if $\nu>\mu$ it has a zero at $x=\infty$ of order % $\nu-\mu$. Basic properties of rational functions are described in books % of complex analysis such as [Ahlfors 1953, Henrici 1974, Markushevich % 1985]. %% % These representations highlight the nonlinearity of rational functions, % but a different perspective is suggested when we represent them by _partial % fractions._ (Partial fractions were the subject of Jacobi's PhD % thesis , and an excellent general reference is Chapter 7 % of [Henrici 1974].) In the simplest situation, consider % $$r(x) = \sum_{k=1}^n {c_k\over x-\xi_k}, \eqno (23.3)$$ % where $\{\xi_k\}$ are distinct real or complex numbers. For any % coefficients $\{c_k\}$, this is a rational function of type $(n-1,n)$. % The number $c_k$ is the _residue_ of $r$ at $\xi_k$. This % representation highlights the linear aspects of rational functions. For % example, whereas computing the integral of $r$ written in the form $p/q$ % looks daunting, in the representation (23.3) we have simply % $$\int^x r(s)ds = C + \sum_{k=1}^n c_k\log (x-\xi_k). \eqno (23.4)$$ % In applications, it is interesting how often a formula like this turns % out to be instrumental in making a rational function useful. %% % The partial fraction form (23.3) does not apply to all rational % functions. One limitation is that it always represents a rational % function of exact type $(\mu,\nu)$ with $\mu< \nu$. Another is that it % does not represent all functions of this kind, since it cannot account % for poles of multiplicity greater than $1$. The following theorem gives a % partial fraction representation for the general case. %% % % {\bf Theorem 23.1. Partial fraction representation.} % {\em Given $m,n\ge 0$, let $r\in {\cal R}_{mn}$ be arbitrary. Then $r$ has % a unique representation in the form % $$r(x) = p_0(x) + \sum_{k=1}^\mu p_k((x-\xi_k)^{-1}), \eqno (23.5)$$ % where $p_0$ is a polynomial of exact degree $\nu_0$ for some % $\nu_0\le m$ {\rm (}unless $p=0)$ and $\{p_k\}$, $1\le k \le \mu$, are polynomials of exact % degrees $\nu_k\ge 1$ with $p_k(0)=0$ and} % $\sum_{k=1}^\mu \nu_k \le n$. % %% % _Proof._ % See Theorem 4.4h of [Henrici 1974]. % $~\hbox{\vrule width 2.5pt depth 2.5 pt height 3.5 pt}$ %% % The function $p_0$ is the _polynomial part_ of $r$, and % $p_k((x-\xi_k)^{-1})$ is its _principal part_ at $\xi_k$. %% % This is all we shall say for the moment about the mathematics of rational % functions. Let us now turn to the main subject of this chapter, the % discussion of why these functions are useful in approximation theory and % approximation practice. %% % The right place to start is with a cautionary observation. Rational % functions are not always better than polynomials. Indeed, consider the % most basic of all situations, in which a fuction $f$ is analytic in a % $\rho$-ellipse $E_\rho$ for some $\rho> 1$. For such a function, by % Theorem 8.2, polynomial approximations will converge at the rate $O(\kern % 1pt\rho^{-n})$. It turns out that a typical convergence rate for type % $(n,n)$ rational functions is $O(\kern 1pt\rho^{-2n})$. So, doubling % the number of parameters to be determined sometimes just approximately % doubles the convergence rate. (In fact, sometimes it does not increase % the convergence rate at all [Szabados 1970].) For applications of this % kind, rational functions may outperform polynomials, but often it is by a % rather modest factor. %% % For example, here are a pair of curves showing $\|f-p_{2n}^*\|$ (dots) % and $\|f-r_{nn}^*\|$ (stars) as functions of $n$ for $f(x) = \exp(-x^4)$, % where $p_{2n}^*$ and $r_{nn}^*$ are the best approximations to $f$ in % ${\cal P}_{2n}$ and ${\cal R}_{nn}$, respectively. (We shall discuss % rational best approximation in the next chapter.) Both curves decrease % geometrically, and there is not much difference between them. (The % rational approximations here should in principle be computed with % |remez|, but Chebfun's rational remez algorithm is currently not robust % enough, so |cf| is used instead.) %% % \vskip -2em x = chebfun('x'); f = exp(-x.^4); warning off nn = 0:20; errp = []; errr = []; for n = nn p2n = remez(f,2*n); errp = [errp norm(f-p2n,inf)]; [p,q,foo] = cf(f,n,n); rnn = p./q; errr = [errr norm(f-rnn,inf)]; end clf, semilogy(nn,errp,'.-','markersize',12), grid on, ylim([1e-16 10]) hold on, semilogy(nn,errr,'h-r','markersize',4), FS = 'fontsize'; text(10.5,2e-8,'E_{2n,0}',FS,10,'color','b') text(9,1e-11,'E_{n,n}',FS,10,'color','r'), xlabel n title(['Convergence of polynomial and rational '... 'best approxs to exp(-x^4) on [-1,1]'],FS,9) %% % \vskip 1pt %% % What makes rational functions important is that, in contrast to this % example, there are many problems where one wants to operate near % singularities, or on unbounded domains. For these problems, rational % approximations may converge much faster than polynomials. For example, % here is an experiment like the last one, but with $f(x) = |x|$. For this % function, a type $(n,n)$ rational approximant with $n=150$ gives 16-digit % accuracy, whereas polynomial approximants would need $n=10^{15}$ to do so % well. (Again this code should in principle use |remez| but cannot, % so known best approximation errors are hardwired into the code.) %% % \vskip -2em f = abs(x); xx = linspace(-1,1,1000); nn = 0:50; errp = []; errr = [.5 4.37e-2 8.50e-3 2.28e-3 7.37e-4 2.69e-4 1.07e-4 ... 4.60e-5 2.09e-5 9.89e-6 4.88e-6 2.49e-6 1.30e-6 ... 6.3*exp(-pi*sqrt(26:2:max(nn)))]; errr = kron(errr,[1 1]); errr(end) = []; for n = nn p2n = remez(f,2*n); errp = [errp norm(f(xx)-p2n(xx),inf)]; end hold off, semilogy(nn,errp,'.-','markersize',12), grid on hold on, semilogy(nn,errr,'h-r','markersize',4) text(37,3e-4,'E_{2n,0}',FS,10,'color','b') text(21,2e-7,'E_{n,n}',FS,10,'color','r'), xlabel n title(['Convergence of polynomial and rational '... 'best approxs to |x| on [-1,1]'],FS,9) %% % \vskip 1pt %% % The approximation of $|x|$ by rational functions is one of the two % famous problems'' to be considered in Chapter 25. Half a century ago % Donald Newman proved that whereas polynomial approximants to $|x|$ % converge just at the rate $O(n^{-1})$, for rational approximants the rate % is $\exp(-C \sqrt n\kern 1pt)$ with $C>1$ [Newman 1964]. This result % rigorously established the possibility of an exponential difference in % effectiveness of the two types of approximations. %% % The rest of this chapter is devoted to an outline of twelve applications % in which rational approximations are useful. In most of these examples, % there is a singularity or unbounded domain in the picture. The exceptions % are applications #1 and #8, where rational functions outperform % polynomials less decisively. %% % _1. Elementary and special functions._ % Classically, approximation theory brings to mind the problem of designing % subroutines for computers to evaluate elementary functions, like $\sin(x)$, % and special functions, like Airy or Bessel functions. For some of % these applications, especially when the number of digits of accuracy % required is known in advance, rational approximations prove to be the % best choice. A classic project in this line is the SPECFUN software % package [Cody 1993], descendant of the earlier FUNPACK [Cody 1975], which % uses rational best approximations to evaluate Bessel functions, error % functions, gamma functions and exponential integrals to 18 digits of % accuracy. For many years a driving force behind these software products % and an expert on the matter of practical rational approximations was % W. J. Cody at the Argonne National Laboratory; Cody's version of the % rational Remez algorithm is described in [Cody, Fraser & Hart 1968]. For % a presentation of some of the state of the art early in the 21st century, % see [Muller 2006]. %% % _2. Digital filters._ In electrical engineering, the construction of % low-pass, high-pass, and other digital filters often involves % approximation of functions with jumps. (For these problems the % approximation domain is usually the unit circle in the complex plane.) % Jumps amount to singularities on or near the domain of approximation, and % Theorem 8.3 implies that polynomials have no chance of rapid convergence % for such functions. As Newman's theorem would lead us to expect, rational % approximations sometimes do much better. Engineers use the term FIR % (Finite Impulse Response) for polynomial filters and IIR (Infinite % Impulse Response) in the rational case [Oppenheim, Schafer & Buck 1999]. %% % _3. Convergence acceleration for sequences and series._ The mathematical % sciences are full of problems of extrapolation. For example, one might % be interested in $\lim_{h\to 0} f(h)$, where $f(h)$ is a quantity % computed numerically on a grid of spacing $h$. For such a problem, $f$ % is often analytic at $h=0$, in which case _Richardson extrapolation,_ % based on interpolating the data by a polynomial, may be beautifully % effective. On the other hand, suppose we want to evaluate $\lim_{n\to % \infty} a_n$ for a sequence $\{a_n\}$. We can regard this problem too as % $\lim_{h\to 0} f(h)$ with the definition $f(1/n) = a_n$, but now, in many % applications, $f(h)$ will not be analytic at $h=0$ and Richardson % extrapolation will be ineffective. The more powerful extrapolation % methods that have been developed for such problems, such as Aitken % extrapolation and the epsilon algorithm, are mostly based on rational % approximations. See Chapter 28. %% % _4. Determination of poles._ % Suppose a function $f$ is analytic on $[a,b]$ and has some real or % complex poles nearby whose positions and residues are of interest. % Classic examples of such problems arise in the study of phase transitions % in condensed matter physics. If we approximate $f$ by polynomials on % $[a,b]$, then by Theorem 8.3, the convergence fails outside a % $\rho$-ellipse of analyticity, so not much information about poles can be % obtained. If we approximate by rational functions, exponential % convergence to some of the poles can often be achieved. Specifically, a % good strategy is to consider the poles of $r_{mn}$ for moderate values of % $n$, where $r_{mn}$ is a rational approximant to $f$ obtained by % $\hbox{Pad\'e}$ or $\hbox{Chebyshev--Pad\'e}$ approximation or rational % interpolation or least-squares. See Chapters 26--28. %% % _5. Analytic continuation._ % If $f$ is analytic on $[a,b]$, then in many applications it can be % analytically continued, in theory, to the rest of the complex plane, % apart from exceptional points and curves in the form of poles, other % singularities, and branch cuts. Computing such continuations % numerically, however, is a difficult problem. One could try % approximating $f$ by a polynomial, but this approach will be useless % outside the largest Bernstein ellipse in which $f$ is analytic. Rational % functions, by contrast, may be effective for continuation much further % out. Again see Chapter 28. %% % _6. Eigenvalues and eigenvectors of matrices._ % Suppose we want to compute an eigenvector of a matrix $A$. One approach, % the _power method,_ is to pick a starting vector $x$ and compute % $\lim_{n\to \infty}$ $A^n x$, but the convergence of this % polynomial-based idea is very slow in general. A much faster method, % _inverse iteration_, is based on rational approximations: find an % approximation $\mu$ to some eigenvalue $\lambda$ and compute % $\lim_{n\to\infty}$ $(A-\mu I\kern.7pt)^{-n} x$. The convergence gets % faster the closer $\mu$ is to the singularity $\lambda$, and exploitation % of this effect leads to the spectacularly effective QR algorithm for % matrix eigenvalues and eigenvectors [Francis 1961]. Experts in numerical % linear algebra do not usually think about rational approximations when % discussing inverse iteration or the QR algorithm, but such approximations % come explicitly to the fore in the analysis of extensions such as % shift-and-invert Arnoldi or rational Krylov iteration $\hbox{[G\"uttel}$ % 2010]. %% % _7. Model reduction and optimal control._ A major topic in numerical % linear algebra and control theory is the approximation of complex % input-output systems by simpler ones for more efficient computation. Via % the Laplace transform, problems of this kind (in the case of continuous % as opposed to discrete time) can in many cases be reduced to problems of % approximation on the imaginary axis in the complex plane. The unbounded % domain makes rational approximations a natural choice, and in fortunate % cases, a system with hundreds of thousands of degrees of freedom may be % reduced to a model with just dozens or hundreds. One set of methods for % such problems goes by the name of $H_\infty$ approximation, based on % results by Adamyan, Arov and Krein  and Glover  that are % related to CF approximation (Chapter 20). For more % information see [Antoulas 2005, Zhou, Doyle & Glover 1996, % Embree & Sorensen 2012]. %% % _8. Exponential of a matrix._ A famous paper in numerical analysis is % Nineteen dubious ways to compute the exponential of a matrix'', by % Moler and Van Loan in 1978, reprinted in expanded form 25 years later % [Moler & Van Loan 2003]. These authors compared many algorithms for % computing $e^A$ and reached the conclusion that the most effective was a % scaling-and-squaring method based on $\hbox{Pad\'e}$ approximation [Ward % 1977]. Here, first $A$ is scaled so that its norm is on the order of % $1$. Then $e^A$ is approximated by $r(A)$, where $r$ is a type $(n,n)$ % $\hbox{Pad\'e}$ approximant to $e^x$ (Chapter 27). This is an example % where rational approximations outperform polynomials not decisively but % by a more or less constant factor. This approach is used by the matrix % exponential program |expm| in Matlab, which for many years was based on % type $(6,6)$ $\hbox{Pad\'e}$ approximation. A more careful analysis of % the scaling-and-squaring algorithm was later provided by Higham , % who concluded that a better choice was type $(13,13)$, and the |expm| % code was adjusted accordingly in Matlab Version 8. In [Higham & Al-Mohy % 2010, Appendix] the authors conclude that $\hbox{Pad\'e}$ approximants % are up to $23\%$ more efficient than Taylor polynomials in this % application. %% % _9. Numerical solution of stiff PDEs._ A particularly important set of % problems related to matrix exponentials are derived from partial % differential equations. The starting point of such applications is the % Laplace operator $\Delta$ on a spatial domain $\Omega$ with Dirichlet % boundary conditions, which has an infinite set of negative real % eigenvalues diverging to $-\infty$. To solve the heat equation $\partial % u/\partial t = \Delta u$ numerically on $\Omega$ with initial data % $u(x,0) = u_0$, one would like to be able to compute the matrix % exponential product $e^{tA} v_0$, where $A$ is a matrix discretization of % $\Delta$ and $v_0$ is a discretization of $u_0$. The wide range of % eigenvalues makes such a problem stiff'', posing challenges for % numerical methods. One method for coping with stiffness is to find a % rational function $r(x)$ that approximates $e^x$ accurately on % $(-\infty,0\kern .5pt]$, hence in particular at all of the eigenvalues of % $A$, and then to compute $r(tA)v_0$. Polynomials cannot approximate a % bounded function on an infinite interval, but rational functions can. % This problem of rational approximation of $e^x$ on $(-\infty,0\kern % .5pt]$ goes back to Cody, Meinardus and Varga , whose 1/9 % conjecture'', eventually settled by Gonchar and Rakhmanov , is the % other famous problem considered in Chapter 25. Generalizations have % become important in scientific computing in recent years in the design of % _exponential integrators_ for the fast numerical solution of stiff % nonlinear ordinary and partial differential equations [Hochbruck & % Ostermann 2010, Kassam & Trefethen 2005, Schmelzer & Trefethen 2007]. %% % _10. Quadrature formulas._ As we have seen in Chapter 19, a quadrature % formula approximates an integral $I = \int_a^b f(x) dx$ by a finite % linear combination $I_n = \sum_{k=0}^n w_k f(x_k)$. If the weights $w_k$ % are interpreted as residues of a rational function $r(x)$ with poles at % the nodes $x_k$, then by estimation of a Cauchy integral over a contour % $\Gamma$ enclosing $[a,b]$ in the complex plane, one can show that the % error $I-I_n$ is bounded in terms of the size of $f$ in the region % enclosed by $\Gamma$ times the error in approximation of the analytic % function $\log((x+1)/(x-1))$ by $r$ over the same region [Takahasi & Mori % 1971]. So every quadrature formula is connected with a rational % approximation problem. In fact, Gauss's original derivation of the % $(n+1)$-point Gauss quadrature formula on $[-1,1]$ was based on exactly % this connection: he used type $(n,n+1)$ $\hbox{Pad\'e}$ approximation of % $\log((x+1)/(x-1))$ at $x=\infty$ [Gauss 1814]. %% % _11. Adaptive spectral methods for PDEs._ The barycentric interpolation % formula has the form of a rational function that reduces to a polynomial % for a special choice of weights (Chapter 5). Regardless of the choice of % weights, however, one still gets an interpolant, and in some applications % there is no compelling reason to force the interpolant to be a % polynomial. This opens up the possibility of much more flexible rational % interpolants, which have the particular advantage of not being so % sensitive to the distribution of the interpolation points. These ideas % originated with Salzer  and Schneider and Werner , building % on earlier work going as far back as Jacobi , and were later % developed by Berrut , and Floater and Hormann . For ordinary % and partial differential equations, they form the basis of adaptive % spectral methods for solving problems whose solutions have singularities % close to the region of approximation [Berrut, Baltensperger & Mittelmann 2005, % Tee & Trefethen 2006, Hale & Tee 2009]. %% % _12. One-way wave equations._ Our final application became well known in % the 1970s and 1980s [Halpern & Trefethen 1988]. The usual wave equation % permits energy propagation in all directions, but there are applications % where one would like to restrict to half the permitted angles, a % $180^\circ$ range. For example, this idea is useful in underwater % acoustics [Tappert 1977], in geophysical migration [Claerbout 1985], and % in the design of absorbing boundary conditions for numerical simulations % [Lindman 1975, Engquist & Majda 1977]. How can we define a system that % behaves like $u_{tt} = u_{xx} + u_{yy}$ for leftgoing waves, say, with % negative $x$-component of velocity, while not propagating rightgoing % waves? (The subscripts represent partial derivatives.) A Fourier % transform shows that the dispersion relation of such a system should be % $\xi = \omega \sqrt{1-s^2}$, where $s = \eta/\omega$ and $\omega, \xi, % \eta$ are the dual variables to $t,x,y$. Only the positive branch of the % square root should be present, making this system a _pseudodifferential % operator._ However, a rational approximation $\sqrt{1-s^2} \approx r(s)$ % simplifies this to a differential equation. For example, the type $(2,2)$ % $\hbox{Pad\'e}$ approximation $r(s) = (1-{3\over 4} s^2)/(1-{1\over 4} % s^2)$ leads to the PDE $u_{xtt} - {1\over 4} u_{xyy} = u_{ttt} - {3\over % 4} u_{tyy}$, sometimes known as the $\hbox{}45^\circ$ equation'' % because it has high accuracy approximately for angles up to $45^\circ$. % In this application, rational functions are superior to polynomials both % because of higher accuracy in view of the singularities at $s=\pm 1$, and % because polynomial approximations lead to PDEs that are ill-posed % [Trefethen & Halpern 1986]. %% % We have just seen a list of twelve applications. In concluding this % chapter I would like to consider what light these may shed on the biggest % question of all, namely, what is the use of approximation theory? %% % % To see some possible views, let us go back to 1901. That was the year of % Runge's landmark paper (Chapter 13), whose title was\footnote{Title: % \"Uber empirische Funktionen und die Interpolation zwischen % \"aquidistanten Ordinaten.'' First sentence: Die Abh\"angigkeit % zwischen zwei messbaren Gr\"ossen kann, strenge genommen, durch % Beobachtung \"uberhaupt nicht gefunden werden.''} % %% % \em On empirical functions and interpolation between % equidistant ordinates.'' %% % In reading this today, one is struck by the word empirical''. The % empirical theme is echoed in the opening sentence: %% % % \em The relationship between two measurable quantities can, strictly % speaking, not be found by observation. % %% % Runge goes on to mention observations'' six times more in the opening % paragraph. It would seem that his motivation is the processing of % scientific data: interpolation in the traditional sense of evaluating a % function at points lying between those at which it is listed in a table. %% % % The next year, 1902, brought another landmark of approximation theory: % Kirchberger's PhD thesis under Hilbert in G\"ottingen, which included the % first systematic statement and proof of the equioscillation theorem for % polynomial approximation (Theorem 10.1). Here is the first paragraph of % Kirchberger's thesis, as reprinted in the first paragraph of his % published paper a year later , setting forth a clear motivation for % approximation theory. We may imagine that this was probably also % Hilbert's view of the subject.\footnote{Mit dem Begriff der Funktion % ist das Postulat der numerischen Berechnung der Funktionswerte f\"ur % irgendwelche Werte der unabh\"angigen Variabeln gegeben. Da aber die % vier elementaren Spezies der Addition, Subtraktion, Multiplikation und % Division, oder streng genommen nur die erste drei derselben, die einzigen % numerisch ausf\"uhrbaren Rechnungsarten, alle andern aber nur insoweit % durchf\"uhrbar sind, als sie sich auf diese zur\"uckf\"uhren lassen, so % folgt hieraus, dass wir s\"amtliche Funktionen nur insoweit numerisch % beherrschen, als sie sich durch rationale Funktionen ersetzen, d.~h.\ % angen\"ahert darstellen lassen. Hieraus erhellt die gro\ss e Bedeutung % der Ann\"aherungsprobleme f\"ur die gesamte Mathematik und die % ausgezeichnete Stellung, die die Probleme der Ann\"aherung durch % rationale oder ganze rationale Funktionen einnehmen. In der Tat setzt, % wenigstens f\"ur die numerische Berechnung, jede Ann\"aherung durch % andere, z.~B. trigonometrische Funktionen, die ann\"aherungsweise % Ersetzbarkeit dieser Funktionen durch rationale voraus.''} % %% % \em % The notion of a function entails the assumption that a numerical value of % the function can be calculated for any value of the independent variable. % But since the only operations that can really be carried out numerically % are the four elementary operations of addition, subtraction, % multiplication and division, or strictly speaking only the first three of % these, it follows that we are really only masters of more general % functions insofar as we can replace them by rational functions, that is, % represent them approximately. This highlights the great significance of % approximation problems for the whole of mathematics and the special role % of approximation by polynomials and rational functions. Indeed, for % numerical calculation at least, any use of other approximations such as % trigonometric functions presupposes that these can in turn be % approximated by rational functions. % %% % Updated to 2012, we may say that Kirchberger's justification of % approximation theory is all about _machine arithmetic._ Approximation by % polynomials and rational functions is important, he is saying, because % ultimately computers can only carry out polynomial and rational % operations. %% % Both Runge's emphasis on data and Kirchberger's emphasis on arithmetic % capture aspects of approximation theory that remain valid today. In % particular, Kirchberger's paragraph seems a remarkably clear statement of % a justification of approximation theory that in a certain philosophical % sense seems almost unarguable (although the line between primitive'' % operations like $+$ and derived'' ones like $\sin(\cdot)$ is not always % so clear on actual computers, with their multiple levels of hardware, % software and microcode). The same argument is often seen nowadays. %% % Nevertheless, I do not think data analysis or machine arithmetic get at % the heart of why approximation theory is important and interesting. In % fact I don't think Runge's words even capture the truth of why _he_ was % interested in the subject! (He becomes more of a mathematician in the % second half of his paper.) What these observations miss is the importance % of _algorithms._ %% % Let us look again at the list of applications. Kirchberger's motivation % could be said to be on target for #1 and #2 (evaluation of functions, % digital filters), and Runge's for #3, #4, and #5 (extrapolation, % determination of poles, analytic continuation). But the remaining seven % items need to be accounted for in other ways. It is noteworthy that % applications #6 to #9 all involve matrices, sometimes of very large % dimension (eigenvalues and eigenvectors, model reduction, exponentials of % matrices, stiff PDEs). Applications #9 to #12 all involve integrals and % differential equations (stiff PDEs, quadrature, adaptive spectral % methods, one-way wave equations). In most of these problems we seem a % long way from scalars $x$ and $r(x)$: the polynomial and rational % operations are applied to matrices and operators, not just numbers. %% % Chebfun provides another interesting data point (for polynomials rather % than rational functions). Chebfun is built on a century of developments % in polynomial interpolation and approximation, and it makes it possible % to work with univariate functions numerically in almost unlimited ways. A % particularly important Chebfun capability is finding roots of a function % $f(x)$, which enables many further operations like computing extrema, % absolute values, and 1-norms. Chebfun finds the roots by the algorithm % proposed by Good  and Boyd  and described in Chapter 18: % approximate $f$ by polynomial interpolants, then find roots of the % polynomials by computing eigenvalues of colleague matrices. This is as % powerful an application of approximation theory as one could ask for, but % it has little to do with data analysis or machine arithmetic. %% % % Why are polynomial and rational approximations useful? Not because % $r(x)$ is easier to evaluate than $\exp(x)$, but because $r(A)$ is easier % to evaluate than $\exp(A)$, and $r(\partial/\partial x)$ is easier to % evaluate than $\exp(\partial/\partial x)\kern .5pt$! Not because we can % {\em evaluate} $p(x)$, but because we can {\em find its roots\kern % 1.5pt}! % %% % % \begin{displaymath} % \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl % {\sc Summary of Chapter 23.} % Rational functions are more powerful than polynomials for approximating % functions near singularities or on unbounded domains. This is the % reason for their importance in approximation theory and approximation % practice.\vspace{2pt}}} % \end{displaymath} % %% % \small\smallskip\parskip=2pt % \par % {\bf Exercise 23.1. Examples of partial fractions.} % Express the following functions in partial % fraction form: (a) $x^3/(1-x)$ (b) $x/(x^2-4)$, (c) $x^2/(x^2-4)^2$, % (d) $(1-x^3) / (1+x^2)$. % \par % {\bf Exercise 23.2. Uses of partial fractions.} % (a) Express the function $r(x) = (x(x+1)(x+2))^{-1}$ in % partial fraction form. (b) What is its integral from $1$ to $t$? % (c) What is the sum of the infinite series $r(1) + r(2) + r(3) + \cdots\,$? % \par % {\bf Exercise 23.3. Another infinite series.} % (a) Based on numerical experiments, conjecture a value of the % infinite sum $1/(1\cdot 3\cdot 5) % + 2/(3\cdot 5\cdot 7) % + 3/(5\cdot 7\cdot 9) + \cdots .$ (b) Verify your conjecture % with partial fractions. % \par % {\bf Exercise 23.4. A trigonometric identity.} % Verify the identity $1/(1\cdot 3\cdot 5) % - 1/(7\cdot 9\cdot 11) % + 1/(13\cdot 15\cdot 17) - \cdots = \pi/48.$ % \par % {\bf Exercise 23.5. Polynomial vs.\ rational experiments.} Produce % plots comparing $E_{2n,0}(f)$ and $E_{n,n}(f)$ for the following % functions $f$ defined on $[-1,1]:$ (a) $\log(1+x^2)$, (b) % $\tanh(5x)$, (c) $\exp(x)/(2-x)$. % \par % {\bf Exercise 23.6. Approximation of a gamma function.} % Consider the function $f(x) = \Gamma(x+2)$ on $[-1,1]$, which has % simple poles at $x=-2, -3, \dots.$\ \ Determine analytically the % geometric convergence rates to be expected as $m\to\infty$ for rational % approximants to $f$ of types (a) $(m,0)$, (b) $(m,1)$, (c) $(m,2)$. % \par