numerics - Higher order periodic interpolation (curve fitting)

I have a list of points in 3D, and I want to get a smooth interpolation or curve fit (it is more for illustration) of these points such that the first and second derivatives at the start and end points agree. With

ListInterpolation[pts1, {0, 1}, InterpolationOrder -> 4, PeriodicInterpolation -> True]

even the first derivatives do not agree.

Answer

It does seem that the options PeriodicInterpolation -> True and Method -> "Spline" are incompatible, so I'll give a method for implementing a genuine cubic periodic spline for curves. First, let's talk about parametrizing the curve.

Eugene Lee, in this paper, introduced what is known as centripetal parametrization that can be used when one wants to interpolate across an arbitrary curve in $\mathbb R^n$. Here's a Mathematica implementation of his method:

parametrizeCurve[pts_List, a : (_?NumericQ) : 1/2] /; MatrixQ[pts, NumericQ] := 
           FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]] 

The default setting of the second parameter for parametrizeCurve[] gives the centripetal parametrization. Other popular settings include a == 0 (uniform parametrization) and a == 1 (chord length parametrization).

Now we turn to generating the derivatives needed for periodic cubic spline interpolation. This is done through the solution of an appropriate cyclic tridiagonal system, for which the functions LinearSolve[], SparseArray[], and Band[] come in handy:

periodicSplineSlopes[pts_?MatrixQ] := 
 Module[{n = Length[pts], dy, ha, xa, ya}, {xa, ya} = Transpose[pts];
        ha = {##, #1} & @@ Differences[xa];
        dy = ({##, #1} & @@ Differences[ya])/ha;
        dy = LinearSolve[SparseArray[{Band[{2, 1}] -> Drop[ha, 2], 

                                      Band[{1, 1}] -> ListConvolve[{2, 2}, ha], 
                                      Band[{1, 2}] -> Drop[ha, -2],
                                      {1, n - 1} -> ha[[2]], {n - 1, 1} -> ha[[-2]]}], 
                         3 MapThread[Dot[#1, Reverse[#2]] &,
                                     Partition[#, 2, 1] & /@ {ha, dy}]];
        Prepend[dy, Last[dy]]]

Using Sjoerd's example:

sc = Table[{Sin[t], Cos[t], Cos[t] Sin[t]}, {t, 0, 2 π, π/5}] // N;
tvals = parametrizeCurve[sc]

   {0, 0.102805, 0.196242, 0.303758, 0.397195, 0.5, 0.602805, 0.696242,
    0.803758, 0.897195, 1.}

cmps = Transpose[sc];
slopes = periodicSplineSlopes[Transpose[{tvals, #}]] & /@ cmps;
{f1, f2, f3} = MapThread[Interpolation[Transpose[{List /@ tvals, #1, #2}], 
                                       InterpolationOrder -> 3, Method -> "Hermite",
                                       PeriodicInterpolation -> True] &,
                         {cmps, slopes}];

Plot the space curve:

Show[ParametricPlot3D[{f1[u], f2[u], f3[u]}, {u, 0, 1}],
     Graphics3D[{AbsolutePointSize[6], Point[sc]}]]

Individually plotting the components and their respective derivatives verify that the interpolating functions are $C^2$, as expected of a cubic spline: