function ch07_composite_trapezoidal_simpson
%==========================================================================
% ch07_composite_trapezoidal_simpson.m
%
% Composite trapezoidal and Simpson rules for numerical integration.
%
% The program evaluates
%
%       integral sin(x) dx,   0 <= x <= pi,
%
% for successively finer uniform grids. The exact value is 2.
%
% The example illustrates the different convergence rates of the
% composite trapezoidal and Simpson rules as the number of subintervals
% is doubled.
%
% The Simpson routine is based on code written by Alexander Godunov,
% October 2009. The companion trapezoidal routine and MATLAB presentation
% were prepared for the book companion website, 2026.
%
% A plot of the integrand (or integration domain) is shown before the
% numerical calculation to make the character of the problem visible.
%==========================================================================

a = 0.0;
b = pi;
exact = 2.0;

% === Plot the integrand =================================================
Nplot = 1000;
xplot = linspace(a,b,Nplot+1);
yplot = f(xplot);

figure('Name','Integrand: sin(x)')
plot(xplot,yplot,'LineWidth',1.5)
xlabel('x')
ylabel('f(x)')
title('Integrand for composite trapezoidal and Simpson rules')
grid on
% === End plot ===========================================================

fprintf('Composite trapezoidal and Simpson rules\n');
fprintf('Integral of sin(x) from 0 to pi\n');
fprintf('Exact value          = %15.7e\n\n', exact);

fprintf('        n       Trapezoidal          Simpson\n');
fprintf('------------------------------------------------\n');

n = 2;

for k = 1:16
    trap = composite_trapezoidal(@f, a, b, n);
    simp = composite_simpson(@f, a, b, n);

    fprintf('%9d%19.7e%19.7e\n', n, trap, simp);

    n = 2*n;
end

end


function y = f(x)
%--------------------------------------------------------------------------
% Test function for the integration example.
%
% To use another integrand, replace only the line below and change the
% integration limits in the main function if needed.
%--------------------------------------------------------------------------
y = sin(x);
end


function integral = composite_trapezoidal(fun, a, b, n)
%--------------------------------------------------------------------------
% Composite trapezoidal rule on [a,b].
%
% Input:
%   fun - function to integrate
%   a,b - integration limits
%   n   - number of equal subintervals
%
% Output:
%   integral - numerical approximation to the integral
%--------------------------------------------------------------------------

if n < 1
    error('composite_trapezoidal: n must be positive.');
end

h = (b-a)/n;

sum_value = 0.5*(fun(a) + fun(b));

for i = 1:n-1
    x = a + i*h;
    sum_value = sum_value + fun(x);
end

integral = h*sum_value;

end


function integral = composite_simpson(fun, a, b, n)
%--------------------------------------------------------------------------
% Composite Simpson rule on [a,b].
%
% Input:
%   fun - function to integrate
%   a,b - integration limits
%   n   - number of equal subintervals; n must be even
%
% Output:
%   integral - numerical approximation to the integral
%--------------------------------------------------------------------------

if n < 2 || mod(n,2) ~= 0
    error('composite_simpson: n must be a positive even integer.');
end

h = (b-a)/n;

sum_odd = 0.0;
for i = 1:2:n-1
    x = a + i*h;
    sum_odd = sum_odd + fun(x);
end

sum_even = 0.0;
for i = 2:2:n-2
    x = a + i*h;
    sum_even = sum_even + fun(x);
end

integral = h*(fun(a) + fun(b) + 4.0*sum_odd + 2.0*sum_even)/3.0;

end