f5-12

194 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).INTEGER NMANY,NFEW REAL e(NMANY),d(NFEW),c(NMANY),a,b Economize NMANYpower series coefficients e(1:NMANY) in the range (a, b )into NFEW coefficients d(1:NFEW) . call pcshft((-2.-b-a)/(b-a),(2.-b-a)/(b-a),e,NMANY) call pccheb(e,c,NMANY) ... Here one would normally examine the Chebyshev coefficients c(1:NMANY) to decide how small NFEWcan be. call chebpc(c,d,NFEW) call 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– 10terms 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 inthis 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 rateof convergence of Chebyshev coefficients has nothing to do with the rate of convergence ofpower series coefficients; and it is the formerthat 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 polynomialapproximations. These can be found by the methods of §5.8, butnotby 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: Mathe- matical Association of America), Chapter 12. Arfken, G. 1970, Mathematical Methods for Physicists , 2nd ed. (New York: Academic Press), p. 631. [1] 5.12 Pad´e Approximants APad´e 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. Ifthe rational function is R(x)≡M/summationdisplay k=0akxk 1+N/summationdisplay k=1bkxk(5.12.1 ) then R(x)is said to be a Pad ´e approximant to the series f(x)≡∞/summationdisplay k=0ckxk(5.12.2 ) if R(0) = f(0) ( 5.12.3 ) 5.12Pad´eApproximants 195Sample 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).and also dk dxkR(x)/vextendsingle/vextendsingle/vextendsingle/vextendsingle x=0=dk dxkf(x)/vextendsingle/vextendsingle/vextendsingle/vextendsingle x=0,k =1 ,2,...,M +N (5.12.4 ) Equations (5.12.3) and (5.12.4) furnish M+N+1equations for the unknowns a0,...,a M and b1,...,b N. 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 ofxthat 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 a 0=c0, with the remaining a’s and b’s satisfying N/summationdisplay m=1bmcN−m+k=−cN+k,k =1 ,...,N (5.12.5 ) k/summationdisplay m=0bmck−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 frequentlyclose to singular, so that one should not solve them by the methods of §2.8, but rather by full LUdecomposition. Additionally, it is a good idea to refine the solution by iterative improvement (routine mprovein§2.5) [1]. Oncethe b’sareknown,thenequation(5.12.6)givesanexplicitformulafortheunknown a’s, completing the solution. Pad´e 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/prime(0),f/prime/prime(0), and so on. These are of course the first few coefficients in the power series expansion off(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, Pad ´e 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 inmathematics, 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 9x+1 81x2−49 8748x3+175 78732x4+··· (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 thatthey derive from some side analytic calculation.) Equation (5.12.7) is plotted as the curvelabeled “power series” in Figure 5.12.1. One sees that for x >∼4it is dominated by its largest, quartic, term. We now take the five coefficients in equation (5.12.7) and run them through the routine padelistedbelow. Itreturnsfiverationalcoefficients,three a’sandtwo b’s,foruseinequation (5.12.1) with M=N=2. The curve in thefigure labeled “Pad ´e” plotsthe resulting rational function. Note that both solid curves derive from the samefive original coefficient values. Toevaluate the results,we need Deus exmachina (auseful fellow, when he isavailable) to tell us that equation (5.12.7) is in fact the power series expansion of the function f(x)=[ 7+( 1+ x)4/3]1/3(5.12.8 ) which isplotted asthedotted curve inthefigure. Thisfunctionhas abranch point at x=−1, so its power series is convergent only in the range −1<x< 1. In most of the range shown in the figure, the series is divergent, and the value of its truncation to five terms is 196 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).0246810 02468 1 0f(x) xPadé (5 coefficients) exactpower 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 Pad´e approximant for a sample function f(x). The full power series converges only for x< 1. Note that the Pad ´e approximant maintains accuracy far outside the radius of convergence of the series. rather meaningless. Nevertheless, those five terms, converted to a Pad ´e approximant, give a remarkably good representation of the function up to at least x∼10. Why does this work? Are there not other functions with the same firstfive terms in their power series, but completely different behavior in the range (say) 2<x< 10? Indeed there are. Pad ´e approximation has the uncanny knack of picking the function you had in mindfrom among all the possibilities. Except when it doesn ’t!That is the downside of Pad´e approximation: it is uncontrolled. There is, in general, no way to tell how accurate it is, or how far out in xit can usefully be extended. It is a powerful, but in the end still mysterious, technique. Hereistheroutinethatgets a’sand b’sfromyour c’s. Notethattheroutineisspecialized tothe case M=N,and also that,on output, therationalcoef ficients arearranged inaformat for use with the evaluation routine ratval(§5.3). (Also for consistency with that routine, the array of c’s is passed in double precision.) SUBROUTINE pade(cof,n,resid) INTEGER n,NMAX REAL resid,BIGDOUBLE PRECISION cof(2*n+1) For consistency with ratval. PARAMETER (NMAX=20,BIG=1.E30) Max expected value of n, and a big number. C USES lubksb,ludcmp,mprove Given cof(1:2*n+1) , the leading termsin the power series expansion ofafunction, solve the linear Pad´ e equations to return the coefficients of a diagonal rational function approx- imation to the same function, namely (cof(1) +cof(2) x+ ··· +cof(n+1) xN)/(1 + cof(n+2) x+··· +cof(2*n+1) xN).T h ev a l u e residisthenormoftheresidualvector; a small value indicates a well-converged solution. INTEGER j,k,indx(NMAX) REAL d,rr,rrold,sum,q(NMAX,NMAX),qlu(NMAX,NMAX),x(NMAX), * y(NMAX),z(NMAX) do12j=1,n Set up matrix for solving. x(j)=cof(n+j+1) 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) Solve by LUdecomposition and backsubstitution. 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 coef ficients. 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.