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
3.1 Polynomial Interpolation and Extrapolation 109 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, E1: h=乃 P12 T2: 2=P2 P123 P23 P1234 (3.1.2) T3: 3=B P234 P34 E4: 4=P Neville's algorithm is a recursive way of filling in the numbers in the tableau 83 a column at a time,from left to right.It is based on the relationship between a “daughter'”P and its two“parents,” 容 P(+1).…(i+m)= (x-xi+m)Pi+1)…(i+m-1)+(zi-t)P+1)t+2.6+m) xi-工i+m RECIPES (3.1.3) This recurrence works because the two parents already agree at points +1... Ti+m-1. 二X分 Press. 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,..., 9 N-1)2 Progra Cm.i=Pi...(i+m)-Pi..(i+m-1) Dm,i≡P.(i+m)-P(i+1).(i+m小 (3.1.4) CIENTIFIC Then one can easily derive from(3.1.3)the relations Dm+li =itm+1-)(Cm.it1-Dm.) 工-工i+m+1 9 Cm+=)(Cmit1-Dm.) (3.1.5) Ti-工i+m+1 Numerical 10521 At each level m.the C's and D's are the corrections that make the interpolation one 43106 order higher.The final answer P1...N is equal to the sum of any yi plus a set of C's (outside Recipes and/or D's that form a path through the family tree to the rightmost daughter. Here is a routine for polynomial interpolation or extrapolation from N input points.Note that the input arrays are assumed to be unit-offset.If you have North Software. zero-offset arrays,remember to subtract 1 (see $1.2): #include #include "nrutil.h" void polint(float xa[],float ya[],int n,float x,float y,float *dy) Given arrays xa[1..n]and ya[1..n],and given a value x,this routine returns a value y,and an error estimate dy.If P(x)is the polynomial of degree N-1 such that P(xai)=ya;,i= 1,...,n,then the returned value y P(x). int i,m,ns=1; float den,dif,dift,ho,hp,w;
3.1 Polynomial Interpolation and Extrapolation 109 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). 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” P and 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 each level m, the C’s and D’s are the corrections that make the interpolation one order higher. The final answer P1...N is equal to the sum of any yi plus a set of C’s and/or D’s that form a path through the family tree to the rightmost daughter. Here is a routine for polynomial interpolation or extrapolation from N input points. Note that the input arrays are assumed to be unit-offset. If you have zero-offset arrays, remember to subtract 1 (see §1.2): #include #include "nrutil.h" void polint(float xa[], float ya[], int n, float x, float *y, float *dy) Given arrays xa[1..n] and ya[1..n], and given a value x, this routine returns a value y, and an error estimate dy. If P(x) is the polynomial of degree N − 1 such that P(xai) = yai, i = 1,..., n, then the returned value y = P(x). { int i,m,ns=1; float den,dif,dift,ho,hp,w;
110 Chapter 3.Interpolation and Extrapolation float *c,*d; dif=fabs(x-xa[1]); c=vector(1,n); d=vector(1,n); for(1=1;1<=n;1++){ Here we find the index ns of the closest table entry, if (dift=fabs(x-xa[i]))<dif){ ns=1; dif-dift; c[i]-ya[i]; and initialize the tableau of c's and d's. d[i]=ya[i]; *y=ya[ns--]; This is the initial approximation to y. for (m=1;m<n;m++){ For each column of the tableau, 83g for(i=1;i<=n-m;1+)[ we loop over the current c's and d's and update ho=xa[i]-x; them. hp=xa[i+m]-x; w=c[i+1]-d[1]: if (den=ho-hp)==0.0)nrerror("Error in routine 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; *y+=(*dy=(2*ns<(n-m)？c[ns+1]:d[ns--]); After each column in the tableau is completed,we decide which correction,c or d, America Press. we want to add to our accumulating 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 ns accordingly to keep track of where we are.This route keeps the partial approximations centered (insofar as possible) 9 Programs on the target x.The last dy added is thus the error indication. free_vector(d,1,n); SCIENTIFIC free_vector(c,1,n); to dir Quite often you will want to call polint with the dummy arguments xa and ya replaced by actual arrays with offsets.For example,the construction 192 COMPUTING (ISBN polint (&xx[14],&yy [14],4,x,y,dy)performs 4-point interpolation on the tab- ulated values xx[15..18],yy [15..18].For more on this,see the end of $3.4. Fuurggoglrion Numerical Recipes 10621 43106 CITED REFERENCES AND FURTHER READING: (outside 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 North Software. Dover Publications,New York),$25.2. Stoer.J.,and Bulirsch,R.1980,Introduction to Numerical Analysis(New York:Springer-Verlag). 62.1. Gear,C.W.1971,Numerical Initial Value Problems in Ordinary Differential Equations(Englewood Cliffs,NJ:Prentice-Hall),$6.1
110 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). float *c,*d; dif=fabs(x-xa[1]); c=vector(1,n); d=vector(1,n); for (i=1;i<=n;i++) { Here we find the index ns of the closest table entry, if ( (dift=fabs(x-xa[i])) < dif) { ns=i; dif=dift; } c[i]=ya[i]; and initialize the tableau of c’s and d’s. d[i]=ya[i]; } *y=ya[ns--]; This is the initial approximation to y. for (m=1;m<n;m++) { For each column of the tableau, for (i=1;i<=n-m;i++) { we loop over the current c’s and d’s and update ho=xa[i]-x; them. hp=xa[i+m]-x; w=c[i+1]-d[i]; if ( (den=ho-hp) == 0.0) nrerror("Error in routine 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; } *y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--])); After each column in the tableau is completed, we decide which correction, c or d, we want to add to our accumulating 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 ns accordingly to keep track of where we are. This route keeps the partial approximations centered (insofar as possible) on the target x. The last dy added is thus the error indication. } free_vector(d,1,n); free_vector(c,1,n); } Quite often you will want to call polint with the dummy arguments xa and ya replaced by actual arrays with offsets. For example, the construction polint(&xx[14],&yy[14],4,x,y,dy) performs 4-point interpolation on the tabulated values xx[15..18], yy[15..18]. For more on this, see the end of §3.4. 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), §2.1. Gear, C.W. 1971, Numerical Initial Value Problems in Ordinary Differential Equations (Englewood Cliffs, NJ: Prentice-Hall), §6.1
3.2 Rational Function Interpolation and Extrapolation 111 3.2 Rational Function Interpolation and Extrapolation Some functions are not well approximated by polynomials,but are well 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+1 points (i,i)...(i+m,yi+m).More explicitly,suppose P(x）_P0+p1E+…+Pu 三 R4(+1).-(+m)=Q(西=q0+g1x+…+9 (3.2.1) Since there are u+v+1 unknown p's and g's(go being arbitrary),we must have 豆君蜀州 ICAL m+1=4+y+1 (3.2.2) In specifying a rational function interpolating function,you must give the desired RECIPES order of both the numerator and the denominator. Rational functions are sometimes superior to polynomials,roughly speaking, 9 because oftheir ability to model functions with poles,that is,zeros ofthe denominator of 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()is finite for all finite real z,but has an analytic continuation with poles in the complex z-plane. 9 Such poles can themselves ruin a polynomial approximation,even one restricted to 、经是g∽ real values of z,just as they can ruin the convergence of an infinite power series in z.If you draw a circle in the complex plane around your m tabulated points, then you should not expect polynomial interpolation to be good unless the nearest 61 pole is rather far outside the circle.A rational function approximation,by contrast, will stay"good"as long as it has enough powers ofz in its denominator to 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 SPod毫 Numerica 10621 work.One sometimes constructs a rational function approximation by the criterion that the rational function of equation(3.2.1)itself have a power series expansion that agrees with the first m+1 terms of the power series expansion of the desired function f(x).This is called Pade approximation,and is discussed in $5.12. Bulirsch and Stoer found an algorithm of the Neville type which performs North 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. The Bulirsch-Stoer algorithm produces the so-called diagonal rational function,with the degrees of numerator and denominator equal (if m is even)or with the degree of the denominator larger by one (if m is odd,cf.equation 3.2.2 above).For the derivation of the algorithm,refer to [1].The algorithm is summarized by a recurrence
3.2 Rational Function Interpolation and Extrapolation 111 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). 3.2 Rational Function Interpolation and Extrapolation Some functions are not well approximated by polynomials, but are well approximated by rational functions, that is quotients of polynomials. We denote by Ri(i+1)...(i+m) a rational function passing through the m + 1 points (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 there are µ + ν + 1 unknown p’s and q’s (q0 being arbitrary), we must have 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, because of their ability to model functions with poles, that is, zeros of the denominator of 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 finite real x, 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 in x. If you draw a circle in the complex plane around your m tabulated 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 enough powers of x in its denominator to 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 approximation by the criterion that the rational function of equation (3.2.1) itself have a power series expansion that agrees with the first m + 1 terms of the power series expansion of the desired function f(x). This is called P ade´ approximation, 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. The Bulirsch-Stoer algorithm produces the so-called diagonal rational function, with the degrees of numerator and denominator equal (if m is even) or with the degree of the denominator larger by one (if m is odd, cf. equation 3.2.2 above). For the derivation of the algorithm, refer to [1]. The algorithm is summarized by a recurrence