Why is the cubic Spline Interpolation wiggly in flat areas?
19 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Dominik Hiltbrunner
el 24 de Dic. de 2021
Editada: John D'Errico
el 24 de Dic. de 2021
Dear community
I am using Matlab's griddedInterpolant() to fit 2-D data with non-uniform samples. Please consider the following example:
clc; clear; close all; format shorteng;
Z = [728.8979e-003 636.7566e-003 606.6035e-003 565.2786e-003 550.3203e-003 544.8502e-003 542.6799e-003 541.8001e-003 541.4345e-003
521.0865e-003 462.6134e-003 468.5421e-003 472.2405e-003 475.2979e-003 476.8306e-003 477.0386e-003 477.5314e-003 477.6013e-003
478.5540e-003 431.6547e-003 451.0393e-003 462.2832e-003 466.5187e-003 469.1578e-003 470.0577e-003 470.9129e-003 470.9733e-003
454.4821e-003 422.4798e-003 447.9669e-003 459.6772e-003 465.3374e-003 467.3769e-003 468.7602e-003 469.0878e-003 469.5432e-003
444.1855e-003 420.1184e-003 446.5284e-003 459.0328e-003 464.4620e-003 466.8222e-003 467.8873e-003 468.6215e-003 468.6755e-003
441.7790e-003 418.7854e-003 445.8779e-003 458.3203e-003 463.9560e-003 466.7335e-003 467.9822e-003 468.2045e-003 468.2045e-003
440.6139e-003 418.6750e-003 445.4984e-003 457.8616e-003 463.8469e-003 466.8397e-003 467.4772e-003 467.8873e-003 468.6773e-003
440.1721e-003 418.5457e-003 445.2773e-003 458.1930e-003 463.8622e-003 466.3328e-003 467.3769e-003 468.3087e-003 468.4529e-003
439.5239e-003 418.0573e-003 445.3255e-003 457.9355e-003 463.7742e-003 466.3328e-003 467.3769e-003 468.3087e-003 468.4529e-003];
[X,Y] = ndgrid(logspace(log10(210e-9),log10(100e-6),9));
surf(X,Y,Z,'FaceAlpha',0.5)
set(gca,'XScale','log')
set(gca,'YScale','log')
title('Original')
F = griddedInterpolant(X,Y,Z,'spline');
[XX,YY] = ndgrid(logspace(log10(210e-9),log10(100e-6),3*9));
Z_approx = F(XX,YY);
figure
surf(XX,YY,Z_approx,'EdgeColor','black')
set(gca,'XScale','log')
set(gca,'YScale','log')
title('Approximation')
As you can see, the original data is mostly flat except for small X and Y values. However, when I compute the cubic Splines and use a finer grid to sample data points from the interpolated curve, then this flat area becomes wiggly.
a) Why? Shouldn't the cubic Spline Interpolation be able to come up with lower degree polynomials (constants) at those points? I mean, with degree zero polynomials (constants) it should also be possible to fulfil the C2 smoothness condition.
b) The documentation states that the 'spline' argument uses not-a-knot end conditions. Does Matlab also support other end conditions?
c) Any suggestions on how to get rid of that wiggly surface? Are there better methods? I also tried with multi-variate least-squares but the results aren't really better. Also, due to the nature of where the data comes from I am unable to have uniformly sampled data.
0 comentarios
Respuesta aceptada
John D'Errico
el 24 de Dic. de 2021
Editada: John D'Errico
el 24 de Dic. de 2021
Should not the cubic spline be able to be ABSOLUTELY flat in some areas? NO. Sorry, but that will not be the case. A cubic spline is twice differentiable everywhere. And that forces the spline to oscillate into regions where the data is otherwise flat. A simple example in 1-d will help to illustrate what happens.
x = 1:10;
y = [ones(1,5),zeros(1,5)];
spl = spline(x,y);
fnplt(spl,'r')
hold on
plot(x,y,'go')
grid on
hold off
So the spline oscillates everywhere, even though you might hope that in the intervals [1,5] and [6,10], the function could arguably be absoutely flat, based on the data we have.
This is NOT a question of end conditions. If I used different end conditions, this issue would still arise.
What happens is in order for the function to be twice differentiable EVERYWHERE, (and three times differentiable at all points that are not breakpoints) we get a function that is forced to oscillate.
There are alternatives to a simple spline though. pchip and makima, for example are subtly different interpolants that are not twice differentiable at the data points. They are only once differentiable at those locations. For the particular data I have made up, they will probably produce the same results, but that will not always be the case.
spl2 = makima(x,y);
fnplt(spl2,'b')
hold on
plot(x,y,'go')
grid on
hold off
Now look at the derivatives of these functions.
spldiff = fnder(spl,2);
fnplt(spldiff,'r')
hold on
spl2diff = fnder(spl2,2);
fnplt(spl2diff,'b')
grid on
hold off
Do you see that at the locations 5 and 6, spl2 (in blue) has a second derivative discontinuity? At the same time, the second derivative was always purely continuous for the spline.
You may want for a spline to behave however you wish, but mathematics is cruel. Mathematics follows its own hard set of rules.
Again, end conditions are completely, utterly, absolutely irrelevant here.
If you want to understand what is happening, consider the curve as produced by makima. Just slightly to the left of x==6, and just slightly to the right of x==5, the curve must obviously have a maximal amount of curvature. Correct? Heavily curved fuctinos at those locations cannot suddenly turn into absoutely flat functions right next door, at least not and be twice continuously differentiable.
And while the data I have used here is far more extreme that that which you have, it is also only a 1-d probnlem. Things are always more difficult in 2-d.
Are there better methods available? The simplest answer is to just use what is known as a bilinear interpolant. That is, essentially a lower order spline. In this case, a LINEAR spline in two dimensions. Interp2 produces this, when the 'linear' option is chosen. Would you still complain? Possibly. I can show you cases where such a bilinear interpolant produces artifacts that will confuse most people. Is the bilinear interpolant produced by interp2 better? In some respects, yes, in others, no. You will no longer see a smooth surface, but more of a faceted result. But where your data is flat, the surface produced will now still be flat.
Más respuestas (0)
Ver también
Categorías
Más información sobre Interpolation en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!