function ch07_quanc8
%=======================================================================
% ch07_quanc8.m
%
% Adaptive numerical integration using QUANC8.
%
% The program evaluates
%
%               pi
%               /
%              |  sin(x) dx = 2
%              /
%               0
%
% with the adaptive 8-panel Newton-Cotes QUANC8 algorithm.
%
% QUANC8 repeatedly subdivides intervals and compares successive
% Newton-Cotes estimates to control the integration error.
%
% Based on the QUANC8 routine presented by
% G. E. Forsythe, M. A. Malcolm, and C. B. Moler,
% Computer Methods for Mathematical Computations,
% Prentice-Hall, 1977.
%
% MATLAB version adapted by Alexander Godunov, February 2022.
% Revised and reorganized for the companion website, 2026.
%=======================================================================

clear; clc;

% ---- Problem setup ----------------------------------------------------
a      = 0.0;
b      = pi;
abserr = 0.0;
relerr = 1.0e-8;

exact = 2.0;

% === Plot the integrand =================================================
% Plotting the function before integration is a useful diagnostic step:
% it shows the numerical algorithm what the user already sees visually.
Nplot = 1000;
xplot = linspace(a,b,Nplot+1);
yplot = arrayfun(@f,xplot);

figure('Name','Integrand for QUANC8')
plot(xplot,yplot,'LineWidth',1.5)
xlabel('x')
ylabel('f(x)')
title('Integrand for QUANC8 adaptive integration')
grid on
% === End plot ===========================================================

[result, errest, flag, nfun] = quanc8(@f, a, b, abserr, relerr);

fprintf('QUANC8 adaptive integration\n');
fprintf('Integral of sin(x) from 0 to pi\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       = %12d\n', nfun);
fprintf('Flag                 = %12.5f\n', flag);

end


function y = f(x)
%-----------------------------------------------------------------------
% Function to be integrated.
%
% To use this program for another problem, replace only the line
% defining y below and change the integration limits in the main
% program if needed.
%
% Input:
%   x - integration variable
%
% Output:
%   y - value of the integrand
%-----------------------------------------------------------------------

y = sin(x);

end


function [result, errest, flag, nfun] = quanc8(fun, a, b, abserr, relerr)
%-----------------------------------------------------------------------
% Adaptive integration using the 8-panel Newton-Cotes QUANC8 algorithm.
%
% Input:
%   fun    - function handle for the integrand
%   a      - lower integration limit
%   b      - upper integration limit; b may be less than a
%   abserr - requested absolute error tolerance, abserr >= 0
%   relerr - requested relative error tolerance, relerr >= 0
%
% Output:
%   result - numerical approximation to the integral
%   errest - estimated magnitude of the numerical error
%   flag   - reliability indicator
%   nfun   - number of function evaluations
%
% Reliability flag:
%   flag = 0 indicates that the requested tolerance was probably met.
%   A nonzero value indicates that one or more intervals did not
%   converge normally or that the function-evaluation limit was
%   approached.
%
% Method:
%   QUANC8 is an automatic adaptive quadrature routine based on an
%   8-panel Newton-Cotes formula.  The interval is subdivided where
%   the difference between successive estimates is too large.
%
%   The local error estimate is based on
%
%       |Q_new - Q_old| / 1023.
%
%   The routine reuses previously computed function values whenever
%   an interval is subdivided.
%
% Based on the QUANC8 routine presented by
% G. E. Forsythe, M. A. Malcolm, and C. B. Moler,
% Computer Methods for Mathematical Computations,
% Prentice-Hall, 1977.
%
% MATLAB version adapted by Alexander Godunov, February 2022.
%-----------------------------------------------------------------------

% Initialize outputs so the zero-length interval case is well defined.
result = 0.0;
errest = 0.0;
flag   = 0.0;
nfun   = 0;

if abserr < 0.0 || relerr < 0.0
    error('quanc8: error tolerances must be non-negative.');
end

if a == b
    return
end

% ---- Storage ----------------------------------------------------------
qright = zeros(32,1);
fval   = zeros(17,1);
x      = zeros(17,1);
fsave  = zeros(31,31);
xsave  = zeros(31,31);

% ---- General initialization ------------------------------------------
levmin = 1;
levmax = 30;
levout = 6;
nomax  = 5000;

% The algorithm changes strategy when the number of function
% evaluations approaches this limit.
nofin = nomax - 8*(levmax - levout + 128);

% Newton-Cotes coefficients.
w0 =   3956.0 / 14175.0;
w1 =  23552.0 / 14175.0;
w2 =  -3712.0 / 14175.0;
w3 =  41984.0 / 14175.0;
w4 = -18160.0 / 14175.0;

cor11 = 0.0;
area  = 0.0;

% ---- Initialize first interval ---------------------------------------
lev   = 0;
nim   = 1;
x0    = a;
x(16) = b;
qprev = 0.0;
f0    = fun(x0);

stone = (b-a)/16.0;

x(8)  = (x0    + x(16))/2.0;
x(4)  = (x0    + x(8))/2.0;
x(12) = (x(8)  + x(16))/2.0;
x(2)  = (x0    + x(4))/2.0;
x(6)  = (x(4)  + x(8))/2.0;
x(10) = (x(8)  + x(12))/2.0;
x(14) = (x(12) + x(16))/2.0;

for j = 2:2:16
    fval(j) = fun(x(j));
end

nfun = 9;

% ---- Main adaptive loop ----------------------------------------------
while nfun <= nomax

    % Complete the 17-point grid on the current interval.
    x(1) = (x0 + x(2))/2.0;
    fval(1) = fun(x(1));

    for j = 3:2:15
        x(j) = (x(j-1) + x(j+1))/2.0;
        fval(j) = fun(x(j));
    end

    nfun = nfun + 8;

    step = (x(16)-x0)/16.0;

    % Newton-Cotes estimates on the left and right halves.
    qleft = ...
        ( w0*(f0      + fval(8)) ...
        + w1*(fval(1) + fval(7)) ...
        + w2*(fval(2) + fval(6)) ...
        + w3*(fval(3) + fval(5)) ...
        + w4*fval(4) )*step;

    qright(lev+1) = ...
        ( w0*(fval(8)  + fval(16)) ...
        + w1*(fval(9)  + fval(15)) ...
        + w2*(fval(10) + fval(14)) ...
        + w3*(fval(11) + fval(13)) ...
        + w4*fval(12) )*step;

    qnow  = qleft + qright(lev+1);
    qdiff = qnow - qprev;
    area  = area + qdiff;

    % ---- Local convergence test --------------------------------------
    esterr = abs(qdiff)/1023.0;
    tolerr = max(abserr, relerr*abs(area));
    tolerr = tolerr*(step/stone);

    if lev < levmin
        key = 1;
    elseif lev >= levmax
        key = 2;
    elseif nfun > nofin
        key = 3;
    elseif esterr <= tolerr
        key = 4;
    else
        key = 1;
    end

    switch key

        case 1
            % No convergence: subdivide the current interval.
            nim = 2*nim;
            lev = lev + 1;

            % Save the right half for later.
            for i = 1:8
                fsave(i,lev) = fval(i+8);
                xsave(i,lev) = x(i+8);
            end

            % Continue immediately with the left half.
            qprev = qleft;

            for i = 1:8
                j = -i;
                fval(2*j+18) = fval(j+9);
                x(2*j+18)    = x(j+9);
            end

            continue

        case 2
            % Maximum subdivision level reached.
            flag = flag + 1.0;

        case 3
            % Function-evaluation limit is being approached.
            nofin  = 2*nofin;
            levmax = levout;
            flag   = flag + (b-x0)/(b-a);

        case 4
            % Current interval satisfies the convergence test.

    end

    % ---- Accept current interval -------------------------------------
    result = result + qnow;
    errest = errest + esterr;
    cor11  = cor11 + qdiff/1023.0;

    % Locate the next interval that still has to be processed.
    while rem(nim,2) ~= 0
        nim = floor(nim/2);
        lev = lev - 1;
    end

    nim = nim + 1;

    if lev <= 0
        break
    end

    % Restore saved data for the next interval.
    qprev = qright(lev);
    x0    = x(16);
    f0    = fval(16);

    for i = 1:8
        fval(2*i) = fsave(i,lev);
        x(2*i)    = xsave(i,lev);
    end

end

% ---- Final correction and error estimate -----------------------------
result = result + cor11;

if errest == 0.0
    return
end

% Ensure that the reported error is representable relative to result.
temp = abs(result) + errest;

while temp == abs(result)
    errest = 2.0*errest;
    temp = abs(result) + errest;
end

end