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 [ ]: