%computes the sequence of approximations to a function f by polynomials of maximal degree j so that the distance between f and the polynomial is minimal in the norm induced by the scalar product (f,g)=int_{-1}^1 f(t)g(t)dt. The polynomials are computed by computing the projections on the orthogonal basis of Legendre polynomials; the integrals are evaluated using a simple trapezoidal rule

%grid size
N=10000;
h=1/N;
%x grid
x=[-1:h:1];
%function we want to approximate
f=1./(1+x.^2*25);
%f=x.^5-3;
plot(x,f)
%order of approximation
m=10;
%compute coefficients
q=0;a=0;
for j=0:m
  %this gives a vector with entries given by the j'th Legendre polynomial at the points x
  p=legendre(j,x);p=p(1,:); 
  %compute coefficients
  %this computes the scalar product between the Legendre polynomial and our function f
  a(j+1)=sum(p.*f)/N-0.5*(p(1)*f(1)-p(N)*f(N))/N;
  %normalize by the norm
  a(j+1)=a(j+1)/2*(2*j+1);
  %compute the next approximation
  q=q+a(j+1)*p;
  %plot both f and approximation p
  figure(1);
  plot(x,f,x,q)
  axis([-1 1 -0.2 1.2])
  drawnow;
%    pause(1);
end
figure(1);
  plot(x,f,x,q)
  axis([-1 1 -0.2 1.2])
  drawnow;
  hold
  
figure(2);
  plot(x,f-q)
  hold



