%
% This code demonstrates Gibbs phenomena. That is, interpolation
% using a polynomial of higher degree does not necessarily yield
% a better result.
%
clear
figure(1);
set(gcf,'position',[50,50,600,600])
clf reset
%
% We first plot the function f = 1/(1+x^2).
%
x0=[-5:0.01:5];
f0=1./(1+x0.^2);
n0=size(x0,2);
plot(x0,f0,'k-','linewidth',2.0)
hold on
%
% interpolation polynomial of degree 4 through (4+1) points
%
p=4;
dx=10/p;
x=[-5:dx:5];
f=1./(1+x.^2);
yp=zeros(1,n0);
for i=1:p+1,
 prd=ones(1,n0);
 for j=1:p+1,
  if j~=i
   prd=prd.*(x(j)-x0)/(x(j)-x(i));
  end
 end
 yp=yp+f(i)*prd;
end
y4=yp;
%
% interpolation polynomial of degree 8 through (8+1) points
%
p=8;
dx=10/p;
x=[-5:dx:5];
f=1./(1+x.^2);
yp=zeros(1,n0);
for i=1:p+1,
 prd=ones(1,n0);
 for j=1:p+1,
  if j~=i
   prd=prd.*(x(j)-x0)/(x(j)-x(i));
  end
 end
 yp=yp+f(i)*prd;
end
y8=yp;
%
% interpolation polynomial of degree 12 through (12+1) points
%
p=12;
dx=10/p;
x=[-5:dx:5];
f=1./(1+x.^2);
yp=zeros(1,n0);
for i=1:p+1,
 prd=ones(1,n0);
 for j=1:p+1,
  if j~=i
   prd=prd.*(x(j)-x0)/(x(j)-x(i));
  end
 end
 yp=yp+f(i)*prd;
end
y12=yp;
%
plot(x0, y4,'k--','linewidth',1.0)
plot(x0, y8,'b--','linewidth',1.0)
plot(x0, y12,'r-','linewidth',1.0)
%
axis([-5.0,5.0,-4.0,1.5])
set(gca,'fontsize',14)
set(gca,'xtick',[-4:2:4])
set(gca,'ytick',[-4:1:1])
xlabel('x')
ylabel('f')
title('Gibbs phenomena')
h1=legend('f(x)=1/(1+x^2)','p_{4}(x) through 5 points',...
'p_{8}(x) through 9 points','p_{12}(x) through 13 points',0);
set(h1,'fontsize',12)
%
%
