3
$\begingroup$

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]&

plot

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

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:

NIntegrate`GaussRuleData[n, MachinePrecision]                 (* n points & weights *)
NIntegrate`GaussBerntsenEspelidRuleData[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@NIntegrate`GaussRuleData[deg + 1, MachinePrecision];
   xj = Rescale[xj, {0, 1}, {-1, 1}];
   wj = 2 wj;
   fj = Developer`ToPackedArray[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]

Mathematica graphics


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]

Mathematica graphics

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

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}]

Mathematica graphics

$\endgroup$
1
$\begingroup$

(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

real and interpolated arclength derivative

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.