Determine the cubic spline from four points without using built-in matlab functions?

My problem...

He has given us this "formula from our notes...

My function is as follows right now...

function [ yy ] = Kable(x,y,xx)
%  xx is a vector of x-values to be interpolated and 
%  yy is an output vector of interpolated values corresponding to xx
%  also plots data points and cubic spline interpolation
%   x = [x1 x2 x3 x4], y = [y1 y2 y3 y4]
A = [x(1)^3 x(1)^2 x(1) x(1)^0 0 0 0 0;...
     x(2)^3 x(2)^2 x(2) x(2)^0 0 0 0 0;...
     0 0 0 0 x(2)^3 x(2)^2 x(2) x(2)^0;...
     x(3)^3 x(3)^2 x(3) x(3)^0 0 0 0 0;...
     0 0 0 0 x(3)^3 x(3)^2 x(3) x(3)^0;...
     0 0 0 0 x(4)^3 x(4)^2 x(4) x(4)^0;...
     3*x(2)^2 2*x(2) 1 0 -3*x(2)^2 -2*x(2) -1 0;...
     3*x(3)^2 2*x(3) 1 0 -3*x(3)^2 -2*x(3) -1 0;...
     3*x(2) 1 0 0 -3*x(2) -1 0 0;...
     3*x(3) 1 0 0 -3*x(3) -1 0 0;...
     3*x(1)^2 2*x(1) 1 0 0 0 0 0;...
     0 0 0 0 3*x(4)^2 2*x(4) 1 0];
 disp(A)
d = [y(1) y(2) y(2) y(3) y(3) y(4) 0 0 0 0 0 0]';
yy = A\d;
end

My problem is I don't understand how to extend the matrix to 4 points given... I am pretty sure that since it is a cubic spline there will still be only 8 columns, but with 4 points there would be 12 rows. I understand the coding part I believe but I do not understand the mathematics to extend this towards this problem.

I would sincerely appreciate any help.

[John D'Errico]

I am surely the person to answer our question, and you did make an effort (which seems to be a rare thing these days) so you deserve a good response. I'll see what I can do after dinner.

One quick question though - I see the words zero-derivative clamping. I assume this refers to the second derivative at the end points, which is fairly standard for a cubic spline. Is that assumption correct? You could, for example, choose to clamp the first derivative at the end points instead. Your answer would make things slightly easier for me, since although I could dig through those equations to see what is implied there for boundary conditions, it would save me some time.

[MR MISTER]

Thanks a lot. I appreciate the quick response.

[MR MISTER]

I believe we clamp the first derivative in this problem (though you are correct about the second derivative). Here is another snippet from my notes...

If you would like more from my notes as to how the professor attained the first six equations, let me know.

Answer 1

[John D'Errico]

That is sufficient information, thanks. So let me see if I can help. (Best of course, is to not write your own spline code, but to use code written by someone who does that for a living. In fact, you might even find some spline code on the File Exchange with my name on it. :) At the same time, it is ALWAYS good to learn how to write something like a spline tool, at least if you are interested.) Ok,...

Define two vectors, (x,y) pairs, so 4 points.

x = [x1 x2 x3 x4], y = [y1 y2 y3 y4]

One thing I recommend STRONGLY. Learn about the use of semi-colons at the end of a line. It suppresses junk output to the command line.

Look at each of those equations in the matrix equation given to you by your instructor. I'll explain what they mean, in context of the spline as you are building it.

So, first that matrix equation relates to a 2 segment spline. So a pair of cubic polynomials, passing through 3 points, connected at the intermediate point.

Given two cubic segments, we would have 8 coefficients in the representation as it is defined, thus 4 coefficients for each cubic segment.

When we multiply that first row of A times the coefficient vector, we will get

x(1)^3*a_13 + x(1)^2*a_12 + x(1)*a_11 + 1*a_10 = y(1)

So just expand the dot product for that first row. What does that equation mean? It states that the first cubic segment, evaluates at the point x(1), yields y(1). Essentially, it requires that the spline interpolate the first data point. So it says that f(x(1))==y(1), for the FIRST cubic segment.

Second row?

x(2)^3*a_13 + x(2)^2*a_12 + x(2)*a_11 + 1*a_10 = y(1)

These are the coefficients of the first cubic segment, used to evaluate the function at the second data point. So this row enforces f(x(2))==y(2).

Rows 3 and 4 are similar, but you can see they apply to the SECOND cubic segment.

x(2)^3*a_23 + x(2)^2*a_22 + x(2)*a_21 + 1*a_20 = y(2)
x(3)^3*a_23 + x(3)^2*a_22 + x(3)*a_21 + 1*a_20 = y(3)

Again, these are interpolation equations, in the sense that they force the second segment to pass through the points (x(2),y(2)) and (x(3),y(3)).

Note that that constraint for point number 2 is applied to both segments. You can think of that as continuity. The two segments are tied together, so the function is continuous at that point.

Had you a 4th point in the curve, then you would have 3 cubic segments, and therefore two more rows similar to the first set of rows in A. Effectively, all of those rows were talking about the function values of the spline at the data points.

What do lines 5,6,7, and 8 do? For this, we need to look at the derivative of those cubic segments. What is the first derivative of the first segment (I'll call it f_1) evaluated at some arbitrary point u? It is just a cubic polynomial

f_1(u) = u^3*a_23 + u^2*a_22 + u*a_21 + a_20

so the derivative is just

f_1'(u) = 3*u^2*a_13 + 2*u*a_12 + a_11

Similarly, the first derivative of f_2 at u is just:

f_2'(u) = 3*u^2*a_23 + 2*u*a_22 + a_21

How about the second derivatives of f_1 and f_2?

f_1''(u) = 6*u*a_13 + 2*a_12
f_2''(u) = 6*u*a_23 + 2*a_22

Ok, so in light of all that knowledge, what does row 5 tell us? Look at it in light of what I've just said. Row 5 is just:

f_1'(x2) - f_2'(x2) = 0

Just as rows 2 and 3 implied continuity of the spline at x(2) as well as forcing the curve to pass through the point (x(2),y(2)), row 5 implies that the spline is differentiable across that point, with a continuous first derivative there. Effectively row 5 says, "While I don't know what that derivative is at x(2), I do know the derivative is continuous across that knot."

Row 6 is similar, but it applies to the second derivative. (Ok, your instructor divided by 2 in that line. Irrelevant.) Row 6 tells us that

f_1''(x2) - f_2''(x2) = 0

So it enforces second derivative continuity at x(2).

Finally, look at rows 7 and 8. EXPAND THOSE LINES AS DOT PRODUCTS, as we did above. Row 7 enforces the value of the first derivative of f_1, thus f_1', evaluated at the point x(1). Likewise, row 8 enforces the value of f_2', evaluated at x(3).

So rows 7 and 8 are the boundary conditions, those zero-clamp conditions that were given to you.

I honestly don't want to write the code for you, because I think you will gain more by writing it yourself. I'd rather offer suggestions as to what you need to have though.

In your code, for a 4 knot spline written out in the above form...

You will have 12 unknowns, 4 unknown coefficients for each of the 3 cubic segments.

You need 6 rows that correspond to the values of the function at each point. Thus...

f_1(x(1)) == y(1)
f_1(x(2)) == y(2)

f_2(x(2)) == y(2)
f_2(x(3)) == y(3)

f_3(x(3)) == y(3)
f_3(x(4)) == y(4)

Next, you need to enforce first AND second derivative continuity at the interior knot points, x(2) and x(3). So we would have 4 more rows in A, corresponding to:

f_1'(x2) - f_2'(x2) = 0
f_1''(x2) - f_2''(x2) = 0

f_2'(x3) - f_3'(x3) = 0
f_2''(x3) - f_3''(x3) = 0

Finally, the zero-clamp first derivative end-conditions.

f_1'(x1) = 0
f_3'(x4) = 0

In total, your matrix should have 12 rows, and 12 columns, since you have 12 coefficients to solve for.

The resulting matrix A will be theoretically non-singular as long as x(1), x(2), x(3), and x(4) are distinct.

So, having said all of this, this formulation of a cubic spline, is NOT what I would like to see done. It may suffer from numerical problems because those polynomials may get nasty. But that may well be the subject of your next assignment, and since I've now written a LOT, time to stop writing. :) But feel free to ask a question if I've not been clear in what I said above. Or, post what you write, and I can look it over to see if you got it right. I'll always try to answer as long as I see a spark of interest in a poster, as long as I see you making an effort.

[MR MISTER]

Hey thanks a lot. I think I do understand it theoretically, it was just the implementation that threw me off. I will show you my code when it is completed, or if I am still having trouble I will show you the updated code. Thanks for taking the time to address my problem, your answer is currently very helpful.

[MR MISTER]

Hey, I believe it works! Thanks for your help.

function [ yy ] = Kable(x,y,xx)
%  xx is a vector of x-values to be interpolated and 
%  yy is an output vector of interpolated values corresponding to xx
%  also plots data points and cubic spline interpolation
%  x = [x1 x2 x3 x4], y = [y1 y2 y3 y4]
%  plots the function generated cubic spline 
%  plots the user input data points 
%  plots built-in matlab cubic spline interpolater for comparison

A = [x(1)^3 x(1)^2 x(1) x(1)^0 0 0 0 0 0 0 0 0;...
     x(2)^3 x(2)^2 x(2) x(2)^0 0 0 0 0 0 0 0 0;...
     0 0 0 0 x(2)^3 x(2)^2 x(2) x(2)^0 0 0 0 0;...
     0 0 0 0 x(3)^3 x(3)^2 x(3) x(3)^0 0 0 0 0;...
     0 0 0 0 0 0 0 0 x(3)^3 x(3)^2 x(3) x(3)^0;...
     0 0 0 0 0 0 0 0 x(4)^3 x(4)^2 x(4) x(4)^0;...
     3*x(2)^2 2*x(2) 1 0 -3*x(2)^2 -2*x(2) -1 0 0 0 0 0 ;...
     0 0 0 0 3*x(3)^2 2*x(3) 1 0 -3*x(3)^2 -2*x(3) -1 0;...
     0 0 0 0 3*x(2) 1 0 0 -3*x(2) -1 0 0 ;...
     0 0 0 0 3*x(3) 1 0 0 -3*x(3) -1 0 0 ;...
     3*x(1)^2 2*x(1) 1 0 0 0 0 0 0 0 0 0;...
     0 0 0 0 0 0 0 0 3*x(4)^2 2*x(4) 1 0 ];
d = [y(1) y(2) y(2) y(3) y(3) y(4) 0 0 0 0 0 0]';
yy = A\d;

z1=linspace(x(1),x(2));
z2=linspace(x(2),x(3));
z3=linspace(x(3),x(4));
s1 = yy(1).*z1.^3 + yy(2).*z1.^2 + yy(3).*z1 + yy(4);
s2 = yy(5).*z2.^3 + yy(6).*z2.^2 + yy(7).*z2 + yy(8);
s3 = yy(9).*z3.^3 + yy(10).*z3.^2 + yy(11).*z3 + yy(12);

s=[s1 s2 s3];
z=[z1 z2 z3];
plot(z,s),hold
plot(x,y,'or'),axis([1 4 -2.5 2])
zz=spline(x,y,xx);
plot(xx,zz);

To compare, I plotted them at the bottom of the function (while the problem doesn't ask for the plot in the function, I thought it would be nice to compare for a user and also for you to see how it looks).

Here is the plot...

• red circle == data points
• blue line == function generated cubic spline
• red line == built-in matlab function spline()

%

[John D'Errico]

I'll just reiterate my comment about this not actually being a good way to implement a spline. Why not?

Suppose we use the above code to fit a spline through the points in your test problem with the x vector as:

x = [0 1 2 3]

Now, change the x vector to be

x = [10000000 10000001 10000002 10000003]

while leaving the vector of y values unchanged. In theory, the shape of the spline in both cases will be identical. In practice there will be serious numerical problems.

What will happen is just trying to evaluate those cubic polynomials for such large values of x will be a problem. As well, the resulting coefficients will be obscenely nasty in the second case, because the linear algebra will blow up.

This is why one would typically wish to formulate a spline in a different way, that would not be susceptible to numerical problems.

[per isakson]

"not fully sure how to write code on this site" &nbsp

1. Select your code. Selected code is displayed white over a blue background.
2. Click the button, {}Code, which will intend every line two position

or start every line of code by two spaces.

[MR MISTER]

Any hint on how one would go about creating a more effective method of determining the spline for all values of x and y? The problem does state to make it work for any four points but mine is failing for sure

[John D'Errico]

Oh, I see you posted your code. Good job for doing the plot. A picture is worth a thousand words here. It tells me, without even checking the code you wrote, that your spline appears ALMOST correct.

Here is the check I did, where I used my own spline code to do the same.

xx = linspace(x(1),x(4),1000);
yy = spline(x,y,xx);

slm = slmengine(x,y,'knots',x,'plot','on','leftslope',0,'rightslope',0);

hold on
plot(xx,yy,'-b')

As you can see (besides the fact that I sleepily swapped the colors in my plot from yours) the correct fit should have a significant dip in the curve to the right of x(2).

I could have made essentially the same plot like this. (I think csape is in the curve fitting toolbox)

pp = csape(x,y,'clamped',[0 0]);
yy0 = ppval(pp,xx);
plot(xx,yy,'-b')
hold on
plot(xx,yy0,'-r')
plot(x,y,'go')

So I checked the curve that you built, and you got the second derivative continuity equations wrong. VERY close though.

It looks like a copy/paste mistake on your part. This is a common error to make when writing code. Copy a similar line that you wrote, then paste it in, but forget to edit the pasted-in line to make it correct. Yeah, like I've never done that before. :)

Here are the lines you had:

0 0 0 0 3*x(2) 1 0 0 -3*x(2) -1 0 0 ;...
0 0 0 0 3*x(3) 1 0 0 -3*x(3) -1 0 0 ;...

As you can see, they were identical. What you should have had was this:

3*x(2) 1 0 0 -3*x(2) -1 0 0 0 0 0 0 ;...
0 0 0 0 3*x(3) 1 0 0 -3*x(3) -1 0 0 ;...

So that first line now should enforce second derivative continuity across x(2), and the second one as you wrote it should be correct for x(3).

[John D'Errico]

I'm responding separately to your followup question. Once you fix the lines I pointed out that were in error, it should work better of course.

Your code will have problems if you provided ANY set of 4 points. For example, this set will fail, because they lie at the 4 corners of a square.

x = [0 1 1 0];
y = [0 0 1 1];

As such, they do not represent y as a single valued function of x. So you cannot use ANY set of points. To start with, your code would require that the vector x be monotonically increasing. So, all(diff(x)>=tol) must at least be true, for some reasonable value of a tolerance tol.

As well, if you merely added a large number to every member of the vector x, your code will fail due to numerical problems. A fix for that is to work by always shifting your data, thus internally subtracting off any large offset before any computations are done. You would add that offset back in when you do any plot of course.

The road to a more perfect code is a long and winding road. One good start point is the absolute classic "A Practical Guide to Splines", by Carl de Boor. I'm sure there are others since. Personally, I still love to use a Hermite form for splines, as that is what I taught myself almost 40 years ago.

An important point is that your code as it is written will be difficult to extend to very many points, simply because that form has every condition explicitly written out.

So to extend your code to handle bigger sets, one might start with loops to build the respective equations. Or, you could formulate the spline in a different form, where the necessary matrix reduces to a tridiagonal system of equations. Of course, that matrix is easily constructed using loops, or you should be able to do it in only a few heavily vectorized lines.

In fact, there are many ways to build a spline fit, all of which will work, if care is applied in the process. I won't claim that any way is the best.