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.
