function rich_extr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Richardson extrapolation on the example of 
% calulcation \int_{a}^{b} f(x)dx using the composite trapezoidal rule.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=0; b=0.99; % the interval [a,b]
% Calculate the exact value of the integral
%exact=(-cos(3*b)+cos(3*a))/3;
exact=asin(b)-asin(a);
% Calculate the approximation by trapezoidal rule with step sizes h(=(a-b)/10), h/2, h/10
S=comp_trapez(a,b,100);
Shalf=comp_trapez(a,b,200);
Stenth=comp_trapez(a,b,1000);
Stenthhalf=comp_trapez(a,b,2000);
%compute Richardson extrapolation
M=Shalf+(Shalf-S)/3;
Mtenth=Stenthhalf+(Stenthhalf-Stenth)/3;
% compare results: 
fprintf(1,'Trapezoidal rule, number of nodes %4i, Integral: %13.10f Error: %13.10f \n', 100, S, exact-S);
fprintf(1,'Trapezoidal rule, number of nodes %4i, Integral: %13.10f Error: %13.10f \n', 200, Shalf, exact-Shalf);
fprintf(1,'Richardson extrp, number of nodes %4i, Integral: %13.10f Error: %13.10f \n', 200, M, exact-M);
fprintf(1,'Trapezoidal rule, number of nodes %4i, Integral: %13.10f Error: %13.10f \n', 1000, Stenth, exact-Stenth);
fprintf(1,'Trapezoidal rule, number of nodes %4i, Integral: %13.10f Error: %13.10f \n', 2000, Stenthhalf, exact-Stenthhalf);
fprintf(1,'Richardson extrp, number of nodes %4i, Integral: %13.10f Error: %13.10f \n \n \n', 2000, Mtenth, exact-Mtenth);

% the composite trapezoidal rule:
% comp_trapez(a--left end,b--right end,n -- number of nodes)
function S=comp_trapez(a,b,n);
h=(b-a)/n;   %the step size
%%%%% THE Composite Trapezoidal RULE %%%%%
%%%%% TO CALCULATE \int_{a}^{b} f(x)dx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S=f(a);
for i=1:1:n-1
  S=S+2*f(a+i*h);
end;
S=S+f(b);
S=(h/2)*S;

function y=f(x);
%y=sin(3*x);
y=1/sqrt(1-x.^2);