# FunctionInterpolation over an open interval

I'm trying to obtain an arc-length parametrization of a spline curve as per https://mathematica.stackexchange.com/a/8456. To do so, I need to calculate partial derivatives of such curve and then perform a FunctionInterpolation. The spline curve is continuous over its whole domain, but not differentiable at its knots (left-side derivatives are different from the right-side ones). Therefore, I'd need to create several FunctionInterpolations, each defined on an open interval between two consecutive knots.

Is that possible? And if so, how? If not, how can I overcome this limitation?

Here is the actual data I'm working on: curve = BSplineCurve[{{198.2063059205538,402.121623269958},{191.2621031932776,243.2810494404048},{350.2070352491573,246.3653323645185},{307.7702781407043,406.7480444184375},{624.1169877141547,382.0737810255268},{367.1817351705904,154.607913753287},{408.0753377407309,28.15231386462074},{192.0336877673071,89.837972346897}}, SplineDegree -> 3, SplineKnots -> {0,0,0,0,1.,2.,3.,4.,5,5,5,5}, SplineWeights -> {1,1,1,1,1,1,1,1}];

cp = curve[[1]];
dg = curve[[2]] // Values;
kn = curve[[3]] // Values;
wg = curve[[4]] // Values;

x[t_]:=Total[Table[BSplineBasis[{dg,kn},j,t],{j,0,Length[kn]-dg-2}]*wg*cp[[All,1]]]/Table[BSplineBasis[{dg,kn},j,t],{j,0,Length[kn]-dg-2}].wg;
y[t_]:=Total[Table[BSplineBasis[{dg,kn},j,t],{j,0,Length[kn]-dg-2}]*wg*cp[[All,2]]]/Table[BSplineBasis[{dg,kn},j,t],{j,0,Length[kn]-dg-2}].wg;
f[t_]:=Sqrt[(x'[t])^2+(y'[t])^2];


You can see that the interpolating function has spikes (that I'm trying to avoid) at the knot values:

FunctionInterpolation[f[t],{t,Min[kn],Max[kn]}]//Plot[#[g]-f[g],{g,Min[kn],Max[kn]},PlotRange->Full]&


• @MichaelE2 Sorry, typo, fixed (zzz was supposed to be f) – Marco Feb 18 '18 at 19:39

I've been waiting for a good problem to use barycentric interpolation at the Gauss-Legendre points. It has excellent properties for analytic functions. Here we have a piecewise analytic function, so it should be good once we have divided the interval up. For some reason, I was having trouble getting a good result until I hit on the idea of using NIntegrate to subdivide the interval.

The formula for barycentric interpolation is $$f_{IF}={\sum_j \displaystyle {\lambda_j \over x-x_j} \, f(x_j) \over \sum_j \displaystyle {\lambda_j \over x-x_j}}$$ where the weights $\lambda_j$ are given by $$\lambda_j = \prod_{k \ne j} {1 \over x_j-x_k}\,.$$ When the nodes $x_j$ are the zeros of the Legendre polynomial $P_{n+1}$, the weights are given by $$\lambda_j = (-1)^j \sqrt{\left(1-x_j^2\right)\,w_j}$$ where the $w_j$ are the Gauss rule integration weights. The nodes $x_j$ and weights $w_j$ are available in Mathematica in both of the following:

NIntegrateGaussRuleData[n, MachinePrecision]                 (* n points & weights *)
NIntegrateGaussBerntsenEspelidRuleData[n, MachinePrecision]  (* 2n+1 points & weights *)


Here then functions for constructing the Gauss-Legendre barycentric interpolation of degree deg of a function f over an interval {a, b}:

ClearAll[baryInterp, gaussLegendreInterpolation];
(* construct barycentric interpolation at the Gauss-Legendre nodes *)
gaussLegendreInterpolation[deg_Integer?Positive, f_, {a_, b_}] :=
Module[{xj, fj, lj, wj},
{xj, wj} = Most@NIntegrateGaussRuleData[deg + 1, MachinePrecision];
xj = Rescale[xj, {0, 1}, {-1, 1}];
wj = 2 wj;
fj = DeveloperToPackedArray[f /@ Rescale[xj, {-1, 1}, {a, b}],
Real];
lj = (-1)^Range[0, deg] Sqrt[(1 - xj^2) wj];
baryInterp[xj, fj, lj, {a, b}]
];

(* generic barycentric interpolation *)
baryInterp[xj_List, fj_List, lj_, {a_, b_}][x_?NumericQ] :=
Module[{xx, j, c},
xx = Rescale[x, {a, b}, {-1, 1}];
j = First@Nearest[xj -> "Index", xx];
If[xx == xj[[j]],
fj[[j]],
c = lj/(xx - xj);
c.fj/Total[c]
]
];


Let's apply it to the OP's problem. The Gauss-Legendre interpolation converges quickly, like Gaussian quadrature. For a function that is analytic on an interval, it pays to pick a moderately high degree.

nNodes = 21; (* must be odd *)
knots = Cases[f[t], BSplineBasis[{_, u_}, __] :> u, Infinity] // N //
Flatten // DeleteDuplicates
(*  {0., 1., 2., 3., 4., 5.}  *)

Clear[df, t];
df[t0_?NumericQ] := f'[t0];
NIntegrate[
df[t], Evaluate@Flatten@{t, knots},
Method -> {"GaussBerntsenEspelidRule", "Points" -> (nNodes + 1)/2},
PrecisionGoal -> 12,
IntegrationMonitor :> ((intervals = #@ "Boundaries" & /@ #) &)
]; // AbsoluteTiming
(*  {1.36377, Null}  *)

Flatten[intervals, 1]  (* intervals used by NIntegrate to approx. int(f) *)
(*
{{1.5, 2.}, {0.5, 1.}, {0., 0.5}, {2., 2.5}, {2.75, 3.}, {3.5, 4.},
{4.25, 4.5}, {2.5, 2.75}, {3., 3.25}, {4., 4.25}, {3.25, 3.5},
{4.5, 4.75}, {1., 1.5}, {4.75, 5.}}
*)

(* construct a Piecwise function of interpolations *)
ClearAll[fIF];
fIF[t_?NumericQ] = Piecewise[
{gaussLegendreInterpolation[nNodes, f, #][t],
First@# <= t < Last@#} & /@ Flatten[intervals, 1]
];

(* Plot the relative error. Quite often it is accurate to nearly machine precision *)
Plot[{RealExponent@$MachineEpsilon, (fIF[t] - f[t])/f[t] // RealExponent}, {t, 0, 5}, PlotPoints -> 200, MaxRecursion -> 2, PlotRange -> {-18, 0}, Filling -> {1 -> Bottom}, Frame -> True]  Update Here's a way to convert the above to a built-in form of InterpolatingFunction. It is based on chebInterpolation. ClearAll[chebdeval, chebInterpolation, chebSeries]; (* https://mathematica.stackexchange.com/questions/13645/small-issue-with-chebyshev-derivative-approximation/13681#13681 *) chebdeval[coeff_List, x_] := Module[{m = Length[coeff], d, e, f, j, u, v, w}, w = v = 0; f = e = 0; Do[d = 2 (x e + v) - f; f = e; e = d; u = 2 x v - w + coeff[[j]]; w = v; v = u;, {j, m, 2, -1}]; {x v - w + First[coeff], x e + v - f}]; (* Constructs a piecewise InterpolatingFunction, * whose interpolating units are Chebyshev series c1 *) (* data0 = {{x0,x1},c1},..} *) chebInterpolation::usage = "chebInterpolation[data:{{{x0,x1},c1},..}] constructs a piecewise \ InterpolatingFunction, whose interpolating units are Chebyshev series \ c1 on the interval {x0, x1}"; chebInterpolation[data0 : {{{_, _}, _List} ..}] := Module[{data = Sort@data0, domain1, coeffs1, domain, grid, ngrid, coeffs, order, y0, yp0}, domain1 = data[[1, 1]]; coeffs1 = data[[1, 2]]; domain = List @@ Interval @@ data[[All, 1]]; {y0, yp0} = (Last@domain - First@domain)/2 chebdeval[coeffs1, First@domain]; grid = Union @@ data[[All, 1]]; ngrid = Length@grid; coeffs = PadRight[data[[All, 2]], Automatic, 0.]; order = Length[coeffs[[1]]] - 1; InterpolatingFunction[ domain, {5, 1, order, {ngrid}, {order}, 0, 0, 0, 0, Automatic, {}, {}, False}, {grid}, (* if[[3]] *) {{y0, yp0}} ~Join~ coeffs, (* if[[4]] *) {{{{1}}~Join~Partition[Range@ngrid, 2, 1] ~Join~ {{ngrid - 1, ngrid}}, {Automatic} ~Join~ ConstantArray[ChebyshevT, ngrid]}} ] /; Length[domain] == 1 ]; chebSeries[f_, {a_, b_}, n_] := Module[{x, y, c}, x = Rescale[Sin[Pi/2*Range[N@n, -N@n, -2]/n], {-1, 1}, {a, b}]; y = f /@ x; (* function values at Chebyshev points *) c = Sqrt[2/n] FourierDCT[y, 1]; (* get coeffs from values *) c[[{1, -1}]] /= 2; (* adjust first & last coeffs *) c ];  OP's example: foo = {#, chebSeries[gaussLegendreInterpolation[nNodes, f, #], #, nNodes]} & /@ Flatten[intervals, 1] // chebInterpolation  The error is nearly the same: Plot[{RealExponent@$MachineEpsilon, (foo[t] - f[t])/f[t] //
RealExponent},
{t, 0, 5},
PlotPoints -> 200, MaxRecursion -> 2, PlotRange -> {-18, 0},
Filling -> {1 -> Bottom}, Frame -> True]


• There is a built-in, but undocumented routine for barycentric interpolation. In your notation, one would call it as f = StatisticsLibraryBarycentricInterpolation[xj, fj, Weights -> lj]; f[x] – J. M. will be back soon Mar 19 '18 at 4:11
• @J.M. Thanks. And f[x, 1] yields the first derivative. There's good, general stuff in StatisticsLibrary that should accessible to users at the System  level. – Michael E2 Mar 19 '18 at 11:57

Perhaps the easiest way is to use NDSolve, with an event to restart integration at the knots:

Clear[ff, t];
df[t0_?NumericQ] := f'[t0];
fIF = NDSolveValue[
{ff'[t] == df[t], ff[0] == f[0],
WhenEvent[Mod[t, 1] == 0, "RestartIntegration"]},
ff, {t, 0, 5},
(*Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 15},*)
(*PrecisionGoal -> 10,*)
InterpolationOrder -> All];


The optional options result in a higher quality interpolation, but it takes over a minute.

Plot of the relative error:

Plot[(fIF[t] - f[t])/f[t] // RealExponent, {t, 0, 5},
PlotPoints -> 100, MaxRecursion -> 2, PlotRange -> {-20, 0}]


(This is too long for a comment.)

Note that all your weights are set to 1, and so you do not need to construct the full NURBS expression (recall that the B-spline basis functions form a partition of unity); thus, the polynomial case will suffice:

cp = {{198.2063059205538, 402.121623269958}, {191.2621031932776, 243.2810494404048},
{350.2070352491573, 246.3653323645185}, {307.7702781407043, 406.7480444184375},
{624.1169877141547, 382.0737810255268}, {367.1817351705904, 154.607913753287},
{408.0753377407309, 28.15231386462074}, {192.0336877673071, 89.837972346897}};
d = 3;
knots = {0, 0, 0, 0, 1, 2, 3, 4, 5, 5, 5, 5};
wts = {1, 1, 1, 1, 1, 1, 1, 1};

x[t_] = Table[BSplineBasis[{d, knots}, j, t], {j, 0, Length[knots] - d - 2}].cp[[All, 1]];
y[t_] = Table[BSplineBasis[{d, knots}, j, t], {j, 0, Length[knots] - d - 2}].cp[[All, 2]];


Having gotten that out of the way, here is another proposal: let Plot[] do the sampling for you, and then perform a Hermite interpolation of the arclength derivative. Here's one way to go about it:

f[t_] := Sqrt[x'[t]^2 + y'[t]^2]

(* sampling; adjust PlotPoints or MaxRecursion as needed *)
pts = First[Cases[Normal[Plot[f[t], {t, 0, 5}]], Line[l_] :> l, ∞]];
(* ensure endpoints are present *)
{xmin, xmax} = MinMax[knots];
pts = DeleteDuplicatesBy[{{xmin, f[xmin]}} ~Join~ pts ~Join~ {{xmax, f[xmax]}}, First]

(* build Hermite interpolant *)
ff = Interpolation[{{#1}, #2, f'[#1]} & @@@ pts, Method -> "Hermite"];


Compare:

{Plot[{f[t], ff[t]}, {t, 0, 5}, PlotStyle -> {AbsoluteThickness[4], AbsoluteThickness[1]}],
Plot[f[t] - ff[t], {t, 0, 5}, PlotRange -> All]} // GraphicsRow