Relative error when computing derivatives via FFT

Question

I want to compute a discrete derivative via the FFT. This amounts to multiplication by the wave number in Fourier space, as detailed in the stack exchange answer here. When I increase the resolution, I expected that the relative error should decrease. However, it seems that it grows (linearly) with the resolution.

Can anybody explain what is going on? I will be using this to take many derivatives in a PDE solver, so I would like as much accuracy in the computation as possible. Here is some Matlab code I wrote to test it:

close all; clear all;

index = 1;
resolutions = 2.^(2:16);

a = 0;
b = 2*pi;

for N = resolutions
    dx = (b-a)/N;
    x = a + dx*(0:N-1);
    k = 2*pi/(b-a)*[0:N/2-1, 0, -N/2+1:-1];

    f     = sin(x);
    f_der = cos(x);

    fft_f = fft(f);
    f_approx = real(ifft(fft_f));
    error_fcn(index) = norm(f - f_approx)/norm(f);

    fft_f_der = 1i*k.*fft_f;
    f_der_approx = real(ifft(fft_f_der));
    error_der(index) = norm(f_der_approx - f_der)/norm(f_der);

    index = index + 1;
end

subplot(1,2,1);
loglog(resolutions,error_fcn);
title('fft error');

subplot(1,2,2);
loglog(resolutions,error_der);
title('derivative fft error');

Is the FFT just a bad method to use for this, or am I doing something wrong? Any help would be much appreciated.

EDIT: By request, the error plots are now attached.

The error in the derivative starts out near machine error, but then gets all the way up to 10^(-11) when N = 2^16 = 65536. On the other hand, I just got a decreasing error when I used:

f = x.^2.*(2*pi - x).^2
f_der = 4*x.*(2*pi - x).*(pi - x)

(that is, a function with zero derivative at 0 and 2*pi). Maybe the error is coming from a small jump in the sine function at the endpoint?

Answer 1

Your plots help rule out any obvious bugs like incorrect signs, thank you for posting them.

That is simply roundoff error. One way to check is to compute the same thing with higher precision. Also note that for $n$ terms the error is about $O(n\epsilon_\mathrm{mach})$.

In general, the error will be the sum of discretization error that goes as $O(n^{-\alpha})$ and roundoff error, which goes as $O(n \epsilon_\mathrm{mach})$. The former is decreasing with $n$, the latter increasing, so for $n$ above some value, the roundoff error will dominate. In this case, for $\sin x$, the discretization error is exactly zero because it has a finite Fourier series, so you only see the roundoff error, and it is increasing with $n$.