Re [Ψ] is an even function with damped oscillations. Interestingly, for small τ, oscillations appear to be
more persistent and the function is still non-negligibly small at ζ 1 = ±300 and ζ2 ∈ (—300, 300). For large
τ, the integrand has a pronounced peak at (0, 0) and oscillations quickly fade out. The function approaches
zero at different rates in each direction. As predicted, the rate of convergence is slower in ζ2 direction.
Therefore, in FFT, the cutoff points may be chosen (and in Fortran programs are chosen) such that:
∣a1∣ < ∣a2∣ and ∣δ1∣ < ∣½∣. Moreover, for smaller τ it is desirable to have relatively large cutoffs (in absolute
value). For larger τ, ∣α1∣ , ∣0,2∣, ∣δ1∣ , ∣½∣ may be much smaller, but ∆1 and Δ2 must be as tiny as possible to
“capture” the peak.
It should be evident that numerical integration of the integrand function is unlikely to succeed. My
experiments with different numerical integration routines reveal the following. Fast 5-degree polynomial
methods produce very innacurate results and large negative values for many choices of sγ and -v√. Quadrature
methods (Gauss-Kronrod with adaptive integration) usually fail to reach desired precision when the cutoff
points are large in absolute value, even if the number of abscissae is high.8 Romberg integration routines are
reasonably accurate, but accuracy comes at a cost of very long time spans to compute integrals with desired
precision.
6.4. P.D.F.
Numerical integration of с г<8|“7'+^2vτ) . ψ (ζ 1,ζ2) to recover the p.d.f. of (sγ,vγ)' under P is impractical.
Accurate evaluation of the integral for just one point (sγ,vγ)' takes a long time, and such evaluations must
be done (iteratively) for all points of a reasonably fine grid of (sγ,vγ)' to implement recursion (12).
Contrastingly, FFT allows to recover f (s^ ,vγ ) on the whole grid quickly in one step. While being
economical in terms of computational time, FFT places high demand on computer memory.
It is desirable to have relatively large (in absolute value) cutoff points and small step sizes, Δ1 and Δ2,
in (ζ 1,ζ2)' grid. Since Δ1 = blj^°1, Δ2 = ^ʌ-ʃ2, the dimensions of the Fourier coefficient matrix, 2N1 and
N2, are necessarily large.9 Worse, FFT algorithms require N1 and N2 to be powers of 2. This leaves little
leeway in choosing the step sizes Δ1 and Δ2, once the cutoffs are set to values beyond which the integrand is
negligibly small. Note that the step sizes in (sγ ,vτ )' grid are b 2^α and b 'jj∖1, respectively; bounds of this
grid are: s⅞lin = v^in = 0, and s⅞lax = 2",..'∖. ',v^ax = 2⅞^2o~2υ .
Picking a finer grid of (ζ 1, ζ2)' implies at least doubling the size of the Fourier coefficient matrix. Modern
32-bit computers theoretically cannot address arrays with more than 232 — 1 elements. In practice, Fortran
8A standard recommendation for oscillatory functions is to set this number to the maximum allowed by the algorithm.
9The Fourier coefficients, Φj1,J2, are complex numbers. The bivariate discrete complex FFT routine I use to program the
transformation stores an Ni × N array of complex numbers as a 2Ni X N array of real numbers: real parts are in odd-numbered
rows, imaginary parts are in even-numbered rows.
15