Spline Interpolation with derivative condition for knots

I have been using "spline" function in MATLAB to generate splines, although it is only possible to enforce derivative conditions at the first and final knot of the spline like this,

x = -4:4;
y = [0 .15 1.12 2.36 2.36 1.46 .49 .06 0];
cs = spline(x,[0 y 0]);

How can derivative conditions for each knot (x and y) be specified while evaluating the spline?

Any help will be greatly appreciated. :)

Answer 1

[John D'Errico]

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.)

[Pranjal Biswas]

Thanks, now I understand the limitation in that case w.r.t. cubic splines.

Although when I read documentation about pchip, I dont see how first derivative condition can be enforced (I have sampled derivative values for each knot that I want the piecewise polynomial to have).

"It chooses a set of derivatives at the breask to enforce local monotonicity instead."

In this how can I give my own set of derivative values for pchip to enforce?

Also. w.r.t. mkpp, it uses coeffs as well along with breaks, moreover the documentation doesnt mention a way to enforce first derivative condition on the breaks. Am I missing out on anything?

[John D'Errico]

I did not say it would be trivial in either case. :)

The simplest solution is to use my SLM toolbox. (I'll add a link to it on the FEX.)

For a simple example. Lets see how I would do it, as well as the cost of the cubic interpolant that results. For a function, I'll use exp(x), since we know the derivative at those points too.

x = linspace(0,1,5)';
y = sin(x);
dydx = cos(x);
slm.form = 'slm';
slm.degree = 3;
slm.knots = x;
slm.coef = [y,dydx];
pp = slm2pp(slm)

Now, suppose I use pchip to interpolate those points?

pp = pchip(x,y);
pp2 = fnder(pp,2);
fnplt(pp2)

Looking at the second derivative, we see jump discontinuities.

This is the second derivative of the resulting interpolant. In some cases, you can actually see those second derivative jumps in the shape of the curve. I had one client who found them to be objectionable in some cases. In fact of course, we expect the second derivitive to look roughly like -sin(x).

I'll use one of the utilities in my SLM toolbox to create a Hermite form spline, then convert it to a pp form using slm2pp.

slm.form = 'slm';
slm.degree = 3;
slm.knots = x;
slm.coef = [y,dydx];
pp = slm2pp(slm);

fnplt(pp)
hold on
plot(x,y,'ro')

The interpolation looks reasonable. It has exactly the speified first derivatives at those points.

pp2 = fnder(pp,2);
fnplt(pp2)
hold on
plot(x,-sin(x),'r*')

Which looks reasonably continuous, and it even looks like -sin(x). So the solution using slm2pp is easy enough. Let me push things just a bit now.

dydx = [1 -1 1 -1 1]';
slm.coef = [y,dydx];
pp = slm2pp(slm);
fnplt(pp)
fnplt(pp)
hold on
plot(x,y,'ro')

It clearly interpolates the data, and has the desired slopes at each data point. Looks like crap, but then I was going for that look. ;-)

I could give you code to do a quintic spline too. But as long as the chosen derivitives are reasonable, then you don't really need it. (And I'm being lazy right now.)

Find SLM here:

https://www.mathworks.com/matlabcentral/fileexchange/24443-slm-shape-language-modeling

[Pranjal Biswas]

Thanks a lot for the extensive explanation, that works for me. Also is there any good resource that you can suggest to read all about diiferent types of curves in literature to get a better understanding.

[Pranjal Biswas]

Although the toolbox works, I was trying to work without using any external toolboxes. For that case I was thinking of working with quartic splines to impose an extra condition for the slope at the knots as follows,

1) C0 continuity: f1(x) = f2(x)

2) C1 continuity: f1'(x) = f2'(x)

3) C2 continuity: f1''(x) = f2''(x)

4) Slope condition at knot: f1'(x) = a, f2'(x) = a

5) End conditions for f1 and f2

Therefore 5 conditions for each segment of the spline for 5 coefficients. This should result in on a linear system like, Ma=y, thus a=M^(-1)y.

I can code this for a specific number of points, but is there something already which can do this for n points?