function ch05_polynomial_interpolation
%==========================================================================
% General polynomial interpolation
%
% Function values f(x) are calculated at n base points. A local
% polynomial of selected order is then used to interpolate the
% function at nint points across the interval. The program also
% calculates the average absolute interpolation error.
%
% The interpolating polynomial is evaluated in Lagrange form.
% For order = 1, 2, 3, ... the routine uses 2, 3, 4, ...
% neighboring data points, respectively.
%
% Original implementation by Alexander Godunov, January 2010.
% MATLAB version revised for the companion website, 2026.
%==========================================================================

n = 11;       % base points for interpolation
nint = 21;    % compute interpolation in nint points

xmin = 0.0;
xmax = 2.0;

xi = zeros(1,n);
yi = zeros(1,n);

% Step 1: generate xi and yi from f(x), xmin, xmax, and n
step = (xmax-xmin)/(n-1);

for i = 1:n
    xi(i) = xmin + step*(i-1);
    yi(i) = f(xi(i));
end

% Step 2: choose the interpolation order
% order = 1 -> linear
% order = 2 -> quadratic
% order = 3 -> cubic
order = 3;

fprintf('    General Polynomial Interpolation\n');
fprintf('     order of interpolation = %2d\n',order);

% Step 3: evaluate the interpolant at nint points
errav = 0.0;
step = (xmax-xmin)/(nint-1);

xplot = zeros(1,nint);
ysplot = zeros(1,nint);

fprintf('           x    interpolated       error\n');

for i = 1:nint
    x = xmin + step*(i-1);
    y = f(x);

    ys = polynomial_interp(x,xi,yi,order+1);

    error = ys-y;

    fprintf('%12.5f%16.5f%12.5f\n',x,ys,error);

    % Step 4: calculate the average absolute interpolation error
    errav = errav + abs(y-ys)/nint;

    xplot(i) = x;
    ysplot(i) = ys;
end

fprintf('%28s%12.5f\n','Average error',errav);

% Plot the exact function, interpolation points, and interpolated curve.
xfine = linspace(xmin,xmax,400);
yfine = arrayfun(@f,xfine);

% Evaluate the interpolant on a fine grid for a smooth plotted curve.
yinterp_fine = zeros(size(xfine));

for i = 1:length(xfine)
    yinterp_fine(i) = polynomial_interp(xfine(i),xi,yi,order+1);
end

figure;
plot(xfine,yfine,'-','LineWidth',1.2);
hold on;
plot(xfine,yinterp_fine,'--','LineWidth',1.2);
plot(xi,yi,'o','MarkerSize',6,'MarkerFaceColor','auto');
hold off;

xlabel('x');
ylabel('f(x)');
title(sprintf('Polynomial interpolation, order = %d',order));
legend('Exact function','Interpolated curve','Interpolation points', ...
       'Location','best');
grid on;

end


function y = f(x)
%--------------------------------------------------------------------------
% Function used to generate the interpolation data.
%--------------------------------------------------------------------------

y = sin(x);

end


function y = polynomial_interp(xx,xi,yi,npts)
%==========================================================================
% Local polynomial interpolation of selected order.
%
% The interpolating polynomial is evaluated in Lagrange form using
% npts neighboring data points. The interpolation order is npts-1.
%
% Input:
%   xx   = abscissa at which the interpolation is evaluated
%   xi   = array of data abscissas
%   yi   = array of data ordinates
%   npts = number of points used for interpolation
%          (interpolation order = npts-1)
%
% Output:
%   y = interpolated value at xx
%
% Comments:
%   The routine works for both equally and unequally spaced xi.
%   A binary search locates the interval containing xx, and a local
%   group of neighboring data points is selected for interpolation.
%   If xx lies outside the data interval, the nearest boundary value
%   is returned.
%
% Original implementation by Alexander Godunov, January 2010.
% MATLAB version revised for the companion website, 2026.
%==========================================================================

ni = length(xi);

% Check the number of interpolation points
nuse = npts;

if nuse > ni
    nuse = ni;
end

% If xx is outside the xi interval, return a boundary value
if xx <= xi(1)
    y = yi(1);
    return
end

if xx >= xi(ni)
    y = yi(ni);
    return
end

% Binary search to find i such that xi(i) < xx < xi(i+1)
i = 1;
j = ni;

while j > i+1
    k = floor((i+j)/2);

    if xx < xi(k)
        j = k;
    else
        i = k;
    end
end

% Shift i so that xx lies near the middle of the selected data points
i = i + 1 - floor(nuse/2);

% Keep the interpolation points inside the available data range
if i < 1
    i = 1;
end

if i + nuse - 1 > ni
    i = ni-nuse+1;
end

% Evaluate the interpolating polynomial in Lagrange form
y = 0.0;

for js = i:i+nuse-1
    lambda = 1.0;

    for jl = i:i+nuse-1
        if jl ~= js
            lambda = lambda*(xx-xi(jl))/(xi(js)-xi(jl));
        end
    end

    y = y + yi(js)*lambda;
end

end