% Continuous least squares approximation of sqrt(x) by polynomials for x in [0,1].
% Construction uses the monomial basis, giving a numerically treacherous Hilbert matrix.

 N = input('polynomial degree? ');

 A = zeros(N+1,N+1);
 f  = zeros(N+1,1);
 fe1 = zeros(N+1,1);
 fe2 = zeros(N+1,1);

 integrand = inline('sqrt(x).*polyval(p,x)','x','p');


 for j=0:N, 
    for k=0:N
       A(j+1,k+1) = 1/(j+k+1);
    end
    f(j+1)   = (1/(j+1.5));                  % integral of sqrt(x)*x^j from 0 to 1
    fe1(j+1) = (1/(j+1.5))*(1+1e-8*randn);   % f(j+1), but with a small error
    fe2(j+1) = quad(integrand,0,1,1e-8,[],[1;zeros(j,1)]); % f(j+1) computed using numerical integration

 end

 p  = A\f; 
 pe1 = A\fe1;
 pe2 = A\fe2;

 figure(1), clf
 xx = linspace(0,1,500);
 plot(xx,sqrt(xx), 'b-','linewidth',2), hold on
 plot(xx,polyval(flipud(p),xx), 'k--','linewidth',2)
 plot(xx,polyval(flipud(pe1),xx), 'r--','linewidth',2)
 plot(xx,polyval(flipud(pe2),xx), 'm--','linewidth',2)
 legend('f(x) = sqrt(x)', 'poly approx, exact <f,x^k>', 'poly approx, perturbed <f,x^k>', 'poly approx, numerically integrated <f,x^k>',2)
 xlabel('x','fontsize',20)
 format long
 fprintf('Polynomial coefficients, monomial basis\n')
 fprintf('exact <f,x^k>, left; perturbed <f,x^k>, middle; numerically integrated <f,x^k> (right)\n')
 disp([p pe1 pe2])