echo off;
% NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
% To accompany the text:
% NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
% Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
% This free software is complements of the author.

% Algorithm 7.3 (Recursive Trapezoidal Rule).
% Section 7.3, Recursive Rules and Romberg Integration, Page 378
echo on; clc; format long; hold off; clear

% This program implements the recursive trapezoidal rule.

% The function f(x) is integrated over the interval [a,b].

% Define and store the function f(x) in the M-file f.m
% function y = f(x)
% y = 1 + exp(-x).*sin(4*x);

delete f.m
diary f.m; disp('function y = f(x)');...
disp('y = 1 + exp(-x).*sin(4*x);');...
diary off;

f(0); % Remark. f.m rctrap.m grpoly.m are used for Algorithm 7.3

pause % Press any key to continue.

clc;, clg;
% Plot f(x) over the interval [a,b].

clc; clg;
a = 0;
b = 1;
h = (b-a)/200;
x0 = a:h:b;
y0 = f(x0);
c = 0;
d = 2;
axis([a b c d]);...
plot(x0,y0,'g');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
title('The curve y = f(x) over [a,b].');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

clc;

% The recursive trapezoidal rule is applied over [a,b].

% Place the endpoints of [a,b] in a and b.

% Place the number of times to bisect [a,b] in N.

a = 0;
b = 1;
n = 6;

T = rctrap('f',a,b,n);

pause % Press any key to continue.

clc;
for i=0:n, J(i+1) = 2^i; end
values = [J;T];
Mx1 = 'The recursive trapezoidal rule results: ';
Mx2 = ' # subintervals trapezoidal rule';
clc,echo off,diary output,...
disp(''),disp(Mx1),disp(''),disp(Mx2),disp(values'),diary off,echo on

pause % Press any key to see a nice graph y = f(x).

clg;
c = 0;
d = 2;
axis([a b c d]);...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
hold on;...
for i=1:(n+1),
if J(i) <= 16,
grpoly(a,b,J(i),1);
hold on;...
end
end;...
xlabel('x');...
ylabel('y');...
title('The recursive trapezoidal rule.');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

Mx1 = 'The recursive trapezoidal rule was applied with ';
Mx2 = [Mx1,num2str(n),' recursions.'];
Mx3 = 'The approximate value of the integral is:';
clc,echo off,diary output,disp(''),disp(Mx2),...
disp(''),disp(Mx3),disp(T(n+1)),diary off, echo on