I'm trying to calculate the derivative of a function using FFT in Matlab. I've coded the following function to calculate it:
function dfdx = derivative(x,psi_0,n)
Nx = max(size(psi_0));
k = 2*pi/(range(x))*[0:Nx/2-1 0 -Nx/2+1:-1];
dfdx = ifft(((1i*k).^n).*fft(psi_0));
end
For higher derivatives (let's say for the $20$th derivative) it doesn't work. More precisely, the plot of the $20$th derivative doesn't make sense.
Any idea what causes this, and how this can be corrected?
Let's suppose you have a function $f$ in $L^{2}$, and for simplicity, let's suppose it's a single-variable function defined over the whole real line. If its Fourier transform $\mathcal{F}(f)$ is $\hat{f}$, then $\mathcal{F}(f')(\xi) = 2\pi i \xi\hat{f}(\xi)$.
If you're calculating the $n$th derivative, iterating the process above yields $\mathcal{F}(f^{(n)})(\xi) = (2 \pi i \xi)^{n}\hat{f}(\xi)$, suggesting that the calculation of higher derivatives becomes ill-conditioned. Assuming $\hat{f}$ is a well-conditioned function at $\xi$, the large coefficient in the derivative expression implies that small perturbations in $\xi$ perturb values of the Fourier transform of the $n$th derivative significantly for large $n$.
For instance, if $n = 20$, the magnitude of the constant pre-factor will be $(2\pi)^{20}$, or roughly $10^{16}$. With such a large constant factor, it would be surprising if you obtained meaningful results in floating-point arithmetic. You could compare this approach to calculating derivatives numerically through finite-difference approximations; calculating a 20th derivative with this approach would yield garbage in most cases.
If you used arbitrary precision arithmetic, you would probably see more accurate results. As $n$ increases, you'll require more digits of precision in your input to obtain an output accurate to at least 8 significant figures (this cutoff is arbitrary, but presumably, you want many accurate significant figures).
