f5-13

5.13RationalChebyshevApproximation 197Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).y(j)=x(j) do11k=1,n q(j,k)=cof(j-k+n+1) qlu(j,k)=q(j,k) enddo 11 enddo 12 call ludcmp(qlu,n,NMAX,indx,d) Solveby LUdecomposition andbacksubstitution. call lubksb(qlu,n,NMAX,indx,x)rr=BIG 1 continue Important to use iterative improvement, since the Pad´e equations tend to be ill-conditioned. rrold=rr do 13j=1,n z(j)=x(j) enddo 13 call mprove(q,qlu,n,NMAX,indx,y,x) rr=0. do14j=1,n Calculate residual. rr=rr+(z(j)-x(j))**2 enddo 14 if(rr.lt.rrold)goto 1 If it is no longer improving, call it quits. resid=sqrt(rrold) do16k=1,n Calculate the remaining coefficients. sum=cof(k+1) do15j=1,k sum=sum-z(j)*cof(k-j+1) enddo 15 y(k)=sum enddo 16 Copy answers to output. do17j=1,n cof(j+1)=y(j) cof(j+n+1)=-z(j) enddo 17 returnEND CITED REFERENCES AND FURTHER READING: Ralston,A.andWilf,H.S.1960, MathematicalMethodsforDigitalComputers (NewYork:Wiley), p. 14. Cuyt, A., and Wuytack, L. 1987, Nonlinear Methods in Numerical Analysis (Amsterdam: North- Holland), Chapter 2. Graves-Morris, P.R. 1979, in Pad´e Approximation and Its Applications , Lecture Notes in Mathe- matics, vol. 765, L. Wuytack, ed. (Berlin: Springer-Verlag). [1] 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 ableto achieve substantially higher accuracy than the optimal polynomial approximation with thesame number of coefficients. This must be weighed against the fact that finding a rationalfunction approximation is not as straightforward as finding a polynomial approximation,which, as we saw, could be done elegantly via Chebyshev polynomials. 198 Chapter5. EvaluationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Let the desired rational function R(x)have numerator of degree mand denominator of degree k. Then we have R(x)≡p0+p1x+···+pmxm 1+q1x+···+qkxk≈f(x)fora≤x≤b (5.13.1 ) Theunknown quantitiesthatweneedtofindare p0,...,p mandq1,...,q k,thatis, m+k+1 quantities in all. Let r(x)denote the deviation of R(x)from f(x), and let rdenote its maximum absolute value, r(x)≡R(x)−f(x) r≡max a≤x≤b|r(x)| (5.13.2 ) The ideal minimaxsolution would be that choice of p’s and q’s that minimizes r. Obviously there issomeminimax solution, since ris bounded below by zero. How can we find it, or a reasonable approximation to it? Afirsthintisfurnishedbythefollowingfundamentaltheorem: If R(x)isnondegenerate (has no common polynomial factors in numerator and denominator), then there is a uniquechoice of p’s and q’s that minimizes r; for this choice, r(x)hasm+k+2extrema in a≤x≤b,all of magnitude rand 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 wesawthatthe errortermofan nthorder approximation, with n+1Chebyshev coefficients, was generally dominated by the first neglected Chebyshev term, namely Tn+1, which itself hasn+2extrema of equal magnitude and alternating sign. So, here, the number of rational coefficients, m+k+1,plays thesame role ofthenumber ofpolynomial coefficients, n+1. A different way to see why r(x)should have m+k+2extrema is to note that R(x) can bemade exactlyequal to f(x)atany m+k+1points 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+1linear equations for the unknown p’s and q’s, which can be solved by standard methods (e.g., LUdecomposition). If we choose the xi’s to all be in the interval (a, b), then there will generically be an extremum between each chosen xiand xi+1, plus also extrema where the function goes out of the interval at aandb, for a total ofm+k+2extrema. 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)andR(xi)equal at the points xi, one can instead force the residual r(xi)to any desired values yiby solving the linear equations p0+p1xi+···+pmxm i=[f(xi)−yi](1 + q1xi+···+qkxk i) i=1,2,...,m +k+1(5.13.4 ) In fact, if the xi’s are chosen to be the extrema (not the zeros) of the minimax solution, then the equations satisfied will be p0+p1xi+···+pmxm i=[f(xi)±r](1 + q1xi+···+qkxk i) i=1,2,...,m +k+2(5.13.5 ) wherethe ±alternates forthealternatingextrema. Noticethatequation (5.13.5)issatisfiedat m+k+2extrema,whileequation (5.13.4) was satisfiedonly at m+k+1arbitrarypoints. How can this be? The answer is that rin equation (5.13.5) is an additional unknown, so that the number of both equations and unknowns is m+k+2. True, the set is mildly nonlinear (inr),but ingeneral itis stillperfectly soluble by methods that we willdevelop in Chapter 9. We thus see that, given only the locations of the extrema of the minimax rational function, we can solve for its coefficients and maximum deviation. Additional theorems, 5.13RationalChebyshevApproximation 199Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).R(x) − f(x)2 × 10−6 10−6 0 −1 × 10−6 −2 × 10−6 0 .5 1 1.5 2 2.5 3 xm = k = 4 f(x) = cos(x)/(1 + ex) 0 < x < π Figure 5.13.1. Solid curves show deviations r(x)forfive successive iterations of the routine ratlsq for an arbitrary test problem. The algorithm does not converge to exactly the minimax solution (shown as the dotted curve). But, after one iteration, the discrepancy is a small fraction of the last signi ficant bit of accuracy. leading up to the so-called Remes algorithms [1], tell how to converge to these locations by an iterative process. For example, here is a (slightly simpli fied) statement of Remes’ Second Algorithm : (1) Find an initial rational function with m+k+2extrema xi(not having equal deviation). (2) Solve equation (5.13.5) for new rational coef ficients and r. (3) Evaluate the resulting R(x)tofind its actual extrema (which will not be the same as the guessed values). (4) Replace each guessed value with the nearest actual extremum of the same sign. (5) Goback to step 2 and iterate to convergence. Under a broad set of assumptions, this method willconverge. Ralston [1]fillsinthenecessary details,including how to findtheinitialsetof xi’s. Up to this point, our discussion has been textbook-standard. We now reveal ourselves as heretics. We don ’t much like the elegant Remes algorithm. Its two nested iterations (on rin the nonlinear set 5.13.5, and on the new sets of xi’s) arefinicky and require a lot of special logic for degenerate cases. Even more heretical, we doubt that compulsive searchingfor theexactly best , equal deviation, approximation is worth the effort —except perhaps for those few people in the world whose business it is to find optimal approximations that get built into compilers and microchips. Whenweuserationalfunctionapproximation, thegoalisusuallymuchmorepragmatic: Inside some inner loop we are evaluating some function a zillion times, and we want tospeed up its evaluation. Almost never do we need this function to the last bit of machineaccuracy. Suppose (heresy!) we use an approximation whose error has m+k+2extrema whose deviations differ by a factor of 2. The theorems on which the Remes algorithmsare based guarantee that the perfect minimax solution will have extrema somewhere withinthis factor of 2 range –forcing down the higher extrema will cause the lower ones to rise, until all are equal. So our “sloppy”approximation is in fact within a fraction of a least significant bit of the minimax one. That is good enough for us, especially when we have available a very robust method forfinding the so-called “sloppy”approximation. Such a method is the least-squares solution of overdetermined linear equations by singular value decomposition ( §2.6 and §15.4). We proceed as follows: First, solve (in the least-squares sense) equation (5.13.3), not just form+k+1values of x i, but for a signi ficantly larger number of xi’s, spaced approximately 200 Chapter5. EvaluationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).like the zeros of a high-order Chebyshev polynomial. This gives an initial guess for R(x). Second, tabulate the resulting deviations, find the mean absolute deviation, call it r, and then solve (again in the least-squares sense) equation (5.13.5) with rfixed and the ±chosen to be the sign of the observed deviation ateach point xi. Third,repeat the second step a few times. You can spot some Remes orthodoxy lurking in our algorithm: The equations we solve are trying to bring the deviations not to zero, but rather to plus-or-minus some consistentvalue. However, we dispense with keeping track of actual extrema; and we solve only linearequations at each stage. One additional trick is to solve a weighted least-squares problem, where the weights are chosen to beat down the largest deviations fastest. Here is a program implementing these ideas. Notice that the only calls to the function fnoccurintheinitial fillingofthetable fs. Youcouldeasilymodifythecodetodothis filling outsideoftheroutine. Itisnotevennecessarythatyourabscissas xsbeexactlytheonesthatwe use, though the quality of the fit will deteriorate if you do not have several abscissas between each extremum of the (underlying) minimax solution. Notice that the rational coef ficients are output in a format suitable for evaluation by the routine ratvalin§5.3. SUBROUTINE ratlsq(fn,a,b,mm,kk,cof,dev) INTEGER kk,mm,NPFAC,MAXC,MAXP,MAXITDOUBLE PRECISION a,b,dev,cof(mm+kk+1),fn,PIO2,BIG PARAMETER (NPFAC=8,MAXC=20,MAXP=NPFAC*MAXC+1, * MAXIT=5,PIO2=3.141592653589793D0/2.D0,BIG=1.D30) EXTERNAL fn C USES fn,ratval,dsvbksb,dsvdcmp DOUBLE PRECISION versionsof svdcmp,svbksb. Returns in cof(1:mm+kk+1) the coefficients of a rational function approximation to the function fnintheinterval (a,b). Inputquantities mmandkkspecifytheorderofthenumer- atoranddenominator, respectively. Themaximumabsolutedeviationoftheapproximation (insofar as is known) is returned as dev. INTEGER i,it,j,ncof,nptDOUBLE PRECISION devmax,e,hth,pow,sum,bb(MAXP),coff(MAXC),ee(MAXP), * fs(MAXP),u(MAXP,MAXC),v(MAXC,MAXC),w(MAXC),wt(MAXP),xs(MAXP), * ratval ncof=mm+kk+1npt=NPFAC*ncof Number of points where function is evaluated, i.e., fineness of the mesh. dev=BIG do 11i=1,npt Fill arrays with mesh abscissas and function values. if (i.lt.npt/2) then Ateachend,useformulathatminimizesroundoffsensitivity. hth=PIO2*(i-1)/(npt-1.d0) xs(i)=a+(b-a)*sin(hth)**2 else hth=PIO2*(npt-i)/(npt-1.d0) xs(i)=b-(b-a)*sin(hth)**2 endiffs(i)=fn(xs(i))wt(i)=1.d0 Inlateriterations wewilladjusttheseweights to combatthe largest deviations. ee(i)=1.d0 enddo 11 e=0.d0 do17it=1,MAXIT Loop over iterations. do14i=1,npt Set up the “design matrix” for the least-squares fit. pow=wt(i)bb(i)=pow*(fs(i)+sign(e,ee(i))) Key idea here: Fit to fn(x)+ ewhere thedeviationispositive,to fn(x)−e where it is negative. Then eis sup- posed to become an approximation to the equal-ripple deviation.do 12j=1,mm+1 u(i,j)=powpow=pow*xs(i) enddo 12 pow=-bb(i) do13j=mm+2,ncof pow=pow*xs(i) u(i,j)=pow enddo 13 enddo 14 call dsvdcmp(u,npt,ncof,MAXP,MAXC,w,v) Singular Value Decomposition. 5.14EvaluationofFunctionsbyPathIntegration 201Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Inespeciallysingularordifficultcases,onemighthereeditthesingularvalues w(1:ncof), replacing small values by zero. call dsvbksb(u,w,v,npt,ncof,MAXP,MAXC,bb,coff) devmax=0.d0sum=0.d0 do 15j=1,npt Tabulate the deviations and revise the weights. ee(j)=ratval(xs(j),coff,mm,kk)-fs(j)wt(j)=abs(ee(j)) Use weighting to emphasize most deviant points. sum=sum+wt(j) if(wt(j).gt.devmax)devmax=wt(j) enddo 15 e=sum/npt Update eto be the mean absolute deviation. if (devmax.le.dev) then Save only the best coefficient set found. do16j=1,ncof cof(j)=coff(j) enddo 16 dev=devmax endifwrite (*,10) it,devmax enddo 17 return 10 FORMAT (1x,’ratlsq iteration=’,i2,’ max error=’,1pe10.3) END Figure 5.13.1 shows the discrepancies for the firstfive iterations of ratlsqwhen it is applied to find the m=k=4rationalfit to the function f(x)=c o s x/(1 + ex)in the interval (0,π). One sees that after the first iteration, the results are virtually as good as the minimax solution. The iterations do not converge in the order that the figure suggests: In fact, it is the second iteration that is best (has smallest maximum deviation). The routineratlsqaccordingly returns the best of its iterations, not necessarily the last one; there is no advantage in doing more than five iterations. CITED REFERENCES AND FURTHER READING: Ralston,A.andWilf,H.S.1960, MathematicalMethodsforDigitalComputers (NewYork:Wiley), Chapter 13. [1] 5.14 Evaluation of Functions by Path Integration In computer programming,the technique of choice is not necessarily the most efficient,orelegant,orfastestexecutingone. Instead,itmaybetheonethatis quick to implement, general, and easy to check. One sometimes needs only a few, or a few thousand, evaluations of a special function, perhaps a complex valued function of a complex variable, that has manydifferentparameters,or asymptoticregimes,or both. Use of the usual tricks (series, continued fractions, rational function approximations, recurrence relations, and so forth) may result in a patchwork program with tests and branches to differentformulas. While such a programmaybe highlyef ficient in execution,it is oftennot the shortest way to the answer from a standing start. A different technique of considerable generality is direct integration of a function’sd efining differential equation –an ab initio integration for each desired