Plot Natural Cubic Spline

Does Matlab have a built in code to plot a 'Natural' Cubic spline? I was using pp = spline(... , but have realised this is not a 'natural spline. I have also tried pp = csape(x,y,'second'), where x and y are my nodes, but im not sure this is right either.

 Respuesta aceptada

John D'Errico
John D'Errico el 8 de Mzo. de 2018
Editada: John D'Errico el 8 de Mzo. de 2018

2 votos

You were close, in that csape can do it. (My SLM toolbox also has a natural cubic as an option, but there is no need to go that far.) The proper choice in csape is:
pp = csape(x,y,'variational');
The flag 'variational' comes from the calculus of variations, I assume. The calculus of variations can also be seen to call those corresponding end conditions the "natural" boundary conditions.
'variational' : set end second derivatives equal to zero
(ignoring given end condition values if any)
A natural cubic spline is the one that has zero second derivatives at the ends.
As a check, does it produce what I expect? So here a function that has decidedly NON-zero second derivatives at the ends.
x = linspace(0,2*pi,10);
y = cos(x);
pp = csape(x,y,'variational');
pp2 = fnder(pp,2);
fnplt(pp2)
grid on
Yet as you see, the second derivatives at the ends are forced to zero. This is the classic example that I always used when teaching about splines, as to why a natural cubic spline might often be a poor choice.
So if we look at the fit that arises, compare the natural cubic to that produced by spline.
fnplt(pp,[0 1])
hold on
pps = spline(x,y);
fnplt(pps,[0 1])
In blue is the natural cubic, whereas the green curve is the result of spline, which uses not-a-knot end conditions, generally a safer choice. If we remember these curves are an approximation to cos(x), the green curve is clearly much better.
So while you can find a natural cubic in MATLAB, you need to look deeply enough, and know what you should be looking for. That is in general a good thing in this regard.

4 comentarios

JJ
JJ el 8 de Mzo. de 2018
Many thanks. One quick question? Does this matter that I am plotting between points and my y is not a function?
John D'Errico
John D'Errico el 8 de Mzo. de 2018
Editada: John D'Errico el 8 de Mzo. de 2018
Why would it matter in the least?
Examples that use functions are common, because you know what to expect as in the examples I gave. But this will produce an entirely valid natural cubic spline:
pp = csape(linspace(0,1,10),rand(1,10),'variational');
It is complete garbage, in the sense that there is no valid reason to interpolate random numbers. But the result is indeed a natural cubic spline. You can easily prove that for yourself. I won't even plot the result, since I hope you can guess what to expect.
Splines are utterly commonly applied to all sorts of data that did not arise from a function with no error in it. A risk of course is when the spline begins to exhibit artifacts. You need to know what to expect then, and how to deal with it.
x = 0:10;
y = double(x > 5);
pp = csape(x,y,'variationa;');
fnplt(pp)
hold on
plot(x,y,'ro')
So a common problem with splines, and a good reason why spline-like tools like pchip were invented.
If your question about y not being a function implies that the relation you want to plot is not a single valued function, then indeed it DOES matter. For example, if your curve arises as points around the perimeter of a circle, then the classic spline tools are NOT appropriate.
t = linspace(0,2*pi,20);
x = cos(t);
y = sin(t);
Do NOT use spline or its immediate cousins to interpolate those points. You can use tools like cscvn, or my own interparc, found on the file exchange. Both will produce viable results.
JJ
JJ el 8 de Mzo. de 2018
Yes, thank you. I have it now. Plotting the default one along side the Natural one helped me see what was going on. Thanks again
Joaquín Vidal Peña
Joaquín Vidal Peña el 15 de Nov. de 2022
BRO MANY THANKS I WAS LOOKING FOR THIS
I wanted to corroborate my work of Cubic Natural Spline by hand, and the command 'spline' didn't seem to be the same.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Preguntada:

JJ
el 8 de Mzo. de 2018

Comentada:

el 15 de Nov. de 2022

Community Treasure Hunt

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

Start Hunting!

Translated by