numpy - Fast Fourier Transform performance in Python

Question

So I wrote a short Python program to estimate the accuracy of the Python's FFT method.

import numpy as np
import matplotlib.pyplot as plt

#Aufgabe 1
x0=0
a=2.5
k0=3
X=np.linspace(-4,4,100)
timestep=0.1
k=np.fft.fftfreq(X.size,d=timestep)
psi_analytical=[(2/(np.pi*a**2))**(1/4)*np.exp(-((i-x0)**2)/a**2)*np.exp(1j*k0*(i-x0)) for i in X]
psi_tilde_numerical=np.fft.fft(psi_analytical)
psi_tilde_analytical=[(2/(np.pi*a**2))**(1/4)*(a/2)*np.exp(-(a*(i-k0))**2/4)*np.exp(-1j*i*x0) for i in k]
psi_numerical=np.fft.ifft(psi_tilde_analytical)


#plt.plot(k,np.abs(psi_tilde_numerical),label='numerical psi tilde')
#plt.plot(k,np.abs(psi_tilde_analytical),'--',color='tab:orange', label='analytical psi tilde')

plt.plot(X,np.abs(psi_analytical),label='analytical psi, real')
plt.plot(X,np.abs(psi_numerical),'--',color='tab:orange',label='numerical psi, real')
plt.legend()
plt.show()

The analytical function is as follows:

To my surprise, the numerical and analytical functions are totally different. However, I'm not sure why this is the case.

The normalisation constant N is (2/(np.pi*a**2))**(1/4)

Answer 1

Doing a bit more research, I think I might have an answer for you.

Discrete Fourier Transform (DFT), which is computed efficiently using the Fast Fourier Transform algorithm (FFT), operates on discrete time domain signals. The Fourier Transform (FT) operates on function in continuous time domain.

DFT will approximate the FT under certain condition. One of those conditions is that the signal has to be band limited. This means that the FT for the function has to be zero for all frequencies above a certain frequency threshold α and the DFT has to have a sample rate that is at least 2*α This goes back to Nyquist-Shannon sampling theorem.

In your case, you are trying to approximate a Gaussian function exp(-x2) which is not band limited. This is because as you can see from your formulas, FT of a Gaussian is also a Gaussian. This means that it has negligible but non-zero components for frequencies all the way to infinity. As a result, you won't be able to approximate the FT using a DFT since you would need to have infinite sampling rate.

In conclusion, its important to realize that the DFT and FT are vastly different transforms and thus can not just be compared.