Different results using spline interpolation in Wolfram and MATLAB

Question

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:

Answer 1

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]}]

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

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

---

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.