cubic spline interpolation - mixed boundary conditions possible?

60 次查看（过去 30 天）

SA-W 2024-2-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2081156-cubic-spline-interpolation-mixed-boundary-conditions-possible

评论： SA-W 2024-3-10

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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

SA-W 2024-2-12

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 2024-2-12

编辑：Torsten 2024-2-12

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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Bruno Luong 2024-2-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2081156-cubic-spline-interpolation-mixed-boundary-conditions-possible#answer_1407096

编辑：Bruno Luong 2024-2-12

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.

21 个评论
显示 19更早的评论隐藏 19更早的评论

Bruno Luong 2024-2-24

编辑：Bruno Luong 2024-2-24

"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 2024-2-24

在 MATLAB Online 中打开

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.

Bruno Luong 2024-2-26

编辑：Bruno Luong 2024-2-26

"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 2024-3-10

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

请先登录，再进行评论。

更多回答（2 个）

Matt J 2024-2-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2081156-cubic-spline-interpolation-mixed-boundary-conditions-possible#answer_1407051

This FEX download will let you impose conditions like that,

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

8 个评论
显示 6更早的评论隐藏 6更早的评论

Bruno Luong 2024-2-24

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

Matt J 2024-2-24

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

请先登录，再进行评论。

John D'Errico 2024-2-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2081156-cubic-spline-interpolation-mixed-boundary-conditions-possible#answer_1407061

编辑：John D'Errico 2024-2-12

在 MATLAB Online 中打开

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.

3 个评论
显示 1更早的评论隐藏 1更早的评论

John D'Errico 2024-2-12

编辑：John D'Errico 2024-2-12

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 2024-2-12

在 MATLAB Online 中打开

@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?

请先登录，再进行评论。

请先登录，再回答此问题。