# Why does spectral accuracy of laplacian decrease with sampling size?

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.

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?