function ch07_improper_infinite
%==========================================================================
% ch07_improper_infinite.m
%
% Improper integral with an infinite upper limit.
%
% The program evaluates
%
%              infinity
%                 /
%                |  exp(-x) sin(x) dx = 1/2
%                /
%                0
%
% by transforming [0,infinity) to [0,1):
%
%       x = t/(1-t),        dx = dt/(1-t)^2 .
%
% Gauss-Legendre quadrature is then applied on [0,1]. Since its nodes do
% not include the endpoints, t = 1 is never sampled.
%
% Based on the Fortran implementation by Alexander Godunov.
% 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;
exact = 0.5;
nvalues = [8, 16, 32, 64];

% === Plot the integrands ================================================
% First show the original integrand on a finite window.  The exponential
% damping makes x = 20 sufficient to display the important structure.
xplot = linspace(0.0,20.0,2001);
yplot = f(xplot);

figure('Name','Original improper-integral integrand')
plot(xplot,yplot,'LineWidth',1.5)
xlabel('x')
ylabel('f(x)')
title('Original integrand: exp(-x) sin(x)')
grid on

% Now show the transformed integrand on 0 <= t < 1.  The endpoint t = 1
% corresponds to x = infinity and is intentionally excluded.
tplot = linspace(0.0,0.98,2001);
gplot = g(tplot);

figure('Name','Transformed improper-integral integrand')
plot(tplot,gplot,'LineWidth',1.5)
xlabel('t')
ylabel('g(t)')
title('Transformed integrand after x = t/(1-t)')
grid on
% === End plots ==========================================================

fprintf('Improper integral: infinite interval\n');
fprintf('Integral of exp(-x) sin(x) from 0 to infinity\n');
fprintf('Exact value          = %15.7e\n\n', exact);

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

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

end


function y = f(x)
%--------------------------------------------------------------------------
% Original integrand on the semi-infinite interval.
%
% To use another problem, replace only the line below.
%--------------------------------------------------------------------------
y = exp(-x).*sin(x);
end


function y = g(t)
%--------------------------------------------------------------------------
% Transformed integrand after x = t/(1-t).
%--------------------------------------------------------------------------
one_minus_t = 1.0 - t;
x = t./one_minus_t;
y = f(x)./(one_minus_t.^2);
end


function integral = gauss_legendre(fun, a, b, n)
%--------------------------------------------------------------------------
% n-point Gauss-Legendre quadrature on [a,b].
%--------------------------------------------------------------------------

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)
if 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