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

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

%  parameters
Nterms = 21;  Npts = max(4*Nterms,64);

%  x is a vector of 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,2,3, ...
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));

%  plot:  f(x) vs f_S(s) and error
figure(1);  clf
subplot(2,1,1)
plot(x,Fexact,'r');  hold on
plot(x,Fser,'b');
plot(x,Fser,'k.');
axis([-L L min(Fser)-0.2 max(Fser)+0.2])
title(['\bf fourier series f(x) = x^2:  up to j=' num2str(Nterms) ' terms'],'fontsize',14)
xlabel('\bf x-axis','fontsize',12)
ylabel('\bf partial sum, S_j(x)','fontsize',12)

subplot(2,1,2);
plot(x,Fexact-Fser,'b');  hold on
plot(x,Fexact-Fser,'k.')
axis([-L L min(Fexact-Fser)-0.2 max(Fexact-Fser)+0.2])
title(['\bf fourier series error:  max error = ' num2str(max_error)],'fontsize',14)
xlabel('\bf x-axis','fontsize',12)
ylabel('\bf partial sum, f(x)-S_j(x)','fontsize',12)
