function ch07_adaptive_simpson
%==========================================================================
% ch07_adaptive_simpson.m
%
% Adaptive numerical integration using Simpson's rule.
%
% The routine compares Simpson estimates on one interval and on its two
% equal subintervals. Intervals that do not satisfy the local error test
% are subdivided further.
%
% The algorithm is non-recursive: interval data are stored explicitly
% in arrays that act as a stack.
%
% Example:
%
%             1
%     -----------------
%     1 + 400(x-0.2)^2
%
% integrated from 0 to 1.
%
% Exact value:
%   [atan(16) + atan(4)] / 20
%
% Original adaptive Simpson routine written by Alexander Godunov,
% July 2012.
% MATLAB version prepared for the 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 = 1.0;
eps = 1.0e-8;

exact = (atan(16.0) + atan(4.0))/20.0;

% === Plot the integrand =================================================
% The narrow peak is the main reason this is a useful adaptive example.
Nplot = 2000;
xplot = linspace(a,b,Nplot+1);
yplot = arrayfun(@f,xplot);

figure('Name','Narrow peak integrand')
plot(xplot,yplot,'LineWidth',1.5)
xlabel('x')
ylabel('f(x)')
title('Narrow peak: 1/[1+400(x-0.2)^2]')
grid on
% === End plot ===========================================================

[result, errest, nfun] = adaptive_simpson(@f, a, b, eps);

fprintf('Adaptive Simpson integration\n');
fprintf('Integral of 1/[1+400(x-0.2)^2] from 0 to 1\n');
fprintf('Exact value          = %15.7e\n', exact);
fprintf('Calculated value     = %15.7e\n', result);
fprintf('Estimated error      = %12.4e\n', errest);
fprintf('Function calls       = %10d\n', nfun);

end


function y = f(x)
%--------------------------------------------------------------------------
% Test function with a narrow peak centered near x = 0.2.
%--------------------------------------------------------------------------
y = 1.0/(1.0 + 400.0*(x-0.2)^2);
end


function [result, errest, nfun] = adaptive_simpson(fun, a, b, eps)
%--------------------------------------------------------------------------
% Adaptive non-recursive integration using Simpson's rule.
%
% Input:
%   fun - function to integrate
%   a,b - integration limits
%   eps - requested error tolerance
%
% Output:
%   result - numerical approximation
%   errest - accumulated estimate of the absolute error
%   nfun   - total number of function evaluations
%
% Method:
%   Simpson estimates on the left and right halves are compared with the
%   estimate on the whole interval. The local error estimate is
%
%       |s1+s2-s0|/15.
%
%   Intervals that fail the test are placed on an explicit stack.
%--------------------------------------------------------------------------

im = 32;

tol   = zeros(1,im);
x     = zeros(1,im);
h     = zeros(1,im);
fa    = zeros(1,im);
fm    = zeros(1,im);
fb    = zeros(1,im);
s     = zeros(1,im);
level = zeros(1,im);

result = 0.0;
errest = 0.0;
i = 1;

x(i) = a;
h(i) = (b-a)/2.0;
fa(i) = fun(a);
fm(i) = fun(a+h(i));
fb(i) = fun(b);
tol(i) = 15.0*eps;
level(i) = 1;

s(i) = h(i)*(fa(i) + 4.0*fm(i) + fb(i))/3.0;
nfun = 3;

while i > 0

    f1 = fun(x(i) +       h(i)/2.0);
    f3 = fun(x(i) + 3.0*h(i)/2.0);
    nfun = nfun + 2;

    s1 = h(i)*(fa(i) + 4.0*f1 + fm(i))/6.0;
    s2 = h(i)*(fm(i) + 4.0*f3 + fb(i))/6.0;

    x0   = x(i);
    f0   = fa(i);
    f2   = fm(i);
    f4   = fb(i);
    step = h(i);
    err  = tol(i);
    s0   = s(i);
    deep = level(i);

    i = i - 1;

    if abs(s1+s2-s0) <= err

        result = result + s1 + s2;
        errest = errest + abs(s1+s2-s0)/15.0;

    else

        if deep >= im
            error('adaptive_simpson: maximum subdivision depth reached.');
        end

        % Right subinterval.
        i = i + 1;
        x(i)     = x0 + step;
        fa(i)    = f2;
        fm(i)    = f3;
        fb(i)    = f4;
        h(i)     = step/2.0;
        tol(i)   = err/2.0;
        s(i)     = s2;
        level(i) = deep + 1;

        % Left subinterval.
        i = i + 1;
        x(i)     = x0;
        fa(i)    = f0;
        fm(i)    = f1;
        fb(i)    = f2;
        h(i)     = h(i-1);
        tol(i)   = tol(i-1);
        s(i)     = s1;
        level(i) = level(i-1);

    end
end

end