clear
reset(0)
echo on
format short
%It is interesting to explore the weights associated
%with Romberg Iteration (Repeated Richardson Extrapolation)
%since it, as with any quadrature rule, can be represented
%as the product of a vector of weights and a vector of
%function evaluations.  The node locations are, of course,
%simple, since they are evenly spaced with some 2^n panels.
%The formula for the weights, however, is very complex and
%is most easily generated numerically.  The code for
%producing these weights, from the trapezoidal rule, through
%the Simpson's rule, and on to the most accurate extrapolation
%is given below:
pause

%To accomplish this, we must produce the trapezoidal rule
%estimate resulting from a function which is unity at one
%of the nodes and zero at all the rest.  We will then 
%subject this array to Romberg iteration, and determine
%the weight for that particular non-zero node.  By repeading
%this for all 1+2^n nodes, we may construct a complete set
%of weights.

pause

n=4; %We set the discretization.  There will be 1+2^n points.

total=zeros(n+1,1+2^n);
for i=1:n+1;
 total(n+2-i,[1:2^(i-1):1+2^n])=2^(i-n-1);
 echo off
end
echo on
total(:,1)=total(:,1)/2;
total(:,1+2^n)=total(:,1+2^n)/2;

weights=zeros(n+1,2^n+1);
for kk=1:2^n+1;
 t=zeros(n+1,n+1);
 t(:,1)=total(:,kk);
 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
 weights(:,kk)=t(n+1,:)';
end
echo on

pause

%We may examine these weights.  Each of the rows of the matrix
%weights correspond to progressively more accurate quadrature
%rules: The first is the Trapezoidal rule, the second is
%Simpson's rule, and so forth.  The final row is the one we really
%want.  It is more convenient to look at the transpose.  We also
%multiply by the inverse of the panel width:

weights'*2^n

pause

%We can examine the properties of these weights.  First we examine
%the sum of the weights:

sum(weights')'

%which yeilds a vector of ones - for each quadrature rule all of
%the weights sum to unity.  This makes sense, since it would have
%to be true if integrating constant functions (e.g., f(x)=2, for
%example) is to be correct.  Note that if we integrate over [a,b]
%we would multiply the weights by the width of the interval (b-a).

pause

%Next we look at the norm of the weights.  This is important in
%examining error propagation in the integration of noisy data.
%Assuming that the function evaluation error (data noise) is
%independent, the relevant item of interest is the 2-norm:

norms=zeros(1,n+1);
for i=1:n+1
 norms(i)=norm(weights(i,:));
 echo off
end
echo on

norms

%As you can see, the norm of the TR rule is a bit less than that
%of Simpson's rule, but the norm rapidly approaches a constant.
%Thus there is no real reason why you shouldn't continue the
%extrapolation process, as it does not significantly increase
%the contribution of the random error, and does decrease the
%algorithm error.

pause

%To finish, we shall look at the accuracy of the method again.
%We apply it to two functions, one which is well behaved and
%one which isn't.  First, we look at the integral of exp(x) 
%over the domain [0:2].  Thus:

a=0;
b=2;
x=a+[0:1/2^n:1]'*(b-a);
%which gives the node locations.  The exact integral is just
%exp(2)-1, so we can calculate the error:

format short e

error=(b-a)*weights*exp(x)-(exp(2)-1)

%and we can also plot this up:
pause

figure(1)
semilogy(abs(error))
xlabel('iteration number')
ylabel('error')
set(gca,'FontSize',18)
pause

%In the above example, the error dropped off rapidly.  That is
%because the function was well behaved.  If the function is
%NOT well behaved, the error is much larger.  An example of
%this is the integral of x*log(x) over the domain [0,1]:

a=0;
b=1;
x=a+(b-a)*[0:1/2^n:1]';
f=x.*log(x);
f(1)=0; %Note that we have to "fix" things at x=0.

integ=(b-a)*weights*f
%And now for the error (the exact value is -1/4):

error=integ+1/4

%which doesn't die off at all!  We can also plot this
%up:
pause

figure(1)
hold on
semilogy(abs(error),'g')
hold off
pause

%So in conclusion, if you are going to invoke Simpson's
%rule, you might as well go all the way on the iteration
%procedure.  The algorithm error will be basically zero
%if all derivatives are well behaved in the domain, and
%the 2-norm of the weights doesn't increase significantly.
%If the function is poorly behaved (e.g., the derivatives
%blow up) then there isn't much benefit to continuing the
%iterative procedure.

echo off