%
%  lect11b.m -- code for FFT -- spectral accuracy
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  geometric series example

nC = [];
maxC = [];
m = [2 3 5 7 9 11];

%  double loop for N having simple factorization

for mm = 1:1:6
k = m(mm);

for nn = 1:1:6
L= pi;  N = k*(2^nn);  h = 2*L/N;

%  physical & spectral grids

x = 0:h:2*L-h;
n = -N/2:1:N/2-1;

%  sum formula for geometric series

r = 0.25;  f = (1 - r^2)./(1 + r^2 - 2*r*cos(x));

c  = fftshift(fft(f))/N;
cx = r.^abs(n);

%  data for this N

nC   = [nC ; N];
maxC = [maxC ; max(abs(real(c)-cx))];

end
end

clf
subplot(2,1,1)
semilogy(nC,maxC,'x')
title('\bf geometric series')
xlabel('\bf N-axis');  ylabel('\bf max|c_n^{fft}-c_n^{exact}|')

axis([2 80 10^(-16) 1])

%- - - - - - - 

%  cosine exponential example

nC = [];
maxC = [];

%  double loop for N having simple factorization

for mm = 1:1:6
k = m(mm);

for nn = 1:1:6
L= pi;  N = k*(2^nn);  h = 2*L/N;

%  physical & spectral grids

x = 0:h:2*L-h;
n = -N/2:1:N/2-1;

%  sum formula for geometric series

r = 2.0;  f = exp(r*cos(x));

c  = fftshift(fft(f))/N;
cx = besseli(abs(n),r);

%  data for this N

nC   = [nC ; N];
maxC = [maxC ; max(abs(real(c)-cx))];

end
end

subplot(2,1,2)
semilogy(nC,maxC,'x')
title('\bf cosine exponential')
xlabel('\bf N-axis');  ylabel('\bf max|c_n^{fft}-c_n^{exact}|')

axis([2 80 10^(-16) 1])