Spline cubic with tridiagonal matrix

Question

I wrote this code for my homework on MATLAB about cubic spline interpolation with a tridiagonal matrix. I follow all the steps of the algorithm and I really don't find my error. The code it's ok with second grade functions but when I put, for example, sin(x) the result is not a spline and don't know why because with the other function I have no problems. Can anyone help me to find the error? Thanks

CUBIC SPLINE SCRIPT:

close all
clear all
clf reset

ff = @(x) sin(x);
x = [-2 0 2 4 6 7 8 9 10 11 12];
for i = 1: length(x),
    omega(i) = ff(x(i));
end

n = length(x);
h = zeros(n - 1, 1);
for i = 1: n - 1,
    h(i) = x(i + 1) - x(i);
end

a = zeros(n, 1);
b = zeros(n, 1);
d = zeros(n, 1);
f = zeros(n, 1);

for j = 2: n - 1,
    a(j) = 2*(h(j) + h(j - 1));
    b(j) = h(j - 1);
    d(j) = h(j);
    f(j) = 6 * (omega(j + 1) - omega(j) / h(j)) - (6 * (omega(j) - omega(j - 1)) / h(j - 1));
end

% Starting conditions
a(1) = -2;
f(1) = 0;
a(n) = 4;
f(n) = 0;

% Coefficents
c = tridiag(a, b, d, f);

t = linspace (x(1), x(n), 301);
for k = 1: length(t),
    tk = t(k);
    y(k) = spline_aux(x, omega, c, tk);
    z(k) = ff(tk);
end

plot(t, z, 'b-', 'linewidth', 2)
hold on;
plot(t, y, 'r.', x, omega, 'go')
grid on

TRIADIAGONAL MATRIX:

function [x] = tridiag(a, b, d, f)

n = length(f);
alfa = zeros(n, 1);
beta = zeros(n, 1);

alfa(1) = a(1);
for i = 2: n,
    beta(i) = b(i) / alfa(i - 1);
    fprintf ('  i: %d     beta: %12.8f\n', i, beta(i))
    alfa(i) = a(i) - (beta(i)*d(i - 1));
    fprintf ('  i: %d     alfa: %12.8f\n', i, alfa(i))
end

y(1) = f(1);
for i = 2: n,
    y(i) = f(i) - beta(i)*y(i - 1);
end

x(n) = y(n) / alfa(n);
for i = n - 1: 1,
    x(i) = (y(i) - (d(i)*x(i + 1))) / alfa(i);
end

SPLINE EVALUATION IN tk POINT:

function [s] = spline_aux(x, w, c, tk)

n = length(x);

h = zeros(n - 1, 1);
for i = 1: n - 1,
    h(i) = x(i+1) - x(i);
end

 for i = 1: n - 1,
     if (x(i) <= tk && tk <= x(i+1))
        break
     end
 end

s1 = c(i)*((x(i+1) - tk)^3)/(6*h(i));
s2 = c(i+1)*((tk - x(i))^3)/(6*h(i));
s3 = (w(i)/h(i) - (c(i)*h(i)/6))*(x(i+1) - tk);
s4 = (w(i+1)/h(i) - (c(i+1)*h(i)/6))*(tk - x(i));
s = s1 + s2 + s3 + s4;

Answer 1

Thats because you are not using Matlab's for correctly

in the function function [x] = tridiag(a, b, d, f)

the last for reads for i = n - 1: 1 but that will never execute, you shoudl writte:

for i = n - 1:-1:1

Then works. You should notice that it didt work in ANY previous attempt, not only with sin(x)