Chapter 3. Interpolation and Extrapolation http://www.nr.com or call 3.0 Introduction 11-800-72 We sometimes know the value ofa function f()at a set of points 1,x2,...,N Cambridge (say,withz1<...<N),but we don't have an analytic expression for f(x)that lets NUMERICAL RECIPES IN us calculate its value at an arbitrary point.For example,the f(i)'s might result from server some physical measurement or from long numerical calculation that cannot be cast into a simple functional form.Often the zi's are equally spaced,but not necessarily. compu Press The task now is to estimate f(r)for arbitrary x by,in some sense,drawing a smooth curve through(and perhaps beyond)the z.If the desired r is in between the C:THE ART OF largest and smallest of the xi's,the problem is called interpolation;if is outside 马 that range,it is called extrapolation,which is considerably more hazardous (as many former stock-market analysts can attest). SCIENTIFIC Interpolation and extrapolation schemes must model the function,between or beyond the known points,by some plausible functional form.The form should 6 be 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 used are polynomials($3.1).Rational functions(quotients of polynomials)also turn out to be extremely useful(83.2).Trigonometric functions,sines and cosines,give r Numerical 92y rise to trigonometric interpolation and related Fourier methods,which we defer to Chapters 12 and 13 Recipes There is an extensive mathematical literature devoted to theorems about what sort of functions can be well approximated by which interpolating functions.These theorems are,alas,almost completely useless in day-to-day work:If we know ecipes 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! Interpolation is related to,but distinct from,function approximation.That task North Software. consists of finding an approximate(but easily computable)function to use in place of a more complicated one.In the case of interpolation,you are given the function f America) at points not ofyour own choosing.For the case of function approximation,you are visit website machine allowed to compute the function f at any desired points for the purpose of developing your approximation.We deal with function approximation in Chapter 5. One can easily find pathological functions that make a mockery of any interpo- lation scheme.Consider,for example,the function f回=3x2+-9]+1 (3.0.1) 105
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Chapter 3. Interpolation and Extrapolation 3.0 Introduction We sometimes know the value of a function f(x) at a set of points x1, x2,...,xN (say, with x1 <...<xN ), but we don’t have an analytic expression for f(x)that lets us calculate its value at an arbitrary point. For example, the f(xi)’s might result from some physical measurement or from long numerical calculation that cannot be cast into a simple functional form. Often the xi’s are equally spaced, but not necessarily. The task now is to estimate f(x) for arbitrary x by, in some sense, drawing a smooth curve through (and perhaps beyond) the x i. If the desired x is in between the largest and smallest of the xi’s, the problem is called interpolation; if x is outside that range, it is called extrapolation, which is considerably more hazardous (as many 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 should be 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 used are polynomials (§3.1). Rational functions (quotients of polynomials) also turn out to be extremely useful (§3.2). Trigonometric functions, sines and cosines, give rise to trigonometric 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 functions can be well approximated by which interpolating functions. These theorems 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! Interpolation is related to, but distinct from, function approximation. That task consists of finding an approximate (but easily computable) function to use in place of a more complicated one. In the case of interpolation, you are given the function f at points not of your own choosing. For the case of function approximation, you are allowed to compute the function f at any desired points for the purpose of developing your approximation. We deal with function approximation in Chapter 5. One can easily find pathological functions that make a mockery of any interpolation scheme. Consider, for example, the function f(x)=3x2 + 1 π4 ln (π − x) 2 +1 (3.0.1) 105
106 Chapter 3.Interpolation and Extrapolation which is well-behaved everywhere except at x=m,very mildly singular at =m. and otherwise takes on all positive and negative values.Any interpolation based on the values x 3.13,3.14,3.15,3.16,will assuredly get a very wrong answer for the value z =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 provide 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 smoothness for the function interpolated,but within this framework of presumption,deviations from smoothness can be detected. Conceptually,the interpolation process has two stages:(1)Fit an interpolating function to the data points provided.(2)Evaluate that interpolating function at the target point z. 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(z)directly from the N tabulated values every time one is desired.Most practical schemes start 9 at a nearby point f(i),then add a sequence of(hopefully)decreasing corrections, as information from other f(i)'s is incorporated.The procedure typically takes O(N2)operations.If everything is well behaved,the last correction will be the smallest,and it can be used as an informal (though not rigorous)bound on the error. 9 In the case of polynomial interpolation,it sometimes does happen that the A之白 coefficients of the interpolating polynomial are of interest,even though their use in evaluating the interpolating function should be frowned on.We deal with this eventuality in $3.5. 6 Local interpolation,using a finite number of"nearest-neighbor"points,gives interpolated values f()that do not,in general,have continuous first or higher derivatives.That happens because,as z crosses the tabulated values zi,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 discontinuity in the interpolated function itself at that point.Bad idea!) Numeric 10621 In situations where continuity of derivatives is a concern,one must use the "stiffer"interpolation provided by a so-called spline function.A spline is a polynomial between each pair of table points,but one whose coefficients are 、三GN determined"slightly"nonlocally.The nonlocality is designed to guarantee global smoothness in the interpolated function up to some order of derivative.Cubic splines (83.3)are the most popular.They produce an interpolated function that is continuous North through the second derivative.Splines tend to be stabler than polynomials,with less possibility of wild oscillation between the tabulated points. The number of points(minus one)used in an interpolation scheme is called the order of the interpolation.Increasing the order does not necessarily increase the accuracy,especially in polynomial interpolation.If the added points are distant from the point of interest z,the resulting higher-order polynomial,with its additional 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,adding points close to the desired point usually does help
106 Chapter 3. Interpolation and Extrapolation Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). which is well-behaved everywhere except at x = π, very mildly singular at x = π, and otherwise takes on all positive and negative values. Any interpolation based on the 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 interpolation and extrapolation routine should provide 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 smoothness for the function interpolated, but within this framework of presumption, deviations from smoothness can be detected. Conceptually, the interpolation process 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 from the N tabulated values every time one is desired. Most practical schemes start at a nearby point f(xi), 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, and it can be used as an informal (though not rigorous) bound on the 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 in evaluating 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 x crosses 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 discontinuity in 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 spline function. 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 global smoothness in the interpolated function up to some order of derivative. Cubic splines (§3.3) are the most popular. They produce an interpolated function that is continuous through the second derivative. Splines tend to be stabler than polynomials, with less possibility of wild oscillation between the tabulated points. The number of points (minus one) used in an interpolation scheme is called the order of the interpolation. Increasing the order does not necessarily increase the accuracy, especially in polynomial interpolation. If the added points are distant from the point of interest x, the resulting higher-order polynomial, with its additional 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, adding points close to the desired point usually does help
3.0 Introduction 107 http://www.nr. read able files .com or call granted for 19881992 11-800-872 (including this one) to any /Cambridge →423 74 uae us e University Press. n NUMERICAL RECIPES IN C: server computer, THE 是 (b) strictly proh Programs 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 by a high-order polynomial (dotted line),which is too "stiff,"than by a low-order polynomial (dashed lines).Even some smooth functions,such as exponentials or rational functions,can to dir be badly approximated by high-order polynomials. 1881892 ART OF SCIENTIFIC COMPUTING(ISBN but a finer mesh implies a larger table of values,not always available. Unless there is solid evidence that the interpolating function is close in form to Numerical 10-621 the true function f,it is a good idea to be cautious about high-order interpolation. We enthusiastically endorse interpolations with 3 or 4 points,we are perhaps tolerant .Further reproduction, Recipes 43108 of 5 or 6:but we rarely go higher than that unless there is quite rigorous monitoring of estimated errors. When your table of values contains many more points than the desirable order (outside of interpolation,you must begin each interpolation with a search for the right"local" 首 Software. place in the table.While not strictly a part of the subject of interpolation,this task is important enough(and often enough botched)that we devote 83.4 to its discussion The routines given for interpolation are also routines for extrapolation.An visit website important application,in Chapter 16,is their use in the integration of ordinary machine differential equations.There,considerable care is taken with the monitoring of errors. Otherwise,the dangers of extrapolation cannot be overemphasized:An interpolating function,which is perforce an extrapolating function,will typically go berserk when the argument z is 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
3.0 Introduction 107 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). (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 by a high-order polynomial (dotted line), which is too “stiff,” than by a low-order polynomial (dashed lines). Even some smooth functions, such as exponentials or rational functions, can be badly approximated by high-order polynomials. but a finer mesh implies a larger table of values, not always available. Unless there is solid evidence that the interpolating function is close in form to the true function f, it is a good idea to be cautious about high-order interpolation. We enthusiastically endorse interpolations with 3 or 4 points, we are perhaps tolerant of 5 or 6; but we rarely go higher than that unless there is quite rigorous monitoring of estimated errors. When your table of values contains many more points than the desirable order of interpolation, you must begin each interpolation with a search for the right “local” place in the table. While not strictly a part of the subject of interpolation, this task is important enough (and often enough botched) that we devote §3.4 to 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 is taken with the monitoring of errors. Otherwise, the dangers of extrapolation cannot be overemphasized: An interpolating function, which is perforce an extrapolating function, will typically go berserk when the argument x is 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
108 Chapter 3.Interpolation and Extrapolation f(,y,z).Multidimensional interpolation is often accomplished by a sequence of one-dimensional interpolations.We discuss this in 83.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 by Dover Publications,New York),$25.2. Stoer,J.,and Bulirsch,R.1980,Introduction to Numerical Analysis(New York:Springer-Verlag). Chapter 2. 三 Acton,F.S.1970,Numerica/Methods That Work;1990,corrected edition (Washington:Mathe- matical Association of America),Chapter 3. Kahaner,D.,Moler,C.,and Nash,S.1989,Numerica/Methods and Software (Englewood Cliffs. 心复a NJ:Prentice Hall),Chapter 4. Johnson,L.W.,and Riess,R.D.1982,Numerical Analysis,2nd ed.(Reading,MA:Addison- ICAL Wesley).Chapter 5. Ralston,A.,and Rabinowitz,P.1978,A First Course in Numerical Analysis,2nd ed.(New York: McGraw-Hill),Chapter 3. RECIPES Isaacson,E.,and Keller,H.B.1966,Analysis of Numerica/Methods(New York:Wiley).Chapter 6. 9 Ea 9 3.1 Polynomial Interpolation and Extrapolation 9 Through any two points there is a unique line.Through any three points,a IENTIFIC unique quadratic.Et cetera.The interpolating polynomial of degree N-1 through the N points y1 f(1),y2 =f(2),...,yN f(N)is given explicitly by Lagrange's classical formula, (x-x2)(x-t3)(x-xN) -21--h+2c-- (22-1(2-3)(2-xN (x-x1)(x-x2)(x-xN-1) 2.5兰 Numerica 10.621 十…·十 (EN-1)(EN -)..(N -EN-1)UN (3.1.1) There are N terms,each a polynomial of degree N-1 and each constructed to be (outside zero at all of the i except one,at which it is constructed to be yi. It is not terribly wrong to implement the Lagrange formula straightforwardly, North but it is not terribly right either.The resulting algorithm gives no error estimate,and it is also somewhat awkward to program.A much better algorithm(for constructing the same,unique,interpolating polynomial)is Neville's algorithm,closely related to and sometimes confused with Aitken's algorithm,the latter now considered obsolete. Let P be the value at z of the unique polynomial of degree zero (i.e., a constant)passing through the point (x1,y1);so P=y1.Likewise define P2,P3,...,P.Now let P12 be the value at z of the unique polynomial of degree one passing through both (1,1)and (z2,y2).Likewise P23,P34,..., PN-1)N.Similarly,for higher-order polynomials,up to P123...N,which is the value of the unique interpolating polynomial through all N points,i.e.,the desired answer
108 Chapter 3. Interpolation and Extrapolation Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover Publications, New York), §25.2. Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), Chapter 2. Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), Chapter 3. Kahaner, D., Moler, C., and Nash, S. 1989, Numerical Methods and Software (Englewood Cliffs, NJ: Prentice Hall), Chapter 4. Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis, 2nd ed. (Reading, MA: AddisonWesley), Chapter 5. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed. (New York: McGraw-Hill), Chapter 3. Isaacson, E., and Keller, H.B. 1966, Analysis of Numerical Methods (New York: Wiley), Chapter 6. 3.1 Polynomial Interpolation and Extrapolation Through any two points there is a unique line. Through any three points, a unique quadratic. Et cetera. The interpolating polynomial of degree N − 1 through the N points y1 = f(x1), y2 = f(x2),...,yN = 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 N terms, each a polynomial of degree N − 1 and each constructed to be zero at all of the xi except one, at which it is constructed to be yi. It is not terribly wrong to implement the Lagrange formula straightforwardly, but it is not terribly right either. The resulting algorithm gives no error estimate, and it is also somewhat awkward to program. A much better algorithm (for constructing the same, unique, interpolating polynomial) is Neville’s algorithm, closely related to and sometimes confused with Aitken’s algorithm, the latter now considered obsolete. Let P1 be the value at x of the unique polynomial of degree zero (i.e., a constant) passing through the point (x1, y1); so P1 = y1. Likewise define P2, P3,...,PN . Now let P12 be the value at x of the unique polynomial of degree one passing through both (x1, y1) and (x2, y2). Likewise P23, P34,..., P(N−1)N . Similarly, for higher-order polynomials, up to P123...N , which is the value of the unique interpolating polynomial through all N points, i.e., the desired answer