Below are the MATLAB scripts to illustrate the composite Trapezoidal rule and the improvement of the convergence by extrapolation.
Save this as comptrap.m:
function t = comptrap(f,a,b,n)
%
% function t = comptrap(f,a,b,n)
%
% Applies the composite trapezoidal rule to approximate the integral
% of the function f over the interval [a,b], using n subintervals.
%
% Example:
% t = comptrap('cos',0,pi/2,100)
%
% Note: compare this to the result of MATLAB built-in command quad:
% quad('cos',0,pi/2)
%
h = (b-a)/n;
t = (feval(f,a) + feval(f,b))/2;
for i = 1:n-1
t = t + feval(f,a+i*h);
end;
t = t*h;
We apply the composite trapezoidal rule "adaptively", running adaptrap:
function t = adaptrap(f,a,b,n)
%
% function t = adaptrap(f,a,b,n)
%
% Returns a vector of n approximations for the definite integral of
% the function f over the interval [a,b], the i-th entry in the vector
% uses 2^i function evaluations.
%
% Example:
% t = adaptrap('cos',0,pi/2,10);
%
t = zeros(n,1);
h = (b-a); % size of subinterval
m = 1; % number of subintervals
t(1) = (feval(f,a) + feval(f,b))*h/2;
for i = 2:n
h = h/2;
for j=0:m-1
t(i) = t(i) + feval(f,a+h+j*2*h);
end;
t(i) = t(i-1)/2 + h*t(i);
m = 2*m;
end;
The order of convergence is O(h^2), which means that we need one million function evaluations, or h = 10^(-6), to get about twelve decimal places correct. This is illustrated below on the computation of cos(x) over the interval [0,pi/2]:
>> t = adaptrap('cos',0,pi/2,20)
t =
0.78539816339745
0.94805944896852
0.98711580097278
0.99678517188617
0.99919668048507
0.99979919432002
0.99994980009210
0.99998745011753
0.99999686253529
0.99999921563419
0.99999980390857
0.99999995097714
0.99999998774429
0.99999999693607
0.99999999923402
0.99999999980851
0.99999999995213
0.99999999998802
0.99999999999699
0.99999999999925
We can accelerate the convergence by extrapolation:
function et = romberg(t)
%
% function et = romberg(t)
%
% Applies extrapolation to the approximations in t,
% t is a sequence of results of the composite Trapezoidal rule.
%
% Example:
% t = adaptrap('cos',0,pi/2,6)
% et = romberg(t);
%
n = size(t,1);
et = zeros(n,n);
et(:,1) = t;
for i = 2:n
m = 2;
for j = 2:i
r = 2^m;
et(i,j) = (r*et(i,j-1) - et(i-1,j-1))/(r - 1);
m = m+2;
end;
end;
The results are quite amazing: we get the integral accurate up to full precision with only 64 subintervals:
>> t = adaptrap('cos',0,pi/2,6)
t =
7.8540e-01
9.4806e-01
9.8712e-01
9.9679e-01
9.9920e-01
9.9980e-01
>> et = romberg(t)
et =
7.8540e-01 0 0 0 0 0
9.4806e-01 1.0023e+00 0 0 0 0
9.8712e-01 1.0001e+00 9.9999e-01 0 0 0
9.9679e-01 1.0000e+00 1.0000e+00 1.0000e+00 0 0
9.9920e-01 1.0000e+00 1.0000e+00 1.0000e+00 1.0000e+00 0
9.9980e-01 1.0000e+00 1.0000e+00 1.0000e+00 1.0000e+00 1.0000e+00
>> et-ones(6,6) % exact value of the integral is 1, here are the errors
ans =
-2.1460e-01 -1.0000e+00 -1.0000e+00 -1.0000e+00 -1.0000e+00 -1.0000e+00
-5.1941e-02 2.2799e-03 -1.0000e+00 -1.0000e+00 -1.0000e+00 -1.0000e+00
-1.2884e-02 1.3458e-04 -8.4345e-06 -1.0000e+00 -1.0000e+00 -1.0000e+00
-3.2148e-03 8.2955e-06 -1.2377e-07 8.1440e-09 -1.0000e+00 -1.0000e+00
-8.0332e-04 5.1668e-07 -1.9046e-09 2.9837e-11 -1.9831e-12 -1.0000e+00
-2.0081e-04 3.2265e-08 -2.9646e-11 1.1480e-13 -1.7764e-15 2.2204e-16
Form the table with errors above we see that the differences between the exponents in two consecutive columns is about 2.