c3-1 (779471)
Текст из файла
108Chapter 3.Interpolation and Extrapolationf(x, y, z). Multidimensional interpolation is often accomplished by a sequence ofone-dimensional interpolations. We discuss this in §3.6.CITED REFERENCES AND FURTHER READING: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) =Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).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.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 = P2P12P123x3 :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;Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).P23110Chapter 3.Interpolation and Extrapolationfloat *c,*d;}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.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).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.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.