cubic spline interpolation - mixed boundary conditions possible?

We need two boundary conditions. Is it possible to specify the first and second derivative at the same boundary point?

In csape and others, I have not seen such mixed boundary conditions. I also think it is not possible as these two equations lead to conflicting entries in the matrix and right-hand-side.

I want to "connect" two interpolation domains ( f1 on x1 between [a,b] and f2 on x2 between [b, c]) together by first creating f1 and using f1'(b) and f2''(b) as boundary conditions for the construction of f2. But I am not sure if this makes sense.

[Torsten]

Why don't you interpolate on [a,c] right from the beginning with the usual boundary conditions on a and c ?

Does f1(b) not equal f2(b) ?

[SA-W]

Does f1(b) not equal f2(b) ?

No.

The y-values of my data are very linear at the beginning ( on [a,b] ) and become non-linear then ( on [b,c] ). One spline on [a,c] is not appropriate to represent the transition region, which closes the circle to my question.

[Torsten]

If there is a jump in the function value at b, it does not make sense to connect the two pieces by conditions on derivative and second derivative. Use two separate splines on [a,b] and [b,c].

[SA-W]

If there is a jump in the function value at b, it does not make sense to connect the two pieces by conditions on derivative and second derivative

Sorry, f1(b) = f2(b), but the derivatives are not the same. So do you see a way to feed f1'(b) and f1''(b) into f2 as BC?

[Torsten]

But if the derivatives of f1 and f2 are different at b, why do you want to feed f1'(b) into f2 as boundary condition ? Use two separate splines on [a, b] and [b, c]. They will automatically incorporate the different slopes of f1 and f2 at b.

[SA-W]

Because, the function is involed in a Newton-method, so it must be continuously differentiable. At the moment, this is not the case because the derivatives at the break do not match.

[Torsten]

All information you give is contradictory - so it doesn't make sense to continue replying.

Either you want to have f1'(b) = f2'(b) to make your Newton-method work. In this case, use one spline over the complete interval [a c].

If you want to incorporate that the derivatives are different, use two usual splines - one over [a, b] and one over [b, c]. They will automatically reflect a jump of f' at b.

Answer 1

[Bruno Luong]

In principle yes. The spline on the second segment (b,c) are piecewise cubic poynomials; starting from b sequentially to each intervals [xi,xj], x0=b < x1 < x2 < ... < xn = c; up to c, each polynomial have 4 requirements

for 0 < i+1 = j <= n :

• f(xi) = yi
• f(xj) = yj
• f'(xi) = ydi
• f''(xi) = yddi

So the cubic polynomial is complelely deterined on the first interval then second, etc... On each interval coefficents can be estimated by using linear algebra with "\" or even directly computed (see comment below).

No need for any toolbox.

Now as John have warned, it can become unstable and the spline might strongly oscillate when you get to the last interval.

[SA-W]

Nice, did not think about that!

Of course, the question arises whether the curve's derivatives constructed with this method look way different than those coming from a single fit over the entire domain [a,c].

[Bruno Luong]

If it is not stable, chance that the curves are way different

b = 0;
c = 1;
n = 5;
fd0 = 0;  % left first derivative
fdd0 = 0; % left second derivative
% random input data
x = linspace(b, c, n+1);
y = randn(size(x));

%% standard spline interpolation
pp = spline(x, y);

% copy pp struct 
pp_bc0 = pp;
fi = y(1);
fdi = fd0;
fddi = fdd0;
% upate coefs field, rexursively from left to right
for i=1:n
    xi = x(i);
    xj = x(i+1);
    h = xj - xi;
    fj = y(i+1);
    ci(4) = fi;
    ci(3) = fdi;
    ci(2) = fddi/2;
    ci(1) = (fj-polyval(ci(2:end),h))/h^3;
    % update left bc for the next iteration
    fi = fj;
    pp_bc0.coefs(i,:) = ci;
    p = polyder(ci);
    fdi = polyval(p, h);
    p = polyder(p);
    fddi = polyval(p, h);
end

% Graphical output
xi = linspace(min(x), max(x), 257);
figure(1)
clf
plot(x, y, 'or', xi, ppval(pp_bc0, xi), 'b', xi, ppval(pp, xi), 'k');
legend('data', 'spline with left bc', 'standard not-a-knot spline')

[Bruno Luong]

Alternativly the for-loop can be achieve using BSFK FEX with

pntcon = struct('p', {1 2},'x', b, 'v', {fd0 fdd0});
opt = struct('KnotRemoval','none','pntcon',pntcon);
pp_bc0 = BSFK(x,y,4,n,x,opt);

[SA-W]

@Bruno Luong

This approach with the for-loop allows me to take the derivative of the spline w.r.t y-values. Right?

[Bruno Luong]

@SA-W Unconstrained spline function f(x) is linear wrt to interpolation value y and in this case also two input boundary values fd0 and fd00,for all x.

I just repeat myself again.

[SA-W]

@Bruno Luong

I have some concern as for the calculation of the derivative df/dy. In particular, how can we compute df/dy(xj) in every loop iteration?

For the interpolating spline f, we know all interpolation values. As for the derivatives df/dy, all we know is the first and second derivative + the value on the left side xi, but the value on the right side xj is not known. Do I miss something here?

[Bruno Luong]

If you know 4 values of the cubic polynomial (2 interpolation data + first and second derivatives on the left side) then you know the function (and its derivatives) everywhere, as show in the calculation inside the loop, the values of fdi and fddi are the derivative on the right side of the current iteration

p = polyder(ci);    
fdi = polyval(p, h);
p = polyder(p);
fddi = polyval(p, h);

[SA-W]

I know. What you say is correct for the interpolationg spline f, but I talk about the sensitivity df/dy. We know the first and second derivative of df/dy, the value on the left. But how do we know the value of df/dy on the right when I do the first loop iteration?

[Bruno Luong]

Your questions [quote "but the value on the right side xj is not known"] is about first and second derivative on the right side, so I understand you want to compute df/dx and d^2f/dx^2 on the right knot, and I give the answer.

WHere as df/dy I repeat myself again The cubic spline function is linear with respect to y.

I try to understand why you make a confusion betwen df/dx d^2f/dx^2 on one hand and df/dy in other hand. It seems it's not the first time you confuse the partial derivatives

[SA-W]

Everything about the construction of f is clear. My question pertains to the derivatives df/dy.

Suppose we have constructed df/dy on all segments of the first domain [a,b]. Now we are at [b,c] and do the first iteration of the loop. To construct a cubic polynomial that describes df/dy on the first segment of [b,c], we know the value of df/dy and the first and second derivative of df/dy by evaluating df/dy at the point b. But these are only three conditions and we also need the value of df/dy at xj as a fourth condition. It is not clear to me where this value comes from. Maybe I misunderstand the construction of df/dy on the domain [b,c] and you can comment again on this.

[Bruno Luong]

"we also need the value of df/dy at xj as a fourth condition.

I din't understand that. If you insist: df/dyi at xj is 1 if i == j, 0 if not, since f is interpolation spline, meaning f(xj) = yj. Period. And why you need it???

Something is not clear in your reasoning, and I don't understand exactly what. All I observe is that many things you wrote did not make much sense to me.

BTW the algorithm itself (with for-loop or not by other mean) is irrelevant to df/dy. You can consider the interpolation cubic-spline function f with 2 bc at b on [b,c] (thus so global f on [a,c]) is a black-box, and compute df/fy by finite difference as you already experiment in other thread. The finite difference result is independent of step sizes and independent of the specific operating data; since unconstrained spline is linear. It ends up as doing spline interpolation with the method you describe ([a,b] and [b,c] with matching derivatives) with every column of eye matrix as input y-data, each function returned is the derivative df/dyi(x)

[SA-W]

n = 10;
x = linspace(1,10,n);
xi = linspace(min(x), max(x), 1e3);
y = rand(length(x));
pp = csape(x,y); % interpolating spline f
ei = eye(length(x));
dfdy = csape(x,ei); % derivative df/dy1i for i=1,...,n
dfdy = ppval(dfdy, xi); 
plot(xi, dfdy');

If you insist: df/dyi at xj is 1 if i == j, 0 if not, since f is interpolation spline, meaning f(xj) = yj. Period. And why you need it???

I do not think this is always true. See the above code, df/dyi at xi is not exactly 1 for some curves.

It ends up as doing spline interpolation with the method you describe ([a,b] and [b,c] with matching derivatives) with every column of eye matrix as input y-data, each function returned is the derivative df/dyi(x)

I guess I understand now what you wanted to tell me: We first construct f on [a,b], then on [b,c] to have f on [a,c]. After that, we construct df/dy using the identify matrix approach. Is that correct?

What I incorrectly thought is that we construct f and df/dy on [a,b]. Then we construct both f and df/dy on [b,c] by using the approach sketched above.

[Torsten]

I do not think this is always true. See the above code, df/dyi at xi is not exactly 1 for some curves.

if you add the command

ppval(dfdy,x)

after the line

dfdy = csape(x,ei); % derivative df/dy1i for i=1,...,n

you will see that you get eye(10).

[Bruno Luong]

"I do not think this is always true. See the above code, df/dyi at xi is not exactly 1 for some curves."

I disagree, the circle marker (value at the knots) are all in y=1 line or intersect to y=0 lines.

Of course, it INTERPOLATES the basis vectors (column of eye)

You seem to look at peak values of the df/dyi and not values at the knots x.

n = 10;
x = linspace(1,10,n);
xi = linspace(min(x), max(x), 1e3);

%y = rand(length(x)); % not really relevant
%pp = csape(x,y); % interpolating spline f % not really relevant

ei = eye(length(x));
dfdy = csape(x,ei); % derivative df/dy1i for i=1,...,n
dfdyv = ppval(dfdy, xi);

plot(xi, dfdyv');
hold on
for i=1:length(x)
    yi = ppval(dfdy, x(i));
    yi = yi(i);
    plot(x(i), yi, 'o')
    xline(x(i))
end
yline(1)

[SA-W]

Ok. I did not pay attention that the maximum is not precisely at the knots but shifted for some curves. Out of curiosity, is there an intuitive explanation why that is the case? E.g., why has the second curve associated with y2 its maximum not at x=2?

[Bruno Luong]

There is nothing that ensures mathematically that the maximum is located at the knot. It is shift due partly of the boundary condition (here not-a-knot) that amplifies the df/fy function NOT symmetrical at the knot, and shift the maximum. I guess natural bc will be better.

[SA-W]

One last question:

At the moment, my software design is that each interpolator (class) only knows about the points and values on the subdomain. E.g. if f1 interpolates on [a,b], f1 only knows the points subset of points x and values y living on [a,b].

To construct the spline f, this works, but for the cuves df/dyi (analytically using eye matrix approach) this can not work, can it?

[Bruno Luong]

Computing df/dyi is no different than construction f. For a given i index, df/dyi is a special spline interpolant "f" with input y(j) = 0 for all j ~= i and y(i)=1.

Once you take this as an evident fact, and your claim; quote "To construct the spline f, this works," you won't ask such odd question.

Now may be you can't use directly eye as input, but may be you have to do for-loop,or reorganize, extend your class etc... This is simply code structuring, code calling, class definition design, I don't see it get into the way of acheiving the desire calculation of df/dy. MATLAB after all is a programming language

[SA-W]

Let me demonstrate my concern:

x = [a,c] = [1,10]
dfdy_all = csape(x,eye(numel(x));
xl = [a,b] = [1,5]
dfdy_left = csape(x,eye(numel(xl));
xr = [b,c] = [5,10]
% not start the loop ...

f is linear w.r.t y, so dfdy_all is the correct solution on the entire domain [1,10]. Clearly, dfdy_all ploted over [a,b] and dfdy_left are not the same curves. So regardless of what we do in the loop to construct dfdy_right, we construct a solution different to dfdy_all. Makes sense?

[Bruno Luong]

"Makes sense? "

Certainly NOT.

I don't understand what you want to demonstrate with your code snippet. There is no output, no comment and then I don't see the logic. For the starter you at least mix-match the boundary conditions

• not-a-knot on [a,c] and
• chaining not_a_knot (?) on [a,b] then[b,c] with C2 condition(s) at b as described in the orignal post.

Those two are NOT the same. Then your code compare apple with orange. Torsen's asked you over and over why don't you do global interpolation on [a,c] intead of chaining [a,b] & [b,c]. So you knew they are NOT the same.

[SA-W]

@Bruno Luong

I appreciate your help. Can you please comment on that question?

https://de.mathworks.com/matlabcentral/answers/2092321-steps-to-convert-spline-from-b-form-to-pp-form-fn2fm?s_tid=srchtitle

I think the question is clear and concise and I will not start a longer discussion.

Answer 2

[Matt J]

This FEX download will let you impose conditions like that,

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

[SA-W]

I think what I want to do (two separate interpolations + imposing C2 contnuity at the break) does not make sense as Torsten pointed out.

Or do you see a way to do it with the slm engine?

[Matt J]

Only you know if it makes sense for you. But the SLM engine will definitely let you do it. Once you do the spline fit on [a,b], you will know its derivatives at a,b. You can then impose the derivatives at b as constraints when you fit [b,c].

[Bruno Luong]

But OP asks for interpolation not fitting

[Matt J]

Fitting and interpolation become equivalent when you put the control points at the data points.

[Bruno Luong]

Not if you provide extra BC (meaning more than two independent).

[Matt J]

I'm not sure why not. You probably have to add more degrees of freedom to the model to force the spline through all the points and satisfy boundary conditions, but SLM claims to be able to do that.

[Bruno Luong]

Of course just be specific, what DOF you add and how you add with SLM

[Matt J]

It appears you use slmset() with the 'xy','xyp' and 'xypp' settings to force the curve to have certain values and derivatives at prescribed points.

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

Answer 3

[John D'Errico]

You generally don't want to do this. That is, fit one spline, then construct a second spline to match the first. The result will not be the optimal one, because the second curve shape is now entirely dependent on the first. The most extreme case of this would be a two segment spline. So you have some data in the first segment. Fit a cubic polynomial through the data for the first segment, ignoring everything to the right of it.. Then try to fit the second segment, while forcing the curve to be C2 across the boundary. That would result in complete crap. We can do it easily enough. For example...

x = linspace(0,2,100)';
y = sin(2*x);
ind1 = x<=1; ind2 = x>=1;

Seg1 = fit(x(ind1),y(ind1),'poly3') % polyfit would also be efficient here

Seg1 = 
     Linear model Poly3:
     Seg1(x) = p1*x^3 + p2*x^2 + p3*x + p4
     Coefficients (with 95% confidence bounds):
       p1 =     -0.6913  (-0.7314, -0.6511)
       p2 =     -0.5272  (-0.5877, -0.4667)
       p3 =       2.125  (2.1, 2.151)
       p4 =   -0.005957  (-0.008857, -0.003058)

% we want the curve to be C2 across the break at x==1
Bc0 = Seg1(1);
[Bc1,Bc2] = differentiate(Seg1,1);

Now, in this silly example of what not to do, we choose to fit a cubic segment to the second half of our data. I'll use lsqlin to do that, to enforce the boundary constraints match.

beq = [Bc0;Bc1;Bc2];
Aeq = [1 1 1 1;3 2 1 0;6 2 0 0];
C = x(ind2).^[3 2 1 0];
D = y(ind2);

coef2 = lsqlin(C,D,[],[],Aeq,beq)

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

coef2 = 4x1

2.1070 -8.9220 10.5201 -2.8042

plot(x,y,'k.',x(ind1),Seg1(x(ind1)),'r-',x(ind2),polyval(coef2,x(ind2)),'b-')

As you can see, the red ucrve fits nicely to the first part of the data, but then a fit to the second part, where the function is constrained to match so it is C2 with the first segment is terrible.

Instead, you want to perform a fit where both sections are fit at the same time, with a constraint that the function is C2 at the join point. A simple least squares spline with the same knots, will fit very nicely.

SIGH. But, CAN you do that? Well, yes, using my SLM toolbox. That part is easy enough. I don't know that you can do it using the MATLAB included tools though. Even so, I would suggest it is better to do the entire modeling part together, rather than in two halves. If not, then you can see crapola result, as I showed above.

[SA-W]

@John D'Errico I totally agree with you that making one fit over the entire domain is the better approach.

However, I wanted to take the shape of my data into account. Consider a piecewise linear function (two segments) where we make the fit independently in each segment. We will see that standard cubic spline interpolators will produce constant first and zero second derivative on each segment. This is clearly not the case if we wanted to fit both segments at once. Makes sense?

[John D'Errico]

NO. You apparently have multiple misconceptions about splines. A cubic spline interpolant will not enforce both constant first AND second derivatives at the end points. That is simply wrong.

Yes, one approach are the natural end conditions, which says the 2nd derivative at the ends is set to zero. This can be derived from a calculus of variations solution for the shape of a flexible beam, thus a mathematical model for a thin flexible beam. A zero endpoint second derivative makes sense there.

However, a more common (and arguably better most of the time) solution is to use what is known as the not-a-knot end condtions. Effectively, the idea is to enforce third derivative continuity at the 2nd and next to last break. Done that way, there is no derivative constraint applied at either boundary.

So what you are saying does not make sense. It merely says you don't fully understand splines.

[SA-W]

@John D'Errico

A cubic spline interpolant will not enforce both constant first AND second derivatives at the end points. That is simply wrong.

x1 = 1:5;
y1 = x1;
figure;
fnplt(fnder(csape(x1,y1), 1)); % perfectly 1
hold on;
fnplt(fnder(csape(x1,y1), 2)); % perfectly 0


x2 = 5:10;
y2 = 2.*(x2-5) + y1(end);
figure;
fnplt(fnder(csape(x2,y2), 1)); % perfectly 2
hold on;
fnplt(fnder(csape(x2,y2), 2)); % perfectly 0

x = [x1, x2];
y = [y1, y2];
figure;
fnplt(fnder(csape(x,y), 1)); 
hold on;
fnplt(fnder(csape(x,y), 2)); 

I do not know what you think was was wrong in my comment. What I said is that a cubic spline interpolator will produce constant first and zero second derivative when applied to a linear data set. Secondly, I said that a cubic spline interpolator applied to a data set composed of two linear regions with different slopes may not be the best choice. The code above demonstrates that. So would you say that the last fit (x,y) is a good choice here?