clear
reset(0)
echo on
%This is an example of three types of quadrature
%rules.  We first use a quadrature rule where f(x)
%is evaluated only on the left edge of each panel.
%We shall evaluate the integral for several dif-
%ferent discretizations to determine the order of
%the rule.  We apply the code for f(x)=exp(x).
pause
%first we set up the interval
a=1;
b=3;
m=5;%This is the number of different discretizations
for i=1:m;
 %next we set the panel size
 n=2^(i-1);
 dx=(b-a)/n;
 x=[a:dx:b];
 %now for the weighting function
 w=ones(1,n);
 w(n+1)=0;
 %and now the result:
 f=exp(x);
 integral(i)=dx * w * f';
 h(i)=dx;
 echo off
end;
echo on
%we calculate the error:
error=integral-(exp(b)-exp(a));
loglog(h,abs(error))
set(gca,'FontSize',18)
xlabel('panel width')
ylabel('quadrature error')
pause
%
%Now we do it again for the trapezoidal rule
for i=1:m,
 %next we set the panel size
 n=2^(i-1);
 dx=(b-a)/n;
 x=[a:dx:b];
 %now for the weighting function
 w=ones(1,n);
 w(n+1)=.5;
 w(1)=.5;
 %and now the result:
 f=exp(x);
 integral(i)=dx * w * f';
 trh(i)=dx;
 echo off
end;
echo on
%we calculate the error:
trerror=integral-(exp(b)-exp(a));
loglog(h,abs(error),trh,abs(trerror))
set(gca,'FontSize',18)
xlabel('panel width')
ylabel('quadrature error')
pause
%And, finally, for Simpson's rule:
for i=1:m,
 %next we set the panel size
 n=2^(i-1);
 dx=(b-a)/n;
 x=[a:dx:b];
 %now for the weighting function
 w=2*ones(1,n)/3.;
 w(2:2:n)=w(2:2:n)+2./3.;
 w(n+1)=1/3.;
 w(1)=1/3.;
 %and now the result:
 f=exp(x);
 integral(i)=dx * w * f';
 sih(i)=dx;
 echo off
end;
echo on
%we calculate the error:
sierror=integral-(exp(b)-exp(a));
%Finally, we compare the errors.
loglog(h,abs(error),trh,abs(trerror),sih,abs(sierror))
set(gca,'FontSize',18)
xlabel('panel width')
ylabel('quadrature error')
pause
%As you can see, the slope of the first
%algorithm was unity, giving an error of O(h),
%the slope of the trapezoidal rule is of O(h^2),
%and the slope of Simpson's rule is O(h^4).
echo off