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 local approximations to the derivative operator, which brings about error. Also, the convergence with decreasing grid spacing is slow.

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 global representation of the derivative operator and is highly accurate for functions with periodic boundary conditions.

For one-dimensional data sets, implementation of equation (4) is straightforward. The kernel is

and the nth discrete derivative of the data set, data, with spacing , reads

(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 DFTDerivative is based on that of Derivative. Re is included to eliminate the small imaginary component that can arise during computation. Although the derivation refers to the interval , the implementation works for any interval with spacing .

For the Gaussian data set,

the results of DFTDerivative[n] for are overlaid upon the exact value of the derivative.

Here, DFTDerivative[9][data][] is compared to the exact value.

These figures show that DFTDerivative[n] is extremely accurate even for high and large .

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 ListCorrelate. ListCorrelate[kernel, list, 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 FD7.

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 arbitrary dimensions, with possibly different grid spacings for each dimension, using Outer.

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 global, giving it improved accuracy over local representations. Moreover, the calculation of higher-order derivatives using the DFT scheme requires little extra computing in contrast to, for example, the finite-difference methods. The code is simple, its structure transparent, and the generalization to arbitrary dimension is straightforward.

Robert Knapp (rknapp@wolfram.com) comments:

In applications it is important to cache the computation of Kernel. If the data consists of machine numbers, improved efficiency is possible if Kernel is a packed array. Doing this takes more work (and some type-detecting code), but it can make a huge difference in performance.

Code for automatic computation of numerical derivatives using an extension of these ideas is planned for incorporation into a future version of Mathematica.