% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% gibbsFSE      GIBBS Phenomenon in Fourier Series (Example).  A DEMO
% ____________________________________________________________________________
% 18.385 (Nonlinear Dynamical Systems).         MIT, Fall 1999. R. R. Rosales.
%
% Consider the 2*pi periodic function given in |x| < pi by 
%        F(x) =  0.5 for        |x| < pi/2
%        F(x) = -0.5 for pi/2 < |x| < pi.
% 
% For this function the Fourier Series is
%
%        F(x) = Sum(n=1:infty) (2/(pi*n))*sin(n*pi/2)*cos(n*x)
%
% This sample script calculates this sum -- with very many terms -- in a small
% neighborhood of the discontinuity at x = pi/2. The answer is then plotted.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n Demo will sum and plot the first nnL terms in the Fourier series')
fprintf('\n for a square wave -- in a small neighborhood of a discontinuity.')
fprintf('\n Pick nnL large --- say,  nnL = 5000 --- to see how the overshoot')
fprintf('\n at the discontinuity does not go down to zero  (as the number of')
fprintf('\n terms goes to infinity).    Note also how the width of the "bad"')
fprintf('\n region goes down to zero as 1/nnL.  In fact, the plotting window')
fprintf('\n scales with 1/nnL,   so you will see that the plots are all very')
fprintf('\n similar for nnL large!    This behavior is known by the name of:')
fprintf('\n -----> GIBBS phenomenon. <----- \n')
%
nnL = input(' Enter the number of terms in series to sum. nnL = ');
%       
dz  = 2*pi/nnL/25;  % Want to have at least 25 points per wavelength.
NN  = 500;          % Width of grid on left side of discontinuity is NN*dz.
%                   % Thus will see about 20 oscillations in the sum.
nnLTEXT = num2str(nnL);
z = (pi/2 - NN*dz):dz:(pi/2 + 25*dz);
u = zeros(size(z));
for nn=1:nnL
    k = nn*pi/2;
    u = u + (1/k)*sin(k)*cos(nn*z);
end
figure
whitebg('k')
plot(z, u, '-r', 'LineWidth', 3)
grid on
axis([(pi/2-NN*dz) (pi/2 + 25*dz) 0 0.6])
xlabel('Jump discontinuity from 0.5 to -0.5 at x = pi/2.', ...
       'FontName', 'HelveticaNarrowBold', 'FontSize', 16)
ylabel('Fourier Sum', 'FontName', 'HelveticaNarrowBold', 'FontSize', 16)
title(['Fourier Sum with ', nnLTEXT, ' terms.'], 'FontName', ...
       'HelveticaNarrowBold', 'FontSize', 18)
%
fprintf('\n DONE. The sum is in the array u. The grid is in the array z.\n')
%
clear NN dz k nn nnL nnLTEXT
%
% EOF