%
%  w03fourier_error.m -- djm -- 25 jan 2008 
%

%  interval (-L < x < +L)
L = pi;

%  data vectors
Nsave = [];
Esave = [];

%  parameter loop
for jj = 4:1:14
    Nterms = 2^jj;  Npts = 128;

    %  plotting points
    dx = 2*L/Npts;  x = -L:dx:L;

    %
    %  NOTE:  ".^" is an example of element-wise arithmetic
    %
    %  exact function
    Fexact = x.^2;

    %  summation loop
    %  j = 0 term
    Fser = (L^2/3)*ones(size(x));
    %  j = 1,3,5 ...
    for j = 1:1:Nterms
        Fser = Fser + (4*L^2*(-1)^j/pi^2/(j^2))*cos((j*pi/L)*x);
    end

    max_error = max(abs(Fexact-Fser));
    
    [Nsave] = [Nsave Nterms];
    [Esave] = [Esave max_error];
end

%  error vs N plots
figure(2);  clf
subplot(2,1,1)
plot(Nsave,Esave,'kx')
title('\bf max error vs N','fontsize',14);
xlabel('\bf N','fontsize',12);
ylabel('\bf max|f(x)-f_S(x)|','fontsize',12);

%  log-log plot is so much more useful
subplot(2,1,2)
plot(log(Nsave),log(Esave),'kx');  hold on
title('\bf log-log error plot','fontsize',14);
xlabel('\bf log(N)','fontsize',12);
ylabel('\bf log(max|f(x)-f_S(x)|)','fontsize',12);

coeff = polyfit(log(Nsave),log(Esave),1);
plot(log(Nsave),coeff(1)*log(Nsave)+coeff(2),'r--')

text(log(Nsave(end))-2,log(Esave(end))+4,['\bf slope = ' num2str(coeff(1))])

slope = coeff(1)
