Effacer les filtres
Effacer les filtres

cubic spline interpolation - mixed boundary conditions possible?

21 vues (au cours des 30 derniers jours)
SA-W
SA-W le 12 Fév 2024
Commenté : SA-W le 10 Mar 2024
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 commentaires
SA-W
SA-W le 12 Fév 2024
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
Torsten le 12 Fév 2024
Modifié(e) : Torsten le 12 Fév 2024
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.

Connectez-vous pour commenter.

Réponse acceptée

Bruno Luong
Bruno Luong le 12 Fév 2024
Modifié(e) : Bruno Luong le 12 Fév 2024
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 commentaires
Bruno Luong
Bruno Luong le 26 Fév 2024
Modifié(e) : Bruno Luong le 26 Fév 2024
"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
SA-W le 10 Mar 2024
I appreciate your help. Can you please comment on that question?
I think the question is clear and concise and I will not start a longer discussion.

Connectez-vous pour commenter.

Plus de réponses (2)

Matt J
Matt J le 12 Fév 2024
This FEX download will let you impose conditions like that,
  8 commentaires
Bruno Luong
Bruno Luong le 24 Fév 2024
Of course just be specific, what DOF you add and how you add with SLM
Matt J
Matt J le 24 Fév 2024
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.

Connectez-vous pour commenter.


John D'Errico
John D'Errico le 12 Fév 2024
Modifié(e) : John D'Errico le 12 Fév 2024
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 commentaires
John D'Errico
John D'Errico le 12 Fév 2024
Modifié(e) : John D'Errico le 12 Fév 2024
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
SA-W le 12 Fév 2024
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?

Connectez-vous pour commenter.

Catégories

En savoir plus sur Spline Construction dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by