% Code to demonstrate computation of continuous least squares approximation.
% Uses MATLAB's built-in codes to compute inner products.

% Use the weight function w(x) = 1 on the interval [-1,1].
% Construct the orthogonal polynomials for this weight, interval.
% These are the Legendre polynomials; one can look up their coefficients
% in mathematical tables.  We input them in MATLAB's standard format
% for polynomials.  (We have normalized the standard Legendre polynomials.)

Leg = [[     0     0       0      0      0      0     1]*sqrt(1/2);  % p_0(x) 
       [     0     0       0      0      0      1     0]*sqrt(3/2);  % p_1(x)
       [     0     0       0      0    3/2      0  -1/2]*sqrt(5/2);  % p_2(x)
       [     0     0       0    5/2      0   -3/2     0]*sqrt(7/2);  % p_3(x)
       [     0     0    35/8      0  -15/4      0   3/8]*sqrt(9/2);  % p_4(x)
       [     0  63/8       0  -35/4      0   15/8     0]*sqrt(11/2); % p_5(x)
       [231/16     0 -315/16      0 105/16      0 -5/16]*sqrt(13/2)];% p_6(x)

% All the necessary integrals have integrands that are the product of our
% target function f(x) = exp(x)*sin(5*x) and some polynomial.
% The following inline function defines this general form of integrand.
% f = inline('sin(pi*x) + 3*exp(-(50*(x-.5)).^2)'); 
 f = inline('sin(pi*x)'); 
 integrand = inline('feval(f,x).*polyval(p,x)','x','f','p');

% We also include a function to evaluate the 2-norm of the error
 errintegrand = inline('(feval(f,x)-polyval(p,x)).*polyval(q,x)','x','f','p','q');

% compute the expansion coefficients for the optimal polynomial approximation
 x = linspace(-1.1,1.1,1000)';
 figure(1),clf
 plot(x,f(x),'b-','linewidth',3), hold on
 axis([-1.1 1.1 -2 5])
 set(gca,'fontsize',20)
 drawnow 
 px = zeros(1,size(Leg,1));

 clear pxplt
 for j=1:size(Leg,1)
     input('press return to continue')
     c(j) = quad(integrand,-1,1,1e-10,[],f,Leg(j,:));
     px = px + c(j)*Leg(j,:);
     fprintf(' c_%d = %10.7f \n', j-1, c(j))
     if exist('pxplt','var'), set(pxplt,'linewidth',1); end
     pxplt = plot(x,polyval(px,x),'r-','linewidth',3);
     fprintf('Test orthogonality of error: <f-p, 1> = %10.7e\n', quad(errintegrand,-1,1,1e-10,[],f,px,[1]))
     title(sprintf('Degree %d Least-Squares Approximation',j),'fontsize',20)
 end