% Matlab code polyfit31.m
% For "Applied Numerical Linear Algebra",  Figure 3.1
% Written by James Demmel, Feb  9, 1997
%                Modified, Jun  2, 1997
%
% Produce plot for figure 3.1 of book
% 
y = (-5:.5:6)';
b = sin(pi*y/5) + y/5;
for i=1:20,
  A = ones(size(y));
  for deg = 1:i, A = [A,y.^deg]; end
  x = A\b;
  e = A*x-b;
  err(i)=norm(e,2);
end
figure(1),

hold off, clf

subplot(121),
cc= ['k','r','c','b'];
% handl=plot(y,b,'k+'),
% set(handl,'MarkerSize',8)
axis([-10 10 -3 3])
axis('square')
title('Original data (o), and polynomial fits')
hold on
clr = 0;
yplot = (-10:.05:13)';
% for i=[2 3 5 18],
for i=[1 3 6 19],
  clr = clr + 1;
  A = ones(size(y));
  for deg = 1:i, A = [A,y.^deg]; end
  x = A\b;
  Aplot = ones(size(yplot));
  for deg = 1:i, Aplot = [Aplot,yplot.^deg]; end
  fit = Aplot*x;
  handl=plot(yplot,fit,cc(clr));
  set(handl,'LineWidth',2);
end
handl=plot(y,b,'ko');
set(handl,'MarkerSize',7);
grid
text(-2,2.7, 'deg=  1'), 
text(-2,2.45,'deg=  3'), 
text(-2,2.2, 'deg=  6'), 
text(-2,1.95,'deg= 19'), 
where = [-4, -2];
handl=plot(where,2.71*[1 1],cc(1));
set(handl,'LineWidth',2);
handl=plot(where,2.46*[1 1],cc(2));
set(handl,'LineWidth',2);
handl=plot(where,2.21*[1 1],cc(3));
set(handl,'LineWidth',2);
handl=plot(where,1.95*[1 1],cc(4));
set(handl,'LineWidth',2);
  
subplot(122),
handl=semilogy((1:20),err,'k+');
set(handl,'MarkerSize',7);
axis([0 21 1e-13 10])
axis('square')
xlabel('Degree of fitted polynomial')
title('Least squares residual norm')
grid