function ch05_runge_chebyshev
%==========================================================================
% Runge phenomenon and Chebyshev nodes
%
% The Runge function f(x) is sampled using two sets of n base points:
% equally spaced points and Chebyshev points. A global polynomial
% interpolant is then evaluated for both cases at nint points across
% the interval. The program also calculates the average absolute
% interpolation error for each choice of nodes.
%
% The interpolating polynomial is evaluated in Lagrange form.
%
% This example illustrates that increasing the polynomial order with
% equally spaced points can produce large oscillations near the ends
% of the interval, while Chebyshev nodes greatly reduce this effect.
%
% Alexander Godunov
% Prepared for the companion website, 2026.
%==========================================================================

n = 13;       % base points for interpolation
nint = 25;    % compute interpolation in nint points

xmin = -1.0;
xmax =  1.0;

xe = zeros(1,n);
ye = zeros(1,n);
xc = zeros(1,n);
yc = zeros(1,n);

% Step 1: generate equally spaced interpolation points
step = (xmax-xmin)/(n-1);

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

% Step 2: generate Chebyshev interpolation points
% The nodes are written in increasing order from left to right.
for i = 1:n
    xc(i) = cos((2*(n-i)+1)*pi/(2*n));
    yc(i) = f(xc(i));
end

% Step 3: evaluate both global interpolants at nint points
errav_e = 0.0;
errav_c = 0.0;
step = (xmax-xmin)/(nint-1);

fprintf('     Runge Function: Equally Spaced and Chebyshev Nodes\n');
fprintf('     number of interpolation points = %2d\n',n);
fprintf('            x   equal poly  equal error    Cheb poly   Cheb error\n');

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

    yse = lagrange_global(x,xe,ye);
    ysc = lagrange_global(x,xc,yc);

    erre = yse-y;
    errc = ysc-y;

    fprintf('%13.6f%13.6f%13.6f%13.6f%13.6f\n', ...
            x,yse,erre,ysc,errc);

    % Step 4: calculate the average absolute interpolation errors
    errav_e = errav_e + abs(erre)/nint;
    errav_c = errav_c + abs(errc)/nint;
end

fprintf('\n Average absolute error:\n');
fprintf('   equally spaced = %12.6f\n',errav_e);
fprintf('   Chebyshev      = %12.6f\n',errav_c);

% Plot the Runge function and both interpolants on a fine grid.
xfine = linspace(xmin,xmax,600);
yfine = arrayfun(@f,xfine);

ye_fine = zeros(size(xfine));
yc_fine = zeros(size(xfine));

for i = 1:length(xfine)
    ye_fine(i) = lagrange_global(xfine(i),xe,ye);
    yc_fine(i) = lagrange_global(xfine(i),xc,yc);
end

figure;
plot(xfine,yfine,'-','LineWidth',1.2);
hold on;
plot(xfine,ye_fine,'--','LineWidth',1.2);
plot(xfine,yc_fine,'-.','LineWidth',1.2);
plot(xe,ye,'o','MarkerSize',5);
plot(xc,yc,'s','MarkerSize',5);
hold off;

xlabel('x');
ylabel('f(x)');
title('Runge function: equally spaced and Chebyshev nodes');
legend('Exact function', ...
       'Equally spaced interpolation', ...
       'Chebyshev interpolation', ...
       'Equally spaced nodes', ...
       'Chebyshev nodes', ...
       'Location','best');
grid on;

end


function y = f(x)
%--------------------------------------------------------------------------
% Function f(x): Runge function
%--------------------------------------------------------------------------

y = 1.0/(1.0 + 25.0*x*x);

end


function y = lagrange_global(xx,xi,yi)
%==========================================================================
% Global polynomial interpolation in Lagrange form.
%
% All data points are used to construct and evaluate the
% interpolating polynomial.
%
% Input:
%   xx = abscissa at which the interpolation is evaluated
%   xi = array of data abscissas
%   yi = array of data ordinates
%
% Output:
%   y = interpolated value at xx
%
% Comments:
%   The routine works for both equally and unequally spaced points.
%   Unlike the local interpolation routine used elsewhere in this
%   chapter, all available points are used here so that the effect
%   of node placement on a high-order global polynomial is visible.
%
% Alexander Godunov
% Prepared for the companion website, 2026.
%==========================================================================

n = length(xi);
y = 0.0;

for i = 1:n
    lambda = 1.0;

    for j = 1:n
        if j ~= i
            lambda = lambda*(xx-xi(j))/(xi(i)-xi(j));
        end
    end

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

end