MATLAB: Deviated results of natural cubic spline function

cubicinterpolationspline

I'm setting up a function to output to evaluate the spline function in certain points based on a given set of abscis points with its function values, but it's giving me very deviated results.
The inputs of my function are a vector of abscis points (x_a) with its function values (f) and also a vector of equidistant interpolation points in the interval (t) from which I want the function values of the spline as the outputs. So I test my code with the following vectors: x_a is a 1×10 double vector of 10 equidistant points in the interval [-2,2], f are the sin(x)-function values of x_a and t is a 1×20 double vector of 20 equidistant points in the interval [-2,2]:
x_a = linspace(-2,2,10);
f = sin(x);
t = linspace(-2,2,20);
This is my function (based on the natural interpolating cubic spline method):
function [ y ] = naturalspline( x_a, f, t )
if length(x_a) ~= length(f) 
  error('vectors x_a and f must have same length'); 
end
n = length(x_a);
d = zeros(n,1); %vector of differences between abscis points
d(1)=x_a(1);
for i = 2:n
  d(i) = x_a(i) - x_a(i-1); 
end
% Coefficientmatrix C (tridiagonal): 
C = zeros(n); 
% Natural Spline boundary conditions: 
C(1,1)= 2*(d(1)+d(2));
C(n,n) = 2*(d(n-1)+d(n));
for j = 2:n-1 
  C(j,j-1) = d(j-1); 
  C(j,j) = 2*(d(j-1)+d(j)); 
  C(j-1,j) = d(j); 
end
% Vector b: 
b = zeros(n,1);
for i = 2:n-1 
  b(i) = (6/d(i))*(f(i+1)-f(i)) - (6/d(i-1))*(f(i)-f(i-1));     
end
% Vector s (second derivatives, solution of tridiagonal system): 
s = C\b;
% Integration constant c1 : 
c1 = zeros(n,1); 
for i = 1:n
  c1(i) = (f(i)/d(i)) - d(i)/(6*s(i));
end
% Integration constant c2: 
c2 = zeros(n,1); 
for i = 1:n-1 
  c2(i) = (f(i)/d(i+1)) - d(i+1)/(6*s(i));
end
%Setting up polynomials for each interval
for i = 1:n-1 
  p{i} = @(x) [s(i+1)*((x-x_a(i))^3)/(6*d(i+1)) - s(i)*((x-x_a(i+1))^3)/(6*d(i+1)) ...
      + c1(i+1)*(x-x_a(i)) + c2(i+1)*(x_a(i+1)-x)]; 
end
% evaluating spline function values of the interpolation points
j = length(t);
y = zeros(j,1);
for i = 1:j
[~,~,k] = histcounts(t(i),x_a); %k decides which interval t(i) is located
y(i) = p{k}(t(i));
end
When I execute this function, it returns me y (20×1 double). But I checked the values and compared to the real sin(x)-values and it's not even close to them. The last 3 y-values are -Inf so that makes it even more doubtful. Somebody familiar with the cubic spline method and sees the fault in my code or used formulas? Every help is appreciated.

Best Answer

The only possible thing I can say to help is that your code is wrong, starting at the beginning.
Even if I assume that you really intended to write
f = sin(x_a);
instead of f=sin(x), LOOK AT EVERY intermediate RESULT. Think about what you see.
What is a natural cubic spline? It has second derivatives that are zero at each end. Now, look at s as computed by your code:
s
s =
  -0.024467
   -0.17127
     1.2292
    0.54368
    0.24562
   -0.22541
   -0.64476
   -0.84511
     -1.264
          0
You claim that s is a vector of the second derivatives at the breaks. With 10 data points, I see that s is a vector of length 10. In fact, I even see that s(10) was zero.
What was s(1)? Is it zero? If not, then are these the second derivatives from a Natural cubic spline?
When you write code, start at the beginning. Work through EVERY line. Test EVERY result to see if it makes sense. This is especially true if you are just learning to write code.
I'm afraid I can't go much deeper, lacking the algorithm you wrote this code from. I can give you one thing, the correct values of the second derivatives at the knots for a natural cubic spline, assuming that f was computed as sin(x_a).
ans =
          0
     1.2642
    0.84458
    0.64674
    0.21801
   -0.21801
   -0.64674
   -0.84458
    -1.2642
          0
(That a natural cubic is actually a poor choice for this function on this support is not terribly relevant. That is something you might learn at some point in time. But it is something worth thinking about. Do the assumed boundary conditions fit well here?)
So, were I you, I would check your equations, CAREFULLY. Look for an error in how you defined C or possibly b, or both.
Are there other errors?I do see something that I would suspect highly.
d = zeros(n,1); %vector of differences between abscis points
d(1)=x_a(1);
for i = 2:n
  d(i) = x_a(i) - x_a(i-1); 
end
Do the first element of the difference vector is x_a(1)? The vector of differences would seem to be correctly obtained simply as
d = diff(x_a);
I would expect d to be of length n-1.
So part of your problem may simply be the garbage introduced in the first element of d.
Again, write code ONE line at a time. Verify EVERY SINGLE LINE. Does it do what you expect it to do? If not, then fix it, before you go to the next line. Do NOT write an untested large block of code, hoping it will all work the first time. It won't work, and you will just need to start at the very beginning anyway, testing and verifying every line anyway.
Related Question