We know that for any real-valued function $f(x,y,z)$ whose Fourier transform is $\mathcal F[f]$, its laplacian can be computed from a spectral interpolant as follows.

$$ \Delta f(x,y,z) \simeq \sum_{k_1=-N/2+1}^{N/2}\sum_{k_2=-N/2+1}^{N/2}\sum_{k_3=-N/2+1}^{N/2} -(k_1^2+k_2^2+k_3^2)\, \mathcal F[f]\, e^{i(k_1x + k_2y + k_3z)} $$

One could easily implement this in Python or C, for instance, and check the spectral accuracy of this scheme for different values of the sampling size, $N$.

from numpy import linspace, arange, concatenate, meshgrid, pi, sin, cos
from numpy.fft import rfftn as fft, irfftn as ifft

def rms_err(N):
    # set up the coordinates
    xj       = linspace(0,2*pi,N+1)[:-1]
    y,x,z    = meshgrid(xj,xj,xj)
    last_bin = arange(0,N//2+1)
    othrbins = concatenate((last_bin, arange(-N//2+1,0)), axis=0)
    k2,k1,k3 = meshgrid(othrbins,othrbins,last_bin)

    # periodic test function and its spectral laplacian
    f = cos(4*z)*sin(x)
    F = ifft( -( k1**2 + k2**2 + k3**2 )*fft(f) )

    return np.sqrt(np.sum((F+17*f)**2))

For double-precision floating point numbers in Python, the error increases with the sampling size.

root mean squared error

For single-precision floating point numbers in C, the error is much worse. For $N=1024$, the error is of the order of $10^{-2}$, which is very bad.

Why is this? And how can I improve the spectral accuracy for large sampling sizes and single-precision floats?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.