3

Consider

data = {{0, 0}, {1, 0}, {2, 0}, {3, 1}, {4, 0}, {5, 0}, {6, 0}};
f = Interpolation[data, InterpolationOrder -> 3, Method -> "Spline"];

Then the second derivative is continuous and piecewise linear as expected but at the end points I get second derivative f''[0] = f''[6] = -2.57143. How does Interpolation decide on that at the end points?

Also, how can I get a natural cubic spline with second derivative = 0 at each end point?

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Ted Ersek
  • 7,124
  • 19
  • 41

2 Answers2

4

It is interesting to compare and spot the differences between the routine in this answer and the routine given below:

naturalSpline[pts_?MatrixQ] := Module[{dy, h, sl, tr},
       h = Differences[pts[[All, 1]]]; dy = Differences[pts[[All, 2]]]/h;
       tr = SparseArray[{Band[{2, 1}] -> Append[Rest[h], 1], 
                         Band[{1, 1}] -> Join[{2}, ListCorrelate[{2, 2}, h], {2}], 
                         Band[{1, 2}] -> Prepend[Most[h], 1]}];
       sl = LinearSolve[tr, Join[{3 dy[[1]]}, 
                                 3 Total[Partition[dy, 2, 1]
                                         Reverse[Partition[h, 2, 1], 2], {2}],
                                 {3 dy[[-1]]}]];
       Interpolation[MapThread[{{#1[[1]]}, #1[[2]], #2} &, {pts, sl}], 
                     InterpolationOrder -> 3, Method -> "Hermite"]]

Test on the OP's data:

data = {{0, 0}, {1, 0}, {2, 0}, {3, 1}, {4, 0}, {5, 0}, {6, 0}};
spl = naturalSpline[data];

{spl''[0], spl''[6]}
   {0, 0}

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

natural cubic spline interpolant

Verify $C^2$ continuity:

Plot[{spl[x], spl'[x], spl''[x]}, {x, 0, 6}, PlotRange -> All]

natural cubic spline and first two derivatives

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
2

I augmented your data a bit to force saddle points at the ends

data = {{0, 0}, {1, 0}, {2, 0}, {3, 1}, {4, 0}, {5, 0}, {6, 0}};
f = Interpolation[data, InterpolationOrder -> 3, Method -> "Spline"];
g = Interpolation[
   Join[{(# // First) - {1, 0}}, #, {(# // Last) + {1, 0}}] &@data, 
   Method -> "Spline", InterpolationOrder -> 3];

enter image description here

Now g''[0]==g''[6]==0

user2757771
  • 829
  • 4
  • 9
  • Will that give the desired result for any data? – Ted Ersek Feb 07 '20 at 01:33
  • This should be a bit more general data = Sort@RandomReal[{0, 10}, {10, 2}]; data' = Join[{2*data[[1]] - data[[2]]}, data, {2*data[[-1]] - data[[-2]]}]; f = Interpolation[data']; Plot[f[x], Prepend[MinMax[data[[All, 1]]], x], PlotRange -> Full] – user2757771 Feb 08 '20 at 21:37