%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)/sqrt{1-t^2}dt. The polynomials are computed by computing the projections on the orthogonal basis of Chebyshev polynomials; the integrals are evaluated using a simple trapezoidal rule


%grid size
pi=acos(0)*2;
N=10000;
h=1/N;
%x grid
x=[-1+h/2:h:1-h/2];
%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
 %the Chebychev polynomial of order j
 p=cos(j*acos(x));
 %the weighted scalar product
 a(j+1)=sum(p.*f./sqrt(1-x.^2))/N;
 %normalization by the norm
 a(j+1)=2*a(j+1)/pi;
 %the norm is different for the 0th polynomial
 if (j==0)
  a(j+1)=a(j+1)/2;
 end
 %now compute the approximating polynomial
 q=q+a(j+1)*p;
 figure(3);
   plot(x,f,x,q)
   axis([-1 1 -0.2 1.2])
   drawnow;
   pause(1);
end
figure(3);
  plot(x,f,x,q)
  axis([-1 1 -0.2 1.2])

figure(4);
  plot(x,f-q)