## Numerical Differentiation Using FourierJing B. WangUniversity of Western Australia Nedlands WA 6009 Australia wang@physics.uwa.edu.au A common numerical task is to evaluate the derivatives of a data set given on a numerical mesh. For example, in simulating quantum wavepacket scattering using the time-dependent Schrödinger's equation, the dominant computational effort is in the evaluation of the second-order spatial derivatives of system wavefunctions [1, 2]. In the five-point finite-difference method (FD5), the second derivative of a data set is given by where is the grid spacing. If higher accuracy is required, the seven-point finite-difference method (FD7) can be applied However, finite-difference methods are based on A more accurate and computationally efficient scheme [3], based on the discrete Fourier transform (DFT), uses the fact that if for a function , then the partial derivatives are It follows that the derivatives can be evaluated by firstly calculating using the DFT and secondly taking the inverse DFT of the product . For data points of a one-dimensional function , separated by uniform spacing , where , the continuous Fourier integral (3) can be approximated using the DFT where . The corresponding inverse DFT is where and . For any sampling interval , the maximum frequency in its Fourier transform is , which is the Nyquist critical frequency [4]. In other words, the valid frequency components are contained in with , while for , where indicates complex conjugation. Consequently, the inverse Fourier transform given by equation (3) is discretized to a good approximation as where for and for . Similarly, for a two-dimensional data set of points, the discrete form of equation (3) is The DFT scheme provides a For one-dimensional data sets, implementation of equation (4) is straightforward. The kernel is and the (Since there must be constant spacing, , between points, one could use the number of points as a parameter instead.) The following code is included to handle derivatives of order zero, required for the generalization to higher dimensions. The syntax of For the Gaussian data set, the results of Here, These figures show that To compare the numerical error in evaluating the second derivatives using the FD5, FD7, and DFT methods, and to examine the
relative convergence of these methods as a function of grid spacing, , we first implement equations (1) and (2) using j`,` ] forms the cyclic correlation in which the jth element of kernel is aligned with each element in list. Since the central element of is , we have for FD5.
Similarly, for When is reasonably large, say , the DFT scheme is several orders of magnitude more accurate than the FD5 and FD7 methods. There is a limitation to the DFT scheme because the Fourier transformation implies periodic boundary conditions in which the last point of the data set needs to be connected smoothly to the starting point. For Gaussian wavepackets, which rapidly approach zero far from the "center" of the wavepacket, such a boundary condition is clearly satisfied since the wavepackets have essentially zero values at the boundaries. This scheme works equally well for functions with periodic boundaries. The code is easily extended to For a two-dimensional Gaussian data set, the difference between the exact derivative, , and the DFT scheme is everywhere less than . The DFT representation of the derivative operator is Robert Knapp (rknapp@wolfram.com) comments: In applications it is important to cache the computation of Code for automatic computation of numerical derivatives using an extension of these ideas is planned for incorporation into
a future version of Copyright © 2002 Wolfram Media, Inc. All rights reserved. |