Deviated results of natural cubic spline function

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 1x10 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 1x20 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 (20x1 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.

[Stephen23]

@Mindanao: you might like to read this discussion about natural splines:

https://www.mathworks.com/matlabcentral/answers/325706-natural-spline-unmkpp-gives-wrong-answer

[Mindanao]

Thanks for the link but it doesn't really help me out with this.

[dpb]

If have the CurveFitting Toolbox, there's a natural spline function cscvn available prepackaged rather than "roll your own".

A little much to ask volunteers here to debug your implementation from scratch, though...

[John D'Errico]

I tried answering the last question on this same problem by Mindanao, but I never got a definition of what the variables were.

One can ask for help on fixing your code. But this code essentially refers to a derivation that we don't have, with variables that have no meaning except for that derivation.

Thus, s is apparently a vector of second derivatives of the spline at the breaks. That I can handle, and I could in theory verify their values. But then we see vectors c1 and c2, defined only as integration constants. Since I don't have the derivation, they are meaningless. So there is little I can really check.

I might point out that the code MAY be fine, but that it is simply being used improperly. Thus see how it was called:

x_a = linspace(-2,2,10);
f = sin(x);

So mindanao defines the variable x_a, but then uses x to compute f.

It is possible that a variable x already existed in the workspace, where x was completely different. That MIGHT explain the garbage results.

Sorry, but while I'd be a good person to answer these questions, I can't try to guess how the code was actually called, or what was intended by c1 and c2. You need to be careful when writing code. Write slowly, carefully. Make sure you know what information was passed in. And if you seriously want help, then you need to explain what the variables mean.

Answer 1

[John D'Errico]

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.