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

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

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)

114 Chapter 3.Interpolation and Extrapolation where A三巧+1~x B三1-A=r- (3.3.2) Tj+1-工j Tj+1-Tj Equations(3.3.1)and(3.3.2)are a special case of the general Lagrange interpolation formula (3.1.1). Since it is (piecewise)linear,equation(3.3.1)has zero second derivative in the interior of each interval,and an undefined,or infinite,second derivative at the abscissas ri.The goal of cubic spline interpolation is to get an interpolation formula that is smooth in the first derivative,and continuous in the second derivative,both within an interval and at its boundaries. Suppose,contrary to fact,that in addition to the tabulated values of y,we also have tabulated values for the function's second derivatives,y",that is,a set 黑生总馆 81 鱼君 of numbers y.Then,within each interval,we can add to the right-hand side of equation (3.3.1)a cubic polynomial whose second derivative varies linearly from a valueon the left to a value on the right.Doing so,we will have the desired RECIPESI continuous second derivative.If we also construct the cubic polynomial to have 2 zero values at and j1,then adding it in will not spoil the agreement with the tabulated functional values yj and y+at the endpoints j and j+1 A little side calculation shows that there is only one way to arrange this construction,namely replacing (3.3.1)by y=A+B+1+C+D+1 (3.3.3) a。会 where A and B are defined in (3.3.2)and C=4-4G+1-PD=B-BG+1-马P (3.3.4) Notice that the dependence on the independent variable in equations(3.3.3)and (3.3.4)is entirely through the linear z-dependence of A and B,and(through A and B)the cubic z-dependence of C and D. We can readily check that y"is in fact the second derivative of the new 度 Numerica 10-521 interpolating polynomial.We take derivatives of equation(3.3.3)with respect to z, using the definitions of A,B,C,D to compute dA/dr,dB/dr,dC/dr,and dD/dx. The result is A =5+1-当-3A2-1 dx tj+1-xj 6-a*1-马财+3B2-1 6—(+1-)g+163.5) for the first derivative,and =A+B Py (3.3.6) for the second derivative.Since A=1 at j,A=0 at j+1,while B is just the other way around,(3.3.6)shows that y"is just the tabulated second derivative,and also that the second derivative will be continuous across(e.g.)the boundary between the two intervals (xj-1,j)and (j,j+1)

114 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). where A ≡ xj+1 − x xj+1 − xj B ≡ 1 − A = x − xj xj+1 − xj (3.3.2) Equations (3.3.1) and (3.3.2) are a special case of the general Lagrange interpolation formula (3.1.1). Since it is (piecewise) linear, equation (3.3.1) has zero second derivative in the interior of each interval, and an undefined, or infinite, second derivative at the abscissas xj . The goal of cubic spline interpolation is to get an interpolation formula that is smooth in the first derivative, and continuous in the second derivative, both within an interval and at its boundaries. Suppose, contrary to fact, that in addition to the tabulated values of y i, we also have tabulated values for the function’s second derivatives, y , that is, a set of numbers y i . Then, within each interval, we can add to the right-hand side of equation (3.3.1) a cubic polynomial whose second derivative varies linearly from a value y j on the left to a value y j+1 on the right. Doing so, we will have the desired continuous second derivative. If we also construct the cubic polynomial to have zero values at xj and xj+1, then adding it in will not spoil the agreement with the tabulated functional values yj and yj+1 at the endpoints xj and xj+1. A little side calculation shows that there is only one way to arrange this construction, namely replacing (3.3.1) by y = Ayj + Byj+1 + Cy j + Dy j+1 (3.3.3) where A and B are defined in (3.3.2) and C ≡ 1 6 (A3 − A)(xj+1 − xj ) 2 D ≡ 1 6 (B3 − B)(xj+1 − xj ) 2 (3.3.4) Notice that the dependence on the independent variable x in equations (3.3.3) and (3.3.4) is entirely through the linear x-dependence of A and B, and (through A and B) the cubic x-dependence of C and D. We can readily check that y is in fact the second derivative of the new interpolating polynomial. We take derivatives of equation (3.3.3) with respect to x, using the definitions of A, B, C, D to compute dA/dx, dB/dx, dC/dx, and dD/dx. The result is dy dx = yj+1 − yj xj+1 − xj − 3A2 − 1 6 (xj+1 − xj )y j + 3B2 − 1 6 (xj+1 − xj )y j+1 (3.3.5) for the first derivative, and d2y dx2 = Ay j + By j+1 (3.3.6) for the second derivative. Since A = 1 at xj , A = 0 at xj+1, while B is just the other way around, (3.3.6) shows that y is just the tabulated second derivative, and also that the second derivative will be continuous across (e.g.) the boundary between the two intervals (xj−1, xj ) and (xj , xj+1)

3.3 Cubic Spline Interpolation 115 The only problem now is that we supposed the y's to be known,when,actually, they are not.However,we have not yet required that the first derivative,computed from equation(3.3.5),be continuous across the boundary between two intervals.The key idea of a cubic spline is to require this continuity and to use it to get equations for the second derivatives y. The required equations are obtained by setting equation(3.3.5)evaluated for x=z;in the interval (j1,j)equal to the same equation evaluated forx=x;but in the interval (j,j+1).With some rearrangement,this gives(for j=2,...,N-1) 。1+1二+。五=二路-为二 81 6 3 x1+1-Ejxj-工j-1 (3.3.7) These are N-2 linear equations in the N unknowns y,i=1,...,N.Therefore ICAL there is a two-parameter family of possible solutions. For a unique solution,we need to specify two further conditions,typically taken as boundary conditions at z1 andrN.The most common ways of doing this are either RECIPES set one or both of y"and y equal to zero,giving the so-called natural 9 cubic spline,which has zero second derivative on one or both of its boundaries,or 器三。学州 set either of y and y to values calculated from equation (3.3.5)so as to make the first derivative of the interpolating function have a specified value on either or both boundaries. One reason that cubic splines are especially practical is that the set ofequations (3.3.7),along with the two additional boundary conditions,are not only linear,but also tridiagonal.Eachy is coupled only to its nearest neighbors atj+1.Therefore, 6 the equations can be solved in O(N)operations by the tridiagonal algorithm(82.4). That algorithm is concise enough to build right into the spline calculational routine This makes the routine not completely transparent as an implementation of(3.3.7), so we encourage you to study it carefully,comparing with tridag(S2.4).Arrays are assumed to be unit-offset.If you have zero-offset arrays,see $1.2. N传是 Numerica 10621 #include "nrutil.h" 4310 void spline(float x[],float y[],int n,float ypi,float ypn,float y2[]) Recipes Given arrays x[1..n]and y[1..n]containing a tabulated function,i.e.,y;f(xi),with x0.99e30) The lower boundary condition is set either to be "nat- y2[1]=u[1]=0.0; ural" else or else to have a specified first derivative y2[1]=-0.5; u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);

3.3 Cubic Spline Interpolation 115 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 only problem now is that we supposed the y i ’s to be known, when, actually, they are not. However, we have not yet required that the first derivative, computed from equation (3.3.5), be continuous across the boundary between two intervals. The key idea of a cubic spline is to require this continuity and to use it to get equations for the second derivatives y i . The required equations are obtained by setting equation (3.3.5) evaluated for x = xj in the interval(xj−1, xj ) equal to the same equation evaluated for x = xj but in the interval(xj , xj+1). With some rearrangement, this gives (for j = 2,...,N −1) xj − xj−1 6 y j−1 + xj+1 − xj−1 3 y j + xj+1 − xj 6 y j+1 = yj+1 − yj xj+1 − xj − yj − yj−1 xj − xj−1 (3.3.7) These are N − 2 linear equations in the N unknowns y i , i = 1,...,N. Therefore there is a two-parameter family of possible solutions. For a unique solution, we need to specify two further conditions, typically taken as boundary conditions at x1 and xN . The most common ways of doing this are either • set one or both of y 1 and y N equal to zero, giving the so-called natural cubic spline, which has zero second derivative on one or both of its boundaries, or • set either of y 1 and y N to values calculated from equation (3.3.5) so as to make the first derivative of the interpolating function have a specified value on either or both boundaries. One reason that cubic splines are especially practical is that the set of equations (3.3.7), along with the two additional boundary conditions, are not only linear, but also tridiagonal. Each y j is coupled only to its nearest neighbors at j ±1. Therefore, the equations can be solved in O(N) operations by the tridiagonal algorithm (§2.4). That algorithm is concise enough to build right into the spline calculational routine. This makes the routine not completely transparent as an implementation of (3.3.7), so we encourage you to study it carefully, comparing with tridag (§2.4). Arrays are assumed to be unit-offset. If you have zero-offset arrays, see §1.2. #include "nrutil.h" void spline(float x[], float y[], int n, float yp1, float ypn, float y2[]) Given arrays x[1..n] and y[1..n] containing a tabulated function, i.e., yi = f(xi), with x1 0.99e30) The lower boundary condition is set either to be “naty2[1]=u[1]=0.0; ural” else { or else to have a specified first derivative. y2[1] = -0.5; u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1); }

116 Chapter 3. Interpolation and Extrapolation for(i=2;10.99e30) The upper boundary condition is set either to be qn=un=0.0; “natura" else or else to have a specified first derivative. qn=0.5; n=(3.0/(x[n]-x[n-1])*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]); y2[n]=(un-qn*u[n-1])/(gn*y2[n-1]+1.0); for(k=n-1;k>=1;k--) This is the backsubstitution loop of the tridiagonal y2[]=y2[k]*y2[k+1]+u[k]; algorithm. free_vector(u,1,n-1); It is important to understand that the program spline is called only once to 。 nted for process an entire tabulated function in arrays xi and yi.Once this has been done, values of the interpolated function for any value of x are obtained by calls (as many as desired)to a separate routine splint (for"spline interpolation"): void splint(float xa[],float ya[],float y2a[],int n,float x,float *y) Press. Given the arrays xa[1..n]and ya[1..n],which tabulate a function (with the xai's in order), and given the array y2a[1..n],which is the output from spline above,and given a value of x,this routine returns a cubic-spline interpolated value y. Programs void nrerror(char error._text☐)； int klo,khi,k; SCIENTIFIC float h,b,a; k1o=1; We will find the right place in the table by means of to dir khisn; bisection.This is optimal if sequential calls to this while (khi-klo 1){ routine are at random values of x.If sequential calls k=(kh1+k1o)>>1; are in order,and closely spaced,one would do better if (xa[k]x)khi=k; to store previous values of klo and khi and test if 1920 COMPUTING (ISBN else klo=k; they remain appropriate on the next call. klo and khi now bracket the input value of x. h=xa[khi]-xa[klo]; 10621 if (h ==0.0)nrerror("Bad xa input to routine splint"); The xa's must be dis- a=(xa[khi]-x)/h; tinct uction Numerical Recipes 43108 b=(x-xa[klo])/h; Cubic spline polynomial is now evaluated. *y=a*ya [klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0; (outside North Software. CITED REFERENCES AND FURTHER READING: De Boor,C.1978,A Practical Guide to Splines (New York:Springer-Verlag). Forsythe,G.E.,Malcolm,M.A.,and Moler,C.B.1977,Computer Methods for Mathematical Computations (Englewood Cliffs,NJ:Prentice-Hall),84.4-4.5. Stoer.J.,and Bulirsch.R.1980.Introduction to Numerical Analysis(New York:Springer-Verlag). 82.4. Ralston,A.,and Rabinowitz.P.1978,A First Course in Numerical Analysis,2nd ed.(New York: McGraw-Hill),83.8

116 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). for (i=2;i 0.99e30) The upper boundary condition is set either to be qn=un=0.0; “natural” else { or else to have a specified first derivative. qn=0.5; un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1])); } y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0); for (k=n-1;k>=1;k--) This is the backsubstitution loop of the tridiagonal y2[k]=y2[k]*y2[k+1]+u[k]; algorithm. free_vector(u,1,n-1); } It is important to understand that the program spline is called only once to process an entire tabulated function in arrays xi and yi. Once this has been done, values of the interpolated function for any value of x are obtained by calls (as many as desired) to a separate routine splint (for “spline interpolation”): void splint(float xa[], float ya[], float y2a[], int n, float x, float *y) Given the arrays xa[1..n] and ya[1..n], which tabulate a function (with the xai’s in order), and given the array y2a[1..n], which is the output from spline above, and given a value of x, this routine returns a cubic-spline interpolated value y. { void nrerror(char error_text[]); int klo,khi,k; float h,b,a; klo=1; We will find the right place in the table by means of bisection. This is optimal if sequential calls to this routine are at random values of x. If sequential calls are in order, and closely spaced, one would do better to store previous values of klo and khi and test if they remain appropriate on the next call. khi=n; while (khi-klo > 1) { k=(khi+klo) >> 1; if (xa[k] > x) khi=k; else klo=k; } klo and khi now bracket the input value of x. h=xa[khi]-xa[klo]; if (h == 0.0) nrerror("Bad xa input to routine splint"); The xa’s must be disa=(xa[khi]-x)/h; tinct. b=(x-xa[klo])/h; Cubic spline polynomial is now evaluated. *y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0; } CITED REFERENCES AND FURTHER READING: De Boor, C. 1978, A Practical Guide to Splines (New York: Springer-Verlag). Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), §§4.4–4.5. Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), §2.4. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed. (New York: McGraw-Hill), §3.8

3.4 How to Search an Ordered Table 117 3.4 How to Search an Ordered Table Suppose that you have decided to use some particular interpolation scheme, such as fourth-order polynomial interpolation,to compute a function f(z)from a set of tabulated x;'s and f;'s.Then you will need a fast way of finding your place in the table of zi's,given some particular value z at which the function evaluation is desired.This problem is not properly one of numerical analysis,but it occurs so often in practice that it would be negligent of us to ignore it. Formally,the problem is this:Given an array ofabscissas xx[j],j=1,2,...n, with the elements either monotonically increasing or monotonically decreasing,and given a number x,find an integer j such that x lies between xx[j]and xx [j+1]. For this task,let us define fictitious array elements xx[o]and xx[n+1]equal to plus or minus infinity (in whichever order is consistent with the monotonicity of the 含 table).Then j will always be between 0 and n,inclusive;a value of 0 indicates 'off-scale"at one end of the table.n indicates off-scale at the other end. In most cases,when all is said and done,it is hard to do better than bisection, RECIPES I which will find the right place in the table in about log2n tries.We already did use 2 bisection in the spline evaluation routine splint of the preceding section,so you might glance back at that.Standing by itself,a bisection routine looks like this: void locate(float xx[],unsigned long n,float x,unsigned long *j) Given an array xx[1..n],and given a value x,returns a value j such that x is between xx[j] and xx[j+1].xx must be monotonic,either increasing or decreasing.j=0 or j=n is returned to indicate that x is out of range. 、9 OF SCIENTIFIC unsigned long ju,jm,jl; int ascnd; to dir j1=0; Initialize lower ju=n+1; and upper limits. COMPUTING ascnd=(xx[n]>xx[1]); h11e(ju-j1>1)( If we are not yet done, jm=(ju+j1)>1; compute a midpoint, (ISBN if (x >xx[im]=ascnd) jl=jm; and replace either the lower limit v@cambri 10621 else ju=jm; or the upper limit,as appropriate Repeat until the test condition is satisfied. 1988-1992 by Numerical Recipes 43108 1f(x=xx[1])*j=1; Then set the output else if(x =xx[n])*j=n-1; else *j=jl; and return. (outside North Software. A unit-offset array xx is assumed.To use locate with a zero-offset array, Ame ying of remember to subtract 1 from the address of xx,and also from the returned value j. Search with Correlated Values Sometimes you will be in the situation of searching a large table many times, and with nearly identical abscissas on consecutive searches.For example,you may be generating a function that is used on the right-hand side of a differential equation:Most differential-equation integrators,as we shall see in Chapter 16,call

3.4 How to Search an Ordered Table 117 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.4 How to Search an Ordered Table Suppose that you have decided to use some particular interpolation scheme, such as fourth-order polynomial interpolation, to compute a function f(x) from a set of tabulated xi’s and fi’s. Then you will need a fast way of finding your place in the table of xi’s, given some particular value x at which the function evaluation is desired. This problem is not properly one of numerical analysis, but it occurs so often in practice that it would be negligent of us to ignore it. Formally, the problem is this: Given an array of abscissas xx[j], j=1, 2, ...,n, with the elements either monotonically increasing or monotonically decreasing, and given a number x, find an integer j such that x lies between xx[j] and xx[j+1]. For this task, let us define fictitious array elements xx[0] and xx[n+1] equal to plus or minus infinity (in whichever order is consistent with the monotonicity of the table). Then j will always be between 0 and n, inclusive; a value of 0 indicates “off-scale” at one end of the table, n indicates off-scale at the other end. In most cases, when all is said and done, it is hard to do better than bisection, which will find the right place in the table in about log 2n tries. We already did use bisection in the spline evaluation routine splint of the preceding section, so you might glance back at that. Standing by itself, a bisection routine looks like this: void locate(float xx[], unsigned long n, float x, unsigned long *j) Given an array xx[1..n], and given a value x, returns a value j such that x is between xx[j] and xx[j+1]. xx must be monotonic, either increasing or decreasing. j=0 or j=n is returned to indicate that x is out of range. { unsigned long ju,jm,jl; int ascnd; jl=0; Initialize lower ju=n+1; and upper limits. ascnd=(xx[n] >= xx[1]); while (ju-jl > 1) { If we are not yet done, jm=(ju+jl) >> 1; compute a midpoint, if (x >= xx[jm] == ascnd) jl=jm; and replace either the lower limit else ju=jm; or the upper limit, as appropriate. } Repeat until the test condition is satisfied. if (x == xx[1]) *j=1; Then set the output else if(x == xx[n]) *j=n-1; else *j=jl; } and return. A unit-offset array xx is assumed. To use locate with a zero-offset array, remember to subtract 1 from the address of xx, and also from the returned value j. Search with Correlated Values Sometimes you will be in the situation of searching a large table many times, and with nearly identical abscissas on consecutive searches. For example, you may be generating a function that is used on the right-hand side of a differential equation: Most differential-equation integrators, as we shall see in Chapter 16, call

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

已到末页，全文结束

点击下载（PDF格式）

浏览记录