606 Chapter 13.Fourier and Spectral Applications Various numerical schemes can be used to solve sparse linear systems of this "hierarchically band diagonal"form.Beylkin,Coifman,and Rokhlin [1]make the interesting observations that (1)the product of two such matrices is itself hierarchically band diagonal(truncating,of course,newly generated elements that are smaller than the predetermined threshold e);and moreover that(2)the product can be formed in order N operations. Fast matrix multiplication makes it possible to find the matrix inverse by Schultz's (or Hotelling's)method,see $2.5. Other schemes are also possible for fast solution ofhierarchically band diagonal forms.For example,one can use the conjugate gradient method,implemented in $2.7 as linbcg. 菲 CITED REFERENCES AND FURTHER READING: Daubechies,I.1992,Wavelets (Philadelphia:S.I.A.M.). Strang.G.1989,S/AM Review,vol.31.pp.614-627. Beylkin,G.,Coifman,R.,and Rokhlin,V.1991,Communications on Pure and Applied Mathe- RECIPES matics,vol.44,Pp.141-183.[1] Daubechies,1.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 Theorem 61 In 86.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 formula derives from numerical application of the sampling theorem (812.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: Numerica 10621 Consider an arbitrary function g(t)and the grid of sampling pointsn=o+nh,where n ranges over the integers and a is a constant that allows an arbitrary shift of the sampling 43106 grid.We then write %g7G g(t)= g(tn)sinc (t-tn)+e(t) (13.11.1) h where sincx =sin x/r.The summation over the sampling points is called the sampling representation of g(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(w)= g(t)edt (13.11.2) vanishes identically forw>h When can sampling representations be used to advantage for the approximate numerical computation of functions?In order that the error term be small,the Fourier transform G(w) 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)
606 Chapter 13. Fourier and Spectral Applications 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). Various numerical schemes can be used to solve sparse linear systems of this “hierarchically band diagonal” form. Beylkin, Coifman, and Rokhlin [1] make the interesting observations that (1) the product of two such matrices is itself hierarchically band diagonal (truncating, of course, newly generated elements that are smaller than the predetermined threshold ); and moreover that (2) the product can be formed in order N operations. Fast matrix multiplication makes it possible to find the matrix inverse by Schultz’s (or Hotelling’s) method, see §2.5. Other schemes are also possible for fast solution of hierarchically band diagonal forms. For example, one can use the conjugate gradient method, implemented in §2.7 as linbcg. CITED 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 Theorem 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 formula derives 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) = ∞ n=−∞ g(tn) sinc π h(t − tn) + e(t) (13.11.1) where sinc x ≡ sin x/x. The summation over the sampling points is called the sampling representation of g(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(ω) = ∞ −∞ g(t)eiωt dt (13.11.2) vanishes identically for |ω| ≥ π/h. When can sampling representations be used to advantage for the approximate 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)
13.11 Numerical Use of the Sampling Theorem 607 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 the function g(t)and its Fourier transform G(w)must rapidly approach zero for large values of their respective arguments. Unfortunately,these two conditions are mutually antagonistic-the Uncertainty Princi- ple in quantum mechanics.There exist strict limits on how rapidly the simultaneous approach to zero can be in both arguments.According to a theorem of Hardy [2],ifg(t)=O(e) asll一and Gw)=O(ew2)asw一o,theng()≡Ce-t2 where C is a constant.This can be interpreted as saying that of all functions the Gaussian is the most rapidly decaying in both t and w,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-2 81 e-13 =∑esimc-)+e0 (13.11.3) n=-00 The error e(t)depends on the parameters h and a as well as on t,but it is sufficient for 3 the present purposes to state the bound, le(t)h. When the summation in(13.11.3)is approximated by one with finite limits,say from Press. No-N to No +N,where No is the integer nearest to-/h,there is a further truncation error.However,if N is chosen so that N>/(2h2),the truncation error in the summation 9 is less than the bound given by (13.11.4),and,since this bound is an overestimate,we shall continue to use it for (13.11.3)as well.The truncated summation gives a remarkably accurate representation for the Gaussian even for moderate values of N.For example, le(t)<5 x 10-5 for h=1/2 and N 7;le(t)<2 x 10-10 for h 1/3 and N =15; IENTIFIC and le(t)<7x 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 many transcendental functions can be expressed as an integral involving the Gaussian,and by substituting (13.11.3)one can often find excellent approximations to the integrals as a sum over elementary functions. Let us consider as an example the function w(z)of the complex variable z=+iy, related to the complex error function by Numerical 色 w(z)=e-=2erfc(-iz) (13.11.5) 43126 having the integral representation w(z)= 1 fe-dt (13.11.6) mi Jc t-z where the contour C extends from-oo to oo,passing below z(see,e.g.,[31).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 1 w(z)≈ he-号1-(-1)”e-w(a-/h (13.11.7) tn-z n=-0 where we now omit the error term.One should note that there is no singularity as zm for some n=m,but a special treatment of the mth term will be required in this case (for example,by power series expansion). An alternative form of equation(13.11.7)can be found by expressing the complex expo- nential in(13.11.7)in terms of trigonometric functions and using the sampling representation
13.11 Numerical Use of the Sampling Theorem 607 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). 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 the function g(t) and its Fourier transform G(ω) must rapidly approach zero for large values of 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 approach to zero can be in both arguments. According to a theorem of Hardy [2], if g(t) = O(e−t2 ) as |t|→∞ and G(ω) = O(e−ω2/4) as |ω|→∞, then g(t) ≡ Ce−t2 , where C is a constant. This can be interpreted as saying that of all functions the Gaussian is the most rapidly decaying in both t and ω, 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 = ∞ n=−∞ e −t2 n sinc π h(t − tn) + e(t) (13.11.3) The error e(t) depends on the parameters h and α as well as on t, but it is sufficient for the present purposes to state the bound, |e(t)| π/h. When the summation in (13.11.3) is approximated by one with finite limits, say from N0 − N to N0 + N, where N0 is the integer nearest to −α/h, there is a further truncation error. However, if N is 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, we shall continue to use it for (13.11.3) as well. The truncated summation gives a remarkably accurate 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 many transcendental functions can be expressed as an integral involving the Gaussian, and by substituting (13.11.3) one can often find excellent approximations to the integrals as a sum over 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−z2 erfc(−iz) (13.11.5) having the integral representation w(z) = 1 πi C e−t2 dt t − z (13.11.6) where the contour C extends 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 ∞ n=−∞ he−t2 n 1 − (−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). 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
608 Chapter 13.Fourier and Spectral Applications (13.11.3)with z replacing t.This yields r(e)≈e2+1he-片1--l)"cas(-2/h (13.11.8) πi tn-z n=-00 This form is particularly useful in obtaining Re w(z)when ly 1.Note that in evaluating (13.11.7)the exponential inside the summation is a constant and needs to be evaluated only once;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 a.Eight interesting choices are:a=0,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 a,the choices involving a complex a are useful only if the imaginary part of z is not too large.This is not the place to catalog all sixteen possible formulas,and we give only two particular cases that show some of the important features. First of all let =0 in equation (13.11.8),which yields, w(2)≈e2 1 Che-(h21-（-1)”cos(rz/ (13.11.9) πi n=-00 nh-z 3 This approximation is good over the entire z-plane.As stated previously,one has to treat the RECIPES I case where one denominator becomes small by expansion in a power series.Formulas for the case a =0 were 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 a z in (13.11.7),which yields w(a)≈e2一2Te（e-nh)2 (13.11.10) T n the sum being over all odd integers(positive and negative).Note that we have made the substitution n-n in the summation.This formula is simpler than(13.11.9)and contains half the number of terms,but its error is worse ify is large.Equation(13.11.10)is the source IENTIFIC of the approximation formula(6.10.3)for Dawson's integral,used in $6.10. CITED REFERENCES AND FURTHER READING: Rybicki,G.B.1989,Computers in Physics,vol.3,no.2,pp.85-87.[1] 1992 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 by Dover Publications,New York).[3] Numerica 10621 Gautschi,W.1970,SIAM Journal on Numerical Analysis,vol.7,pp.187-198.[4] 43108 Armstrong,B.H.,and Nicholls,R.W.1972,Emission,Absorption and Transfer of Radiation in Recipes Heated Atmospheres (New York:Pergamon).[5] (outside Chiarella,C.,and Reichel,A.1968,Mathematics of Computation,vol.22,pp.137-143.[6] Goodwin,E.T.1949,Proceedings of the Cambridge Philosophical Society,vol.45,pp.241-245. North Software. 7]
608 Chapter 13. Fourier and Spectral Applications 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). (13.11.3) with z replacing t. This yields w(z) ≈ e−z2 + 1 πi ∞ n=−∞ he−t2 n 1 − (−1)n cos π(α − z)/h 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 only once; 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 imaginary part of z is not too large. This is not the place to catalog all sixteen possible formulas, and we give only two particular cases that show some of the important features. First of all let α = 0 in equation (13.11.8), which yields, w(z) ≈ e−z2 + 1 πi ∞ n=−∞ he−(nh)2 1 − (−1)n cos(π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 for the case α = 0 were 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 α = z in (13.11.7), which yields w(z) ≈ e−z2 − 2 πi n odd e−(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 → −n in the summation. This formula is simpler than (13.11.9) and contains half the number of terms, but its error is worse if y is large. Equation (13.11.10) is the source of the approximation formula (6.10.3) for Dawson’s integral, used in §6.10. 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 Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover 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., and Reichel, A. 1968, Mathematics of Computation, vol. 22, pp. 137–143. [6] Goodwin, E.T. 1949, Proceedings of the Cambridge Philosophical Society, vol. 45, pp. 241–245. [7]