% Ett enkelt exempel som visar hur man kan beräkna
% en splinefunktion som uppfyller not-a-knot-villkoret.
%
% En riktig kod ser inte ut så här.
%
% Thomas Ericsson, 2007-11-12

% I varje delintervall [t_k, t_k+1], k = 1, ..., n - 1,
% har vi % p_k = a_k t^3 + b_k t^2 + c_k t + d_k.
% Koefficienterna läggs i vektorn x:
% x = [a_1 b_1 c_1 d_1   a_2 b_2 c_2 d_2 ...]^T

% lite data...
n = 7;  % n måste vara >= 5
t = linspace(-3, 3, n);
y = exp(-t.^2);

m = 4 * (n - 1);  % antal koefficienter
b = zeros(m, 1);  % högerled
A = zeros(m);     % matris för systemet

% p_k(t_k) = y_k, p_k(t_k+1) = y_k+1, interpolationsvillkor
cols = [1 2 3 4];
pow  = [3 2 1 0];
row  = 0;         % nummer på ekvation

for k = 1:n-1     % för alla polynom
  row = row + 1;
  A(row, cols) = t(k).^pow;
  b(row) = y(k);

  row = row + 1;
  A(row, cols) = t(k+1).^pow;
  b(row) = y(k+1);

  cols = 4 + cols;
end

% p_k'(t_k) = p_k+1'(t_k), kont. förstaderivator
% p_k' = 3 a_k t^2 + 2 b_k t + c_k
cols  = [1 2 3];
pow   = [2 1 0];
coeff = [3 2 1];
for k = 2:n-1       % för alla inre t_k
  row = row + 1;
  A(row, cols) = coeff .* t(k).^pow;
  cols = 4 + cols;
  A(row, cols) = -coeff .* t(k).^pow;  % obs MINUS
end

% p_k''(t_k) = p_k+1''(t_k), kont. andraderivator
% p_k'' = 6 a_k t + 2 b_k
cols  = [1 2];
pow   = [1 0];
coeff = [6 2];
for k = 2:n-1       % för alla inre t_k
  row = row + 1;
  A(row, cols) = coeff .* t(k).^pow;
  cols = 4 + cols;
  A(row, cols) = -coeff .* t(k).^pow;  % obs MINUS
end

% p_1'''(t_1) = p_2'''(t_1), not-a-knot
row = row + 1;
A(row, 1) =  1;
A(row, 5) = -1;

% p_n-1'''(t_n-1) = p_n'''(t_n-1), not-a-knot
row = row + 1;
A(row, end-7) =  1;
A(row, end-3) = -1;

x = A \ b;  % beräkna koefficienterna

% lite plottar
figure(1)
hold off
plot(t, y, 'ko')
hold on
t_fine = linspace(t(1), t(end));
plot(t_fine, exp(-t_fine.^2), 'k--')  % facit

% evaluera splinefunktionen
for k = 1:n-1
  t_fine = linspace(t(k), t(k+1), 20);
  plot(t_fine, polyval(x(4*(k-1)+1:4*k), t_fine), 'k-')
end

% kontrollera not-a-knot
[x(1:4)'; x(5:8)']
[x(end-7:end-4)'; x(end-3:end)']
norm(x(1:4) - x(5:8))
norm(x(end-7:end-4) - x(end-3:end))

% ovanstående rader blir
% ans =
%     0.0799    0.6210    1.6042    1.3805
%     0.0799    0.6210    1.6042    1.3805
% ans =
%    -0.0799    0.6210   -1.6042    1.3805
%    -0.0799    0.6210   -1.6042    1.3805
% ans = 1.0497e-15
% ans = 1.8201e-15
