Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 49
Текст из файла (страница 49)
Wereach Z0 by m applications of the algorithm (6.69), starting with Zm = Z.The following FORTRAN subprogram implements the algorithm justdescribed.SUBROUTINE FFT ( Zl, Z2, N, INZEE )CONSTRUCTS THE DISCRETE FOURIER TRANSFORM OF Zl (OR Z2) IN THE COOLEYC TUKEY WAY, BUT WITH A TWIST.INTEGER INZEE,N, AFTER,BEFORE,NEXT,NEXTMX,NOW,PRIME(12)COMPLEX Zl(N),Z2(N)C****** I N P U T ******C Z1, Z2 COMPLEX N-VECTORSC N LENGTH OF Zl AND 22C INZEE INTEGER INDICATING WHETHER Z1 OR Z2 IS TO BE TRANSFORMEDC=l , TRANSFORM Z1C=2 , TRANSFORM Z2C****** W O R K A R E A S ******C Z1, Z2 ARE BOTH USED AS WORKARRAYS282APPROXIMATIONC****** 0 U T P U T ******CZ1 OR Z2 CONTAINS THE DESIRED TRANSFORM (IN THE CORRECT ORDER)C INZEE INTEGER INDICATING WHETHER Zl OR Z2 CONTAINS THE TRANSFORM,C= 1 , TRANSFORM IS IN Z1C= 2 , TRANSFORM IS TN Z2C****** M E T H 0 D ******THE INTEGER N IS DIVIDED INTO ITS PRIME FACTORS (UP TO A POINT).CC FOR EACH SUCH FACTOR P , THE P-TRANSFORM OF APPROPRIATE P-SUBVECTORSC OF Zl (OR Z2) IS CALCULATED IN F F T S T P AND STORED IN A SUITC ABLE WAY IN Z2 (OR Z1).
SEE TEXT FOR DETAILS.CDATA NEXTMX,PRIME / 12, 2,3,5,7,11,13,17,19,23,29,31,37 /AFTER = 1BEFORE = NNEXT = 1C10 IF ((BEFORE/PRIME(NEXT))*PRIME(NEXT) .LT. BEFORE) THENNEXT = NEXT + 1IF (NEXT .LE. NEXTMX) THENGO TO 10ELSENOW = BEFOREBEFORE = 1END IFELSENOW = PRIME(NEXT)BEFORE = BEFORE/PRIME(NEXT)END IFCIF (INZEE .EQ. 1) THENCALL FFTSTP( Zl, AFTER, NOW, BEFORE, Z2 )ELSECALL FFTSTP( Z2, AFTER, NOW, BEFORE, Z1 )END IFINZEE = 3 - INZEERETURNIF (BEFORE .EQ. 1)AFTER = AFTER*NOWGO TO 10ENDSUBROUTINE FFTSTP ( ZIN, AFTER, NOW, BEFORE, ZOUT )CALLED IN F F T .CARRIES OUT ONE STEP OF THE DISCRETE FAST FOURIER TRANSFORM.INTEGER AFTER, BEFORE, NOW, IA,IB,IN,JREAL ANGLE,RATIO,TWOPICOMPLEX ZIN(AFTER,BEFORE,NOW),ZOUT(AFTER,NOW,BEFORE),*DATA TWOPI / 6.2831 85307 17958 64769 /ANGLE = TWOPI/FLOAT(NOW*AFTER)OMEGA = CMPLX(COS(ANGLE) ,-SIN(ANGLE))ARG = CMPLX(1.,0.)DO 100 J=l,NOWDO 90 IA=l,AFTERDO 80 IB=1,BEFOREVALUE = ZIN(IA,IE,NOW)DO 70 IN=NOW-1,1,-lVALUE = VALUE*ARG + ZIN(IA,IB,IN)70ZOUT(IA,J,IB) = VALUE80ARG = ARG*OMEGA90100 CONTINUERETURNENDARG,OMEGA,VALUEIf N is the product of m integers,N = P1P2 · · · Pmthen a program like the above makes it possible to compute the transform*6.6FAST FOURIER TRANSFORMS283W = N(P1 + P2 + · · · + Pm )operations (rather than N2 ).
Since, for integers Q, R greater than 1,Q + R < QR unless Q = R = 2, this number W is minimized if everyfactor of N is actually used, except that factors of 2 may be combined to 4without loss. Further,W/N = P1 + · · · + Pmandlog N = log P1 + · · · + log PmsoThis shows W/(N log N) to be a weighted average of the numbersP j/log Pj, j = 1, . . .
, m. It is easy to see that P /log2 P, as a function ofthe integer P, has the minimum value 1.89 . . . at P = 3, and has thevalue 2 at P = 2 and P = 4, and is only 3.0 1 . . . at P = 10. Hence1.89N log2 N < Wwhile, even for factors Pj as big as 10, W is no bigger than 3.02N log, N.Further savings occur in case the data vector z is real, since then(6.70)See (6.44) and (6.65).There are other FFTs available when N is a prime or when N is aproduct of integers which are pairwise relatively prime; see Winograd’sarticle [36].EXERCISES6.6-l Prove directly from the definition (6.64) that (6.70) holds in case z is real.6.6-2 Use FFT (with N = 81, say) to check your answers for Exercises 6.5-1 and 6.5-2. [Thiswill force you to pay close attention to all the details in (6.65)!]for6.6-3 Use FFT to calculate (approximately) the Fourier coefficientsb.
f(x) = sin (π x)a. f(x) = sin 3xusing, e.g., N = 81 or 324 or whatever. Why do the Fourier coefficients for f(x) = sin ( π x )fail to decay rapidly as |j| increases?6.6-4 Tailor the FFT program to the specific case N = 3.4, making whatever savings incalculations and storage you can.6.6-4 Improve FFTSTP by adding special coding as a replacement for the range of the DOloop over IB in case NOW = 2, 3, or 4 (say).6.6-6 Discuss the use of FFT for evaluating the trigonometric sumat the points x j = a + 2πj/(2n + 1), j = 0, .
. . , n, for some fixed a in the interval[0, 2π/(2n + 1)].84APPROXIMATION6.6-7 Make use of FFT to construct the trigonometric polynomial interpolant at the N = 2n+ 1 points xj = 2 π j/N, j = 0, . . . , N - 1, to the square wave f(x) = signum (sin x ), usingN = 35 Then use FFT again to evaluate the interpolant at the 105 points y j + 2 πj /105,j = 0, . . . , 104.
(Hint: Use Exercise 6.6-6.)6.6-8 Use FFT to construct an approximation to the spectrum of a function f(x) whosevalues at the points xj = 2π j/N, j = 0, . . . , N - 1, with N - 128 say, are obtained from a(pseudo-)random number generator giving numbers uniformly distributed between 0 and 1.Compare it with the spectrum of the function considered in Exercise 6.5-2.6.6-9 Using Exercise 6.6-8, discuss how one might use FFT to recover the values f(xj ),j = 0, .
. . , N - 1 of a “smooth” 2 π-periodic function f(x) from given data f(xj) + e j, all j,with ej uniformly distributed noise.6.6-10 Show that(This means that you get back the N-vector z from itsdiscrete Fourier transformby (a) changing all entries of to their complexconjugates, then (b) constructing the discrete Fourier transform of the resulting vector andthen (c) dividing each entry of the resulting vector by N.)6.6-11 Describe how you would use FFT to construct the polynomial interpolant of degree< n at the Chebyshev points (6.18) to given data. (Hint: Construct the interpolant as a linearcombination of the n + 1 Chebyshev polynomials T0 , . . . , Tn , using (6.9).
Subsequentevaluation would, of course, be via the FUNCTION CHEB.)6.7 PIECEWISE-POLYNOMIAL APPROXIMATIONA simple and familiar example of piecewise-polynomial approximation islinear interpolation in a table of values f(xi ), i = 1, . . . , N + 1, wherea = x1 < x2 < . . . < xN+1 = b. Here f(x) is approximated at a point bylocating the interval [xk,xk+1 ] which contains and then takingas the approximation toIn effect, f(x) is approximated over [a,b] bythe “broken line” or piecewise-linear function g 1 (x) (see Fig. 6.9) withbreakpoints x2, . . .
, xN, which interpolates f(x) at x1, . . . , xN+1. It followsfrom Example 2.6, applied to each of the subintervals [xk,xk+1] k =1, . . . , N, thatFor all xprovided that f(x) is twice differentiable on [a,b]. Note that we can makethe interpolation error as small as we wish by makingsmall for all k.Note further that such an increase in interpolation points does not complicate further work with g 1 (x), since g 1 (x) is “locally” a very simplefunction.By using a piecewise-polynomial function gr(x) of degree r > 1 insteadof the piecewise-linear g1(x), we can produce approximations to f(x) whoseerror term contains the (r + 1)st power of maxkhence goes to zerofaster than the error (6.71) for piecewise-linear interpolation as max6.7PIECEWISE-POLYNOMIAL APPROXIMATION285Figure 6.9 Broken-line interpolation.becomes small.
Piecewise-cubic approximation has become particularlypopular. We now discuss several piecewise-cubic interpolation schemes.Let f(x) be a real-valued function defined on some interval [a,b]. Wewish to construct a piecewise-cubic (polynomial) function g 3 (x) whichinterpolates f(x) at the points x1, . .
. , xN+1, wherea = x1 < x2 < · · · < xN+1 = b(6.72)As with piecewise-linear interpolation, we choose the interior interpolationpoints x2, . . . , xN to be the breakpoints for g3(x); that is, on each interval[xi , xi+1 ], we construct g3(x) as a certain cubic polynomial Pi (x), i =1, . . . , N.To facilitate the use of g3(x) in subsequent calculations, we write eachcubic piece Pi (x) of g3(x) asPi (x) = c 1 , i + c 2 , i (x - xi ) + c3,i (x - xi )2 + c4,i ( x - xi )3 (6.73)Once we know the coefficients cj,i , j = 1, . .
. , 4, i = 1, . . . , N, then thefollowing FORTRAN function PCUBIC efficiently evaluates g3(x) for anyparticular point x =REAL FUNCTION PCUBIC ( XBAR, XT, C, N )C RETURNS THE VALUE AT XBAR OF THE PIECEWISE CUBIC FUNCTION ON NC INTERVALS WITH BREAKPOINT SEQUENCE XT AND COEFFICIENTS C .INTEGER N, I,JREAL C(4,N),XBAR,XI(N+l), DXDATA I /l/IF (XBAR .GE. XT(I)) THENDO 10 J=I,NGO TO 30IF (XBAR .LT. XI(J+l))10CONTINUEJ = NELSEDO 20 J=I-1,1,-1GO TO 30IF (XBAR .GE.
XI(J))20CONTINUEJ = 1END IF30 I = JDX = XBAR - XI(I)PCUBIC = C(l,I) + DX*(C(2,1) + DX*(C(3,I) + DX*C(4,1)))RETURNEND286APPROXIMATIONWe now turn to the determination of the piecewise-cubic interpolatingfunction g3(x). Since we wanti = 1, . . . , N + lg (x ) = f(x )3iiwe must havePi (x i ) = f(x i )Note that (6.74) impliesP i (x i + l ) = f(x i + 1 )P i - 1 (x i ) = P i (x i )i = 1, . . . , N(6.74)i = 2, .
. . , Nso that g3(x) is guaranteed to be continuous on [a,b].Recall from Theorem 2.1 or 2.4 that we can always interpolate a givenfunction at four points by a cubic polynomial. So far, each of the cubicpieces Pi (x) is required to interpolate f(x) only at two points. Hence wehave still quite a bit of freedom in choosing the Pi (x). Different interpolation methods differ only in how this freedom is used.In piecewise-cubic Hermite interpolation, one determines Pi (x) so as tointerpolate f(x) at xi , xi , xi+1, xi+1, that is, so that alsoi = 1, . .
. , N (6.75)Pi ´(xi ) = f´(xi )Pi ´(xi+l ) = f´(x i+1 )It then follows from the Newton formula (2.32) that, for i = 1, . . . , N,Pi (x) = f(xi ) + f[xi , xi ](x - xi ) + f[xi , xi , xi+1 ](x - xi ) 2+ f [xi , xj, xj+1 xi+1](x - xi ) 2 (x - xi+1 )Since (x - xi+1 ) = (x - xi ) + (xi - xi+1 ), this givesPi(x) = f(xi ) + f´(xi )(x - xi ) + (f [ xi , xi , xi+1 ] - f [ xi , xi , xi+1, xi+1 ]∆ x i )× (x - xi )2 + f[ xi , xi , xi+l, xi+l ](x - xi ) 3where ∆xi = xi+1 - xi , from which we can read off directly thecoefficients c1,i, c2,i, c3,i, c4,i for Pi (x), Using the abbreviationsf i = f(x i )we getc1,i = fisi = f´(xi )i = 1 .
. . , N + l(6.76)c2,i = sic3,i = f[ xi , xi , xi+1 ] - f[ xi , xi , xi+1, xi+1 ] ∆x i(6.77)6.7PIECEWISE-POLYNOMIAL APPROXIMATION287With fi stored in c1,i and si stored in c2,i, i = 1, . . . , N + 1, the followingFORTRAN subroutine utilizes (6.77) to calculate c3,i, c4,i, i = 1, . . . , N.SUBROUTINE CALCCF ( XI, C, N )INTEGER N, IREAL C(4,N+l) ,XI (N+l) , DIVDFl,DIVDF3,DXC****** I N P U T ******C XI(l), .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
