f3-4

110 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).values ofthe interpolatedfunctionforanyvalue of xare obtainedbycalls (as many as desired) to a separate routine splint(for “splineinterpolation”): SUBROUTINE splint(xa,ya,y2a,n,x,y) INTEGER n REAL x,y,xa(n),y2a(n),ya(n) Given the arrays xa(1:n) andya(1:n) of length n, which tabulate a function (with the xa i’s in order), and given the array y2a(1:n) , which is the output from splineabove, and given a value of x, this routine returns a cubic-spline interpolated value y. INTEGER k,khi,kloREAL a,b,h klo=1 Wewill find the right place in the table by means of bisection. Thisisoptimal ifsequential callsto thisroutine areatrandom values of x. If sequential calls are in order, and closely spaced, one would do better to store previous values of kloand khiand test if they remain appropriate on the next call.khi=n 1 if (khi-klo.gt.1) then k=(khi+klo)/2 if(xa(k).gt.x)then khi=k else klo=k endif goto 1endif klo andkhinow bracket the input value of x. h=xa(khi)-xa(klo) if (h.eq.0.) pause ’bad xa input in splint’ The xa’s must be distinct. a=(xa(khi)-x)/h Cubic spline polynomial is now evaluated. b=(x-xa(klo))/h y=a*ya(klo)+b*ya(khi)+ * ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6. return END CITED REFERENCES AND FURTHER READING: De Boor, C. 1978, A Practical Guide to Splines (New York: Springer-Verlag). Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), §§4.4–4.5. Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §2.4. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), §3.8. 3.4 How to Search an Ordered Table Suppose that you have decided to use some particular interpolation scheme, such as fourth-order polynomial interpolation, to compute a function f(x)from a set of tabulated xi’s and fi’s. Then you will need a fast way of finding your place in the table of xi’s, given some particular value xat which the function evaluation is desired. This problem is not properly one of numerical analysis, but it occurs sooften in practice that it would be negligent of us to ignore it. Formally,theproblemisthis: Givenanarrayofabscissas xx(j),j=1,2, ...,n, with theelements eithermonotonicallyincreasingor monotonicallydecreasing,and givenanumber x,findaninteger jsuchthat xliesbetween xx(j)andxx(j+1).F o r 3.4HowtoSearchanOrderedTable 111Sample 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).this task, let us define fictitious array elements xx(0)andxx(n+1) equal to plus or minus infinity (in whichever order is consistent with the monotonicityof the table).Then jwill always be between 0 and n, inclusive; a returned value of 0 indicates “off-scale” at one end of the table, nindicates off-scale at the other end. In most cases, when all is said and done, it is hard to do better than bisection, which will find the right place in the table in about log 2ntries. We already did use bisection in the spline evaluation routine splintof the preceding section, so you might glanceback at that. Standingby itself, a bisection routinelooks like this: SUBROUTINE locate(xx,n,x,j) INTEGER j,nREAL x,xx(n) Given an array xx(1:n) , and given a value x, returns a value jsuch that xis between xx(j)andxx(j+1) .xx(1:n) must be monotonic, either increasing or decreasing. j=0 orj=nis returned to indicate that xis out of range. INTEGER jl,jm,jujl=0 Initialize lower ju=n+1 and upper limits. 10 if(ju-jl.gt.1)then If we are not yet done, jm=(ju+jl)/2 compute a midpoint, if((xx(n).ge.xx(1)).eqv.(x.ge.xx(jm)))then jl=jm and replace either the lower limit else ju=jm or the upper limit, as appropriate. endif goto 10 Repeat until endif the test condition 10is satisfied. if(x.eq.xx(1))then Then set the output j=1 else if(x.eq.xx(n))then j=n-1 else j=jl endif return and return. END Note the use of the logical equality relation .eqv., which is true when its two logical operands are either both true or both false. This relation allows the routine to work for both monotonically increasing and monotonically decreasing orders of xx(1:n). Search withCorrelated Values Sometimes you will be in the situation of searching a large table many times, and with nearly identical abscissas on consecutive searches. For example, you may be generating a function that is used on the right-hand side of a differential equation: Most differential-equationintegrators, as we shall see in Chapter 16, call for right-hand side evaluations at points that hop back and forth a bit, but whose trend moves slowly in the direction of the integration. In such cases it is wasteful to do a full bisection, ab initio, on each call. The following routine instead starts with a guessed position in the table. It first “hunts,” either up or down, in increments of 1, then 2, then 4, etc., until the desired value isbracketed. Second,it then bisects in the bracketedinterval. At worst, this routineis about a factor of 2 slower than locateabove (if the hunt phase expands to include thewholetable). Atbest,itcanbeafactorof log 2nfasterthan locate,ifthedesired pointisusuallyquiteclosetotheinputguess. Figure3.4.1comparesthetworoutines. 112 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).hunt phase bisection phase1 71 08 14 2232 3832 1 (a) (b)5164 Figure 3.4.1. (a) The routine locatefinds a table entry by bisection. Shown here is the sequence of steps that converge to element 51 in a table of length 64. (b) The routine huntsearches from a previous known position in the table by increasing steps, then converges by bisection. Shown here is a particularly un favorable example, converging to element 32 from element 7. A favorable example would be convergence to an element near 7, such as 9, which would require just three “hops.” SUBROUTINE hunt(xx,n,x,jlo) INTEGER jlo,nREAL x,xx(n) G i v e na na r r a y xx(1:n) , and given a value x, returns a value jlosuch that xis between xx(jlo) andxx(jlo+1) .xx(1:n) must be monotonic, either increasing or decreasing. jlo=0orjlo=nis returned to indicate that xis out of range. jloon input is taken as the initial guess for jloon output. INTEGER inc,jhi,jmLOGICAL ascndascnd=xx(n).ge.xx(1) True if ascending order of table, false otherwise. if(jlo.le.0.or.jlo.gt.n)then Inputguess notuseful. Goimmediatelytobisection. jlo=0jhi=n+1 goto 3 endifinc=1 Set the hunting increment. if(x.ge.xx(jlo).eqv.ascnd)then Hunt up: 1 jhi=jlo+inc if(jhi.gt.n)then Done hunting, since off end of table. jhi=n+1 else if(x.ge.xx(jhi).eqv.ascnd)then Not done hunting, jlo=jhi inc=inc+inc so double the increment goto 1 and try again. endif Done hunting, value bracketed. else Hunt down: jhi=jlo 2 jlo=jhi-inc if(jlo.lt.1)then Done hunting, since off end of table. jlo=0 else if(x.lt.xx(jlo).eqv.ascnd)then Not done hunting, jhi=jlo inc=inc+inc so double the increment goto 2 and try again. endif Done hunting, value bracketed. endif Hunt is done, so begin the final bisection phase: 3 if(jhi-jlo.eq.1)then if(x.eq.xx(n))jlo=n-1 if(x.eq.xx(1))jlo=1 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 coef ficients of that poly- nomial. A valid use of the coef ficients 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,pleasebecertainthatthecoef ficientsarewhatyouneed. Generallythe coefficientsof the interpolatingpolynomialcan be determinedmuchless accurately than 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.