We'll use this notebook to verify some of the analytic results we found in the tutorial problems. Finding Fourier transforms numerically is fast and very simple.
import numpy as np
import matplotlib.pyplot as plt
import scipy.fft as spf
def f(x):
if x > -1.0 and x < 1.0:
return x
else:
return 0.0
def analytic(k):
return 1j*np.sqrt(2/np.pi)*(k*np.cos(k) - np.sin(k))/k**2
f=np.vectorize(f)
xv = np.linspace(-10,10,500)
plt.figure()
plt.plot(xv, f(xv), 'k-')
plt.show()
F = spf.fft(f(xv))
kv = spf.fftfreq(len(xv), xv[1]-xv[0])
kvals = np.linspace(-5,5,500)
plt.figure()
plt.plot(kvals, np.abs(analytic(kvals*2.0*np.pi))/np.max(np.abs(analytic(kvals*2.0*np.pi))), 'k-', lw=3)
plt.plot(kv, np.abs(F)/np.max(np.abs(F)), 'r--')
plt.xlim(-5,5)
plt.show()
gamma = 0.5
k0 = 2.0*np.pi
xv = np.linspace(-1,10,500)
def f(x):
if x > 0:
return np.exp(-gamma*x)*np.cos(k0*x)
else:
return 0.0
f=np.vectorize(f)
def analytic(k):
return (1.0/np.sqrt(2.0*np.pi))*(1j*k + gamma)/((1j*k + gamma)**2 + k0**2)
plt.figure()
plt.plot(xv, f(xv), 'k-')
plt.show()
F = spf.fft(f(xv))
kv = spf.fftfreq(len(xv), xv[1]-xv[0])
kvals = np.linspace(-5,5,500)
plt.figure()
plt.plot(kvals, np.abs(analytic(kvals*2.0*np.pi))/np.max(np.abs(analytic(kvals*2.0*np.pi))), 'k-', lw=3)
plt.plot(kv, np.abs(F)/np.max(np.abs(F)), 'r--')
plt.xlim(-5,5)
plt.show()
