6
$\begingroup$

I use the same data for interpolation in Mathematica and MATLAB, but the result is different.

x={-1.00,-0.96,-0.65,0.10,0.40,1.00};
y={-1.0000,-0.1512,0.3860,0.4802,0.8838,1.0000};
Interpolation[{x,y}//Transpose,Method->"Spline"][-0.3]

result: -0.87332

x=[-1.00,-0.96,-0.65,0.10,0.40,1.00];
y=[-1.0000,-0.1512,0.3860,0.4802,0.8838,1.0000];
splinetx(x,y,-0.3)

result: -0.1957

I tried different InterpolationOrder but still different.

Is the "Spline" same as splinetx?

If not, is there a function in Wolfram like splinetx?


The splinetx function is available here or here.

MATLAB also has a built-in spline function which gives identical results to splinetx.


This is a comparison of the results given by Mathematica and MATLAB:

comparison

$\endgroup$
  • 1
    $\begingroup$ What methods does splinetx use? The short answer to your first question is no. And the short answer to your second question is yes. The methods used in MATLAB or gnu octave and other similar programs can be realized with WL functions. $\endgroup$ – CA Trevillian Mar 17 at 6:27
  • 3
    $\begingroup$ The source code of splinetx can be found here: mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/… $\endgroup$ – xzczd Mar 17 at 7:20
  • 2
    $\begingroup$ MATLAB also has a built-in spline function which also gives -0.1957. I think this is a good question and it is worth trying to understand the source of the differences. $\endgroup$ – Szabolcs Mar 17 at 9:51
  • $\begingroup$ Since we have the source code for the MATLAB version, it could no doubt be reproduced in Mathematica with some effort. It is not particularly interesting. But can we reproduce Mathematica's version? I think it is important to fully understand and to be able to reproduce what Mathematica does ... I don't see enough in the documentation to be able to easily do this. The difference must be coming from the conditions on the endpoints. $\endgroup$ – Szabolcs Mar 17 at 10:09
  • 3
    $\begingroup$ @AlexTrounev That much is clear. What I am looking for is not the "correct" solution (as there isn't any single one) but a precise understanding of what Interpolation does here (I guess the question is what is assumed for the derivatives at the boundaries) $\endgroup$ – Szabolcs Mar 17 at 12:28
8
$\begingroup$

Is the "Spline" same as splinetx?

  • No.

If not, is there a function in Wolfram like splinetx?

  • No.

That was a bit negative. However, it is not too difficult to apply the formulae in this answer and this answer to derive a routine that generates not-a-knot cubic splines (as was astutely observed by CA Trevillian and others in the comments.)

Of course, one can use SparseArray[] + LinearSolve[] to solve the underlying tridiagonal system, so I'll do that in the function below:

notAKnotSpline[pts_?MatrixQ] := Module[{dy, h, p1, p2, sl, s1, s2, tr},
    h = Differences[pts[[All, 1]]]; dy = Differences[pts[[All, 2]]]/h;
    s1 = Total[Take[h, 2]]; s2 = Total[Take[h, -2]];
    p1 = ({3, 2}.Take[h, 2] h[[2]] dy[[1]] + h[[1]]^2 dy[[2]])/s1;
    p2 = (h[[-1]]^2 dy[[-2]] + {2, 3}.Take[h, -2] h[[-2]] dy[[-1]])/s2;
    tr = SparseArray[{Band[{2, 1}] -> Append[Rest[h], s2], 
                      Band[{1, 1}] -> Join[{h[[2]]}, ListCorrelate[{2, 2}, h], {h[[-2]]}], 
                      Band[{1, 2}] -> Prepend[Most[h], s1]}];
    sl = LinearSolve[tr, Join[{p1}, 
                              3 Total[Partition[dy, 2, 1]
                                      Reverse[Partition[h, 2, 1], 2], {2}],
                              {p2}]];
    Interpolation[MapThread[{{#1[[1]]}, #1[[2]], #2} &, {pts, sl}], 
                  InterpolationOrder -> 3, Method -> "Hermite"]]

Try it out on the points in the OP:

pts = {{-1., -1.}, {-0.96, -0.1512}, {-0.65, 0.386},
       {0.1, 0.4802}, {0.4, 0.8838}, {1., 1.}};
spl = notAKnotSpline[pts];

spl[-0.3]
   -0.195695

Plot[spl[x], {x, -1, 1},
     Epilog -> {Directive[AbsolutePointSize[6], ColorData[97, 4]], Point[pts]}]

points with not-a-knot spline interpolant

Demonstrate the $C^2$ property of the cubic spline:

Plot[{spl[x], spl'[x], spl''[x]}, {x, -1, 1}, PlotRange -> {-10, 30}]

not-a-knot spline with first two derivatives


Szabolcs's desire to reproduce the results of Method -> "Spline" is a bit more difficult, because the exact formulae being used are not disclosed publicly. That being said, I was able to reverse-engineer and reproduce it some time ago, so go look at that answer if you want more details.

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ (I actually have a more general spline interpolation function than this, but the code is in a hard disk many kilometers away, so releasing that will have to wait.) $\endgroup$ – J. M.'s technical difficulties Mar 18 at 1:59
  • $\begingroup$ Thanks!It's very helpful. I've read your answers and write a natural cubic spline function. It seems that the Method -> "Spline" in Interpolation is not a natural cubic spline. So what the conditions are actually used for it? Why Mathematica use an uncommon method? $\endgroup$ – srtie Mar 19 at 16:45
  • $\begingroup$ I wrote a natural spline routine in another answer, and the formulae used by Method -> "Spline" are in the answer I linked to in the last paragraph. Why these choices were made, I cannot say. $\endgroup$ – J. M.'s technical difficulties Mar 19 at 21:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.