natural spline, unmkpp gives 'wrong' answer

Question

Uwe Brauer 2017-2-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/325706-natural-spline-unmkpp-gives-wrong-answer

编辑： John D'Errico 2019-4-5

I have the following data

x=[-3 -2 1 4 6]
y=[2 0 3 1 4]

And I want to calculate the coefficient of the corresponding natural(!) Spline which is given by

 $s(x_j) = a_0 + a_1*(x-x_j) + a_2*(x-x_j)^2 + a_3*(x-x_j)^3$

I obtain

 |    |     a_3 |     a_2 |     a_1 |    a_0 | intv    |
 | s0 |  0.5044 |       0 | -2.5044 | 2.0000 | [-3,-2] |
 | s1 | -0.2831 |  1.5132 | -0.9912 |      0 | [-2,1]  |
 | s2 |  0.2217 | -1.0351 |  0.4430 | 3.0000 | [1,4]   |
 | s3 | -0.1601 |  0.9605 |  0.2193 | 1.0000 | [4,6]   |

However when I am using matlab

clear all

 x=[-3 -2 1 4 6]
 y=[2 0 3 1 4]
 pspl=spline(x,[0 y 0])% a natural spline 
 [breaks,coefs,l,k,d] = unmkpp(pspl)

I obtain

 |  2.0595 | -4.0595 |       0 | 2.0000 | 
 | -0.3796 |  2.1190 | -1.9405 |      0 |
 |  0.3003 | -1.2976 |  0.5238 | 3.0000 |
 | -0.5387 |  1.4048 |  0.8452 | 1.0000 |

I don't know how to interpret this, the only terms with coincide are that of the right column which must correspond to the terms $a_j$ but for the rest I don't know.

Any help is appreciated.

Thanks

Uwe Brauer

2 个评论
显示无隐藏无

Stephen23 2017-2-18

@Uwe Brauer: this time I formatted your code for you. In future you can do it yourself by selecting the code text and clicking the {} Code button.

Uwe Brauer 2017-2-18

sorry and thanks that you reformatted the code. I will keep this in mind!

请先登录，再进行评论。

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

Answer 1

John D'Errico 2017-2-18

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/325706-natural-spline-unmkpp-gives-wrong-answer#answer_255346

编辑：John D'Errico 2017-2-18

在 MATLAB Online 中打开

Well, it gives the right answer, to a different problem than the one that you think you posed. :)

So first of all, if you read the help for spline:

"However, if Y contains two more values than X has entries, then the first and last value in Y are used as the endslopes for the cubic spline."

That is NOT the definition of a NATURAL spline, although I can see how one might mistake that easily enough. A natural cubic spline has SECOND derivatives zero at the ends.

Next, the coefficients that you see here:

pspl.coefs
ans =
       2.0595      -4.0595            0            2
     -0.37963        2.119      -1.9405            0
      0.30026      -1.2976      0.52381            3
     -0.53869       1.4048      0.84524            1

are listed with the constant term LAST, so in decreasing order. The first element in each row is the coefficient of the cubic term (as you show it, with the knot in that interval subtracted off).

So, if we differentiate the spline ONCE, then evaluate it at the end points, we see:

fnval(fnder(pspl,1),[-3 6])
ans =
            0  -6.6613e-16

The first derivatives are indeed zero at the ends, although not the second derivatives.

fnval(fnder(pspl,2),[-3 6])
ans =
       -8.119      -3.6548

If your goal is to generate a natural cubic spline, I'll need to think for a second. :) A smoothing spline will suffice, if we crank down on tol.

pspl = spaps(x,y,0)
pspl = 
  struct with fields:
        form: 'B-'
       knots: [-3 -3 -3 -3 -2 1 4 6 6 6 6]
       coefs: [2 1.1652 -2.174 6.1053 -0.99415 2.5731 4]
      number: 7
       order: 4
         dim: 1

It produces a result in B-spline form, but s we see, it is indeed a natural cubic spline.

fnval(fnder(pspl,2),[-3 6])
ans =
  -2.6645e-15  -4.4409e-16

I'm pretty sure there is a direct way to generate a natural cubic spline in pp form, but my mind is drawing a blank at the moment. Drat, I almost forgot about csape.

pp = csape(x,y,'var')
pp = 
  struct with fields:
        form: 'pp'
      breaks: [-3 -2 1 4 6]
       coefs: [4×4 double]
      pieces: 4
       order: 4
         dim: 1
fnval(fnder(pp,2),[-3 6])
ans =
     0     0
fnval(fnder(pp,1),[-3 6])
ans =
      -2.5044       2.1404
pp.coefs
ans =
      0.50439            0      -2.5044            2
     -0.28314       1.5132     -0.99123            0
      0.22173      -1.0351      0.44298            3
     -0.16009      0.96053       0.2193            1

So a cubic spline in pp form, with natural end conditions. Of course, if you were to use my SLM toolbox, a spline with natural end conditions is one of the many options.

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

John D'Errico 2017-2-19

编辑：John D'Errico 2017-2-19

Natural splines are also not that valuable. Yes, they get written about a lot, because they are easy to build, and fairly easy to understand.

Splines are what I call a metaphorical model. The idea is that a cubic spline is a good model for the deflection of a thin flexible beam, that passes through a set of points. If one presumes that the deflection of a thin, flexible beam forms a good approximation for the shape of your curve, even though the physics behind the shape of that curve has NOTHING to do with the shape of a true beam, if that is true, then a spline serves as a metaphor. Such a metaphor allows you to predict the behavior of your process, to the extent that the metaphor is a good approximation to the physical process.

Some metaphors serve as better approximations than others, and the quality of that approximation allows you to trust the predictions of the metaphorical model.

So how good are natural splines as a model for a general curve? Sometimes they are great, other times not as good. The example I like to use when I taught about the use of splines as models is to fit the curves sin(x) and cos(x) over the same domains, that is [0,2*pi]. How well will a natural cubic spline do on these two curves? The problem is a flexible beam uses what are called natural boundary conditions, coming from the variational calculus solution for the shape of that beam. This implies that the second derivatives at the end points are zero. That works great to fit sin(x) over the domain [0,2*pi], since the second derivatives at the ends of the domain are indeed zero. The problem is when you try to fit cos(x) over the same domain, you see the curve fits more poorly.

In real life, most curves that you will fit will not truly have zero curvature at the ends. So a natural cubic spline is actually not that great of a choice for some random curve that will be fit. In general, you will probably get better performance form other boundary conditions, and the BCs that seems to make the most sense for some random spline model are often the not-a-knot end conditions.

The point is that while students are taught about natural cubic splines as a jumping off point, it is actually not that great of a choice for your general spline. To be honest, it is actually a good idea that a natural cubic fit is somewhat hard to obtain, since if it was easy, then novices who want to use a spline fit will gravitate to what they learned about early in their studies. As well, for the novice who knows absolutely nothing, they will be told by their compatriots (who also know little beyond what they learned from others) they will be told to use a natural spline fit.

While you may think this does not happen it really does. Such memes propagate all over numerical methods. For example, people are constantly wanting to use tools like the determinant to test for singularity of a matrix or they want to use inv to solve linear systems of equations. Both of these things are bad things in terms of numerical analysis, but because inv and det are do easy to use, people assume they must be appropriate.

A long argument, but as I said, I don't think it is really that good of a thing to make it easy to get a natural cubic spline interpolant. It is in there, if you know where to look. And if you understand enough about splines to also know why it may arguably be a bad thing, then you can find it without too much effort. (As I showed above, there were two distinct ways I came up with to produce a natural cubic interpolant.)

Michael Reshko 2019-4-4

Why does "order" property in your cubic splines shows 4? Shouldn't cubic be order 3?

John D'Errico 2019-4-5

编辑：John D'Errico 2019-4-5

在 MATLAB Online 中打开

Why? Whats in a name?

pp = spline(rand(1,5),rand(1,5))
pp = 
  struct with fields:
      form: 'pp'
    breaks: [0.16679 0.33499 0.46094 0.51324 0.66433]
     coefs: [4×4 double]
    pieces: 4
     order: 4
       dim: 1

Blame it on the designers of splines in MATLAB. See that the CUBIC spline created is shown to have order 4.

Really, a name can indicate anything you want it to indicate. Many people use the words order and degree to be synonyms. My choice was to be consistent with that chosen in MATLAB. I did not write spline, but I need to use it and its cousins.

https://en.wikipedia.org/wiki/Degree_of_a_polynomial

In there we even see it stated that degree of a univariate polynomial in x would essentially be the highest power of x. As well, it states that order was commonly used in the past as a synonym for degree.

However, here we see a statement about the order of a spline as:

https://en.wikipedia.org/wiki/Order_of_a_polynomial

"the order of a spline, either the degree+1 of the polynomials defining the spline or the number of knot points used to determine it."

And here:

https://en.wikipedia.org/wiki/Spline_(mathematics)

"If the n polynomials Pi each have degree at most m, then the spline is said to be of degree m (or is of order m+1)."

A name is just a name. But order and degree have been disambiguated, no longer twins.

请先登录，再进行评论。

natural spline, unmkpp gives 'wrong' answer

2 个评论
显示无隐藏无

采纳的回答

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

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

natural spline, unmkpp gives 'wrong' answer

2 个评论 显示 无隐藏 无

采纳的回答

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

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

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