% compare accuracy of discrete least squares polynomial approximations
% computed using the normal equations and the QR algorithm.

 m = 35;                     % approximation at m+1 points
 n = 25;                     % using a degree n polynomial
 x = linspace(-.5,.5,m+1)';  % the approximation points

% construct Vandermonde matrix

 A = zeros(m+1,n+1);
 A(:,1) = ones(m+1,1);
 for j = 1:n
    A(:,j+1) = A(:,j).*x;
 end

% construct right hand side vector
 f = sin(x);

% find polynomial coefficients using normal equations formulation

 p_ne = (A'*A)\(A'*f);

% find polynomial coefficients using QR formulation

 [Q,R] = qr(A,0);
 p_qr = R\(Q'*f);

 fprintf('\n\n  norm(p_ne - p_qr) = %10.7e\n\n', norm(p_ne-p_qr));

% plot accuracy of approximation on a larger interval
% to see the innacuracy in the normal equation solution


 figure(1),clf 
 xx = linspace(-2,2,1000);
 plot(xx,sin(xx),'k-','linewidth',2), hold on
 plot(xx,polyval(flipud(p_ne),xx),'r-','linewidth',2)
 plot(xx,polyval(flipud(p_qr),xx),'b-','linewidth',2)
 axis([min(xx) max(xx) -2 2])
 legend('sin(x)', 'normal equation solution', 'QR solution',2)
 set(gca,'fontsize',16)
 xlabel('x')
 title('Least squares approximation to sin(x) with 35 nodes [-.5,.5]')