f13-11

600 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).13.11 Numerical Use of the SamplingTheorem In§6.10 we implemented an approximating formula for Dawson’s integral due to Rybicki. Now that we have become Fourier sophisticates, we can learn that the formuladerives from numerical application of the sampling theorem ( §12.1), normally considered to be 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 n ranges over the integers and αis a constant that allows an arbitrary shift of the sampling grid. We then write g(t)=∞/summationdisplay n=−∞g(tn)s i n cπ h(t−tn)+e(t)( 13.11.1 ) where sincx≡sinx/x. The summation over the sampling points is called the sampling representation ofg(t), and e(t)is its error term. The sampling theorem asserts that the sampling representation is exact, that is, e(t)≡0, if the Fourier transform of g(t), G(ω)=/integraldisplay∞ −∞g(t)eiωtdt (13.11.2 ) vanishes identically for |ω|≥π/h. Whencan sampling representations be used toadvantage fortheapproximate numerical computation 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 summation in (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 are led 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 values of their respective arguments. Unfortunately, thesetwoconditions aremutuallyantagonistic—theUncertaintyPrinci- pleinquantum mechanics. Thereexiststrictlimitsonhow rapidlythesimultaneous approach to zero can be in both arguments. According to a theorem of Hardy [2],i fg(t)=O(e−t2) as|t|→∞andG(ω)=O(e−ω2/4)as|ω|→∞, then g(t)≡Ce−t2, where Cis a constant. This can be interpreted as saying that of all functions the Gaussian is the mostrapidly decaying in both tandω, and in this sense is the “best” function to be expressed numerically as a sampling representation. Let us then write for the Gaussian g(t)=e −t2, e−t2=∞/summationdisplay n=−∞e−t2 nsincπ h(t−tn)+e(t)( 13.11.3 ) The error e(t)depends on the parameters handαas well as on t, but it is sufficient for the 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 the Gaussian 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 from N0−NtoN0+N, where N0is the integer nearest to −α/h, there is a further truncation error. However, if Nis chosen so that N>π / (2h2), the truncation error in the summation is 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 −5forh=1/2andN=7;|e(t)|<2×10−10forh=1/3andN=1 5; and|e(t)|<7×10−18forh=1/4andN=2 5. 13.11NumericalUseoftheSamplingTheorem 601Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).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 by w(z)=e −z2erfc(−iz)( 13.11.5 ) having the integral representation w(z)=1 πi/integraldisplay Ce−t2dt t−z(13.11.6 ) where the contour Cextends from −∞to∞,passing below z(see, e.g., [3]). Many methods exist 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 obtain w(z)≈1 πi∞/summationdisplay n=−∞he−t2 n1−(−1)ne−πi(α−z)/h tn−z(13.11.7 ) where we now omit the error term. One should note that there is no singularity as z→tm for some n=m, but a special treatment of the mth term will be required in this case (for example, by power series expansion). Analternativeformofequation(13.11.7)canbefoundbyexpressing thecomplex expo- nential in (13.11.7) interms of trigonometric functions and using the sampling representation(13.11.3) with zreplacing t. This yields w(z)≈e−z2+1 πi∞/summationdisplay n=−∞he−t2 n1−(−1)ncosπ(α−z)/h tn−z(13.11.8 ) This form is particularly useful in obtaining Re w(z)when|y|/lessmuch 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, orz,plus thevalues obtained by adding h/2toeach ofthese. Sincethe errorbound (13.11.3) assumed areal value of α, thechoices involving a complex αareuseful only ifthe imaginary partof zisnottoo large. Thisisnot theplace tocatalog allsixteen possible formulas,and we give only two particular cases that show some of the important features. First of all let α=0in equation (13.11.8), which yields, w(z)≈e−z2+1 πi∞/summationdisplay n=−∞he−(nh)21−(−1)ncos(πz/h ) nh−z(13.11.9 ) This approximation is good over the entire z-plane. As stated previously, one has to treat the case where one denominator becomes small by expansion in a power series. Formulas forthe case α=0were discussed briefly in [5]. They are similar, but not identical, to formulas derived by Chiarella and Reichel [6], using the method of Goodwin [7]. Next, let α=zin (13.11.7), which yields w(z)≈e−z2−2 πi/summationdisplay nodde−(z−nh)2 n(13.11.10 ) the sum being over all odd integers (positive and negative). Note that we have made the substitution n→−nin the summation. This formula is simpler than (13.11.9) and contains half the number of terms,but its error is worse if yis large. Equation (13.11.10) isthe source of the approximation formula (6.10.3) for Dawson’s integral, used in §6.10. 602 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).CITED REFERENCES AND FURTHER READING: Rybicki, G.B. 1989, Computers in Physics , vol. 3, no. 2, pp. 85–87. [1] Hardy, G.H. 1933, Journal of the London Mathematical Society , vol. 8, pp. 227–231. [2] Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York). [3] Gautschi, W. 1970, SIAM Journal on Numerical Analysis , vol. 7, pp. 187–198. [4] Armstrong, B.H., and Nicholls, R.W. 1972, Emission, Absorption and Transfer of Radiation in Heated Atmospheres (New York: Pergamon). [5] Chiarella, C., andReichel, A. 1968, Mathematics of Computation , vol. 22, pp. 137–143. [6] Goodwin,E.T.1949, ProceedingsoftheCambridgePhilosophicalSociety ,vol.45,pp.241–245. [7]