# 538 Day -2: Trigonometric interpolation and spectral differentiation¶

with careful treatment of Nyquist frequency

DFT of sinusoids:

In [ ]:
np.set_printoptions(linewidth=300)
twopi = 2*pi

n = 10
j = arange(n)
phi = twopi*np.random.rand()
kstar = 1
np.around(  fft(cos(twopi*kstar*j/n +phi))  ,  2)


Trig poly interp and spectral derivative examples:

In [18]:
import numpy as np
from numpy import pi,exp,cos,sin,arange
np.seterr(divide='ignore', invalid='ignore')
from scipy.fftpack import fft,ifft
from pylab import plot,ion,ioff,draw,title,figure,axvline,legend
import matplotlib.pyplot as plt
%config InlineBackend.figure_format='retina'

for example in [6,2]:

#figure(example,figsize=(15,10))
fig, axes = plt.subplots(3,1,figsize=(15,10))

if example==1:
mytitle = 'zigzag'
c = 0.
d = 1.
L = d-c
n = 8
x = 0.5 + (-1)**arange(n)#np.array([1,0,1,0,1,0])+0.5
def myf(tt):
if len(tt)!=len(x): return tt/0  # unplottable
return x
def myfprime(tt):
return tt/0   # unplottable
t = L*np.arange(n)/float(n)
elif example==2:
mytitle = 'distorted sine'
n = 28
c = 0.
d = 1.
L = d-c
t = np.linspace(0.,L,n,endpoint=False)
r = 3
def myf(t): return np.exp(r*np.sin(2*pi*t/L))/20
def myfprime(t): return exp(r*sin(2*pi*t/L))*r*cos(2*pi*t/L)*(2*pi/L)/20
x = myf(t)
xp = myfprime(t)
elif example==3:
mytitle = 'Runge function'
n = 18
c = -1.
d =  1.
L = d-c
t = np.linspace(c,d,n,endpoint=False)
def myf(t): return 1/(1+12.*t*t)
def myfprime(t): return -myf(t)**2 * 24*t
x = myf(t)
xp = myfprime(t)

elif example==4:
mytitle = 'odd sine'
n = 8
c = -1.
d =  1.
L = d-c
t = linspace(c,d,n,endpoint=False)
print( len(t) )
def myf(t): return np.sin(3*pi*(t-c)/L)
x = myf(t)

elif example==5:
mytitle = 'vector [0,1,2,3,4,5,6,7,8,9]'
c = 0.
d = 1.
L = d-c
x = np.array([0,1,2,3,4,5,6,7,8,9],dtype=float)
n = len(x)
t = L*np.arange(n)/float(n)

elif example==6:
mytitle = 'sinusoid'
c = 0.
d = 1.
L = d-c
n = 10
amp = 1
kstar = n//2#-2   # issue if kstar = n/2
phaseshift = 1.
t = np.linspace(c,d,n,endpoint=False)
x = amp*np.cos(2*np.pi*kstar*t/L + phaseshift)
def myf(t): return amp*np.cos(2*np.pi*kstar*t/L + phaseshift)
def myfprime(t): return -amp*2*np.pi*kstar/L*np.sin(2*np.pi*kstar*t/L + phaseshift)

else:
print("Which example do you want me to do?")

y = fft(x)/np.sqrt(n) # division to get Sauer def of DFT

tt = np.linspace(c,d,1000)

p  = np.zeros_like(tt,dtype=complex)  # for interpolant to x
pp = np.zeros_like(tt,dtype=complex)  # for interpolant to x'

really_smart = True
if really_smart:
p = np.zeros_like(tt,dtype=complex)
pp = np.zeros_like(tt,dtype=complex)
assert n%2==0
kvals = range(-n//2,n//2+1)
for k in kvals:
term = y[k%n]/np.sqrt(n)*np.exp(1j*2.*np.pi*k*(tt-c)/L)
if k==-n//2 or k==n//2:
p += term/2
pp += 1j*2.*np.pi*k/L*term/2
else:
p += term
pp += 1j*2.*np.pi*k/L*term
else:
assert(0)

ax = axes[0]
ax.plot([c,d],[0,0],'k',alpha=0.5)
for tj in t: ax.axvline(tj,color='k',alpha=0.25)
ax.plot( t, x, 'o',clip_on=False)
ax.plot(tt,myf(tt)   ,'k',label='f',alpha=0.5)
ax.plot(tt,np.real(p),'r',label='real(p(t))')
ax.plot(tt,np.imag(p),'b',label='imag(p(t))')
ax.set_xlim(c,d)
ax.legend()
ax.set_title(mytitle)

ax = axes[1]
for tj in t: ax.axvline(tj,color='k',alpha=0.25)
ax.plot([c,d],[0,0],'k',alpha=0.5)
ax.plot(tt,myfprime(tt)   ,'c',label="f'",alpha=0.5)
ax.plot(tt,np.real(pp),'m',label="real(p'(t))")
ax.set_xlim(c,d)
ax.set_title(mytitle)
ax.plot(t, (np.roll(x,-1)-np.roll(x,1))/2/(L/n), 'ko', label='error in FD derivative')
ax.legend()

ax = axes[2]
for tj in t: ax.axvline(tj,color='k',alpha=0.25)
ax.plot(tt,np.real(pp)-myfprime(tt)   ,'g',label="error in p'",alpha=0.5)
ax.set_xlim(c,d)
ax.set_title(mytitle)
ax.plot(t, (np.roll(x,-1)-np.roll(x,1))/2/(L/n) - myfprime(t), 'ko', label='error in FD derivative')
ax.legend()


In [ ]:


In [ ]: