clear
reset(0)
echo on
%This is an example of repeated Richardson
%extrapolation.  The algorithm is based on
%using evaluations of a function to create
%a high polynomial degree quadrature rule.
%For example, since we know that the trapezoidal
%rule is O(h^2) in accuracy, we can combine
%a trapezoidal rule with one at twice the
%discretization to get a Simpson's rule.  In
%the same way, we can combine two Simpson's
%rules to get a new rule with error O(h^6), and
%so on.
pause
%The algorithm is as follows:
%T1
%T2  S2
%T4  S4  P4
%T8  S8  P8  Q8
%T16 S16 P16 Q16 R16
%
%and so forth.  The idea is that the one and
%two panel trapezoidal rules T1 and T2 combine
%to form the two panel Simpson's rule S2.  In
%the same manner, the four and eight panel
%Simpson's rules combine to form the eight panel
%rule at the next higher order.  This algorithm
%proceeds until you have gone as far as you can,
%which depends on the number of original function
%evaluations you made.
pause
%Let's implement this technique.  We will go up
%to the fifth level of iteration.  We start by
%constructing the trapezoidal rules.  We will
%operate on the function exp(x) over the interval
%[-1,1].  Thus:
n=5;
a=-1;
b=1;
h=(b-a)/2^n;
fname='exp';

t(1,1)=(feval(fname,a)+feval(fname,b))/2;
for i=2:n+1,
 t(i,1)=t(i-1,1);
 for j=1:2^(i-2),
  t(i,1)=t(i,1)+feval(fname,(a+h*2^(n-i+1)+h*(j-1)*2^(n-i+2)));
  echo off
 end
end
echo on
for i=1:n+1,
 t(i,1)=t(i,1)*h*2^(n+1-i);
echo off
end
echo on
pause
%This algorithm for the trapezoidal rule was
%constructed so that we made the minimum possible
%number of function calls.  Remember that the
%points used in determining the n point rule are
%reused in calculating the 2n point rule.  Thus
%at each step only the "in between" function
%evaluations are needed.  If you study the structure
%of the above expression you will see that this is,
%in fact, accomplished.  No further function 
%evaluations are required for the repeated Richardson
%extrapolation.
pause
%Now that we have constructed the trapezoidal rule
%estimates, let's plot the errors up:
semilogy(abs(t(:,1)-(exp(1)-exp(-1))))
set(gca,'FontSize',18)
%As you can see, the error is cut by a factor of
%four each time the interval is cut in half.
pause
%Let's now do the repeated Richardson extrapolation.
%Each succeeding rule will be of order h^2 higher
%than the one used to calculate it.  Thus we calculate
%the 2n point Simpson's rule from the formula
%S2n = (4*T2n - Tn)/(4-1) and the next higher formula
%by P2n = (16*S2n - Sn)/(16-1) and so forth:

for j=2:n+1
 for i=n+1:-1:j
  t(i,j)=(4^(j-1)*t(i,j-1)-t(i-1,j-1))/(4^(j-1)-1);
  echo off
 end
end
echo on

%and that's it!.  The matrix t is in lower triangular
%form, and contains the trapezoidal rule, Simpson's
%rule, and higher order quadrature estimates.  Let's
%look at this matrix:
pause
t
pause
%As you can see, the estimates rapidly converge to
%the correct value.  It is of more interest to plot up
%the errors of the last row.  This row contains the
%2^n quadrature estimates, the first column being the
%trapezoidal rule result, the second being the Simpson's
%rule result, and so on.  Thus:
pause
semilogy(abs(t(n+1,:)-(exp(1)-exp(-1))))
hold on
semilogy(abs(t(:,1)-(exp(1)-exp(-1))))
set(gca,'FontSize',18)
hold off
%Thus, after the fifth iteration we have a rule which
%is of order h^12, and the error is on the order of the
%roundoff error, at least for this smooth function.
%There is not much point going any further, as we become
%limited by the double precision accuracy of the computer.
pause
%In conclusion, repeated Richardson extrapolation is a
%simple and useful approach to making quadrature estimates
%if you don't happen to have a table of the Gaussian 
%quadrature weights and nodes handy.  The difference between
%the last two elements of the diagonal provide an estimate
%of the quadrature error as well.  A G7K15 Gauss-Kronod pair
%would achieve a somewhat higher level of accuracy with an
%error estimate as well, and with a bit fewer than half the
%function evaluations used here.
echo off