function ch07_gauss_legendre
%==========================================================================
% ch07_gauss_legendre.m
%
% Gauss-Legendre quadrature for numerical integration.
%
% The program evaluates
%
%       integral sin(x) dx,   0 <= x <= pi,
%
% using several n-point Gauss-Legendre rules. The exact value is 2.
%
% The Gauss-Legendre nodes and weights are generated numerically.
% Roots of the Legendre polynomial P_n(x) are found by Newton iteration,
% and symmetry is used so that only half of the roots are computed.
%
% Based on the original 8-point and 16-point Gauss quadrature programs
% written by Alexander Godunov, October 2009.
% 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 = pi;
exact = 2.0;
nvalues = [2, 4, 8, 16];

% === 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 Gauss-Legendre quadrature')
grid on
% === End plot ===========================================================

fprintf('Gauss-Legendre quadrature\n');
fprintf('Integral of sin(x) from 0 to pi\n');
fprintf('Exact value          = %15.7e\n\n', exact);

fprintf('        n       Gauss-Legendre\n');
fprintf('--------------------------------\n');

for n = nvalues
    integral = gauss_legendre(@f, a, b, n);
    fprintf('%9d%21.7e\n', n, integral);
end

end


function y = f(x)
%--------------------------------------------------------------------------
% Test function for the integration example.
%--------------------------------------------------------------------------
y = sin(x);
end


function integral = gauss_legendre(fun, a, b, n)
%--------------------------------------------------------------------------
% n-point Gauss-Legendre quadrature on [a,b].
%
% The nodes are the roots of P_n(x), found by Newton iteration.
% The corresponding weights are calculated from P'_n(x).
% Symmetry is used so that only half of the roots are computed.
%--------------------------------------------------------------------------

tol = 1.0e-14;
max_iter = 100;

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

midpoint = (a+b)/2.0;
halfwidth = (b-a)/2.0;

integral = 0.0;
m = floor((n+1)/2);

for i = 1:m

    z = cos(pi*(i-0.25)/(n+0.5));
    converged = false;

    for iter = 1:max_iter
        [pn, pnm1] = legendre_pair(n, z);

        if n == 1
            dp = 1.0;
        else
            dp = n*(z*pn-pnm1)/(z*z-1.0);
        end

        zold = z;
        z = zold - pn/dp;

        if abs(z-zold) <= tol
            converged = true;
            break;
        end
    end

    if ~converged
        error('gauss_legendre: Newton iteration did not converge.');
    end

    [pn, pnm1] = legendre_pair(n, z);

    if n == 1
        dp = 1.0;
    else
        dp = n*(z*pn-pnm1)/(z*z-1.0);
    end

    weight = 2.0/((1.0-z*z)*dp*dp);

    xleft  = midpoint - halfwidth*z;
    xright = midpoint + halfwidth*z;

    if abs(z) <= tol
        integral = integral + weight*fun(midpoint);
    else
        integral = integral + weight*(fun(xleft) + fun(xright));
    end
end

integral = halfwidth*integral;

end


function [pn, pnm1] = legendre_pair(n, z)
%--------------------------------------------------------------------------
% Return P_n(z) and P_(n-1)(z) using the three-term recurrence.
%--------------------------------------------------------------------------

if n == 0
    pn = 1.0;
    pnm1 = 0.0;
    return;
elseif n == 1
    pn = z;
    pnm1 = 1.0;
    return;
end

p0 = 1.0;
p1 = z;

for j = 2:n
    p2 = ((2.0*j-1.0)*z*p1 - (j-1.0)*p0)/j;
    p0 = p1;
    p1 = p2;
end

pn = p1;
pnm1 = p0;

end