%% 13. Equispaced points, Runge phenomenon ATAPformats %% % So far in this book, we have considered three good methods for % approximating functions by polynomials: Chebyshev interpolation, % Chebyshev projection, and best approximation. Now we shall look at a % catastrophically bad method!---interpolation in equally spaced points. % This procedure is so unreliable that for generations, it has tainted % people's views of the whole subject of polynomial interpolation. The % mathematical tools we will need to understand what is going on are the % Hermite integral formula and potential theory, as discussed in the last % two chapters. %% % As mentioned in Chapter 5, polynomial interpolation was an established % tool by the 19th century. The question of whether or not polynomial % interpolants would converge to an underlying function as $n\to \infty$ % was not given much attention. Presumably many mathematicians would have % supposed that if the function was analytic, the answer would be yes. In % 1884 and 1896, $\hbox{M\'eray}$ published a pair of papers in which he % identified the fact that certain interpolation schemes do not converge % $\hbox{[M\'eray}$ 1884 & 1896]. In the first paper he writes, %% % % \em It is rather astonishing that practical applications have not yet % turned up any cases in which the interpolation is illusory.\footnote{Il % est assez \'etonnant que les hasards de la pratique n'aient encore fait % conna\^itre aucun cas dans lequel l'interpolation soit illusoire.'' By % illusory, M\'eray means nonconvergent.} % % % % M\'eray's derivations had the key idea of making use of the % Hermite integral formula. However, the examples he devised were rather % contrived, and his idio\-syncratically written papers had little impact. It % was Runge in 1901 who made the possibility of divergence famous by % showing that divergence of interpolants occurs in general even for % equispaced points in an real interval and evaluation points in the % interior of that interval [Runge 1901]. % %% % % Runge illustrated his discovery with an example that has become known as % the {\em Runge function}: $1/(1+x^2)$ on $[-5,5]$, or equivalently, % $1/(1+25x^2)$ on $[-1,1]$: % %% % \vskip -2em x = chebfun('x'); f = 1./(1+25*x.^2); %% % We already know from Chapter 8 that there is nothing wrong with this % function for polynomial interpolation in Chebyshev points: $f$ is % analytic, and the polynomial interpolants converge geometrically. Now, % however, let us follow Runge and look at interpolants in equally spaced % points, which we can compute using the Chebfun overload of Matlab's % |interp1| command. %% % Here is what we get with 8 points: %% % \vskip -2em s = linspace(-1,1,8); p = interp1(s,f,domain(-1,1)); hold off, plot(f), hold on, plot(p,'r'), grid on plot(s,p(s),'.r'), axis([-1 1 -1 3]), FS = 'fontsize'; title('Equispaced interpolation of Runge function, 8 points',FS,9) %% % \vskip 1pt %% % Here is the result for 16 points: %% % \vskip -2em s = linspace(-1,1,16); p = interp1(s,f,domain(-1,1)); hold off, plot(f), hold on, plot(p,'r'), grid on plot(s,p(s),'.r'), axis([-1 1 -1 3]) title('Equispaced interpolation of Runge function, 16 points',FS,9) %% % \vskip 1pt %% % Is this going to converge as $n\to\infty$? Things look pretty good in the % middle, but not so good at the edges. Here is the result for 20 points: %% % \vskip -2em s = linspace(-1,1,20); p = interp1(s,f,domain(-1,1)); hold off, plot(f), hold on, plot(p,'r'), grid on plot(s,p(s),'.r'), axis([-1 1 -1 3]) title('Equispaced interpolation of Runge function, 20 points',FS,9) %% % \vskip 1pt %% % What is happening is exponential convergence in the middle of the % interval but exponential divergence near the ends. The next figure shows % the maximum error over $[-1,1]$ as a function of the number of % points. We get a hint of convergence at first, but after $n=10$, % things just get worse. Note the log scale. %% % \vskip -2em ee = []; nn = 2:2:50; for np = nn s = linspace(-1,1,np); p = interp1(s,f,domain(-1,1)); ee = [ee norm(f-p,inf)]; end hold off, semilogy(nn,ee,'.-'), grid on, axis([0 50 5e-2 2e6]) xlabel n+1, title('Divergence as n+1 -> \infty',FS,9) %% % \vskip 1pt %% % By now the reader may have suspected that what is going wrong here can be % understood by looking at potentials, as in the last two chapters. Here is % an adaptation of a code segment from Chapter 11 to plot equipotential % curves for $n+1 = 8$ and $20$. %% % \vspace{-2em} d = domain(-1.5,1.5); xgrid = -1.5:.02:1.5; ygrid = -1:.02:1; [xx,yy] = meshgrid(xgrid,ygrid); zz = xx+1i*yy; for np = [8 20] xj = linspace(-1,1,np); ell = poly(xj,d); hold off, plot(xj,ell(xj),'.k','markersize',8) hold on, ylim([-1.2 1.2]), axis equal ellzz = ell(zz); levels = ((1.25:.25:3)/exp(1)).^np; contour(xx,yy,abs(ellzz),levels,'k') title(['Level curves of |l(x)| for '... int2str(np) ' equispaced points'],FS,9) snapnow end %% % \vspace{1pt} %% % What we see here is that $[-1,1]$ is very far from being a level curve % for equispaced interpolation points. From the last two chapters, we % expect serious consequences of this irregularity. In the second plot % just shown, for example, it is the fourth curve (from inside out) that % approximately touches the endpoints $\pm 1$. For Theorem 11.1 to be of % any use in such a landscape, $f$ will have to be analytic in a region % larger than the football'' enclosed by that curve. Analyticity on % $[-1,1]$ is not enough for convergence; we will need analyticity in a % much bigger region of the complex plane. This is what Runge discovered in % 1901. %% % Following the method of the last chapter, we now consider the limit % $n\to\infty$, where the distribution of interpolation points approaches % the constant measure (12.9), % $$\mu(\tau) = {1\over 2} . \eqno (13.1)$$ % The potential (12.7) associated with this measure is % $$u(s) = -1 + {1\over 2}\hbox{Re} % \left[(s+1)\log(s+1)-(s-1)\log(s-1)\right]. \eqno (13.2)$$ % Here is a code that plots just one level curve of this potential, the one % passing through $\pm 1$, where the value of the potential is $-1 + \log % 2$. %% % \vskip -2em x1 = -1.65:.02:1.65; y1 = -0.7:.02:0.7; [xx,yy] = meshgrid(x1,y1); ss = xx+1i*yy; u = @(s) -1 + 0.5*real((s+1).*log(s+1)-(s-1).*log(s-1)); hold off contour(x1,y1,u(ss),(-1+log(2))*[1 1],'k','linewidth',1.4) set(gca,'xtick',-2:.5:2,'ytick',-.5:.5:.5), grid on ylim([-.9 .9]), axis equal hold on, plot(.5255i,'.k') text(0.05,.63,'0.52552491457i') title(['Runge region for equispaced interpolation ' ... 'in the limit n -> \infty'],FS,9) %% % \vskip 1pt %% % For the interpolants to a function $f$ in equispaced nodes to converge as % $n\to\infty$ for all $x\in[-1,1]$, $f$ must be analytic, not just on % $[-1,1]$, but throughout this _Runge region,_ which crosses the real axis % at $\pm 1$ and the imaginary axis at $\pm 0.52552491457\dots i.$ Runge % reports this number correctly to 4 digits, and writes %% % % \em The curve has somewhat the shape of an ellipse. At the ends of the % long axis, however, our curve is more pointed than an % ellipse.\footnote{Die Kurve$\dots$ hat etwa die Gestalt einer % Ellipse$\dots.$ An den Enden der grossen Achse ist unsere Kurve aber % spitzer als eine Ellipse.''} %% % Here are the values of (13.2) at the endpoints and the middle: % $$u(\pm 1) = -1 + \log 2, \quad u(0) = -1,$$ % and thus % $$\exp( u(\pm 1))= {2\over e}, \quad \exp(u(0)) = {1\over e}.$$ % These numbers indicate that in the limit $n\to\infty$, the endpoints of % an equispaced grid in $[-1,1]$ lie at an average distance $2/e$ from the % other grid points, in the geometric mean sense, whereas the midpoint lies % at an average distance of just $1/e$. As emphasized in the last chapter, % notably in equation (12.6), the effect of such a discrepancy grows % exponentially with $n$. %% % Here are some examples. Equispaced interpolation will converge throughout % $[-1,1]$ for $f(x) = \exp(-x^2)$, which is analytic everywhere, and for % $f(x) = (1+x^2)^{-1}$, which has poles at $\pm i$. On the other hand it % will not converge for $f(x) = (1+16x^2)^{-1}$, which has poles at $\pm % i/4$. It will converge slowly for $f(x) = (1+(x/0.53)^2)^{-1}$, and % diverge slowly for $f(x) = (1+(x/0.52)^2)^{-1}$ (Exercise 13.3). %% % The worst-case rate of divergence is $2^n$, and this rate will always % appear if $f$ is not analytic on $[-1,1]$. To be precise, for such a % function the errors will be of size $O(2^n)$ as $n\to\infty$ but not of % size $O(C^n)$ for any $C<2$. Here, for example, we take $f$ to be a hat % function, with just one derivative of bounded variation. The dots show % errors in Chebyshev interpolation, converging at the rate $O(n^{-1})$ in % keeping with Theorem 7.2, and the crosses show errors in equispaced % interpolation, with a dashed line parallel to $2^n$ for comparison. %% % \vskip -2em f = max(0,1-2*abs(x)); eequi = []; echeb = []; nn = 2:2:60; for n = nn s = linspace(-1,1,n+1); pequi = interp1(s,f,domain(-1,1)); eequi = [eequi norm(f-pequi,inf)]; pcheb = chebfun(f,n+1); echeb = [echeb norm(f-pcheb,inf)]; end hold off, semilogy(nn,2.^(nn-12),'--r') hold on, axis([0 60 1e-4 1e14]), grid on semilogy(nn,eequi,'x-r','markersize',8), semilogy(nn,echeb,'.-b') text(47,3e6,'equispaced','color','r') text(41,0.8,'Chebyshev','color','b') text(32,4e8,'C 2^n','fontsize',12,'color','r') xlabel np, ylabel Error, title('Chebyshev vs. equispaced points',FS,9) %% % \vskip 1pt %% % All of the above remarks about equispaced interpolation concern the ideal % mathematics of the problem. On a computer in floating point arithmetic, % however, a further complication arises: even if convergence ought to take % place in theory, rounding errors will be amplified by $O(2^n)$, causing % divergence in practice. Here, for example, are the errors in equispaced % and Chebyshev interpolation of the entire function $\exp(x)$: %% % \vskip -2em f = exp(x); eequi = []; echeb = []; nn = 2:2:80; for n = nn s = linspace(-1,1,n+1); pequi = interp1(s,f,domain(-1,1)); eequi = [eequi norm(f-pequi,inf)]; pcheb = chebfun(f,n+1); echeb = [echeb norm(f-pcheb,inf)]; end hold off, semilogy(nn,2.^(nn-50),'--r') hold on, axis([0 80 1e-17 1e4]), grid on semilogy(nn,eequi,'x-r','markersize',8), semilogy(nn,echeb,'.-b') text(22,6e-6,'C 2^n','fontsize',12,'color','r') text(42,3e-7,'equispaced','color','r') text(51,6e-14,'Chebyshev','color','b') xlabel np, ylabel Error, title('The effect of rounding errors',FS,9) %% % \vskip 1pt %% % In exact arithmetic we would see convergence for both sets of points, but % on the computer the divergence for equispaced points sets in early. The % rate is cleanly $O(2^n)$ until we reach $O(1)$. Notice that the line of % crosses, if extended backward to $n=0$, meets the $y$ axis at % approximately $10^{-18}$, i.e., a digit or two below machine precision. %% % The $2^n$ divergence of equispaced polynomial interpolants is a % fascinating subject, and we must remind ourselves that one should only go % into so much detail in analyzing a method that should never be used! But % perhaps we should qualify that never'' in one respect. As Runge himself % emphasized, though interpolants in equispaced points do not converge on % the whole interval of interpolation, they may still do very well near the % middle. In the numerical solution of differential equations, for % example, higher order centered difference formulas are successfully used % based on 5 or 7 equally spaced grid points. A less happy example would % be Newton--Cotes quadrature formulas, based on polynomial interpolation % in equally spaced points, where the $O(2^n)$ effect is unavoidable and % causes serious problems for larger $n$ and divergence as $n\to\infty$, as % first proved by $\hbox{P\'olya}$ . We shall discuss quadrature in % Chapter 19. %% % We close this chapter with an observation that highlights the fundamental % nature of the Runge phenomenon and its associated mathematics. Suppose % you want to persuade somebody that it is important to know something about % complex variables, even for dealing with real functions. I still % remember the argument my calculus teacher explained to me: to understand % why the Taylor series for $1/(1+x^2)$ only converges for $-1 < x < 1$, % you need to know that Taylor series converge within disks in the complex % plane, and this function has poles at $\pm i$. %% % Runge's observation is precisely a generalization of this fact to % interpolation points equispaced in an interval rather than all at $x=0$. % The convergence or divergence of polynomial interpolants to a function % $f$ again depends on whether $f$ is analytic in a certain region; the % change is that the region is now not a disk, but elongated. Even the % phenomenon of divergence in floating-point arithmetic for functions whose % interpolants ought'' to converge is a generalization of familiar facts % from real arithmetic. Just try to evaluate $\exp(x)$ on a computer for % $x=-20$ using the Taylor series! %% % % \begin{displaymath} % \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl % {\sc Summary of Chapter 13.} % Polynomial interpolation in equispaced points is exponentially % ill-conditioned: the interpolant $p_n$ may have oscillations near the % edge of the interval nearly $\kern .7pt 2^n$ times larger than the function $f$ % being interpolated, even if $f$ is analytic. In particular, even if $f$ % is analytic and the interpolant is computed exactly without rounding % errors, $p_n$ need not converge to $f$ as $n\to\infty$.\vspace{2pt}}} % \end{displaymath} % %% % \smallskip \parskip=2pt\small % {\bf Exercise 13.1. Three examples.} % Draw plots comparing accuracy of equispaced % and Chebyshev interpolants as functions of $n$ for % $\exp(x^2)$, $\exp(-x^2)$, $\exp(-1/x^2)$. % \par % {\bf Exercise 13.2. Computing geometric means in Chebfun.} % (a) What output is produced by the program below? (b) Why? % \par\verb|x = chebfun('x',[0 1]);| % \hfill\break\verb|f = chebfun(@(y) prod(abs(x-1i*y)),[0.1 1],'vectorize');| % \hfill\break\verb|roots(f-2/exp(1))| % \par % {\bf Exercise 13.3. Borderline convergence.} % The claim was made in the text that equi\-spaced % polynomial interpolants on $[-1,1]$ converge % for $f(x) = (1+(x/0.53)^2)^{-1}$ and diverge for % $f(x) = (1+(x/0.52)^2)^{-1}$. Can you observe this difference % numerically? Run appropriate experiments and discuss the results. % \par % {\bf Exercise 13.4. Approaching the sinc function.} % The sinc function is defined for all $x$ by $S(x) = \sin(\pi x)/(\pi x)$ % (and $S(0) = 1$). % For any $n\ge 1$, an approximation to $S$ is given by % $$S_n = \prod_{k=1}^n (1-x^2/k^2).$$ % Construct $S_{20}$ in Chebfun on the interval $[-20,20]$ and % compare it with $S$. On what interval around $x=0$ do you % find $|S_{20}(x)-S(x)| < 0.1$? (It can be shown that for % every $x$, $\lim_{n\to \infty} S_n(x) = S(x)$.) % \par % {\bf Exercise 13.5. Barycentric weights and ill-conditioning.} % (a) Suppose a function is interpolated by a polynomial in $n+1$ % equispaced points in $[-1,1]$, with $n$ even. From the result of Exercise 5.6, % derive a formula for the ratio of the barycentric weights at % the midpoint $x=0$ and the endpoint $x=1$. % (b) With reference to the barycentric formula (5.11), explain what % this implies about sensitivity of these polynomial interpolants % to perturbations in the data at $x=0$. % \par