c13-11 (Numerical Recipes in C)

606Chapter 13.Fourier and Spectral ApplicationsCITED REFERENCES AND FURTHER READING:Daubechies, I. 1992, Wavelets (Philadelphia: S.I.A.M.).Strang, G. 1989, SIAM Review, vol. 31, pp. 614–627.Beylkin, G., Coifman, R., and Rokhlin, V. 1991, Communications on Pure and Applied Mathematics, vol. 44, pp. 141–183.

[1]Daubechies, I. 1988, Communications on Pure and Applied Mathematics, vol. 41, pp. 909–996.[2]Vaidyanathan, P.P. 1990, Proceedings of the IEEE, vol. 78, pp. 56–93. [3]Mallat, S.G. 1989, IEEE Transactions on Pattern Analysis and Machine Intelligence, vol.

11,pp. 674–693. [4]Freedman, M.H., and Press, W.H. 1992, preprint. [5]13.11 Numerical Use of the Sampling TheoremIn §6.10 we implemented an approximating formula for Dawson’s integral due toRybicki. Now that we have become Fourier sophisticates, we can learn that the formuladerives from numerical application of the sampling theorem (§12.1), normally considered tobe a purely analytic tool. Our discussion is identical to Rybicki [1].For present purposes, the sampling theorem is most conveniently stated as follows:Consider an arbitrary function g(t) and the grid of sampling points tn = α + nh, where nranges over the integers and α is a constant that allows an arbitrary shift of the samplinggrid.

We then writeg(t) =∞Xg(tn ) sincn=−∞π(t − tn ) + e(t)h(13.11.1)where sinc x ≡ sin x/x. The summation over the sampling points is called the samplingrepresentation of g(t), and e(t) is its error term. The sampling theorem asserts that thesampling representation is exact, that is, e(t) ≡ 0, if the Fourier transform of g(t),Z ∞G(ω) =g(t)eiωt dt(13.11.2)−∞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).to near neighbors in its own hierarchy (square blocks along the main diagonal) andnear neighbors in other hierarchies (rectangular blocks off the diagonal).The number of nonnegligible elements in a matrix like that in Figure 13.10.5scales only as N , the linear size of the matrix; as a rough rule of thumb it is about10N log10 (1/), where is the truncation level, e.g., 10−6 . For a 2000 by 2000matrix, then, the matrix is sparse by a factor on the order of 30.Various numerical schemes can be used to solve sparse linear systems of this“hierarchically band diagonal” form. Beylkin, Coifman, and Rokhlin [1] makethe interesting observations that (1) the product of two such matrices is itselfhierarchically band diagonal (truncating, of course, newly generated elements thatare smaller than the predetermined threshold ); and moreover that (2) the productcan be formed in order N operations.Fast matrix multiplication makes it possible to find the matrix inverse bySchultz’s (or Hotelling’s) method, see §2.5.Other schemes are also possible for fast solution of hierarchically band diagonalforms.

For example, one can use the conjugate gradient method, implemented in§2.7 as linbcg.13.11 Numerical Use of the Sampling Theorem6072e−t =∞X2e−tn sincn=−∞π(t − tn ) + e(t)h(13.11.3)The error e(t) depends on the parameters h and α as well as on t, but it is sufficient forthe present purposes to state the bound,|e(t)| < e−(π/2h)2(13.11.4)which can be understood simply as the order of magnitude of the Fourier transform of theGaussian at the point where it “spills over” into the region |ω| > π/h.When the summation in (13.11.3) is approximated by one with finite limits, say fromN0 − N to N0 + N , where N0 is the integer nearest to −α/h, there is a further truncationerror. However, if N is chosen so that N > π/(2h2 ), the truncation error in the summationis less than the bound given by (13.11.4), and, since this bound is an overestimate, weshall continue to use it for (13.11.3) as well.

The truncated summation gives a remarkablyaccurate representation for the Gaussian even for moderate values of N . For example,|e(t)| < 5 × 10−5 for h = 1/2 and N = 7; |e(t)| < 2 × 10−10 for h = 1/3 and N = 15;and |e(t)| < 7 × 10−18 for h = 1/4 and N = 25.One may ask, what is the point of such a numerical representation for the Gaussian,which can be computed so easily and quickly as an exponential? The answer is that manytranscendental functions can be expressed as an integral involving the Gaussian, and bysubstituting (13.11.3) one can often find excellent approximations to the integrals as a sumover elementary functions.Let us consider as an example the function w(z) of the complex variable z = x + iy,related to the complex error function by2w(z) = e−z erfc(−iz)(13.11.5)having the integral representationw(z) =1πiZ2Ce−t dtt−z(13.11.6)where the contour C extends from −∞ to ∞, passing below z (see, e.g., [3]).

Many methodsexist for the evaluation of this function (e.g., [4]). Substituting the sampling representation(13.11.3) into (13.11.6) and performing the resulting elementary contour integrals, we obtainw(z) ≈∞n −πi(α−z)/h2 1 − (−1) e1 Xhe−tnπi n=−∞tn − z(13.11.7)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).vanishes identically for |ω| ≥ π/h.When can sampling representations be used to advantage for the approximate numericalcomputation of functions? In order that the error term be small, the Fourier transform G(ω)must be sufficiently small for |ω| ≥ π/h.

On the other hand, in order for the summationin (13.11.1) to be approximated by a reasonably small number of terms, the function g(t)itself should be very small outside of a fairly limited range of values of t. Thus we areled to two conditions to be satisfied in order that (13.11.1) be useful numerically: Both thefunction g(t) and its Fourier transform G(ω) must rapidly approach zero for large valuesof their respective arguments.Unfortunately, these two conditions are mutually antagonistic — the Uncertainty Principle in quantum mechanics. There exist strict limits on how rapidly the simultaneous approach2to zero can be in both arguments.

According to a theorem of Hardy [2], if g(t) = O(e−t )22as |t| → ∞ and G(ω) = O(e−ω /4 ) as |ω| → ∞, then g(t) ≡ Ce−t , where C is aconstant. This can be interpreted as saying that of all functions the Gaussian is the mostrapidly decaying in both t and ω, and in this sense is the “best” function to be expressednumerically as a sampling representation.2Let us then write for the Gaussian g(t) = e−t ,608Chapter 13.Fourier and Spectral Applicationswhere we now omit the error term. One should note that there is no singularity as z → tmfor some n = m, but a special treatment of the mth term will be required in this case (forexample, by power series expansion).An alternative form of equation (13.11.7) can be found by expressing the complex exponential in (13.11.7) in terms of trigonometric functions and using the sampling representation(13.11.3) with z replacing t.

This yields∞n2 1 − (−1) cos π(α − z)/h1 Xhe−tnπi n=−∞tn − z(13.11.8)This form is particularly useful in obtaining Re w(z) when |y| 1. Note that in evaluating(13.11.7) the exponential inside the summation is a constant and needs to be evaluated onlyonce; a similar comment holds for the cosine in (13.11.8).There are a variety of formulas that can now be derived from either equation (13.11.7)or (13.11.8) by choosing particular values of α.

Eight interesting choices are: α = 0, x, iy,or z, plus the values obtained by adding h/2 to each of these. Since the error bound (13.11.3)assumed a real value of α, the choices involving a complex α are useful only if the imaginarypart of z is not too large. This is not the place to catalog all sixteen possible formulas, and wegive only two particular cases that show some of the important features.First of all let α = 0 in equation (13.11.8), which yields,2w(z) ≈ e−z +∞n2 1 − (−1) cos(πz/h)1 Xhe−(nh)πi n=−∞nh − z(13.11.9)This approximation is good over the entire z-plane.

As stated previously, one has to treat thecase where one denominator becomes small by expansion in a power series. Formulas forthe case α = 0 were discussed briefly in [5]. They are similar, but not identical, to formulasderived by Chiarella and Reichel [6], using the method of Goodwin [7].Next, let α = z in (13.11.7), which yields2w(z) ≈ e−z −22 X e−(z−nh)πin(13.11.10)n oddthe sum being over all odd integers (positive and negative). Note that we have made thesubstitution n → −n in the summation. This formula is simpler than (13.11.9) and containshalf the number of terms, but its error is worse if y is large. Equation (13.11.10) is the sourceof the approximation formula (6.10.3) for Dawson’s integral, used in §6.10.CITED REFERENCES AND FURTHER READING:Rybicki, G.B.