f3-1

102 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).f(x, y, z ). Multidimensional interpolation is often accomplished by a sequence of one-dimensional interpolations. We discuss this in §3.6. CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), §25.2. Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), Chapter 2. Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), Chapter 3. Kahaner,D.,Moler,C.,andNash,S.1989, NumericalMethods andSoftware (EnglewoodCliffs, NJ: Prentice Hall), Chapter 4. Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis , 2nd ed. (Reading, MA: Addison- Wesley), Chapter 5. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), Chapter 3. Isaacson,E.,andKeller,H.B.1966, AnalysisofNumericalMethods (NewYork:Wiley),Chapter6. 3.1 PolynomialInterpolation and Extrapolation Through any two points there is a unique line. Through any three points, a uniquequadratic. Et cetera. The interpolatingpolynomialof degree N−1through the Npoints y1=f(x1),y2=f(x2),...,y N=f(xN)is given explicitly by Lagrange’s classical formula, P(x)=(x−x2)(x−x3)...(x−xN) (x1−x2)(x1−x3)...(x1−xN)y1+(x−x1)(x−x3)...(x−xN) (x2−x1)(x2−x3)...(x2−xN)y2 +··· +(x−x1)(x−x2)...(x−xN−1) (xN−x1)(xN−x2)...(xN−xN−1)yN (3.1.1 ) There are Nterms, each a polynomial of degree N−1and each constructed to be zero at all of the xiexcept one, at which it is constructed to be yi. It is not terribly wrong to implement the Lagrange formula straightforwardly, butitis notterriblyrighteither. Theresultingalgorithmgivesnoerrorestimate,and it is also somewhatawkwardto program. A muchbetteralgorithm(forconstructing thesame,unique,interpolatingpolynomial)is Neville’salgorithm ,closelyrelatedto andsometimesconfusedwith Aitken’salgorithm ,thelatternowconsideredobsolete. Let P1be the value at xof the unique polynomial of degree zero (i.e., a constant) passing through the point (x1,y1);s o P1=y1. Likewise define P2,P 3,...,P N.Now let P12be the value at xof the unique polynomial of degree one passing through both (x1,y1)and (x2,y2). Likewise P23,P 34,..., P(N−1)N. Similarly,forhigher-orderpolynomials,upto P123 ...N,whichisthevalue oftheuniqueinterpolatingpolynomialthroughall Npoints,i.e.,thedesiredanswer. 3.1PolynomialInterpolationandExtrapolation 103Sample 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).The various P’s form a “tableau” with “ancestors” on the left leading to a single “descendant” at the extreme right. For example, with N=4, x1: y1=P1 P12 x2: y2=P2 P123 P23 P1234 x3: y3=P3 P234 P34 x4: y4=P4(3.1.2 ) Neville’s algorithm is a recursive way of filling in the numbers in the tableau a column at a time, from left to right. It is based on the relationship between a “daughter” Pand its two “parents,” Pi(i+1) ...(i+m)=(x−xi+m)Pi(i+1) ...(i+m−1)+(xi−x)P(i+1)( i+2) ...(i+m) xi−xi+m (3.1.3 ) This recurrence works because the two parents already agree at points xi+1... xi+m−1. An improvement on the recurrence (3.1.3) is to keep track of the small differences between parents and daughters, namely to define (for m=1 ,2,..., N−1), Cm,i≡Pi...(i+m)−Pi...(i+m−1) Dm,i≡Pi...(i+m)−P(i+1) ...(i+m).(3.1.4 ) Then one can easily derive from (3.1.3) the relations Dm+1 ,i=(xi+m+1−x)(Cm,i +1−Dm,i) xi−xi+m+1 Cm+1 ,i=(xi−x)(Cm,i +1−Dm,i) xi−xi+m+1(3.1.5 ) At eachlevel m,the C’s and D’s are thecorrectionsthatmaketheinterpolationone orderhigher. The final answer P1...Nis equal to the sum of any yiplus a set of C’s and/or D’s that form a path throughthe family tree to the rightmost daughter. Here is a routine for polynomial interpolation or extrapolation: SUBROUTINE polint(xa,ya,n,x,y,dy) INTEGER n,NMAXREAL dy,x,y,xa(n),ya(n) PARAMETER (NMAX=10) Largest anticipated value of n. Given arrays xaandya,e a c ho fl e n g t h n, and given a value x, this routine returns a value y, and an error estimate dy.I f P(x)is the polynomial of degree N−1such that P(xa i)=yai,i=1 ,..., n, then the returned value y=P(x). INTEGER i,m,ns REAL den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)ns=1 dif=abs(x-xa(1)) 104 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).do11i=1,n Here we find the index nsof the closest table entry, dift=abs(x-xa(i)) if (dift.lt.dif) then ns=idif=dift endif c(i)=ya(i) and initialize the tableau of c’s and d’s. d(i)=ya(i) enddo 11 y=ya(ns) This is the initial approximation to y. ns=ns-1do 13m=1,n-1 For each column of the tableau, do12i=1,n-m we loop over the current c’sand d’sandupdate them. ho=xa(i)-x hp=xa(i+m)-xw=c(i+1)-d(i) den=ho-hp if(den.eq.0.)pause ’failure in polint’ This error can occur only iftwo input xa’s are(to within roundoff) identical. den=w/den d(i)=hp*den Here the c’s and d’s are updated. c(i)=ho*den enddo 12 if (2*ns.lt.n-m)then After each column in the tableau is completed, we decide which correction, cord, we want to add to our accu- mulating value of y, i.e., which path to take through the tableau—forking up or down. We do this in such a wayastotakethemost“straightline”routethroughthe tableautoitsapex,updating nsaccordinglytokeeptrack ofwhereweare. Thisroutekeepsthepartialapproxima- tions centered (insofar as possible) onthe target x.Th e lastdyadded is thus the error indication.dy=c(ns+1) else dy=d(ns) ns=ns-1 endify=y+dy enddo 13 return END Quite often you will want to call polintwith the dummy arguments xa andyareplaced by actual arrays with offsets . For example, the construction call polint(xx(15),yy(15),4,x,y,dy) performs4-point interpolation on the tabulatedvalues xx(15:18) ,yy(15:18) . For moreon this, see the end of §3.4. CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), §25.2. Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §2.1. Gear,C.W.1971, NumericalInitialValueProblemsinOrdinaryDifferentialEquations (Englewood Cliffs, NJ: Prentice-Hall), §6.1. 3.2 Rational Function Interpolation and Extrapolation Some functions are not well approximated by polynomials, but arewell approximated by rational functions, that is quotients of polynomials. We de- note by Ri(i+1) ...(i+m)a rational function passing through the m+1points