f3-0

Sample 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).Chapter 3. Interpolation and Extrapolation 3.0 Introduction Wesometimesknowthevalueofafunction f(x)atasetofpoints x1,x2,...,x N (say,with x1<...<x N),butwedon’thaveananalyticexpressionfor f(x)thatlets uscalculateitsvalueatanarbitrarypoint. Forexample,the f(xi)’smightresultfrom some physical measurement or from long numerical calculation that cannot be castintoasimplefunctionalform. Oftenthe x i’s areequallyspaced,butnotnecessarily. The task now is to estimate f(x)for arbitrary xby, in some sense, drawing a smoothcurvethrough(andperhapsbeyond)the xi. Ifthedesired xisinbetweenthe largest and smallest of the xi’s, the problem is called interpolation ;i fxis outside thatrange,itiscalled extrapolation ,whichisconsiderablymorehazardous(asmany former stock-market analysts can attest). Interpolation and extrapolation schemes must model the function, between or beyond the known points, by some plausible functional form. The form shouldbe sufficiently general so as to be able to approximate large classes of functions which might arise in practice. By far most common among the functional forms usedarepolynomials( §3.1). Rationalfunctions(quotientsofpolynomials)alsoturn out to be extremely useful ( §3.2). Trigonometric functions, sines and cosines, give rise totrigonometric interpolation and related Fourier methods, which we defer to Chapters 12 and 13. There is an extensive mathematical literature devoted to theorems about what sort of functionscan be well approximatedby which interpolatingfunctions. Thesetheorems are, alas, almost completely useless in day-to-day work: If we know enough about our function to apply a theorem of any power, we are usually not in the pitiful state of having to interpolate on a table of its values! Interpolationis related to, but distinct from, functionapproximation . That task consists of finding an approximate (but easily computable) function to use in place ofamorecomplicatedone. Inthecaseofinterpolation,youaregiventhefunction f at pointsnotof your ownchoosing . Forthe case offunctionapproximation,youare allowedtocomputethefunction fatanydesiredpointsforthepurposeofdeveloping your approximation. We deal with function approximationin Chapter 5. Onecaneasilyfindpathologicalfunctionsthatmakea mockeryofanyinterpo- lation scheme. Consider, for example, the function f(x)=3 x 2+1 π4ln/bracketleftbig (π−x)2/bracketrightbig +1 ( 3.0.1 ) 99 100 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).which is well-behaved everywhere except at x=π, very mildly singular at x=π, and otherwise takes on all positive and negativevalues. Any interpolationbased onthe values x=3.13,3.14,3.15,3.16, will assuredly get a very wrong answer for the value x=3.1416, even though a graph plotting those five points looks really quite smooth! (Try it on your calculator.) Because pathologies can lurk anywhere, it is highly desirable that an interpo- lation and extrapolation routine should return an estimate of its own error. Such an error estimate can never be foolproof, of course. We could have a function that, for reasons known only to its maker, takes off wildly and unexpectedly between two tabulated points. Interpolation always presumes some degree of smoothnessfor the function interpolated, but within this framework of presumption, deviations from smoothness can be detected. Conceptually,the interpolationprocess has two stages: (1) Fit an interpolating function to the data points provided. (2) Evaluate that interpolating function at the target point x. However, this two-stage method is generally not the best way to proceed in practice. Typically it is computationally less efficient, and more susceptible to roundoff error, than methods which construct a functional estimate f(x)directly fromthe Ntabulatedvalues everytime oneis desired. Most practical schemes start at a nearby point f(x i), then add a sequence of (hopefully) decreasing corrections, as information from other f(xi)’s is incorporated. The procedure typically takes O(N2)operations. If everything is well behaved, the last correction will be the smallest, andit canbeusedas aninformal(thoughnotrigorous)boundonthe error. In the case of polynomial interpolation, it sometimes does happen that the coefficients of the interpolating polynomial are of interest, even though their use inevaluating the interpolating function should be frowned on. We deal with this eventuality in §3.5. Local interpolation, using a finite number of “nearest-neighbor” points, gives interpolated values f(x)that do not, in general, have continuous first or higher derivatives. That happens because, as xcrosses the tabulated values xi, the interpolation scheme switches which tabulated points are the “local” ones. (If such a switch is allowed to occur anywhere else, then there will be a discontinuityin the interpolated function itself at that point. Bad idea!) In situations where continuity of derivatives is a concern, one must use the “stiffer” interpolation provided by a so-called splinefunction. A spline is a polynomial between each pair of table points, but one whose coefficients are determined “slightly” nonlocally. The nonlocality is designed to guarantee globalsmoothnessintheinterpolatedfunctionuptosomeorderofderivative. Cubicsplines (§3.3)arethemostpopular. Theyproduceaninterpolatedfunctionthatiscontinuous throughthesecondderivative. Splinestendtobestablerthanpolynomials,withlesspossibility of wild oscillation between the tabulated points. The number of points (minus one) used in an interpolation scheme is called theorderof the interpolation. Increasing the order does not necessarily increase the accuracy, especially in polynomial interpolation. If the added points are distant fromthepointofinterest x,theresultinghigher-orderpolynomial,withitsadditional constrained points, tends to oscillate wildly between the tabulated values. This oscillation may have no relation at all to the behavior of the “true” function (see Figure 3.0.1). Of course, addingpoints closeto the desired point usually does help, 3.0Introduction 101Sample 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).(a) (b) Figure 3.0.1. (a) A smooth function (solid line) is more accurately interpolated by a high-order polynomial (shown schematically as dotted line) than by a low-order polynomial (shown as a piecewise linear dashed line). (b) A function with sharp corners or rapidly changing higher derivatives is less accurately approximated byahigh-orderpolynomial (dottedline),whichistoo “stiff,”thanbyalow-order polynomial (dashed lines). Even some smooth functions, such as exponentials or rational functions, canbe badly approximated by high-order polynomials. but afiner mesh implies a larger table of values, not always available. Unless there is solidevidencethat the interpolatingfunctionis close in formto the true function f, it is a good idea to be cautious about high-order interpolation. Weenthusiasticallyendorseinterpolationswith3or4points,weareperhapstolerant of 5 or 6; but we rarelygo higherthanthat unless thereis quite rigorousmonitoring of estimated errors. When your table of values contains many more points than the desirable order ofinterpolation,youmustbegineachinterpolationwithasearchfortheright “local” placeinthetable. Whilenotstrictlyapartofthesubjectofinterpolation,thistaskis importantenough(andoftenenoughbotched)that we devote §3.4to its discussion. The routines given for interpolation are also routines for extrapolation. An important application, in Chapter 16, is their use in the integration of ordinary differential equations. There, considerable care istaken with the monitoring of errors. Otherwise, the dangers of extrapolation cannot be overemphasized: Aninterpolatingfunction,which is perforcean extrapolatingfunction,will typically go berserk when the argument xis outside the range of tabulated values by more than the typical spacing of tabulated points. Interpolation can be done in more than one dimension, e.g., for a function 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 theNpoints 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. LetP1be 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 de fine 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.