Extract function from a parametric interpolating function

Question

I've read the answers of

• How to make use of an interpolating function outside of Mathematica?,

• Interpreting the interpolating function and saving data to plot with external program,

• Extracting the function from InterpolatingFunction object.

I have not managed to extrapolate them to my problem. The difference is that in my case, the interpolating function is of Output dimension 2: it's a parametric function (several values for one $x$).

This is an example code:

pts = {{0.504`, 2.79`}, {0.519`, 2.7773}, {0.5349`, 
   2.7642`}, {0.5504`, 2.7515`}, {0.5666`, 2.7398}, {0.5858`, 
   2.7341`}, {0.5914`, 2.8566}, {0.5917`, 2.8364}, {0.5918`, 
   2.8766}, {0.5924`, 2.8164}, {0.5933`, 2.7963`}, {0.5935`, 
   2.7525`}, {0.5935`, 2.8966`}, {0.5940, 2.7763`}, {0.5974`, 2.9163`}}
path = First@FindCurvePath@Standardize@pts;
curve = Interpolation[#] &@MapIndexed[{#2[[1]], #1} &, pts[[path]]];
Show[ParametricPlot[curve[t], {t, 1, Length[path]}, 
  AspectRatio -> 3/4, PlotStyle -> Red], ListPlot[pts]]

In the end, I want to export the function curve in Asymptote, which uses a syntax similar to C. For this, I plan to manipulate strings to replace the Piecewise syntax into if(x<3) {return x*.... ;}. Although it is a bit tedious, it should work. My problem is accessing the Interpolatingfunction curve, as explained above.

Concretely, I would like to have the piecewise polynomials expressions and their boundaries. Typically: pos={{0.5,0.52},{0.52,0.54}...} and expr={2.3+0.1*x-2.6*x^2, ...}. The reason why I can't use the answers of the questions above is that InterpolationPolynomial does not apply to a list of points, but to a indexed list of points (MapIndexed[{#2[[1]], #1} &, pts[[path]]]).

Answer 1

The way curve is built is by using a parameter $t$ corresponding to indexes in list: $(t_i,(x_i,y_i))_{i=1,\dots,n}$ where $t_i=i$ and $x_i=x(t_i)$, $y_i=y(t_i)$.

Hence, curve[i] give the coordinates of point $i$.

What I did is reconstruct two third order polynomials (func) on each interval $[i,i+1]$ and identify their 4 coefficents by solving the system of equation consisting in equating func and curve in i and i+1, as well as the derivatives. This has to be done twice on each interval, once for $x(t)$, once for $y(t)$. Note that I forced the method in Interpolation to be Spline so that curve is not made of Hermite polynomials.

pts = {{0.504`, 2.79`}, {0.519`, 2.7773}, {0.5349`, 
   2.7642`}, {0.5504`, 2.7515`}, {0.5666`, 2.7398}, {0.5858`, 
   2.7341`}, {0.5914`, 2.8566}, {0.5917`, 2.8364}, {0.5918`, 
   2.8766}, {0.5924`, 2.8164}, {0.5933`, 2.7963`}, {0.5935`, 
   2.7525`}, {0.5935`, 2.8966`}, {0.5940, 2.7763`}, {0.5974`, 2.9163`}}
path = First@FindCurvePath@Standardize@pts;
curve = Interpolation[#, Method -> "Spline"] &@
   MapIndexed[{#2[[1]], #1} &, pts[[path]]];
Show[ParametricPlot[curve[t], {t, 1, Length[path]}, 
  AspectRatio -> 3/4, PlotStyle -> Red], ListPlot[pts]]

func[t_] = a0 + a1*t + a2*t^2 + a3*t^3
output = Table[
   Block[{}, 
    eqns = Table[{func[i] == curve[i][[j]], 
       func'[i] == curve'[i][[j]], func[i + 1] == curve[i + 1][[j]], 
       func'[i + 1] == curve'[i + 1][[j]]}, {j, 1, 2}];
    solX[x_] = func[x] /. Solve[eqns[[1]], {a0, a1, a2, a3}][[1]];
    solY[x_] = func[x] /. Solve[eqns[[2]], {a0, a1, a2, a3}][[1]];
    {{curve[i][[1]], curve[i + 1][[1]]}, {solX[t], solY[t]}}], {i, 1, 
    Length[path]}];

Show@{ParametricPlot[curve[t], {t, 1, Length[path]}, 
   AspectRatio -> 3/4, PlotStyle -> {Red, Thickness[0.005]}], 
  Table[ParametricPlot[{output[[i, 2, 1]], output[[i, 2, 2]]}, {t, i, 
     i + 1}], {i, 1, Length[path]}]}

The two plots perfetctly match (in fact, the difference is typically of $10^{-15}$).

Conclusion From what I understand, there is no way to directly extract information from InterpolatingFunction. Instead, it is possible to reconstruct a piecewise solution by using the InterpolatingFunction and its derivatives.

Answer 2

From the documentation of InterpolatingFunction, section "properties and relations":

InterpolatingFunction does a Piecewise polynomial interpolation

and then comes the code to get those polynomials.

After running the first two lines of your code I do (just for the x-coordinates):

pts = First /@ pts[[path]];
pts = {#, pts[[#]]} & /@ Range@Length@pts;
xfunc = Interpolation[pts];

and get the InterpolatingFunction you want to describe with polynomials.

Then to get the polynomial description:

pf[t_] := 
 Piecewise@
   (({InterpolatingPolynomial[pts[[# ;; # + 3]], t], # <= t < (# + 1)} & /@
     Range[Length[pts] - 4])~Join~
       {{InterpolatingPolynomial[pts[[-4 ;; -1]], t], t >= Length[pts] - 3}})

pf[t] returns the Piecewise form.

Caveat: your answer gives much better precision. A plot of xfunc[t] - pf[t] shows discrepancies as high as 0.001, i.e. almost 1% of the range spanned by the points.

Answer 3

How I'd do it:

pts = {{0.504, 2.79}, {0.519, 2.7773}, {0.5349, 2.7642}, {0.5504, 2.7515},
       {0.5666, 2.7398}, {0.5858, 2.7341}, {0.5914, 2.8566}, {0.5917, 2.8364},
       {0.5918, 2.8766}, {0.5924, 2.8164}, {0.5933, 2.7963}, {0.5935, 2.7525},
       {0.5935, 2.8966}, {0.594, 2.7763}, {0.5974, 2.9163}};

(* reordering *)
pts = pts[[First[FindCurvePath[Standardize[pts]]]]];

(* Lee's algorithm *)
parametrize[pts_List, a : (_?NumericQ) : 1/2] := 
            FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]] /; 
            MatrixQ[pts, NumericQ]

tvals = parametrize[pts]; (* parameter values *)

n = Length[pts]; m = 3; (* order of spline *)
knots = Join[ConstantArray[tvals[[1]], m + 1], 
             If[m + 2 <= n, MovingAverage[ArrayPad[tvals, -1], m], {}], 
             ConstantArray[tvals[[-1]], m + 1]];

(* spline control points *)
cp = LinearSolve[Outer[BSplineBasis[{m, knots}, #2, #1] &,
                       tvals, Range[0, n - 1], 1], pts];

(* B-splines in explicit piecewise polynomial form *)
{fc[t_], gc[t_]} = 
    Table[BSplineBasis[{m, knots}, j, t] // PiecewiseExpand, {j, 0, n - 1}].cp //
    Simplify

I've omitted the result of that last one, but you can see the piecewise nature for yourself if you execute the previous lines. Now, a plot:

ParametricPlot[{fc[t], gc[t]}, {t, 0, 1}, 
               AspectRatio -> 1/GoldenRatio, Frame -> True]

(Some of you who are familiar with my previous answers will recognize the pieces I used.)