function ch07_adaptive_gauss_kronrod
%==========================================================================
% ch07_adaptive_gauss_kronrod.m
%
% Adaptive numerical integration using embedded Gauss-Kronrod rules.
%
% Two quadrature pairs are included:
%   1. Gauss 7-point / Kronrod 15-point (G7-K15)
%   2. Gauss 10-point / Kronrod 21-point (G10-K21)
%
% The same non-recursive adaptive subdivision algorithm is used with
% either pair.
%
% Example:
%   Integral of 1/[1+400(x-0.2)^2] from 0 to 1.
%
% Exact value:
%   [atan(16) + atan(4)] / 20
%
% Based on integration codes 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 shows where adaptive subdivision is most useful.
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 for adaptive Gauss-Kronrod integration')
grid on
% === End plot ===========================================================

fprintf('Adaptive Gauss-Kronrod integration\n');
fprintf('Integral of 1/[1+400(x-0.2)^2] from 0 to 1\n');
fprintf('Exact value          = %18.10e\n\n', exact);

[result, errest, nfun] = GKinteg(@f, a, b, eps, 1);
fprintf('G7-K15 result        = %18.10e\n', result);
fprintf('Estimated error      = %12.4e\n', errest);
fprintf('Function calls       = %10d\n\n', nfun);

[result, errest, nfun] = GKinteg(@f, a, b, eps, 2);
fprintf('G10-K21 result       = %18.10e\n', result);
fprintf('Estimated error      = %12.4e\n', errest);
fprintf('Function calls       = %10d\n', nfun);

end


function y = f(x)
%--------------------------------------------------------------------------
% Narrow-peak test function.
%--------------------------------------------------------------------------
y = 1.0/(1.0 + 400.0*(x-0.2)^2);
end


function [g7, k15] = G7K15(fun, a, b)
%--------------------------------------------------------------------------
% Gauss 7-point / Kronrod 15-point pair on one interval [a,b].
%--------------------------------------------------------------------------

x = [ ...
    0.000000000000000, ...
    0.207784955007898, ...
    0.405845151377397, ...
    0.586087235467691, ...
    0.741531185599394, ...
    0.864864423359769, ...
    0.949107912342758, ...
    0.991455371120813];

wg = [ ...
    0.417959183673469, ...
    0.381830050505119, ...
    0.279705391489277, ...
    0.129484966168870];

wk = [ ...
    0.209482141084728, ...
    0.204432940075298, ...
    0.190350578064785, ...
    0.169004726639267, ...
    0.140653259715525, ...
    0.104790010322250, ...
    0.063092092629979, ...
    0.022935322010529];

h = (b-a)/2.0;
c = (b+a)/2.0;

f0 = fun(c);
fl = zeros(1,7);
fr = zeros(1,7);

for i = 1:7
    fl(i) = fun(c - h*x(i+1));
    fr(i) = fun(c + h*x(i+1));
end

g7 = wg(1)*f0;
for i = 1:3
    g7 = g7 + wg(i+1)*(fl(2*i) + fr(2*i));
end
g7 = h*g7;

k15 = wk(1)*f0;
for i = 1:7
    k15 = k15 + wk(i+1)*(fl(i) + fr(i));
end
k15 = h*k15;

end


function [g10, k21] = G10K21(fun, a, b)
%--------------------------------------------------------------------------
% Gauss 10-point / Kronrod 21-point pair on one interval [a,b].
%--------------------------------------------------------------------------

x = [ ...
    0.0000000000000000000000000, ...
    0.1488743389816312108848260, ...
    0.2943928627014601981311266, ...
    0.4333953941292471907992659, ...
    0.5627571346686046833390001, ...
    0.6794095682990244062343274, ...
    0.7808177265864168970637176, ...
    0.8650633666889845107320967, ...
    0.9301574913557082260012072, ...
    0.9739065285171717200779640, ...
    0.9956571630258080807355273];

wk = [ ...
    0.1494455540029169056649365, ...
    0.1477391049013384913748415, ...
    0.1427759385770600807970943, ...
    0.1347092173114733259280540, ...
    0.1234919762620658510779581, ...
    0.1093871588022976418992106, ...
    0.0931254545836976055350655, ...
    0.0750396748109199527670431, ...
    0.0547558965743519960313813, ...
    0.0325581623079647274788190, ...
    0.0116946388673718742780644];

wg = [ ...
    0.2955242247147528701738930, ...
    0.2692667193099963550912269, ...
    0.2190863625159820439955349, ...
    0.1494513491505805931457763, ...
    0.0666713443086881375935688];

h = (b-a)/2.0;
c = (b+a)/2.0;

f0 = fun(c);

fl = zeros(1,11);
fr = zeros(1,11);

for i = 2:11
    fl(i) = fun(c - h*x(i));
    fr(i) = fun(c + h*x(i));
end

g10 = 0.0;
for i = 1:5
    g10 = g10 + wg(i)*(fl(2*i) + fr(2*i));
end
g10 = h*g10;

k21 = wk(1)*f0;
for i = 2:11
    k21 = k21 + wk(i)*(fl(i) + fr(i));
end
k21 = h*k21;

end


function [result, errest, nfun] = GKinteg(fun, a, b, eps, key)
%--------------------------------------------------------------------------
% Adaptive non-recursive Gauss-Kronrod integration.
%
% key = 1  -> G7-K15
% key = 2  -> G10-K21
%
% The current interval is accepted if
%
%       (200*|K-G|)^1.5 <= local tolerance.
%
% Otherwise it is divided into two equal subintervals and both are
% placed on an explicit stack.
%--------------------------------------------------------------------------

im = 32;

tol   = zeros(1,im);
x     = zeros(1,im);
h     = zeros(1,im);
level = zeros(1,im);

result = 0.0;
errest = 0.0;
nfun = 0;

i = 1;
x(i) = a;
h(i) = b-a;
tol(i) = eps;
level(i) = 1;

while i > 0

    if key == 1
        [gauss, kronrod] = G7K15(fun, x(i), x(i)+h(i));
        nfun = nfun + 15;
    elseif key == 2
        [gauss, kronrod] = G10K21(fun, x(i), x(i)+h(i));
        nfun = nfun + 21;
    else
        error('GKinteg: key must be 1 or 2.');
    end

    x0 = x(i);
    step = h(i);
    err = tol(i);
    deep = level(i);

    i = i - 1;

    difference = abs(kronrod-gauss);

    if (200.0*difference)^1.5 <= err

        result = result + kronrod;
        errest = errest + difference;

    else

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

        % Right subinterval.
        i = i + 1;
        h(i) = step/2.0;
        tol(i) = err/2.0;
        x(i) = x0 + h(i);
        level(i) = deep + 1;

        % Left subinterval.
        i = i + 1;
        x(i) = x0;
        h(i) = h(i-1);
        tol(i) = tol(i-1);
        level(i) = level(i-1);

    end
end

end