Arbitrary precision spline interpolation

The OP linked to an answer of mine for interpolating over general point sets; for constructing a single interpolating function, a slight modification of my procedure is needed. (In particular, you don't need centripetal or chord-length parametrization in this case.)

Let's start with some data:

data = {{0, 0}, {1/10, 3/10}, {1/2, 3/5}, {1, -1/5}, {2, 3}, {3, -6/5}};

To review the problems with using the built-in Interpolation[], let's try building an interpolating quartic spline from data:

m = 4; (* spline order *)
sp1 = Interpolation[data, InterpolationOrder -> m, Method -> "Spline"];

Now, evaluate:

sp1[3/2]
   0.26753591659167364

Herein lies the problem: even though the contents of data and the argument fed to sp1 are exact, the output is a machine precision number. Yeesh!

So, we go back to basics and build the interpolating spline ourselves from BSplineBasis[]. First, we generate the knots. The clamped condition is the most appropriate for interpolation, so we generate the appropriate knot sequence like so:

n = Length[data];
{xa, ya} = Transpose[data]
knots = Join[ConstantArray[xa[[1]], m + 1], 
             If[m + 2 <= n, MovingAverage[ArrayPad[xa, -1], m], {}], 
             ConstantArray[xa[[-1]], m + 1]];

To generate the control points for the B-spline, we construct the relevant linear system and solve it with LinearSolve[]:

cp = LinearSolve[Outer[BSplineBasis[{m, knots}, #2, #1] &,
                       xa, Range[0, Length[data] - 1]], ya];

Note that all the entries of cp are exact.

At this juncture, you might think that we can use BSplineFunction[] on the control points and the knots, but alas, the function also only evaluates at machine precision. We are left with no choice but to use BSplineBasis[] once more:

sp2[x_] = cp.Table[BSplineBasis[{m, knots}, j - 1, x], {j, n}];

Try it out:

sp2[3/2]
   25251843/94386740

N[%]
   0.2675359165916738

As a verification that we were able to reproduce the control points used internally by Interpolation[], let's do a comparison:

sp1bs = Cases[sp1, _BSplineFunction, ∞][[1]];
sp1bs["ControlPoints"]
   {0., 0.777830630128766, 1.3470257035045392, -6.670028631670085, 13.100332135107115, -1.2}

N[cp]
   {0., 0.777830630128766, 1.347025703504539, -6.6700286316700845, 13.100332135107113, -1.2}

Certainly, I can't end this post without at least one picture, so:

Plot[{sp1[x], sp2[x]}, {x, 0, 3}]

Plot[sp1[x] - sp2[x], {x, 0, 3}, PlotRange -> All]

I leave the task of bundling everything in a routine as an exercise. If you want to learn more, Piegl and Tiller's The NURBS Book is the canonical reference.

Here is my implementation of quadric spline arbitrary precision interpolation which is defined exactly as in Interpolation with options Method->"Spline", InterpolationOrder->2.

Theoretical background

Quadric spline interpolation for n datapoints is defined as a Piecewise function which consists of n-2 parabolas which are spliced together in the middles of successive datapoints (with exceptions for first and last intervals) by two conditions: successive parabolas at these middle points must be equal and their first derivatives must be equal too. So we have 2(n-3) such conditions.

Let us designate jth datapoint as {x[j],y[j]} and define jth parabola as

s[j_, xx_] = y[j + 1] + b[j] (xx - x[j + 1]) + c[j] (xx - x[j + 1])^2; 

This definition automatically makes jth parabola equal to y[j + 1] at x[j + 1]. We have to add also 2 boundary conditions:

s[1, x[1]] == y[1]
s[n - 2, x[n]] == y[n]

And 2(n-3) splicing conditions:

Table[{
  s[j, (x[j + 1] + x[j + 2])/2] == s[j + 1, (x[j + 1] + x[j + 2])/2], 
  ds[j, (x[j + 1] + x[j + 2])/2] == ds[j + 1, (x[j + 1] + x[j + 2])/2]
  }, {j, 1, n - 3}]

where ds is first derivative:

ds[j_, xx_] = D[s[j, xx], xx];

Now we combine everything together and convert into matrix form:

eqs = Flatten[{
    s[1, x[1]] == y[1],
    Table[{
      s[j, (x[j + 1] + x[j + 2])/2] == s[j + 1, (x[j + 1] + x[j + 2])/2],
      ds[j, (x[j + 1] + x[j + 2])/2] == ds[j + 1, (x[j + 1] + x[j + 2])/2]
      }, {j, 1, n - 3}],
    s[n - 2, x[n]] == y[n]}];
arrs = Simplify@
   Normal@CoefficientArrays[eqs, Flatten@Array[{b[#], c[#]} &, n - 2]];
MatrixForm /@ (4 arrs)

Looking at the matrices it is easy to see that they have periodic structure and contain only elements of the forms x[j+1]-x[j] and y[j+1]-y[j] with some numerical coefficients. This periodic structures can be expressed as SparceArrays.

Implementation

Assuming that data contains an array of {x[j],y[j]}, we define

Δx[i_] := Subtract @@ data[[{i + 1, i}, 1]]; 
Δy[i_] := Subtract @@ data[[{i + 1, i}, 2]];

Two matrices in the arrs can be expressed as follows:

bB = SparseArray[{1 -> -4 Δy[1], -1 -> 4 Δy[n - 1], i_ /; OddQ[i] :> 0, i_ /; EvenQ[i] :> 4 Δy[i/2 + 1]}, 2 (n - 2)];
mM = SparseArray[
   Join[{{1, 1} -> -4 Δx[1], {1, 2} -> 4 Δx[1]^2, {-1, -1} -> 4 Δx[n - 1]^2, {-1, -2} -> 4 Δx[n - 1]}, 
    Table[Band[{2 i, 2 i - 1}] -> {{2 Δx[i + 1], Δx[i + 1]^2, 2 Δx[i + 1], -Δx[i + 1]^2}, {4, 4 Δx[i + 1], -4, 4 Δx[i + 1]}}, {i, 1, n - 3}]], {2 (n - 2), 2 (n - 2)}];

Now for obtaining coefficients b[i], c[i] we can use LinearSolve:

bcM = LinearSolve[mM, bB];

Grouping coefficients with identical indexes together:

bcM = Partition[bcM, 2];

The intervals at which the parabolas are defined:

intervList = 
  Partition[
   Join[{data[[1, 1]]}, 
    MovingAverage[data[[2 ;; -2, 1]], 2], {data[[-1, 1]]}], 2, 1];

Now we compile all the spline data in one object:

splineData = Table[{data[[i + 1]], bcM[[i]], intervList[[i]]}, {i, n - 2}];

splineData contains everything that is needed to construct the spline. It contains elements of the form:

{{x[i + 1], y[i + 1]}, {b[i], c[i]}, {xmin, xmax}}

It is very handy to use splineData for exact numerical integration (see here).

To produce the explicit piecewise function we define the constructor (HornerForm is already applied and gives 27% speedup):

makeSpline[splineData_List, x_Symbol] := 
 Piecewise[
  Append[Table[{d[[1, 2]] - d[[1, 1]] d[[2, 1]] + 
      d[[1, 1]]^2 d[[2, 2]] + x (d[[2, 1]] + x d[[2, 2]] - 
         2 d[[1, 1]] d[[2, 2]]), #1 <= x <= #2 & @@ d[[3]]}, {d, splineData}],
         {Indeterminate, True}]]

It can be used as follows:

spline[\[FormalX]_] = makeSpline[splineData, \[FormalX]];

Compiling the code for creation of splineData into one function:

ClearAll[toSplineData, makeSpline];
Options[toSplineData] = {Method -> Automatic};
toSplineData[data_, OptionsPattern[]] /; MatrixQ[data, NumberQ] := 
  Module[{Δx, Δy, bB, mM, bcM, intervList,
     n = Length[data], dataS = Sort[data]},
   Δx[i_] := 
    Subtract @@ dataS[[{i + 1, i}, 1]]; Δy[i_] := 
    Subtract @@ dataS[[{i + 1, i}, 2]];
   bB = SparseArray[{1 -> -4 Δy[1], -1 -> 
       4 Δy[n - 1], i_ /; OddQ[i] :> 0, 
      i_ /; EvenQ[i] :> 4 Δy[i/2 + 1]}, 2 (n - 2)];
   mM = SparseArray[
     Join[{{1, 1} -> -4 Δx[1], {1, 2} -> 
        4 Δx[1]^2, {-1, -1} -> 
        4 Δx[n - 1]^2, {-1, -2} -> 
        4 Δx[n - 1]}, 
      Table[Band[{2 i, 
          2 i - 1}] -> {{2 Δx[i + 1], Δx[
           i + 1]^2, 
          2 Δx[i + 1], -Δx[i + 1]^2}, {4, 
          4 Δx[i + 1], -4, 
          4 Δx[i + 1]}}, {i, 1, n - 3}]], {2 (n - 2), 
      2 (n - 2)}];
   bcM = LinearSolve[mM, bB, Method -> OptionValue[Method]];
   bcM = Partition[bcM, 2];
   intervList = 
    Partition[
     Join[{dataS[[1, 1]]}, 
      MovingAverage[dataS[[2 ;; -2, 1]], 2], {dataS[[-1, 1]]}], 2, 1];
   Table[{dataS[[i + 1]], bcM[[i]], intervList[[i]]}, {i, n - 2}]];

Now to construct the spline from data we can evaluate

spline[\[FormalX]_] = makeSpline[toSplineData[data], \[FormalX]];

Any suggestions and improvements are welcome!

What is so special about quadric spline interpolation?

Quadric spline interpolation with splicing points in the middle of successive data points is locally invariant and weakly depends on boundary condition. (But the same is not true for quadric spline interpolation with splicings in data points.) It gives much lesser artifacts than "usual" cubic spline interpolation. Moreover, it seems that splines of higher degree do not give any advantages over such quadric splines. As opposed to cubic spline, data points are not special points at which polynomials are spliced, and it makes much easier exact integration in these points.

The following is a shameful plug of J. M.'s answer that you already linked to show how to get the parametrization by using BSplineBasis[] with arbitrary precision and packed into a function:

splineInterp[data_, order_, prec_] := 
  Module[{parametrizeCurve, tvals, bas, ctrlpts, knots},
    parametrizeCurve[pts_ /; MatrixQ[pts, NumericQ], a : (_?NumericQ) : 1/2] := 
     FoldList[
       Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]];
    tvals = parametrizeCurve[data];
    (*knots for interpolating B-spline*)
    knots = 
      Join[
        ConstantArray[0, order + 1],
        MovingAverage[ArrayPad[tvals, -1], order],
        ConstantArray[1, order + 1]];
    (*basis function matrix*)
    bas = 
      Table[
        BSplineBasis[{order, knots}, j - 1, N[tvals[[i]], prec]], 
        {i, Length[data]}, {j, Length[data]}];
    ctrlpts = LinearSolve[bas, data];
    Return[
      Sum[
        ctrlpts[[i + 1]] BSplineBasis[{order, knots}, i, #], 
        {i, 0, Length[testData] - 1}] &]
  ]

Usage

SeedRandom[42];
testData = RandomReal[{0, 1}, {10, 2}, WorkingPrecision -> prec];
f = splineInterp[testData, 5, 50];
f[1/10]
ParametricPlot[
  f[t], {t, 0, 1}, 
  Axes -> None, 
  Epilog -> 
    {Directive[Green, PointSize[Large]], Point[testData]}, Frame -> True
]


(*
  {0.06025513038208479326777055464103084790562316, \
   0.5927955887343273319037301352154452480202510}
*)

Please refer to J. M.'s answer for the details.

---

Explanation

The centripetal distribution:

The function parametrizeCurve complete this algorithm

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