%% CHEBFUN GUIDE 4: CHEBFUN AND APPROXIMATION THEORY % Lloyd N. Trefethen, November 2009, revised February 2011 %% 4.1 Chebyshev series and interpolants % Chebfun is founded on the mathematical % subject of approximation theory, and in particular, on % Chebyshev series and interpolants. % Conversely, it provides a simple environment in % which to demonstrate these approximants and other approximation ideas. %% % The history of "Chebyshev technology" goes back to % the 19th century Russian mathematician P. L. Chebyshev (1821-1894) % and his mathematical descendants such as Zolotarev and % Bernstein (1880-1968). These men realized that % just as Fourier series provide an efficient way to represent % a smooth periodic function, series of Chebyshev polynomials can % do the same for a smooth nonperiodic function. A number % of excellent textbooks and monographs have been published on approximation % theory, including [Davis 1963], [Cheney 1966], [Meinardus 1967], [Lorentz 1986], % and Powell [Powell, 1981], % and in addition there are books devoted entirely to Chebyshev polynomials, % including [Rivlin 1974] and [Mason & Handscomb 2003]. % A Chebfun-based book on approximation theory and its % computational applications is in preparation [Trefethen 2012]. %% % From the dates of publication above it will be clear that % approximation theory flourished in the early computer era, and % in the 1950s and 1960s a number of numerical methods were developed based % on Chebyshev polynomials by Lanczos [Lanczos 1957], Fox [Fox & Parker 1966], % Clenshaw, Elliott, Mason, Good, and others. The Fast Fourier Transform % came in 1965 and Salzer's barycentric interpolation formula for % Chebyshev points in 1972 [Salzer 1972]. Then in the 1970s Orszag % and Gottlieb introduced spectral methods, based on the % application of Chebyshev and Fourier technology to the solution of % PDEs. The subject grew rapidly, and it is in the context of spectral % methods that Chebyshev techniques are best known today % [Boyd 2001], [Trefethen 2000], [Canuto et al. 2006/7]. %% % We must be clear about terminology. We shall rarely use the % term *Chebyshev approximation*, for that expression refers specifically to % an approximation that is optimal in the minimax sense. % Chebyshev approximations are fascinating, and in Section 4.6 we % shall see that Chebfun makes it easy to compute them, but % the core of Chebfun is built on the different techniques of polynomial % interpolation in Chebyshev points and expansion in Chebyshev % polynomials. These approximations are not quite optimal, but % they are nearly optimal and much easier to compute. %% % By *Chebyshev points* we shall mean as usual the set of points in % [-1,1] defined by % % x(j) = -cos(j pi/N), 0 <= j <= N, % % where N is an integer >= 1. (If N=0, we take x(0)=0.) % Through any data values f(j) at these points there is a unique % polynomial interpolant p(x) of degree <= N, which we call % the *Chebyshev interpolant*. % In particular, if the data are f(j) = (-1)^(n-j), then p(x) is % T_N, the Nth Chebyshev polynomial, which can also be defined % by the formula T_N(x) = cos(N acos(x)). In Chebfun, % the command "chebpoly(N)" returns a chebfun corresponding to T_N, and % "poly" returns coefficients in the monomial basis 1,x,x^2,.... % Thus we can print the coefficients of the first few Chebyshev % polynomials like this: for N = 0:8 disp(poly(chebpoly(N))) end %% % Note that that output of "poly" follows the pattern for Matlab's standard "poly" % command: it is a row vector, and the high-order coefficients come first. % Thus, for example, the fourth row above tells us that T_3(x) = 4x^3 - 3x. %% % Here are plots of T_2, T_3, T_15, and T_50. subplot(2,2,1), plot(chebpoly(2)) subplot(2,2,2), plot(chebpoly(3)) subplot(2,2,3), plot(chebpoly(15)) subplot(2,2,4), plot(chebpoly(50)) %% % A *Chebyshev series* is an expansion % % f(x) = SUM_{k=0}^infty a_k T_k(x), % % and the a_k are known as *Chebyshev coefficients*. So long % as f is continuous and at least a little bit smooth (Lipschitz % continuity is enough), it has a unique expansion of % this form, which converges absolutely and % uniformly, and the coefficients are given by the integral % % a_k = (2/pi) INT_{-1}^1 f(x) T_k(x) dx / sqrt(1-x^2). % % except that for k=0, the constant changes from 2/pi to 1/pi. % One way to approximate a function is to form % the polynomials obtained by truncating its Chebyshev expansion, % % f_N(x) = SUM_{k=0}^N a_k T_k(x). %% % This isn't quite what Chebfun does, however, since it % does not compute exact Chebyshev coefficients. Instead Chebfun is based % on Chebyshev interpolants, which can also be regarded as % finite series in Chebyshev polynomials for some coefficients c_k: % % p_N(x) = SUM_{k=0}^N c_k T_k(x). % % Each coefficient c_k will converge to a_k as N->infty (apart from the % effects of rounding errors), but for finite N, c_k and a_k are different. % The system actually stores a function by its values at the Chebyshev % points rather than its Chebyshev coefficients, but this hardly matters % to the user, and both representations are exploited for various % purposes internally in the system. %% 4.2 chebpoly and poly % % Throughout this section, we must be sure to use the "factory" setting splitting off %% % since our purpose is to explore global Chebyshev interpolants % without any division into subintervals. %% % We have just seen that the command chebpoly(N) returns a chebfun % corresponding to the Chebyshev polynomial T_N. Conversely, if f is a % chebfun, then chebpoly(f) is the vector of its Chebyshev coefficients. % For example, here are the Chebyshev coefficients of x^3: x = chebfun('x'); c = chebpoly(x.^3) %% % Like "poly", "chebpoly" returns a row vector with the high-order coefficients first. % Thus this computation reveals the identity x^3 = (1/4)*T_3(x) + (3/4)*T_1(x). %% % If we apply chebpoly to a function that is not % "really" a polynomial, we will usually get a vector whose first % entry (i.e., highest order) is just above machine precision. % This reflects the adaptive nature of the Chebfun % constructor, which always seeks to use a minimal number of points. chebpoly(sin(x)) %% % Of course, machine precision is defined relative to the scale of % the function: chebpoly(1e100*sin(x)) %% % By using "poly" we can print the coefficients of such a chebfun % in the monomial basis. Here for example are the coefficients % of the Chebyshev interpolant of exp(x) compared with the % Taylor series coefficients: cchebfun = flipud(poly(exp(x))'); ctaylor = 1./gamma(1:length(cchebfun))'; disp(' chebfun Taylor') disp([cchebfun ctaylor]) %% % The fact that these differ is not an indication of an % error in the Chebfun approximation. On the contrary, the Chebfun % coefficients do a better job of approximating than the truncated % Taylor series. If f were a function like 1/(1+25x^2), the % Taylor series would not converge at all. %% 4.3 chebfun(...,N) and the Gibbs phenomenon %% % We can examine the approximation qualities of Chebyshev interpolants % by means of a command of the form "chebfun(...,N)". When an integer % N is specified in this manner, it indicates that a Chebyshev interpolant % is to be constructed of precisely length N rather than by the usual % adaptive process. %% % Let us begin with a function that cannot be well approximated % by polynomials, the step function sign(x). To start with we % interpolate it in 10 or 20 points, taking N to be even % to avoid including a value 0 at the middle of the step. f = chebfun('sign(x)',10); subplot(1,2,1), plot(f,'.-'), grid on f = chebfun('sign(x)',20); subplot(1,2,2), plot(f,'.-'), grid on %% % There is an overshoot problem here, known as the Gibbs phenomenon, % that does not go away as N -> infty. % We can zoom in on the overshoot region by resetting the axes: subplot(1,2,1), axis([0 .4 .5 1.5]) subplot(1,2,2), axis([0 .2 .5 1.5]) %% % Here are analogous results with N=100 and 1000. f = chebfun('sign(x)',100); subplot(1,2,1), plot(f,'.-'), grid on, axis([0 .04 .5 1.5]) f = chebfun('sign(x)',1000); subplot(1,2,2), plot(f,'.-'), grid on, axis([0 .004 .5 1.5]) %% % The second plot is jagged, not because there is anything wrong with % the underlying chebfun but because we have zoomed in very closely on the result % of a "plot" command. One way to get it right is by using Chebfun's { } feature, % which enables one to construct a chebfun corresponding to a subinterval. % (Another way would be to use the 'interval' flag in the "plot" command.) clf subplot(1,2,2), plot(f,'.'), grid on, axis([0 .004 .5 1.5]) hold on, plot(f{0,.004}) %% % What is the amplitude of the Gibbs overshoot for Chebyshev % interpolation of a step function? We can find out by using "max": for N = 2.^(1:8) gibbs = max(chebfun('sign(x)',N)); fprintf('%5d %13.8f\n', N, gibbs) end %% % This gets a bit slow for larger N, but knowing that the maximum occurs % around x = 3/N, we can speed it up by repeating the { } trick above. for N = 2.^(4:12) f = chebfun('sign(x)',N); fprintf('%5d %13.8f\n', N, max(f{0,5/N})) end %% % The overshoot converges to a number 1.282283455775.... % [Helmberg & Wagner 1997]. %% 4.4 Smoothness and rate of convergence % The most basic principle in approximation theory is this: % the smoother the function, the faster the convergence as N -> infty. % What this means for Chebfun is that so long % as a function is twice continuously differentiable, it can usually % be approximated to machine precision for a workable value of N, even % without subdivision of the interval. %% % After the step function, a function with "one more derivative" of smoothness % would be the absolute value. Here if we interpolate in N points, the % errors decrease at the rate O(N^(-1)). For example: clf f10 = chebfun('abs(x)',10); subplot(1,2,1), plot(f10,'.-'), grid on f20 = chebfun('abs(x)',20); subplot(1,2,2), plot(f20,'.-'), grid on %% % Chebfun has no difficulty computing interpolants of much higher order: f100 = chebfun('abs(x)',100); subplot(1,2,1), plot(f100), grid on f1000 = chebfun('abs(x)',1000); subplot(1,2,2), plot(f1000), grid on %% % Such plots look good to the eye, but they do not achieve machine precision. % We can confirm this by using "splitting on" % to compute a true absolute value and then measuring some norms. fexact = chebfun('abs(x)','splitting','on'); err10 = norm(f10-fexact,inf) err100 = norm(f100-fexact,inf) err1000 = norm(f1000-fexact,inf) %% % Notice the clean linear decrease of the error as N increases. %% % If f is a bit smoother, polynomial approximation to machine precision % becomes practical: length(chebfun('abs(x).*x')) length(chebfun('abs(x).*x.^2')) length(chebfun('abs(x).*x.^3')) length(chebfun('abs(x).*x.^4')) %% % Of course, these particular functions would be easily % approximated by piecewise smooth chebfuns. %% % It is interesting to plot convergence as a function of N. Here % is an example from [Battles & Trefethen 2004] involving the % next function from the sequence above. s = 'abs(x).^5'; exact = chebfun(s,'splitting','off'); NN = 1:100; e = []; for N = NN e(N) = norm(chebfun(s,N)-exact); end clf subplot(1,2,1) loglog(e), ylim([1e-10 10]), title('loglog scale') hold on, loglog(NN.^(-5),'--r'), grid on text(6,4e-7,'N^{-5}','color','r','fontsize',16) subplot(1,2,2) semilogy(e), ylim([1e-10 10]), grid on, title('semilog scale') hold on, semilogy(NN.^(-5),'--r'), grid on %% % The figure reveals very clean convergence at the rate % N^(-5). According to Theorem 2 of the next section, this % happens because f has a fifth derivative of % bounded variation. %% % Here is an example of a smoother function, one that is in % fact analytic. According to Theorem 3 of the next section, if % f is analytic, its Chebyshev interpolants converge geometrically. % In this example we take f to be the Runge function, for which % interpolants in equally spaced points would not converge at all % (in fact they diverge exponentially -- see Section 4.7). %% s = '1./(1+25*x.^2)'; exact = chebfun(s); for N = NN e(N) = norm(chebfun(s,N)-exact); end clf, subplot(1,2,1) loglog(e), ylim([1e-10 10]), grid on, title('loglog scale') c = 1/5 + sqrt(1+1/25); hold on, loglog(c.^(-NN),'--r'), grid on subplot(1,2,2) semilogy(e), ylim([1e-10 10]), title('semilog scale') hold on, semilogy(c.^(-NN),'--r'), grid on text(45,1e-3,'C^{-N}','color','r','fontsize',16) %% % This time the convergence is equally clean but quite different in % nature. Now the straight line appears on the semilogy axes rather % than the loglog axes, revealing the geometric convergence. %% 4.5 Five theorems % % The mathematics of Chebfun can be captured in five % theorems about interpolants in Chebyshev points. The first three % can be found in [Battles & Trefethen 2004], and all will be % discussed in [Trefethen 2012]. Let f be a % continuous function on [-1,1], and let p denote its interpolant % in N Chebyshev points and p* its best degree N approximation with % respect to the maximum norm *||* *||*. %% % The first theorem asserts that Chebyshev interpolants are "near-best" % [Ehlich & Zeller 1966]. %% % *THEOREM 1.* *||* f - p *||* <= (2 + (2/pi)log(N)) *||* f - p* *||*. %% % This theorem implies that even if N is as large as 100,000, one can lose % no more than one digit by using p instead of p*. Whereas Chebfun % will readily compute such a p, it is unlikely that anybody has ever computed % a nontrivial p* for a value of N so large. %% % The next theorem asserts that if f is k times differentiable, roughly % speaking, then the Chebyshev interpolants converge at the algebraic rate % 1/N^k [Mastroianni & Szabados 1995]. %% % *THEOREM 2*. If Let f, f', ..., f^(k-1) be absolutely continuous for some % k >=1, and let f^(k) be a function of bounded variation. Then % *||* f - p *||* = O(N^(-k)) as N -> infty. %% % Smoother than this would be a C-infty function, i.e. infinitely differentiable, % and smoother still would be a function analytic on [-1,1], i.e., % one whose Taylor series at each point of [-1,1] converges at least in a small % neighborhood of that point. In such a case the convergence is geometric. % The essence of the following theorem is due to Bernstein in 1912, though % I do not know where an explicit statement first appeared in print. %% % *THEOREM 3*. If f is analytic and bounded in the % "Bernstein ellipse" of foci 1 and -1 with % semimajor and semiminor axis lengths summing to r, then % *||* f - p *||* = O(r^(-N)) as N -> infty. %% % More precisely, if abs(f(z)) <= M in the ellipse, then the bound on the right % can be taken as 4Mr^(-n)/(r-1). %% % The next theorem asserts that Chebyshev interpolants can be computed % by the barycentric formula [Salzer 1972]. The notation SUM" denotes the % sum from k=0 to k=N with both terms k=0 and k=N multiplied by 1/2. %% % *THEOREM 4*. p(x) = SUM" (-1)^k f(x_k)/(x-x_k) / SUM" (-1)^k/(x-x_k). %% % See [Berrut & Trefethen 2005] for information about barycentric interpolation. %% % The final theorem asserts that the barycentric formula has no difficulty with % rounding errors. Our "theorem" is really just a placeholder; see % [Higham 2004] for a precise statement and proof. Earlier work on this % subject appeared in [Rack & Reimer 1982]. %% % *THEOREM 5*. The barycentric formula of Theorem 4 is numerically stable. %% % This stability result may seem surprising when one notes that for % x close to x_k, the barycentric formula involves divisions by numbers that % are nearly zero. Nevertheless it is provably stable. If x is exactly equal % to some x_k, then one bypasses the formula and returns the exact value % p(x) = f(x_k). %% 4.6 Best approximations and the Remez algorithm % For practical computations, it is rarely worth the trouble % to compute a best (minimax) approximation rather than simply a % Chebyshev interpolant. Nevertheless best approximations are % a beautiful and well-established idea, and it is certainly interesting % to be able to compute them. Chebfun makes this % possible with the command "remez", named after Evgeny Remez, % who devised the standard algorithm for computing these approximations in 1934. % This capability is due to Ricardo Pachon; % see [Pachon & Trefethen 2009]. %% % For example, here is a function on the interval [0,4] together with % its best approximation by a polynomial of degree 20: f = chebfun('sqrt(abs(x-3))',[0,4],'splitting','on'); p = remez(f,20); clf, plot(f,'b',p,'r'), grid on %% % A plot of the error curve (p-f)(x) shows that it % equioscillates between 20+2 = 22 alternating % extreme values. Note that a second output argument from remez returns % the error as well as the polynomial. [p,err] = remez(f,20); plot(f-p,'m'), hold on plot([0 4],err*[1 1],'--k'), plot([0 4],-err*[1 1],'--k') %% % Let's add the error curve for the degree 20 (i.e. 21-point) % Chebyshev interpolant to the same plot: pinterp = chebfun(f,[0,4],21); plot(f-pinterp,'b') %% % Notice that although the best approximation has a smaller % maximum error, it is a worse approximation for almost all x. %% % Chebfun "remez" command can compute certain rational best approximants % too, though it is somewhat fragile. % If your function is smooth, a possibly more robust approach to computing % best approximations is Caratheodory-Fejer approximation, implemented % in the code "cf" due to Joris Van Deun [Van Deun & Trefethen 2011]. For example: f = chebfun('exp(x)'); [p,q] = cf(f,5,5); r = p./q; err = norm(f-r,inf); clf, plot(f-r,'c'), hold on plot([-1 1],err*[1 1],'--k'), plot([-1 1],-err*[1 1],'--k') ylim(2e-13*[-1 1]) %% % CF approximation often comes close to optimal for non-smooth % functions too, provided you specify a fourth argument to tell the % system on how fine a Chebyshev grid to sample: f = abs(x-.3); [p,q,r_handle,lam] = cf(f,5,5,300); clf, plot(f-p./q,'c'), hold on plot([-1 1],lam*[1 1],'--k'), plot([-1 1],-lam*[1 1],'--k') %% 4.7 The Runge phenomenon % Chebfun is based on polynomial interpolants in Chebyshev % points, not equispaced points. It has been known for over a century % that the latter choice is disastrous, even for interpolation of % smooth functions [Runge 1901]. One should never use equispaced polynomial interpolants % for practical work (unless you will only need the result near the % center of the interval of interpolation), but like best % approximations, they are certainly interesting. %% % In Chebfun, we can compute them with the overloaded % "interp1" command. For example, here is an analytic function and % its equispaced interpolant of degree 9: f = tanh(10*x); s = linspace(-1,1,10); p = interp1(s,f); hold off plot(f), hold on, plot(p,'r'), grid on, plot(s,p(s),'.r') %% % Perhaps this doesn't look too bad, but here is what happens % for degree 19. Note the vertical scale. s = linspace(-1,1,20); p = interp1(s,f); hold off plot(f), hold on, plot(p,'r'), grid on, plot(s,p(s),'.r') %% % Approximation experts will know that one of the tools used % in analyzing effects like this is known as the Lebesgue % function associated with a given set of interpolation points. % Chebfun has a command "lebesgue" for computing these functions. % The problem with interpolation in 20 equispaced points is reflected % in a Lebesgue function of size 10^4 -- note the semilog scale: clf, semilogy(lebesgue(s)) %% % For 40 points it is much worse: semilogy(lebesgue(linspace(-1,1,40))) %% % As the degree increases, polynomial interpolants in equispaced % points diverge exponentially, and no other method of approximation based % on equispaced data can % completely get around this problem [Platte, Trefethen and Kuijlaars 2011]. %% 4.8 Rational approximations % Chebfun contains four different programs, at present, for % computing rational approximants to a function f. We say that % a rational function is of type (m,n) if it can be written as % a quotient of one polynomial of degree at most m and another of degree % at most n. %% % To illustrate the possibilities, consider the function f = chebfun('tanh(pi*x/2) + x/20',[-10,10]); length(f) plot(f) %% % We can use the command "chebpade", developed by Ricardo Pachon, % to compute a Chebyshev-Pade approximant, defined by the % condition that the Chebyshev series of p/q should match that of f as % far as possible [Baker and Graves-Morris 1996]. (This is the so-called "Clenshaw-Lord" Chebyshev-Pade % approximation; if the flag 'maehly' is specified the code alternatively % computes the linearized variation known as the "Maehly" approximation.) % Chebyshev-Pade approximation is the analogue for functions defined on an interval % of Pade approximation for functions defined in a neighborhood of a point. [p,q] = chebpade(f,40,4); r = p./q; %% % The functions f and r match to about 8 digits: norm(f-r) plot(f-r,'r') %% % Mathematically, f has poles in the % complex plane at +-i, +-3i, +5i, and so on. % We can obtain approximations to these values by looking at the roots of q: roots(q,'complex') %% % A similar but perhaps faster and more robust approach to rational interpolation % is encoded in the command "ratinterp", which computes a type (m,n) interpolant % through m+n+1 Chebyshev points (or, optionally, a different set of points). % This capability was developed by Ricardo Pachon and Pedro Gonnet [Pachon % & Gonnet 2010]. The results are similar: [p,q] = ratinterp(f,40,4); r = p./q; norm(f-r) plot(f-r,'m') %% % Again the poles are not bad: roots(q,'complex') %% % The third and fourth options for rational approximation, as mentioned in % Section 4.6, are best approximants computed by "remez" and % Caratheodory-Fejer approximants computed by "cf" [Trefethen & Gutknecht 1983, % Van Deun & Trefethen 2011]. % As mentioned in Section 4.6, CF approximants often agree % with best approximations to machine precision if f is smooth. % We explore the same function yet again, and this time obtain % an equioscillating error curve: [p,q] = cf(f,40,4); r = p./q; norm(f-r) plot(f-r,'c') %% % And the poles: roots(q,'complex') %% % It is tempting to vary parameters and functions to explore more poles % and what accuracy can be obtained for them. But rational approximation % and analytic continuation are % very big subjects and we shall resist the temptation. %% 4.9 References % % % [Baker and Graves-Morris 1996] G. A. Baker, Jr. and P. Graves-Morris, % Pade Approximants, 2nd ed., Cambridge U. Press, 1996. % % [Battles & Trefethen 2004] Z. Battles and L. N. Trefethen, % "An extension of Matlab to continuous functions and % operators", SIAM Journal on Scientific Computing 25 (2004), % 1743-1770. % % [Berrut & Trefethen 2005] J.-P. Berrut and L. N. Trefethen, % "Barycentric Lagrange interpolation", SIAM Review 46 (2004), % 501-517. % % [Boyd 2001] J. P. Boyd, Chebyshev and Fourier Spectral Methods, % 2nd ed., Dover, 2001. % % [Canuto et al. 2006/7] C. Canuto, M. Y. Hussaini, A. Quarteroni % and T. A. Zang, Spectral Methods, 2 vols., Springer, 2006 and 2007. % % [Cheney 1966] E. W. Cheney, Introduction to Approximation % Theory, McGraw-Hill 1966 and AMS/Chelsea, 1999. % % [Davis 1963] P. J. Davis, Interpolation and Approximation, % Blaisdell, 1963 and Dover, 1975. % % [Ehlich & Zeller 1966] H. Ehlich and K. Zeller, "Auswertung der % Normen von Interpolationsoperatoren," Math. Annalen 164 (1966), 105-112. % % [Fox & Parker 1966] L. Fox and I. B. Parker, % Chebyshev Polynomials in Numerical Analysis, Oxford U. Press, 1968. % % [Helmberg & Wagner 1997] G. Helmberg & P. Wagner, % "Manipulating Gibbs' phenomenon for % Fourier interpolation," Journal of Approximation Theory 89 % (1997), 308-320. % % [Higham 2004] N. J. Higham, "The numerical stability % of barycentric Lagrange interpolation", IMA Journal of % Numerical Analysis 24 (2004), 547-556. % % [Lanczos 1956] C. Lanczos, Applied Analysis, Prentice-Hall, % 1956 and Dover, 1988. % % [Lorentz 1986] G. G. Lorentz, The Approximation of Functions, % American Mathematical Society, 1986. % % [Mason & Handscomb 2003] J. C. Mason and D. C. Handscomb, % Chebyshev Polynomials, CRC Press, 2003. % % [Mastroianni & Szabados 1995] G. Mastroianni and J. Szabados, % "Jackson order of approximation by Lagrange interpolation," % Acta. Math. Hungar. 69 (1995), 73-82. % % [Meinardus 1967] G. Meinardus, Approximation of Functions: % Theory and Numerical Methods, Springer, 1967. % % [Pachon & Gonnet 2010] R. Pachon and P. Gonnet, manuscript % in preparation. % % [Pachon & Trefethen 2009] R. Pachon and L. N. Trefethen, % "Barycentric-Remez algorithms for best polynomial approximation in the % chebfun system", BIT Numerical Mathematics 49 (2009), ?-?. % % [Platte, Trefethen & Kuijlaars 2011] R. P. Platte, L. N. Trefethen % and A. B. J. Kuijlaars, "Impossibility of fast stable approximation of analytic % functions from equispaced samples", SIAM Review, to appear. % % [Powell 1981] M. J. D. Powell, Approximation Theory and Methods, % Cambridge University Press, 1981. % % [Rack & Reimer 1982] H.-J. Rack and M. Reimer, "The numerical stability % of evaluation schemes for polynomials based on the Lagrange % interpolation form", BIT Numerical Mathematics 22 (1982), 101-107. % % [Rivlin 1974] T. J. Rivlin, The Chebyshev Polynomials, Wiley, 1974 and 1990. % % [Runge 1901] C. Runge, "Ueber empirische Funktionen und die Interpolation % zwischen aequidistanten Ordinaten", Zeitschrift fuer Mathematik % und Physik 46 (1901), 224-243. % % [Salzer 1972] H. E. Salzer, "Lagrangian interpolation at the % Chebyshev points cos(nu pi/n), nu = 0(1)n; some unnoted % advantages", Computer Journal 15 (1972),156-159. % % [Trefethen 2000] L. N. Trefethen, Spectral Methods in Matlab, SIAM, 2000. % % [Trefethen 2010] L. N. Trefethen, Approximation Theory and Approximation % Practice, book in preparation. % % [Trefethen & Gutknecht 1983] L. N. Trefethen and M. H. Gutknecht, "The % Caratheodory-Fejer method for real rational approximation", % SIAM Journal on Numerical Analysis 20 (1983), 420-436. % % [Van Deun & Trefethen 2011] J. Van Deun and L. N. Trefethen, % A robust implementation of the Caratheodory-Fejer method, BIT Numerical % Mathematics, submitted.