200 Chapter 5.Evaluation of Functions #define NFEW.. #define NMANY . float *c,*d,te,a,bi Economize NMANY power series coefficients e[0..NMANY-1]in the range (a,b)into NFEW coefficients d[o..NFEW-1]. c=vector(0,NMANY-1); d=vector(0,NFEW-1); e=vector(0,NMANY-1); pcshft((-2.0-b-a)/(b-a),(2.0-b-a)/(b-a),e,MANY); pccheb(e,c,NMANY); Here one would normally examine the Chebyshev coefficients c[0..NMANY-1]to decide how small NFEW can be. chebpc(c.d.NFEW): 鱼 18881892 pcshft(a,b,d,NFEW); 三%O州 100 In our example,by the way,the 8th through 10th Chebyshev coefficients turn out to be on the order of-7 x 107 3 x 107',and-9 x 107,so reasonable truncations (for single precision calculations)are somewhere in this range,yielding a polynomial with 8- from NUMERICAL RECIPES I 10 terms instead of the original 13. Replacing a 13-term polynomial with a (say)10-term polynomial without any loss of (Nort server accuracy-that does seem to be getting something for nothing.Is there some magic in this technique?Not really.The 13-term polynomial defined a function f(x).Equivalent to computer, America e University Press. 令 THE economizing the series,we could instead have evaluated f()at enough points to construct its Chebyshev approximation in the interval of interest,by the methods of $5.8.We would ART have obtained just the same lower-order polynomial.The principal lesson is that the rate of convergence of Chebyshev coefficients has nothing to do with the rate of convergence of 9 Programs power series coefficients;and it is the former that dictates the number of terms needed in a polynomial approximation.A function might have a divergent power series in some region of interest,but if the function itself is well-behaved,it will have perfectly good polynomial approximations.These can be found by the methods of $5.8,but not by economization of 、三马h馆 series.There is slightly less to economization of series than meets the eye. to dir OF SCIENTIFIC COMPUTING(ISBN CITED REFERENCES AND FURTHER READING: 1988-1992 Acton,F.S.1970,Numerica/Methods That Work;1990,corrected edition (Washington:Mathe- matical Association of America),Chapter 12. Arfken,G.1970,Mathematical Methods for Physicists,2nd ed.(New York:Academic Press), p.631.[1) Numerical Recipes 10-621 43108 (outside 5.12 Pade Approximants Software. A Pade approximant,so called,is that rational function (of a specified order)whose power series expansion agrees with a given power series to the highest possible order.If the rational function is M R(x)三 k=0 N (5.12.1)
200 Chapter 5. Evaluation of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). #define NFEW .. #define NMANY .. float *c,*d,*e,a,b; Economize NMANY power series coefficients e[0..NMANY-1] in the range (a, b) into NFEW coefficients d[0..NFEW-1]. c=vector(0,NMANY-1); d=vector(0,NFEW-1); e=vector(0,NMANY-1); pcshft((-2.0-b-a)/(b-a),(2.0-b-a)/(b-a),e,NMANY); pccheb(e,c,NMANY); ... Here one would normally examine the Chebyshev coefficients c[0..NMANY-1] to decide how small NFEW can be. chebpc(c,d,NFEW); pcshft(a,b,d,NFEW); In our example, by the way, the 8th through 10th Chebyshev coefficients turn out to be on the order of −7 × 10−6, 3 × 10−7, and −9 × 10−9, so reasonable truncations (for single precision calculations) are somewhere in this range, yielding a polynomial with 8 – 10 terms instead of the original 13. Replacing a 13-term polynomial with a (say) 10-term polynomial without any loss of accuracy — that does seem to be getting something for nothing. Is there some magic in this technique? Not really. The 13-term polynomial defined a function f(x). Equivalent to economizing the series, we could instead have evaluated f(x) at enough points to construct its Chebyshev approximation in the interval of interest, by the methods of §5.8. We would have obtained just the same lower-order polynomial. The principal lesson is that the rate of convergence of Chebyshev coefficients has nothing to do with the rate of convergence of power series coefficients; and it is the former that dictates the number of terms needed in a polynomial approximation. A function might have a divergent power series in some region of interest, but if the function itself is well-behaved, it will have perfectly good polynomial approximations. These can be found by the methods of §5.8, but not by economization of series. There is slightly less to economization of series than meets the eye. CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), Chapter 12. Arfken, G. 1970, Mathematical Methods for Physicists, 2nd ed. (New York: Academic Press), p. 631. [1] 5.12 Pade Approximants ´ A Pade approximant ´ , so called, is that rational function (of a specified order) whose power series expansion agrees with a given power series to the highest possible order. If the rational function is R(x) ≡ M k=0 akxk 1 + N k=1 bkxk (5.12.1)
5.12 Pade Approximants 201 then R(x)is said to be a Pade approximant to the series f(x)≡ckxk (5.12.2) =0 if R(0)=f(0) (5.12.3) and also d (=) k=1,2,.,M+N (5.12.4) 三0 x=0 Equations (5.12.3)and (5.12.4)furnish M+N+1 equations for the unknowns ao,...,aM 81 and b1,...,bN.The easiest way to see what these equations are is to equate (5.12.1)and (5.12.2),multiply both by the denominator of equation (5.12.1),and equate all powers of z that have either a's or b's in their coefficients.If we consider only the special case of a diagonal rational approximation,M=N (cf.$3.2),then we have ao co,with the remaining a's and b's satisfying 3 N bmew-mtt=-e+, k=1,..,W (5.12.5) RECIPES m三1 令 》∑bmCk-m=ak, k=1,,N (5.12.6) m=0 X Press. (note,in equation 5.12.1,that bo =1).To solve these,start with equations (5.12.5).which are a set of linear equations for all the unknown b's.Although the set is in the form of a Toeplitz matrix (compare equation 2.8.8),experience shows that the equations are frequently close to singular,so that one should not solve them by the methods of 82.8,but rather by full LU decomposition.Additionally,it is a good idea to refine the solution by iterative improvement (routine mprove in $2.5)[1]. SCIENTIFIC Once the b's are known,then equation (5.12.6)gives an explicit formula for the unknown a's,completing the solution. 可 Pade approximants are typically used when there is some unknown underlying function f().We suppose that you are able somehow to compute,perhaps by laborious analytic expansions,the values of f(x)and a few of its derivatives at x=0:f(0),f(0),f"(0), and so on.These are of course the first few coefficients in the power series expansion of f(r);but they are not necessarily getting small,and you have no idea where (or whether) the power series is convergent. Numerical 10-621 By contrast with techniques like Chebyshev approximation (85.8)or economization of power series (85.11)that only condense the information that you already know about a function,Pade approximants can give you genuinely new information about your function's uctio 4310 values.It is sometimes quite mysterious how well this can work.(Like other mysteries in Recipes mathematics,it relates to analyricity.)An example will illustrate. Imagine that,by extraordinary labors,you have ground out the first five terms in the power series expansion of an unknown function f(r), North )≈2+1z 1 49 7872+. 175 (5.12.7) (It is not really necessary that you know the coefficients in exact rational form-numerical values are just as good.We here write them as rationals to give you the impression that they derive from some side analytic calculation.)Equation (5.12.7)is plotted as the curve labeled "power series"in Figure 5.12.1.One sees that for4 it is dominated by its largest,quartic,term We now take the five coefficients in equation(5.12.7)and run them through the routine pade listed below.It returns five rational coefficients,three a's and two b's,for use in equation (5.12.1)with M=N =2.The curve in the figure labeled "Pade"plots the resulting rational function.Note that both solid curves derive from the same five original coefficient values
5.12 Pad´e Approximants 201 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). then R(x) is said to be a Pade approximant to the series ´ f(x) ≡ ∞ k=0 ckxk (5.12.2) if R(0) = f(0) (5.12.3) and also dk dxk R(x) x=0 = dk dxk f(x) x=0 , k = 1, 2,...,M + N (5.12.4) Equations (5.12.3) and (5.12.4) furnish M + N + 1 equations for the unknowns a0,...,aM and b1,...,bN . The easiest way to see what these equations are is to equate (5.12.1) and (5.12.2), multiply both by the denominator of equation (5.12.1), and equate all powers of x that have either a’s or b’s in their coefficients. If we consider only the special case of a diagonal rational approximation, M = N (cf. §3.2), then we have a0 = c0, with the remaining a’s and b’s satisfying N m=1 bmcN−m+k = −cN+k, k = 1,...,N (5.12.5) k m=0 bmck−m = ak, k = 1,...,N (5.12.6) (note, in equation 5.12.1, that b0 = 1). To solve these, start with equations (5.12.5), which are a set of linear equations for all the unknown b’s. Although the set is in the form of a Toeplitz matrix (compare equation 2.8.8), experience shows that the equations are frequently close to singular, so that one should not solve them by the methods of §2.8, but rather by full LU decomposition. Additionally, it is a good idea to refine the solution by iterative improvement (routine mprove in §2.5) [1]. Once the b’s are known, then equation (5.12.6) gives an explicit formula for the unknown a’s, completing the solution. Pade approximants are typically used when there is some unknown underlying function ´ f(x). We suppose that you are able somehow to compute, perhaps by laborious analytic expansions, the values of f(x) and a few of its derivatives at x = 0: f(0), f (0), f(0), and so on. These are of course the first few coefficients in the power series expansion of f(x); but they are not necessarily getting small, and you have no idea where (or whether) the power series is convergent. By contrast with techniques like Chebyshev approximation (§5.8) or economization of power series (§5.11) that only condense the information that you already know about a function, Pade approximants can give you genuinely new information about your function’s ´ values. It is sometimes quite mysterious how well this can work. (Like other mysteries in mathematics, it relates to analyticity.) An example will illustrate. Imagine that, by extraordinary labors, you have ground out the first five terms in the power series expansion of an unknown function f(x), f(x) ≈ 2 + 1 9 x + 1 81x2 − 49 8748 x3 + 175 78732 x4 + ··· (5.12.7) (It is not really necessary that you know the coefficients in exact rational form — numerical values are just as good. We here write them as rationals to give you the impression that they derive from some side analytic calculation.) Equation (5.12.7) is plotted as the curve labeled “power series” in Figure 5.12.1. One sees that for x >∼ 4 it is dominated by its largest, quartic, term. We now take the five coefficients in equation (5.12.7) and run them through the routine pade listed below. It returns five rational coefficients, three a’s and two b’s, for use in equation (5.12.1) with M = N = 2. The curve in the figure labeled “Pade” plots the resulting rational ´ function. Note that both solid curves derive from the same five original coefficient values.
202 Chapter 5. Evaluation of Functions 10 fx)=[7+(1+x)43]3 6 power series(5 terms) 83g Pade (5 coefficients) granted for 18881992 1600 (including this one) from NUMERICAL RECIPES IN exact /Cambridge (North 0 2 4 6 8 10 America computer e University Press. THE ART Figure 5.12.1.The five-term power series expansion and the derived five-coefficient Pad approximant for a sample function f().The full power series converges only for x #include "nrutil.h" #define BIG 1.0e30 void pade(double cof,int n,float *resid) Given cof [0..2*n],the leading terms in the power series expansion of a function,solve the linear Pade equations to return the coefficients of a diagonal rational function approximation to the same function,namely (cof [o]cof [1]x+...+cof [n]zN)/(1+cof [n+1]x+...+
202 Chapter 5. Evaluation of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 0 2 4 6 8 10 0 2 4 6 8 10 f(x) x Padé (5 coefficients) exact power series (5 terms) f(x) = [7 + (1 + x)4/3]1/3 Figure 5.12.1. The five-term power series expansion and the derived five-coefficient Pade approximant ´ for a sample function f(x). The full power series converges only for x #include "nrutil.h" #define BIG 1.0e30 void pade(double cof[], int n, float *resid) Given cof[0..2*n], the leading terms in the power series expansion of a function, solve the linear Pad´e equations to return the coefficients of a diagonal rational function approximation to the same function, namely (cof[0] + cof[1]x + ··· + cof[n]xN )/(1 + cof[n+1]x + ··· +
5.12 Pade Approximants 203 cof [2*n]).The value resid is the norm of the residual vector;a small value indicates a well-converged solution.Note that cof is double precision for consistency with ratval. void lubksb(float **a,int n,int *indx,float b[]); void ludcmp(float **a,int n,int *indx,float *d); void mprove(float **a,float **alud,int n,int indx[],float b[], f1oatx☐)； int j,k,*indx; float d,rr,rrold,sum,**q,**qlu,*x,*y,*Z; indx=ivector(1,n); q=matrix(1,n,1,n) http://www.nr read able files Permission is qlu=matrix(1,n,1,n); x=vector(1,n); y=vector(1,n); 839 z-vector(1,n); (including this one) granted for 19881992 for (j=1;j<=n;j++){ Set up matrix for solving. y[j]=x[j]=cof[n+j]; 11-600 for(k=1;k<=n;k++){ q[j][k]=cof [j-k+n]; qlu[j][k]=q[j][k] n NUMERICAL RECIPES IN ludcmp(qlu,n,indx,&d); Solve by LU decomposition and backsubstitu- lubksb(qlu,n,indx,x); (North tion. rr-BIG: America to any server computer. Cambridge University Press. do Important to use iterative improvement,since tusers to make one paper THE rrold-rr; the Pade equations tend to be ill-conditioned. ART for (j=1;j<=m;j++)z[j]=x[j]; mprove(q,qlu,n,indx,y,x); for(rr=0.0,j=1;j<=n;j++) Calculate residual. rr +SQR(z[j]-x[j]); while (rr rrold); If it is no longer improving,call it quits. *resid=sqrt(rrold); 1CIYP ic Programs for (k=1;k<=n;k++){ Calculate the remaining coefficients. for (sum=cof [k],j=1;j<=k;j++) sum -z[j]*cof[k-j]; to dir y[k]=sum; 2 Copy answers to output. for(j=1;j<=n;j++){ OF SCIENTIFIC COMPUTING(ISBN cof [j]=y[j]; cof [i+n]=-z[i]i free_vector(z,1,n); free_vector(y,1,n); free_vector(x,1,n); v@cambridge.org free_matrix(qlu,1,n,1,n); 1988-1992 by Numerical Recipes 10-:6211 43108 free_matrix(q,1,n,1,n); free_ivector(indx,1,n); (outside North Software. CITED REFERENCES AND FURTHER READING: Ralston,A.and Wilf,H.S.1960,Mathematical Methods for Digital Computers (New York:Wiley), visit website p.14. machine Cuyt,A.,and Wuytack,L.1987,Nonlinear Methods in Numerical Analysis (Amsterdam:North- Holland).Chapter 2. Graves-Morris,P.R.1979,in Pade Approximation and Its Applications,Lecture Notes in Mathe- matics,vol.765,L.Wuytack,ed.(Berlin:Springer-Verlag).[1]
5.12 Pad´e Approximants 203 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). cof[2*n]xN ). The value resid is the norm of the residual vector; a small value indicates a well-converged solution. Note that cof is double precision for consistency with ratval. { void lubksb(float **a, int n, int *indx, float b[]); void ludcmp(float **a, int n, int *indx, float *d); void mprove(float **a, float **alud, int n, int indx[], float b[], float x[]); int j,k,*indx; float d,rr,rrold,sum,**q,**qlu,*x,*y,*z; indx=ivector(1,n); q=matrix(1,n,1,n); qlu=matrix(1,n,1,n); x=vector(1,n); y=vector(1,n); z=vector(1,n); for (j=1;j<=n;j++) { Set up matrix for solving. y[j]=x[j]=cof[n+j]; for (k=1;k<=n;k++) { q[j][k]=cof[j-k+n]; qlu[j][k]=q[j][k]; } } ludcmp(qlu,n,indx,&d); Solve by LU decomposition and backsubstitulubksb(qlu,n,indx,x); tion. rr=BIG; do { Important to use iterative improvement, since rrold=rr; the Pad´e equations tend to be ill-conditioned. for (j=1;j<=n;j++) z[j]=x[j]; mprove(q,qlu,n,indx,y,x); for (rr=0.0,j=1;j<=n;j++) Calculate residual. rr += SQR(z[j]-x[j]); } while (rr < rrold); If it is no longer improving, call it quits. *resid=sqrt(rrold); for (k=1;k<=n;k++) { Calculate the remaining coefficients. for (sum=cof[k],j=1;j<=k;j++) sum -= z[j]*cof[k-j]; y[k]=sum; } Copy answers to output. for (j=1;j<=n;j++) { cof[j]=y[j]; cof[j+n] = -z[j]; } free_vector(z,1,n); free_vector(y,1,n); free_vector(x,1,n); free_matrix(qlu,1,n,1,n); free_matrix(q,1,n,1,n); free_ivector(indx,1,n); } CITED REFERENCES AND FURTHER READING: Ralston, A. and Wilf, H.S. 1960, Mathematical Methods for Digital Computers (New York: Wiley), p. 14. Cuyt, A., and Wuytack, L. 1987, Nonlinear Methods in Numerical Analysis (Amsterdam: NorthHolland), Chapter 2. Graves-Morris, P.R. 1979, in Pad´e Approximation and Its Applications, Lecture Notes in Mathematics, vol. 765, L. Wuytack, ed. (Berlin: Springer-Verlag). [1]
204 Chapter 5.Evaluation of Functions 5.13 Rational Chebyshev Approximation In 85.8 and $5.10 we learned how to find good polynomial approximations to a given function f(r)in a given interval a sz<b.Here,we want to generalize the task to find good approximations that are rational functions(see $5.3).The reason for doing so is that, for some functions and some intervals,the optimal rational function approximation is able to achieve substantially higher accuracy than the optimal polynomial approximation with the same number of coefficients.This must be weighed against the fact that finding a rational function approximation is not as straightforward as finding a polynomial approximation, which,as we saw,could be done elegantly via Chebyshev polynomials. Let the desired rational function R(x)have numerator of degree m and denominator 8 of degree k.Then we have R(z)=Popi++pmzm 1+9hx+…+9xk≈f() fora≤x≤b (5.13.1) 区 The unknown quantities that we need to find are po,...pm and g,.q&,that is,m++1 quantities in all.Let r(x)denote the deviation of R(x)from f(x),and let r denote its maximum absolute value, RECIPES ra=R(o)-f)r=ar川 (5.13.2) 9 The ideal minimax solution would be that choice of p's and g's that minimizes r.Obviously Press. there is some minimax solution,since r is bounded below by zero.How can we find it,or a reasonable approximation to it? A first hint is furnished by the following fundamental theorem:If R(x)is nondegenerate 9 (has no common polynomial factors in numerator and denominator),then there is a unique choice of p's and g's that minimizes r;for this choice,r(r)has m+k+2 extrema in a <<b,all of magnitude r and with alternating sign.(We have omitted some technical SCIENTIFIC assumptions in this theorem.See Ralston[1]for a precise statement.)We thus learn that the situation with rational functions is quite analogous to that for minimax polynomials:In 85.8 we saw that the error term of an nth order approximation,with n+1 Chebyshev coefficients, was generally dominated by the first neglected Chebyshev term,namely T,which itself has n+2 extrema of equal magnitude and alternating sign.So,here,the number of rational coefficients,m+k+1,plays the same role of the number of polynomial coefficients,n+1. A different way to see why r()should have m+k:+2 extrema is to note that R(x) can be made exactly equal to f(x)at any m+k+1 points xi.Multiplying equation (5.13.1) by its denominator gives the equations Numerica 10.621 p0+p1x+·+pmx"=f(x)1+1x+…+qkx) 43126 (5.13.3) Recipes i=1,2,..,m+k+1 (outside This is a set of m+k+1 linear equations for the unknown p's and g's,which can be Software. solved by standard methods (e.g.,LU decomposition).If we choose the i's to all be in the interval (a,b),then there will generically be an extremum between each chosen zi and i+1,plus also extrema where the function goes out of the interval at a and b,for a total of m+k +2 extrema.For arbitrary xi's,the extrema will not have the same magnitude. The theorem says that,for one particular choice of i's,the magnitudes can be beaten down to the identical,minimal,value of r. Instead of making f(i)and R()equal at the points one can instead force the residual r(i)to any desired values y by solving the linear equations p0+p1xi+·十pmx"=[U(x)-(1+q1x+·+qkx) (5.13.4) i=1,2,..,m+k+1
204 Chapter 5. Evaluation of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 5.13 Rational Chebyshev Approximation In §5.8 and §5.10 we learned how to find good polynomial approximations to a given function f(x) in a given interval a ≤ x ≤ b. Here, we want to generalize the task to find good approximations that are rational functions (see §5.3). The reason for doing so is that, for some functions and some intervals, the optimal rational function approximation is able to achieve substantially higher accuracy than the optimal polynomial approximation with the same number of coefficients. This must be weighed against the fact that finding a rational function approximation is not as straightforward as finding a polynomial approximation, which, as we saw, could be done elegantly via Chebyshev polynomials. Let the desired rational function R(x) have numerator of degree m and denominator of degree k. Then we have R(x) ≡ p0 + p1x + ··· + pmxm 1 + q1x + ··· + qkxk ≈ f(x) for a ≤ x ≤ b (5.13.1) The unknown quantities that we need to find are p0,...,pm and q1,...,qk, that is, m +k + 1 quantities in all. Let r(x) denote the deviation of R(x) from f(x), and let r denote its maximum absolute value, r(x) ≡ R(x) − f(x) r ≡ max a≤x≤b |r(x)| (5.13.2) The ideal minimax solution would be that choice of p’s and q’s that minimizes r. Obviously there is some minimax solution, since r is bounded below by zero. How can we find it, or a reasonable approximation to it? A first hint is furnished by the following fundamental theorem: If R(x)is nondegenerate (has no common polynomial factors in numerator and denominator), then there is a unique choice of p’s and q’s that minimizes r; for this choice, r(x) has m + k + 2 extrema in a ≤ x ≤ b, all of magnitude r and with alternating sign. (We have omitted some technical assumptions in this theorem. See Ralston [1] for a precise statement.) We thus learn that the situation with rational functions is quite analogous to that for minimax polynomials: In §5.8 we saw that the error term of an nth order approximation, with n + 1 Chebyshev coefficients, was generally dominated by the first neglected Chebyshev term, namely Tn+1, which itself has n + 2 extrema of equal magnitude and alternating sign. So, here, the number of rational coefficients, m + k + 1, plays the same role of the number of polynomial coefficients, n + 1. A different way to see why r(x) should have m + k + 2 extrema is to note that R(x) can be made exactly equal to f(x) at any m + k + 1 points xi. Multiplying equation (5.13.1) by its denominator gives the equations p0 + p1xi + ··· + pmxm i = f(xi)(1 + q1xi + ··· + qkxk i ) i = 1, 2,...,m + k + 1 (5.13.3) This is a set of m + k + 1 linear equations for the unknown p’s and q’s, which can be solved by standard methods (e.g., LU decomposition). If we choose the xi’s to all be in the interval (a, b), then there will generically be an extremum between each chosen xi and xi+1, plus also extrema where the function goes out of the interval at a and b, for a total of m + k + 2 extrema. For arbitrary xi’s, the extrema will not have the same magnitude. The theorem says that, for one particular choice of xi’s, the magnitudes can be beaten down to the identical, minimal, value of r. Instead of making f(xi) and R(xi) equal at the points xi, one can instead force the residual r(xi) to any desired values yi by solving the linear equations p0 + p1xi + ··· + pmxm i = [f(xi) − yi](1 + q1xi + ··· + qkxk i ) i = 1, 2,...,m + k + 1 (5.13.4)