f3-5

3.5CoefficientsoftheInterpolatingPolynomial 113Sample 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).return endif jm=(jhi+jlo)/2 if(x.ge.xx(jm).eqv.ascnd)then jlo=jm else jhi=jm endifgoto 3 END After the Hunt The problem: Routines locateandhuntreturn an index jsuch that your desired value lies between table entries xx(j)andxx(j+1), where xx(1:n) is the full length of the table. But, to obtain an m-point interpolated value using a routine likepolint(§3.1) or ratint(§3.2), you need to supply much shorter xxandyy arrays, of length m. How do you make the connection? The solution: Calculate k=min(max(j-(m-1)/2,1),n+1-m) This expression produces the index of the leftmost member of an m-point set of points centered (insofar as possible) between jandj+1, but bounded by 1 at the left and nat the right. FORTRAN then lets you call the interpolation routine with array addresses offset by k, e.g., call polint(xx(k),yy(k),m, ...) CITED REFERENCES AND FURTHER READING: Knuth,D.E. 1973, SortingandSearching ,v ol.3of TheArtofComputerProgramming (Reading, MA: Addison-Wesley), §6.2.1. 3.5 Coefficients of the Interpolating Polynomial Occasionallyyoumaywishtoknownotthevalueoftheinterpolatingpolynomial that passes through a (small!) number of points, but the coefficients of that poly- nomial. A valid use of the coefficients might be, for example, to compute simultaneousinterpolatedvaluesofthefunctionandofseveralofitsderivatives(see §5.3), or to convolve a segment of the tabulated function with some other function, where the moments of that other function (i.e., its convolution with powers of x) are known analytically. However,pleasebecertainthatthecoefficientsarewhatyouneed. Generallythe coefficientsof the interpolatingpolynomialcan be determinedmuchless accuratelythan its value at a desired abscissa. Thereforeit is not a good idea to determine the coefficients only for use in calculating interpolating values. Values thus calculated willnotpassexactlythroughthetabulatedpoints,forexample,whilevaluescomputed by the routines in §3.1–§3.3 will pass exactly through such points. 114 Chapter3. InterpolationandExtrapolationSample 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).Also,youshouldnotmistaketheinterpolatingpolynomial(anditscoefficients) for its cousin, the best fitpolynomial through a data set. Fitting is a smoothing process, since the number of fitted coefficients is typically much less than the number of data points. Therefore, fitted coefficients can be accurately and stably determined even in the presence of statistical errors in the tabulated values. (See§14.8.) Interpolation, where the number of coefficients and number of tabulated pointsareequal,takesthetabulatedvaluesasperfect. Iftheyinfactcontainstatistical errors, these can be magnified into oscillations of the interpolating polynomial in between the tabulated points. As before, we take the tabulated points to be y i≡y(xi). If the interpolating polynomial is written as y=c1+c2x+c3x2+··· +cNxN−1(3.5.1 ) then the ci’s are required to satisfy the linear equation  1 x 1 x2 1··· xN−1 1 1 x2 x2 2··· xN−1 2 ............ 1 x N x2 N··· xN−1 N · c 1 c2 ... cN = y 1 y2 ... yN (3.5.2 ) This is a Vandermonde matrix , as described in §2.8. One could in principle solve equation(3.5.2)bystandardtechniquesforlinearequationsgenerally( §2.3);however the special method that was derived in §2.8 is more efficient by a large factor, of order N, so it is much better. Remember that Vandermonde systems can be quite ill-conditioned. In such a case,nonumerical method is going to give a very accurate answer. Such cases do not, please note, imply any difficulty in finding interpolated valuesby the methods of§3.1, but only difficulty in finding coefficients . Like the routine in §2.8, the following is due to G.B. Rybicki. SUBROUTINE polcoe(x,y,n,cof) INTEGER n,NMAXREAL cof(n),x(n),y(n)PARAMETER (NMAX=15) Largestanticipatedvalueof n. Givenarrays x(1:n)andy(1:n)containingatabulatedfunction yi=f(xi),thisroutine returnsanarrayofcoefficients cof(1:n),suchthat yi=/summationtext jcof jxj−1 i. INTEGER i,j,k REAL b,ff,phi,s(NMAX) do11i=1,n s(i)=0.cof(i)=0. enddo 11 s(n)=-x(1) do13i=2,n Coefficients siofthemasterpolynomial P(x)arefound byrecurrence. do12j=n+1-i,n-1 s(j)=s(j)-x(i)*s(j+1) enddo 12 s(n)=s(n)-x(i) enddo 13 do16j=1,n phi=n do14k=n-1,1,-1 Thequantity phi =/producttext j/negationslash=k(xj−xk)isfoundasaderiva- tiveof P(xj). 3.5CoefficientsoftheInterpolatingPolynomial 115Sample 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).phi=k*s(k+1)+x(j)*phi enddo 14 ff=y(j)/phib=1. CoefficientsofpolynomialsineachtermoftheLagrange formulaarefoundbysyntheticdivisionof P(x)by (x−x j).T h es o l u t i o n ckisaccumulated.do15k=n,1,-1 cof(k)=cof(k)+b*ff b=s(k)+x(j)*b enddo 15 enddo 16 returnEND Another Method Another technique is to make use of the function value interpolation routine already given ( polint §3.1). If we interpolate (or extrapolate) to find the value of the interpolating polynomial at x=0, then this value will evidently be c1.N o w we can subtract c1fromthe yi’s and divideeach by its corresponding xi. Throwing out one point (the one with smallest xiis a good candidate), we can repeat the procedure to find c2, and so on. It is not instantly obvious that this procedure is stable, but we have generally found it to be somewhat morestable than the routine immediately preceding. This method is of order N3, while the preceding one was of order N2. You will find, however, that neither works very well for large N, because of the intrinsic ill-condition of the Vandermonde problem. In single precision, Nu pt o8o r1 0i s satisfactory; about double this in double precision. SUBROUTINE polcof(xa,ya,n,cof) INTEGER n,NMAXREAL cof(n),xa(n),ya(n)PARAMETER (NMAX=15) Largestanticipatedvalueof n. C USES polint Givenarrays xa(1:n)andya(1:n)oflength ncontainingatabulatedfunction yai= f(xa i),thisroutinereturnsanarrayofcoefficients cof(1:n),alsooflength n,suchthat yai=/summationtext jcof jxaj−1 i. INTEGER i,j,kREAL dy,xmin,x(NMAX),y(NMAX) do 11j=1,n x(j)=xa(j)y(j)=ya(j) enddo 11 do14j=1,n call polint(x,y,n+1-j,0.,cof(j),dy) Thisisthepolynomialinterpolationrou- tineof §3.1. Weextrapolateto x= 0.xmin=1.e38 k=0 do12i=1,n+1-j Findtheremaining xiofsmallestabso- lutevalue, if (abs(x(i)).lt.xmin)then xmin=abs(x(i)) k=i endifif(x(i).ne.0.)y(i)=(y(i)-cof(j))/x(i) (meanwhilereducingalltheterms) enddo 12 do13i=k+1,n+1-j andeliminateit. y(i-1)=y(i)x(i-1)=x(i) enddo 13 enddo 14 returnEND 116 Chapter3. InterpolationandExtrapolationSample 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).If the point x=0is not in (or at least close to) the range of the tabulated xi’s, thenthecoefficientsoftheinterpolatingpolynomialwillingeneralbecomeverylarge.However, the real “information content” of the coefficients is in small differencesfrom the “translation-induced” large values. This is one cause of ill-conditioning, resulting in loss of significance and poorly determined coefficients. You should consider redefiningthe origin of the problem,to put x=0in a sensible place. Another pathologyis that, if too high a degree of interpolationis attempted on a smooth function, the interpolating polynomial will attempt to use its high-degree coefficients,incombinationswithlargeandalmostpreciselycancelingcombinations, to match the tabulated values down to the last possible epsilon of accuracy. This effect is the same as the intrinsic tendencyof the interpolatingpolynomialvalues tooscillate (wildly) between its constrained points, and would be present even if the machine’s floating precision were infinitely good. The above routines polcoeand polcofhave slightly different sensitivities to the pathologies that can occur. Are you still quite certain that using the coefficients is a good idea? CITED REFERENCES AND FURTHER READING: Isaacson, E., and Keller, H.B. 1966, Analysis of Numerical Methods (New York: Wiley), §5.2. 3.6 Interpolation in Two or More Dimensions In multidimensional interpolation, we seek an estimate of y(x1,x2,...,x n) from an n-dimensional grid of tabulated values yand none-dimensional vec- tors giving the tabulated values of each of the independent variables x1,x2,..., xn. We will not here consider the problem of interpolating on a mesh that is not Cartesian, i.e., has tabulated function values at “random” points in n-dimensional space rather than at the vertices of a rectangular array. For clarity, we will consider explicitly only the case of two dimensions, the cases of three or more dimensions being analogous in every way. In two dimensions, we imagine that we are given a matrix of functionalvalues ya(j,k), where jvaries from 1 to m, and kvaries from 1 to n. We are also given an array x1aof length m, and an array x2aof length n. The relation of these input quantities to an underlying function y(x1,x2)is ya(j,k) =y(x1a(j) ,x2a(k) )( 3.6.1 ) We want to estimate, by interpolation, the function yat some untabulated point (x1,x2). An important concept is that of the grid square in which the point (x1,x2) falls, that is, the four tabulated points that surround the desired interior point. For convenience, we will number these points from 1 to 4, counterclockwise startingfrom the lower left (see Figure 3.6.1). More precisely, if x1a(j) ≤x 1≤x1a(j+1) x2a(k) ≤x2≤x2a(k+1)(3.6.2 )