《数字信号处理》教学参考资料（Numerical Recipes in C，The Art of Scientific Computing Second Edition）Chapter 03.2 Interpolation and Extrapolation 3.2 Rational Function Interpolation and Extrapolation

团购合买资源类别：文库，文档格式：PDF，文档页数：3，文件大小：52.17KB

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

112 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). relation 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) x−xi x−xi+m 1 − R(i+1)...(i+m)−Ri...(i+m−1) R(i+1)...(i+m)−R(i+1)...(i+m−1) − 1 (3.2.3) This recurrence generates the rational functions through m + 1 points from the ones throughm and (the term R(i+1)...(i+m−1) in equation 3.2.3) m−1 points. It is started with Ri = yi (3.2.4) and with R ≡ [Ri(i+1)...(i+m) with m = −1] = 0 (3.2.5) Now, exactly as in equations (3.1.4) and (3.1.5) above, we can convert the recurrence (3.2.3) to one involving only the small differences Cm,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 relation Cm+1,i − Dm+1,i = Cm,i+1 − Dm,i (3.2.7) which is useful in proving the recurrences Dm+1,i = Cm,i+1(Cm,i+1 − Dm,i) x−xi x−xi+m+1 Dm,i − Cm,i+1 Cm+1,i = x−xi x−xi+m+1 Dm,i(Cm,i+1 − Dm,i) x−xi x−xi+m+1 Dm,i − Cm,i+1 (3.2.8) This recurrence is implemented in the following function, whose use is analogous in every way to polint in §3.1. Note again that unit-offset input arrays are assumed (§1.2). #include #include "nrutil.h" #define TINY 1.0e-25 A 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 of y 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 Interpolation 113 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). c=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 the requested 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. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), §2.2. [1] Gear, C.W. 1971, Numerical Initial Value Problems in Ordinary Differential Equations (Englewood Cliffs, NJ: Prentice-Hall), §6.2. Cuyt, A., and Wuytack, L. 1987, Nonlinear Methods in Numerical Analysis (Amsterdam: NorthHolland), Chapter 3. 3.3 Cubic Spline Interpolation Given a tabulated function yi = y(xi), i = 1...N, focus attention on one particular interval, between xj and xj+1. Linear interpolation in that interval gives the interpolation formula y = Ayj + Byj+1 (3.3.1)

点击进入文档下载页（PDF格式）

已到末页，全文结束

点击下载（PDF格式）

浏览记录