function ch05_cubic_spline
%==========================================================================
% Cubic spline interpolation
%
% Function values f(x) are calculated at n base points. The spline
% coefficients are then computed, and the spline is evaluated at
% nint points across the interval. The program also calculates the
% average absolute interpolation error.
%
% Origin:
%   Based on the SPLINE and SEVAL routines presented by
%   G. E. Forsythe, M. A. Malcolm, and C. B. Moler,
%   Computer Methods for Mathematical Computations,
%   Prentice-Hall, 1977.
%
% This MATLAB version is based on the Fortran implementation by
% Alexander Godunov.
% 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: calculate spline coefficients
[b,c,d] = spline_coefficients(xi,yi);

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

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

fprintf('        x        spline      error\n');

for i = 1:nint
    x = xmin + step*(i-1);
    y = f(x);
    ys = spline_eval(x,xi,yi,b,c,d);
    error = ys-y;

    fprintf('%12.5f%12.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('           Average error%12.5f\n',errav);

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

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

xlabel('x');
ylabel('f(x)');
title('Cubic spline interpolation');
legend('Exact function','Cubic spline','Interpolation points', ...
       'Location','best');
grid on;

end


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

y = sin(x);

end


function [b,c,d] = spline_coefficients(x,y)
%==========================================================================
% Calculate the coefficients b(i), c(i), and d(i), i=1,2,...,n,
% for cubic spline interpolation
%
%   s(x) = y(i) + b(i)*(x-x(i)) + c(i)*(x-x(i))^2
%                + d(i)*(x-x(i))^3
%
% for x(i) <= x <= x(i+1).
%
% Input:
%   x = array of data abscissas in strictly increasing order
%   y = array of data ordinates
%
% Output:
%   b, c, d = arrays of spline coefficients
%
% Origin:
%   Based on the SPLINE routine presented by
%   G. E. Forsythe, M. A. Malcolm, and C. B. Moler,
%   Computer Methods for Mathematical Computations,
%   Prentice-Hall, 1977.
%
% This MATLAB version is based on the Fortran implementation by
% Alexander Godunov.
%==========================================================================

n = length(x);
b = zeros(1,n);
c = zeros(1,n);
d = zeros(1,n);

gap = n-1;

% Check input
if n < 2
    return
end

if n < 3
    b(1) = (y(2)-y(1))/(x(2)-x(1));   % linear interpolation
    c(1) = 0.0;
    d(1) = 0.0;
    b(2) = b(1);
    c(2) = 0.0;
    d(2) = 0.0;
    return
end

%
% Step 1: preparation
%
d(1) = x(2)-x(1);
c(2) = (y(2)-y(1))/d(1);

for i = 2:gap
    d(i) = x(i+1)-x(i);
    b(i) = 2.0*(d(i-1)+d(i));
    c(i+1) = (y(i+1)-y(i))/d(i);
    c(i) = c(i+1)-c(i);
end

%
% Step 2: end conditions
%
b(1) = -d(1);
b(n) = -d(n-1);
c(1) = 0.0;
c(n) = 0.0;

if n ~= 3
    c(1) = c(3)/(x(4)-x(2)) - c(2)/(x(3)-x(1));
    c(n) = c(n-1)/(x(n)-x(n-2)) - c(n-2)/(x(n-1)-x(n-3));

    c(1) = c(1)*d(1)^2/(x(4)-x(1));
    c(n) = -c(n)*d(n-1)^2/(x(n)-x(n-3));
end

%
% Step 3: forward elimination
%
for i = 2:n
    h = d(i-1)/b(i-1);
    b(i) = b(i)-h*d(i-1);
    c(i) = c(i)-h*c(i-1);
end

%
% Step 4: back substitution
%
c(n) = c(n)/b(n);

for j = 1:gap
    i = n-j;
    c(i) = (c(i)-d(i)*c(i+1))/b(i);
end

%
% Step 5: compute spline coefficients
%
b(n) = (y(n)-y(gap))/d(gap) + d(gap)*(c(gap)+2.0*c(n));

for i = 1:gap
    b(i) = (y(i+1)-y(i))/d(i) - d(i)*(c(i+1)+2.0*c(i));
    d(i) = (c(i+1)-c(i))/d(i);
    c(i) = 3.0*c(i);
end

c(n) = 3.0*c(n);
d(n) = d(n-1);

end


function s = spline_eval(u,x,y,b,c,d)
%==========================================================================
% Evaluate the cubic spline at point u:
%
%   spline_eval = y(i) + b(i)*(u-x(i)) + c(i)*(u-x(i))^2
%                       + d(i)*(u-x(i))^3
%
% where x(i) <= u <= x(i+1).
%
% Input:
%   u       = abscissa at which the spline is evaluated
%   x, y    = arrays of given data points
%   b, c, d = spline coefficients computed by spline_coefficients
%
% Output:
%   s = interpolated value at u
%
% Origin:
%   Based on the SEVAL routine presented by
%   G. E. Forsythe, M. A. Malcolm, and C. B. Moler,
%   Computer Methods for Mathematical Computations,
%   Prentice-Hall, 1977.
%
% This MATLAB version is based on the Fortran implementation by
% Alexander Godunov.
%==========================================================================

n = length(x);

% If u is outside the x interval, return the nearest boundary value.
if u <= x(1)
    s = y(1);
    return
end

if u >= x(n)
    s = y(n);
    return
end

%
% Binary search for i such that x(i) <= u <= x(i+1)
%
i = 1;
j = n+1;

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

    if u < x(k)
        j = k;
    else
        i = k;
    end
end

%
% Evaluate spline interpolation
%
dx = u-x(i);
s = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i)));

end