f3-2

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’s and d’s and update 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 if two 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 way as to take the most “straight line” route through the tableau to its apex, updating nsaccordingly to keep track of where we are. This route keeps the partial approxima- tions centered (insofar as possible)on the target x.T h 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 3.2RationalFunctionInterpolationandExtrapolation 105Sample 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).(xi,yi)...(xi+m,yi+m). More explicitly, suppose Ri(i+1)...(i+m)=Pµ(x) Qν(x)=p0+p1x+···+pµxµ q0+q1x+···+qνxν(3.2.1 ) Since thereare µ+ν+1unknown p’s and q’s (q0beingarbitrary),we musthave m+1= µ+ν+1 ( 3.2.2 ) In specifying a rational function interpolating function, you must give the desired order of both the numerator and the denominator. Rational functions are sometimes superior to polynomials, roughly speaking, becauseoftheirabilitytomodelfunctionswithpoles,thatis,zerosofthedenominatorof equation (3.2.1). These poles might occur for real values of x, if the function to be interpolated itself has poles. More often, the function f(x)is finite for all finiterealx, but has an analytic continuation with poles in the complex x-plane. Such poles can themselves ruin a polynomial approximation, even one restricted to real values of x, just as they can ruin the convergence of an infinite power series inx. If you draw a circle in the complex plane around your mtabulated points, then you should not expect polynomial interpolation to be good unless the nearest pole is rather far outside the circle. A rational function approximation,by contrast,will stay “good”as long as it has enoughpowers of xin its denominatorto account for (cancel) any nearby poles. For the interpolation problem, a rational function is constructed so as to go through a chosen set of tabulated functional values. However, we should also mention in passing that rational function approximations can be used in analytic work. One sometimes constructs a rational function approximationby the criterion that the rational function of equation (3.2.1) itself have a power series expansion that agrees with the first m+1terms of the power series expansion of the desired function f(x). This is called Pad ´eapproximation ,and is discussed in §5.12. Bulirsch and Stoer found an algorithm of the Neville type which performs rational function extrapolation on tabulated data. A tableau like that of equation (3.1.2) is constructed column by column, leading to a result and an error estimate.TheBulirsch-Stoeralgorithmproducestheso-called diagonalrationalfunction,with the degrees of numerator and denominator equal (if mis even) or with the degree of the denominator larger by one (if mis odd, cf. equation 3.2.2 above). For the derivationofthealgorithm,referto [1]. Thealgorithmissummarizedbyarecurrence relation exactly analogousto equation (3.1.3)for polynomialapproximation: Ri(i+1)...(i+m)=R(i+1)...(i+m) +R(i+1)...(i+m)−Ri...(i+m−1)/parenleftBig x−xi x−xi+m/parenrightBig/parenleftBig 1−R(i+1) ...(i+m)−Ri... (i+m−1) R(i+1) ...(i+m)−R(i+1) ...(i+m−1)/parenrightBig −1 (3.2.3 ) This recurrence generates the rational functions through m+1points from the ones through mand (the term R(i+1)...(i+m−1)in equation 3.2.3) m−1points. It is started with Ri=yi (3.2.4 ) 106 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).andwith R≡[Ri(i+1)...(i+m)with m=−1] = 0 ( 3.2.5 ) Now, exactly as in equations (3.1.4) and (3.1.5) above, we can convert the recurrence (3.2.3) to one involving only the small differences Cm,i≡Ri...(i+m)−Ri...(i+m−1) Dm,i≡Ri...(i+m)−R(i+1)...(i+m)(3.2.6 ) Note that these satisfy the relation Cm+1,i−Dm+1,i=Cm,i +1−Dm,i (3.2.7 ) which is useful in proving the recurrences Dm+1,i=Cm,i +1(Cm,i +1−Dm,i)/parenleftBig x−xi x−xi+m+1/parenrightBig Dm,i−Cm,i +1 Cm+1,i=/parenleftBig x−xi x−xi+m+1/parenrightBig Dm,i(Cm,i +1−Dm,i) /parenleftBig x−xi x−xi+m+1/parenrightBig Dm,i−Cm,i +1(3.2.8 ) This recurrenceis implementedin thefollowingsubroutine,whoseuse is analogous in every way to polintin§3.1. SUBROUTINE ratint(xa,ya,n,x,y,dy) INTEGER n,NMAX REAL dy,x,y,xa(n),ya(n),TINY PARAMETER (NMAX=10,TINY=1.e-25) Largest expected value of n, and a small number. Given arrays xaandya,e a c ho fl e n g t h n, and given a value of x, this routine returns a value of yand an accuracy estimate dy. The value returned is that of the diagonal rational function, evaluated at x, which passes through the npoints (xa i,yai),i=1 ...n. INTEGER i,m,nsREAL dd,h,hh,t,w,c(NMAX),d(NMAX)ns=1 hh=abs(x-xa(1)) do 11i=1,n h=abs(x-xa(i)) if (h.eq.0.)then y=ya(i)dy=0.0return else if (h.lt.hh) then ns=ihh=h endif c(i)=ya(i)d(i)=ya(i)+TINY The TINYpart is needed to prevent a rare zero-over- zero condition. enddo 11 y=ya(ns) ns=ns-1do 13m=1,n-1 do12i=1,n-m 3.3CubicSplineInterpolation 107Sample 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).w=c(i+1)-d(i) h=xa(i+m)-x h will never be zero, since this was tested in the ini- tializing loop. t=(xa(i)-x)*d(i)/h dd=t-c(i+1)if(dd.eq.0.)pause ’failure in ratint’ This error condition indicates that the interpolating function has a pole at the re- quested value of x. dd=w/ddd(i)=c(i+1)*dd c(i)=t*dd enddo 12 if (2*ns.lt.n-m)then dy=c(ns+1) else dy=d(ns)ns=ns-1 endif y=y+dy enddo 13 return END CITED REFERENCES AND FURTHER READING: Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §2.2. [1] Gear,C.W.1971, NumericalInitialValueProblemsinOrdinaryDifferentialEquations (Englewood Cliffs, NJ: Prentice-Hall), §6.2. Cuyt, A., and Wuytack, L. 1987, Nonlinear Methods in Numerical Analysis (Amsterdam: North- Holland), Chapter 3. 3.3 Cubic Spline Interpolation Given a tabulated function yi=y(xi),i=1...N, focus attention on one particular interval, between xjandxj+1. Linear interpolation in that interval gives the interpolation formula y=Ayj+Byj+1 (3.3.1 ) where A≡xj+1−x xj+1−xjB≡1−A=x−xj xj+1−xj(3.3.2 ) Equations(3.3.1)and(3.3.2)areaspecialcaseofthegeneralLagrangeinterpolation formula (3.1.1). Since it is (piecewise) linear, equation (3.3.1) has zero second derivative in the interior of each interval, and an undefined, or infinite, second derivative at the abscissas xj. Thegoalofcubicsplineinterpolationistogetaninterpolationformula that is smooth in the first derivative, and continuous in the second derivative, both within an interval and at its boundaries. Suppose, contrary to fact, that in addition to the tabulated values of yi,w e also have tabulated values for the function’s second derivatives, y/prime/prime, that is, a set