f3-3

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 xjand xj+1. Linear interpolation in that interval gives the interpolation formula y=Ay j+By j+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 108 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).of numbers y/prime/prime i. Then, within each interval, we can add to the right-hand side of equation (3.3.1) a cubic polynomial whose second derivative varies linearly from avalue y /prime/prime jon the left to a value y/prime/prime j+1on the right. Doingso, we will havethe desired continuous second derivative. If we also construct the cubic polynomial to have zerovaluesatxjand xj+1, then adding it in will not spoil the agreement with the tabulated functional values yjand yj+1at the endpoints xjand xj+1. A little side calculation shows that there is only one way to arrange this construction, namely replacing (3.3.1) by y=Ay j+By j+1+Cy/prime/prime j+Dy/prime/prime j+1 (3.3.3 ) where Aand Bare defined in (3.3.2) and C≡1 6(A3−A)(xj+1−xj)2D≡1 6(B3−B)(xj+1−xj)2(3.3.4 ) Notice that the dependence on the independent variable xin equations (3.3.3) and (3.3.4)is entirely throughthe linear x-dependenceof Aand B, and (through Aand B) the cubic x-dependence of Cand D. We can readily check that y/prime/primeis in fact the second derivative of the new interpolating polynomial. We take derivatives of equation (3.3.3) with respect to x, usingthedefinitionsof A, B, C, D tocompute dA/dx, dB/dx, dC/dx ,and dD/dx. The result is dy dx=yj+1−yj xj+1−xj−3A2−1 6(xj+1−xj)y/prime/prime j+3B2−1 6(xj+1−xj)y/prime/prime j+1 (3.3.5 ) for the first derivative, and d2y dx2=Ay/prime/prime j+By/prime/prime j+1 (3.3.6 ) for the second derivative. Since A=1atxj,A=0atxj+1, while Bis just the other way around, (3.3.6) shows that y/prime/primeis just the tabulated second derivative, and alsothatthesecondderivativewillbecontinuousacross(e.g.) theboundarybetween the two intervals (xj−1,x j)and (xj,x j+1). Theonlyproblemnowisthatwesupposedthe y/prime/prime i’stobeknown,when,actually, they are not. However, we have not yet required that the firstderivative, computed fromequation(3.3.5),becontinuousacrosstheboundarybetweentwointervals. The key idea of a cubic spline is to require this continuity and to use it to get equationsfor the second derivatives y /prime/prime i. The required equations are obtained by setting equation (3.3.5) evaluated for x=xjintheinterval (xj−1,x j)equaltothesameequationevaluatedfor x=xjbut intheinterval (xj,x j+1). Withsomerearrangement,thisgives(for j=2 ,...,N −1) xj−xj−1 6y/prime/prime j−1+xj+1−xj−1 3y/prime/prime j+xj+1−xj 6y/prime/prime j+1=yj+1−yj xj+1−xj−yj−yj−1 xj−xj−1 (3.3.7 ) These are N−2linear equations in the Nunknowns y/prime/prime i,i=1 ,...,N. Therefore there is a two-parameter family of possible solutions. Forauniquesolution,weneedtospecifytwofurtherconditions,typicallytaken asboundaryconditionsat x1and xN. Themostcommonwaysofdoingthisareeither 3.3CubicSplineInterpolation 109Sample 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).•set one or both of y/prime/prime 1and y/prime/prime Nequal to zero, giving the so-called natural cubic spline , which has zero second derivative on one or both of its boundaries, or •set either of y/prime/prime 1and y/prime/prime Nto values calculated from equation (3.3.5) so as to make the first derivative of the interpolating function have a specifiedvalue on either or both boundaries. Onereasonthatcubicsplinesareespeciallypracticalisthatthesetofequations (3.3.7), along with the two additional boundary conditions, are not only linear, but alsotridiagonal . Each y /prime/prime jiscoupledonlytoitsnearestneighborsat j±1. Therefore, the equationscan be solvedin O(N)operationsby the tridiagonalalgorithm( §2.4). That algorithmis concise enoughto build right into the spline calculational routine. This makes the routine not completely transparent as an implementation of (3.3.7), so we encourageyou to study it carefully, comparingwith tridag(§2.4). SUBROUTINE spline(x,y,n,yp1,ypn,y2) INTEGER n,NMAX REAL yp1,ypn,x(n),y(n),y2(n) PARAMETER (NMAX=500) Given arrays x(1:n) andy(1:n) containing a tabulated function, i.e., yi=f(xi),w i t h x1<x2<. . .< xN, and given values yp1 andypn for the first derivative of the inter- polating function at points 1 and n, respectively, this routine returns an array y2(1:n) of length nwhich contains the second derivatives of the interpolating function at the tabulated points xi.I fyp1 and/or ypn are equal to 1×1030or larger, the routine is signaled to set the corresponding boundary condition for a natural spline, with zero second derivative on that boundary.Parameter: NMAX is the largest anticipated value of n. INTEGER i,k REAL p,qn,sig,un,u(NMAX) if (yp1.gt..99e30) then The lower boundary condition is set either to be “natural” y2(1)=0. u(1)=0. else or else to have a specified first derivative. y2(1)=-0.5u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) endif do 11i=2,n-1 This is the decomposition loop of the tridiagonal algorithm. y2and uare used for temporary storage of the decomposed factors.sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) p=sig*y2(i-1)+2. y2(i)=(sig-1.)/p u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) * /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p enddo 11 if (ypn.gt..99e30) then The upper boundary condition is set either to be “natural” qn=0. un=0. else or else to have a specified first derivative. qn=0.5un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) endif y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)do 12k=n-1,1,-1 This is the backsubstitution loop of the tridiago- nal algorithm. y2(k)=y2(k)*y2(k+1)+u(k) enddo 12 return END It is important to understand that the program splineis called only onceto process an entire tabulated function in arrays xiandyi. Once this has been done, 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 spline above, and given a value of x, this routine returns a cubic-spline interpolated value y. INTEGER k,khi,kloREAL a,b,h klo=1 We will find the right place in the table by means of bisection. This is optimal if sequential calls to this routine are at random values of x. If sequential calls are in order, and closely spaced, one would do better to store previous values of klo and khi and 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 andkhi now 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