MATLAB: Natural Cubic Spline boundary conditions

cubic splinesinterpolationnatural cubic splinesnumerical analysis

I try to interpolate the function

with given datas. I am using Matlab function 'csape' but, I am not sure this is a Natural Cubic Spline Interpolation or just Cubic Spline.

I just found out that Natural Cubic Splines are the one with the condition that is Secondary Derivative is equal to 0. However, I am not sure that to take 2nd derivative and rearrange my code.

I have this code written by me and want to ensure that is accurate and second derivatives are 0. What should I add to find out?

%Data
x= [-2 6 2.4 7.8];
%Function
y=cos(x);
%Interpolation
cubic=csape(x,y);
%Plotting Data and function outputs
figure;
plot(x,y,'ro');
hold on;
%Plotting interpolation
fnplt(cubic,'b',2);
%Needs
grid on;
xlabel('x');
ylabel('f(x)');
axis([-3 9 -2 1.5]);
title('Natural Cubic Spline Interpolation');
legend('Data Points','Cubic Spline');

Best Answer

The name "natural" applied to the boundary conditions on a spline refer to the use of calculus of variations in the derivation of the spline equations. Why do I say this?

There are two boundary conditions options in csape that allow you to create a natural cubic spline in csape. "variational" explicitly creates the zero second derivative BCs. Again, the word "variational" should be the flag there. I suppose however if you have never used the calculus of variations, you might not notice the reference. For myself, I did not encounter that area of mathematics until I was in grad school.

Can you verify that the default for csape is NOT a natural cubic spline? Lets try csape with the default first.

x = linspace(0,2*pi,5);
%Function
y = cos(x);
spl1 = csape(x,y);
% evaluate the second derivative at the end points
fnval(fnder(spl1,2),[0,2*pi])
ans =
     -0.74302     -0.74302

Note that the second derivative of cos(x) at 0 and 2*pi would be expected to be NON-zero. In fact, the true second derivative of cos(x) at those points is -1. So if the spline was a perfect approximation to the curve, we would have gotten -1. In context, -0.74302 is not really that far off, and for a curve based on only 5 points, I wold say that is as good as I could possibly have expected. Therefore, fitting a natural cubic spline through these points with those end conditions would be a poor choice, because the natural cubic explicitly forces the second derivatives to be zero at those points. It is totally counter to what the data implies at those points.

But I expect this is why you were assigned this assignment, to learn why other choices of boundary conditions generally tend to work better. (It is exactly the thing I had been known to do when I taught about splines in the past.) Clearly the default for csape is NOT a natural spline.

Now, try it using the variational boundary conditions.

spl2 = csape(x,y,'variational');
fnval(fnder(spl2,2),[0,2*pi])
ans =
  -1.3878e-16   1.6653e-16

So the second derivative is now (effectively) zero at the ends, as we specified.

fnplt(spl1)
hold on
fnplt(spl2)
plot(x,y,'o')
legend('Default BCs','Natural BCs','Data')

The green curve is the natural spline. Note that the true underlying cosine function would actually have a zero first derivative at the ends of the interval, not a zero second derivative.

As such a natural spline would have been a relatively poor choice for this particular set of points. Regardless, now you know how to choose the boundary conditions you wanted to achieve. We could also have done it as:

spl2 = csape(x,y,'second',[0 0]);
fnval(fnder(spl2,2),[0,2*pi])
ans =
  -1.3878e-16   1.6653e-16

When I teach about the use of cubic splines, I like to explain that splines really fall into a class of models I like to call metaphorical models. If I look for the definition of metaphor, I find in Grammarly:

"A metaphor is a figure of speech that describes an object or action in a way that isn’t literally true, but helps explain an idea or make a comparison.".

We often build mathematical models for some process where we lack a model based on physical principles. In those cases, such a model can be thought of as a metaphor. It helps us to understand the process under study, even though the model is not directly applicable to that process. For example, Shakespeare, in Romeo & Juliet has Romeo in the balcony scene describing Juliet...

But soft! What light through yonder window breaks?

It is the east, and Juliet is the sun.

Arise, fair sun, and kill the envious moon,

Who is already sick and pale with grief,

That thou, her maid, art far more fair than she.

Be not her maid since she is envious.

Her vestal livery is but sick and green,

And none but fools do wear it.

William Shakespeare, Romeo and Juliet

In there we see Romeo describing Juliet, comparing her to the sun, a bright shining light. While we know that Juliet is not actually a star, with surface temperature of many thousands of degrees, the metaphor teaches us how Romeo sees her.

In the same sense, it is not uncommon to see other mathematical models used in a metaphorical context. For example, there are very good differential equation models used to predict the spread of viruses, thus pandemics. However, I can also point out instances where those same models were used in a metaphorical sense to predict the early sales of consumer appliances, even washing machines. The idea is that just as you catch a virus by being in contact with a person who has that virus, you might also buy a washing machine or a camera because your neighbor just bought one, and they told you how much they liked it.

So we use mathematical metaphors all the time when building models. We use a mathematical model of a real world process that we do understand, to help us to predict and understand the behavior of a real world process where we lack a direct model.

Even the name spline is a reference to the use of long thin strips of wood used to create smooth curves for those who would create the shape of the hull of a ship. You can even see such pieces of wood used by those who build wooden boats. For example, the fellows who are building Arabella use physical splines extensively in the lofting process:

https://www.youtube.com/channel/UCAiDWnTP0WB1xCp6uuUo0VA/videos

https://www.youtube.com/watch?v=RjKIl83qI-s

You can also see Leo using splines when lofting the lines for Tally Ho.

https://www.youtube.com/watch?v=gfuPTLSBQfY

I hope that gives you some background on splines when used in shipbuilding. One characteristic of such a curve is the ends of such a chape as created by a long thin piece of wood, is the ends of that curve will always have zero second derivative! That is, if you force a uniform, long thin piece of wood through a set of points, at each end of the curve it creates, there will be no force applied at the ends of the curve to bend it. Lacking a bending moment at the ends, the curve will approach a straight line at each end of the curve.

That behavior is perfectly "natural" for a piece of wood, when forced to pass through a sequence of points. And a spline indeed does a VERY nice job of fitting many such curves. However the natural cubic spline metaphor fails just a bit when we try to use a spline with natural end conditions and force it through points from a function where we know the end second derivatives are not in fact zero. As you see in the plot above, the green curve has the shape a long thin piece of wood would be expected to have when passed through those points. However, it is the wrong fundamental shape for that curve. The difference is not terrible, but it is also not terribly accurate for cos(x).

Now, let me redo the same modeling process for sin(x) over the same domain.

x = linspace(0,2*pi,5);
ysin = sin(x);
spl1 = csape(x,ysin);
spl2 = csape(x,ysin,'variational');
fnplt(spl1)
hold on
plot(x,ysin,'ro')
fnplt(spl2)

In this case, the red curve is a very good approximation to the real shape of sin(x). It looks far closer to my eyes to what I would want, than is the curve in blue. Here you can see the end point second derivatives for the default spline interpolant, compared to the natural spline second derivatives.

fnval(fnder(spl1,2),[0,2*pi]) % default spline f"
ans =
      -1.6211       1.6211
fnval(fnder(spl2,2),[0,2*pi]) % natural cubic f"
ans =
     0     0

How well does the natural cubic predict the first derivatives at the end points? For sin(x), the true curve would have them exactly 1.

syms u
subs(diff(sin(u),u),u,0)
ans =
1
fnval(fnder(spl1,1),[0,2*pi]) % Default csape f' estimates
ans =
       1.6977       1.6977
       
fnval(fnder(spl2,1),[0,2*pi]) % Natural cubic f' estimates
ans =
      0.95493      0.95493

Those end point FIRST derivatives are now almost spot on for the natural cubic spline.

So when the end conditions make perfect sense in context, we truly did get a much better approximation to the desired curve.

In the same way in literature, a good metaphor allows us to accurately understand and predict some aspect about what it describes. What I hope you may take away from this is to choose your metaphors carefully. Now if only all of us could be as good at our jobs as William Shakespeare was at his job. But that is a horse of a different color.

Related Solutions

MATLAB: Plot Natural Cubic Spline

You were close, in that csape can do it. (My SLM toolbox also has a natural cubic as an option, but there is no need to go that far.) The proper choice in csape is:

pp = csape(x,y,'variational');

The flag 'variational' comes from the calculus of variations, I assume. The calculus of variations can also be seen to call those corresponding end conditions the "natural" boundary conditions.

   'variational' : set end second derivatives equal to zero
                    (ignoring given end condition values if any)

A natural cubic spline is the one that has zero second derivatives at the ends.

As a check, does it produce what I expect? So here a function that has decidedly NON-zero second derivatives at the ends.

x = linspace(0,2*pi,10);
y = cos(x);
pp = csape(x,y,'variational');
pp2 = fnder(pp,2);
fnplt(pp2)
grid on

Yet as you see, the second derivatives at the ends are forced to zero. This is the classic example that I always used when teaching about splines, as to why a natural cubic spline might often be a poor choice.

So if we look at the fit that arises, compare the natural cubic to that produced by spline.

fnplt(pp,[0 1])
hold on
pps = spline(x,y);
fnplt(pps,[0 1])

In blue is the natural cubic, whereas the green curve is the result of spline, which uses not-a-knot end conditions, generally a safer choice. If we remember these curves are an approximation to cos(x), the green curve is clearly much better.

So while you can find a natural cubic in MATLAB, you need to look deeply enough, and know what you should be looking for. That is in general a good thing in this regard.

MATLAB: Spline Interpolation with derivative condition for knots

Technically, this is impossible to accompish, at least in terms of a cubic spline.

A cubic spline is twice continuously differentiable across the breaks. For n data points, we require the spline to pass through the data points, and be continuous across the breaks, as well as have a continuous derivative, but also to have a continuous second derivative. All of this leaves 2 degrees of freedom, allowing you to choose exactly TWO end conditions to resolve the ambiguity. Typically, the best choice is arguably the not-a-knot end conditions, however, a not uncommon choice is the natural end conditions, thus zero second derivatives at the ends. (That is usally not as good a choice as not-a-knot though.)

However, you want to enfore a spline that passes through n points, as well as specifying the first derivative at each knot? The cost will be that you can no longer enforce second derivative continuity. You will no longer have a spline.

As such, the curve will be less smooth than is a classical spline. There are tools in MATLAB that build interpolants of this general form. For example, pchip is one. A pchip interpolant is not twice differentiable across the breaks. It chooses a set of derivatives at the breask to enforce local monotonicity instead.

So you essentially would have two options. The first is to create an interpolant (like pchip) that is piecewise cubic, but lacks second derivative continuity. Or, you can build a quintic (degree 5) spline, which has sufficient flexibility to allow you to choose the first derivatives, and still be smooth. (In either case, I'd probably use mkpp in the end.)

Best Answer

Related Solutions

MATLAB: Plot Natural Cubic Spline

MATLAB: Spline Interpolation with derivative condition for knots

Related Question