%matplotlib inline
%config InlineBackend.figure_format='svg'
from numpy import pi,zeros,arange,exp,sin,abs,linspace,dot,outer,real,mod,log
from numpy import sqrt,sum
from numpy.fft import fft
import matplotlib.pyplot as plt
The following function computes the trigonometric interpolant and plots it.
# Wave numbers are arranged as k=[0,1,...,N/2,-N/2+1,-N/2,...,-1]
def fourier_interp(N,f,ne=500,fig=True):
if mod(N,2) != 0:
print("N must be even")
return
h = 2*pi/N; x = h*arange(0,N);
v = f(x);
v_hat = fft(v)
k = zeros(N)
n = int(N/2)
k[0:n+1] = arange(0,n+1)
k[n+1:] = arange(-n+1,0,1)
xx = linspace(0.0,2*pi,ne)
vf = real(dot(exp(1j*outer(xx,k)), v_hat)/N)
ve = f(xx)
# Plot interpolant and exact function
if fig:
plt.plot(x,v,'o',xx,vf,xx,ve)
plt.legend(('Data','Fourier','Exact'))
plt.title("N = "+str(N))
errori = abs(vf-ve).max()
error2 = sqrt(h*sum((vf-ve)**2))
print("Error (max,L2) = ",errori,error2)
return errori, error2
f1 = lambda x: exp(sin(x))
fourier_interp(24,f1);
Error (max,L2) = 8.01581023779363e-14 4.559011236168057e-13
g = lambda x: sin(x/2)
e1i,e12 = fourier_interp(24,g)
Error (max,L2) = 0.024794983262922784 0.06907606978750452
e2i,e22 = fourier_interp(48,g,fig=False)
print("Rate (max,L2) = ", log(e1i/e2i)/log(2), log(e12/e22)/log(2))
Error (max,L2) = 0.012409450043334486 0.017665433009611636 Rate (max,L2) = 0.9986090713584319 1.9672568878077048
def trihat(x):
f = 0*x
for i in range(len(x)):
if x[i] < 0.5*pi or x[i] > 1.5*pi:
f[i] = 0.0
elif x[i] >= 0.5*pi and x[i] <= pi:
f[i] = 2*x[i]/pi - 1
else:
f[i] = 3 - 2*x[i]/pi
return f
e1i,e12 = fourier_interp(24,trihat)
Error (max,L2) = 0.03195085708656353 0.1114394822310562
e2i,e22 = fourier_interp(48,trihat,fig=False)
print("Rate (max,L2) = ", log(e1i/e2i)/log(2), log(e12/e22)/log(2))
Error (max,L2) = 0.015806700636915583 0.028104406730080144 Rate (max,L2) = 1.0153183695992953 1.9873921942542216
f2 = lambda x: (abs(x-pi) < 0.5*pi)
e1i,e12 = fourier_interp(24,f2)
Error (max,L2) = 0.9915112319641801 2.0366611449817986
e2i,e22 = fourier_interp(48,f2)
print("Rate (max,L2) = ", log(e1i/e2i)/log(2), log(e12/e22)/log(2))
Error (max,L2) = 0.9828401566916884 1.041560145620768 Rate (max,L2) = 0.0126723113206373 0.9674598167513287