% A short script file to generate a piecewise cubic spline approximation
% of the positive portion of Runge's function f(x) = 1 / (1 + 25 x^2).
%
% It would be much quicker to code (and probably more efficient to run) if
% we just used Matlab's spline command.  But this script file demonstrates
% the construction of the matrix system which is solved for the spline
% coefficients.
%
% We pretend that we know only that the underlying function is symmetric and
% differentiable.  Therefore, to find the missing two equations (beyond the
% basic interpolation and continuity equations) we use a zero derivative at
% the left side of the interval (x = 0) and not-a-knot at the right side.
%
% Also note that the indexing in Matlab starts at 1, so the indices will
% be different that those in the notes.
%
% This is a script file, so there is no "function" line.  Also, it has no
% internal scope or local variables, but rather executes in the scope of the
% caller.
%
% Ian Mitchell for CPSC 303, January 26, 2006.

% Data.
xs = [ 0; 0.1; 0.25; 0.5; 1.0 ];
ys = 1 ./ (1 + 25 * xs.^2);
hs = diff(xs);

% System.
A = [ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
      0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
      0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
      0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
      1, 0, 0, 0, hs(1), 0, 0, 0, hs(1)^2, 0, 0, 0, hs(1)^3, 0, 0, 0;
      0, 1, 0, 0, 0, hs(2), 0, 0, 0, hs(2)^2, 0, 0, 0, hs(2)^3, 0, 0;
      0, 0, 1, 0, 0, 0, hs(3), 0, 0, 0, hs(3)^2, 0, 0, 0, hs(3)^3, 0;
      0, 0, 0, 1, 0, 0, 0, hs(4), 0, 0, 0, hs(4)^2, 0, 0, 0, hs(4)^3;
      0, 0, 0, 0, 1, -1, 0, 0, 2 * hs(1), 0, 0, 0, 3 * hs(1)^2, 0, 0,0;
      0, 0, 0, 0, 0, 1, -1, 0, 0, 2 * hs(2), 0, 0, 0, 3 * hs(2)^2, 0,0;
      0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 2 * hs(3), 0, 0, 0, 3 * hs(3)^2,0;
      0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 6 * hs(1), 0, 0, 0;
      0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 6 * hs(2), 0, 0;
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 6 * hs(3), 0;
      0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, -6 ];
u = [ ys(1); ys(2); ys(3); ys(4); ys(2); ys(3); ys(4); ys(5);
      0; 0; 0; 0; 0; 0; 0; 0 ];

% Compute coefficients.
w = A \ u;

% Our basis function in inline form.
poly = inline('a + b * (x - xi) + c * (x - xi).^2 + d * (x - xi).^3', ...
              'x', 'xi', 'a', 'b', 'c', 'd');

% Evaluate it.  Check out help find for more information on this very
% useful function.
many_xs = linspace(min(xs), max(xs), 201)';
many_ys = zeros(size(many_xs));

% Not particularly efficient to loop through all the data points.
interval = 1;
for i = 1 : length(many_xs)
  x_value = many_xs(i);
  if(x_value > xs(interval + 1))
    interval = interval + 1;
  end
  if(interval == length(xs))
    warning('Cannot extrapolate');
    break;
  end
  
  many_ys(i) = feval(poly, x_value, xs(interval), w(interval), ...
                     w(interval + 4), w(interval + 8), w(interval + 12));
end

% Actual function.
many_fs = 1 ./ (1 + 25 * many_xs.^2);

h_interp = plot(many_xs, many_ys, 'b-');
hold on;
h_func = plot(many_xs, many_fs, 'r--');
h_data = plot(xs, ys, 'b*');
set([ h_interp(1), h_func(1) ], 'LineWidth', 2);
xlabel('x');  ylabel('y');
title('Spline interpolant for Runges function -- Ian Mitchell for CPSC 303')
legend([ h_interp(1), h_func(1) ], { 'Cubic Spline', 'Runge function' });