Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 31
Текст из файла (страница 31)
1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover 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 ExtrapolationThrough any two points there is a unique line.
Through any three points, aunique quadratic. Et cetera. The interpolating polynomial of degree N − 1 throughthe N points y1 = f(x1 ), y2 = f(x2 ), . . . , yN = f(xN ) is given explicitly byLagrange’s classical formula,(x − x2 )(x − x3 )...(x − xN )(x − x1 )(x − x3 )...(x − xN )y1 +y2(x1 − x2 )(x1 − x3 )...(x1 − xN )(x2 − x1 )(x2 − x3 )...(x2 − xN )(x − x1 )(x − x2 )...(x − xN−1 )yN+···+(xN − x1 )(xN − x2 )...(xN − xN−1 )(3.1.1)There are N terms, each a polynomial of degree N − 1 and each constructed to bezero 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, andit is also somewhat awkward to program. A much better algorithm (for constructingthe same, unique, interpolating polynomial) is Neville’s algorithm, closely related toand 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 defineP2 , P3 , . . . , PN .
Now let P12 be the value at x of the unique polynomial ofdegree 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 valueof the unique interpolating polynomial through all N points, i.e., the desired answer.P (x) =3.1 Polynomial Interpolation and Extrapolation109The 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 = P1x2 :y2 = P2P12P123P23x3 :y3 = P3x4 :y4 = P4P1234(3.1.2)P234P34Neville’s algorithm is a recursive way of filling in the numbers in the tableaua column at a time, from left to right.
It is based on the relationship between a“daughter” P and its two “parents,”(x − xi+m )Pi(i+1)...(i+m−1) + (xi − x)P(i+1)(i+2)...(i+m)xi − xi+m(3.1.3)Pi(i+1)...(i+m) =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 smalldifferences 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(xi+m+1 − x)(Cm,i+1 − Dm,i )xi − xi+m+1(xi − x)(Cm,i+1 − Dm,i )=xi − xi+m+1Dm+1,i =Cm+1,i(3.1.5)At each level m, the C’s and D’s are the corrections that make the interpolation oneorder higher. The final answer P1...N is equal to the sum of any yi plus a set of C’sand/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 inputpoints.
Note that the input arrays are assumed to be unit-offset. If you havezero-offset arrays, remember to subtract 1 (see §1.2):#include <math.h>#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, andan 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;110Chapter 3.Interpolation and Extrapolationfloat *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 updateho=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 thetableau—forking up or down.
We do this in such a way as to take the most “straightline” route through the tableau to its apex, updating ns accordingly to keep track ofwhere 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 xaand ya replaced by actual arrays with offsets.
For example, the constructionpolint(&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 byDover 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 (EnglewoodCliffs, NJ: Prentice-Hall), §6.1.3.2 Rational Function Interpolation and Extrapolation1113.2 Rational Function Interpolation andExtrapolationSome functions are not well approximated by polynomials, but are wellapproximated 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, supposeRi(i+1)...(i+m) =p 0 + p1 x + · · · + pµ x µPµ (x)=Qν (x)q 0 + q 1 x + · · · + q ν xν(3.2.1)Since there are µ + ν + 1 unknown p’s and q’s (q0 being arbitrary), we must havem+1 = µ+ν +1(3.2.2)In specifying a rational function interpolating function, you must give the desiredorder 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 denominatorof equation (3.2.1). These poles might occur for real values of x, if the functionto be interpolated itself has poles.
More often, the function f(x) is finite for allfinite 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 toreal values of x, just as they can ruin the convergence of an infinite power seriesin 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 nearestpole 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 accountfor (cancel) any nearby poles.For the interpolation problem, a rational function is constructed so as to gothrough a chosen set of tabulated functional values. However, we should alsomention in passing that rational function approximations can be used in analyticwork. One sometimes constructs a rational function approximation by the criterionthat the rational function of equation (3.2.1) itself have a power series expansionthat agrees with the first m + 1 terms of the power series expansion of the desiredfunction f(x). This is called P adé approximation, and is discussed in §5.12.Bulirsch and Stoer found an algorithm of the Neville type which performsrational 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, withthe degrees of numerator and denominator equal (if m is even) or with the degreeof the denominator larger by one (if m is odd, cf. equation 3.2.2 above). For thederivation of the algorithm, refer to [1]. The algorithm is summarized by a recurrence112Chapter 3.Interpolation and Extrapolationrelation exactly analogous to equation (3.1.3) for polynomial approximation:Ri(i+1)...(i+m) = R(i+1)...(i+m)+#R(i+1)...(i+m) − Ri...(i+m−1)$#$R(i+1)...(i+m) −Ri...(i+m−1)x−xi1−−1x−xi+mR−R(i+1)...(i+m)(i+1)...(i+m−1)(3.2.3)This recurrence generates the rational functions through m + 1 points from the onesthrough m and (the term R(i+1)...(i+m−1) in equation 3.2.3) m−1 points. It is startedwith(3.2.4)Ri = yiand with(3.2.5)R ≡ [Ri(i+1)...(i+m) with m = −1] = 0Now, exactly as in equations (3.1.4) and (3.1.5) above, we can convert therecurrence (3.2.3) to one involving only the small differencesCm,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 relationCm+1,i − Dm+1,i = Cm,i+1 − Dm,i(3.2.7)which is useful in proving the recurrencesDm+1,i = ##Cm+1,i =Cm,i+1 (Cm,i+1 − Dm,i )$x−xix−xi+m+1 Dm,i − Cm,i+1x−xix−xi+m+1$Dm,i (Cm,i+1 − Dm,i )#$x−xix−xi+m+1 Dm,i − Cm,i+1(3.2.8)This recurrence is implemented in the following function, whose use is analogousin every way to polint in §3.1.
Note again that unit-offset input arrays areassumed (§1.2).#include <math.h>#include "nrutil.h"#define TINY 1.0e-25A small number.#define FREERETURN {free_vector(d,1,n);free_vector(c,1,n);return;}void ratint(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 of x, this routine returns a value ofy and an accuracy estimate dy. The value returned is that of the diagonal rational function,evaluated at x, which passes through the n points (xai , yai ), i = 1...n.{int m,i,ns=1;float w,t,hh,h,dd,*c,*d;3.3 Cubic Spline Interpolation113c=vector(1,n);d=vector(1,n);hh=fabs(x-xa[1]);for (i=1;i<=n;i++) {h=fabs(x-xa[i]);if (h == 0.0) {*y=ya[i];*dy=0.0;FREERETURN} else if (h < hh) {ns=i;hh=h;}c[i]=ya[i];d[i]=ya[i]+TINY;The TINY part is needed to prevent a rare zero-over-zero}condition.*y=ya[ns--];for (m=1;m<n;m++) {for (i=1;i<=n-m;i++) {w=c[i+1]-d[i];h=xa[i+m]-x;h will never be zero, since this was tested in the initialt=(xa[i]-x)*d[i]/h;izing loop.dd=t-c[i+1];if (dd == 0.0) nrerror("Error in routine ratint");This error condition indicates that the interpolating function has a pole at therequested value of x.dd=w/dd;d[i]=c[i+1]*dd;c[i]=t*dd;}*y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));}FREERETURN}CITED REFERENCES AND FURTHER READING:Stoer, J., and Bulirsch, R.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
