% This program analyzes convergence of a particular Fourier series

% Function:
% f(x) = 1-2*|x-0.5| on [0,1]
% Fourier coefficients
% F[2n] = 0 , F[2n+1] = -2/[pi^2 (2n+1)^2] , F[0]=0.5

% Construct finite sums
%digits(100);
N = 200;
N0 = 20;
stepx = 1e-6;
x = [0:stepx:1]; % row vector
f = 1-2*abs(x-0.5);
PartialSum = zeros(N+1,length(x));
L2error = zeros(N+1,1);
for k=0:N
    if k == 0
        PartialSum(k+1,:) = 0.5*ones(1,length(x));
    elseif mod(k,2) == 0
        % Don't do anything: F[2p ] =0
        PartialSum(k+1,:)=PartialSum(k,:);
    else
        Fk = -2/(pi*pi*k*k);
        Fminusk = -2/(pi*pi*(-k)*(-k));
        PartialSum(k+1,:) = PartialSum(k,:) + Fk*exp(2*pi*i*k*x) + Fminusk * exp(-2*pi*i*k*x);
    end
    L2error(k+1) = stepx*norm(f-PartialSum(k+1,:))^2;
end


% Plot results
figure(1)
plot(x,real(PartialSum(N0+1,:)),'r',x,f,'g')
legend('PartialSum','f')
title(sprintf('N = %d',N0))

figure(2)
plot(x,imag(PartialSum(N0+1,:)))
title(sprintf('Imaginary Part for N = %d',N0))

figure(3)
for k = 0:N0
    plot(x,real(PartialSum(k+1,:)))
    hold on
end

% Error analysis
figure(4)
plot(0:N,L2error)
xlabel('N'),ylabel('||f-sum_k^N c_k e_k ||_2^2')

figure(5)
loglog(0:N,L2error)
xlabel('N'),ylabel('||f-sum_k^N c_k e_k ||_2^2')

figure(6)
plot(log(1:N),log(L2error(2:N+1)))
xlabel('log(N)'),ylabel('log ||f-sum_k^N c_k e_k ||_2^2')
grid on

    