# Extract function from a parametric interpolating function

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]]]).

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.

• Upvote for a very precise solution. I'm not sure, what causes InterpolatingPolynomial to stumble on your data, when it clearly works very well in the doc's example and on symbolic data. Commented Sep 25, 2015 at 16:17

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.

• I first thought the discrepancy was due to xfunc being calculated with Hermite polynomials while pf with splines, but that's apparently not the case. I wonder why the discrepancies are so high. Commented Sep 25, 2015 at 16:24

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.)

• The output {fc[t_],gc[t_]} could even be cleaned by removing the redundant ponctual values (t == 0.786984 and 0.786984 < t < 0.858021 chanded into 0.786984 < =t < 0.858021. I'll add it as comment when I find time to do it. Commented Oct 8, 2015 at 15:06
• Actually FullSimplify[] does a better job, but it takes longer. Commented Oct 8, 2015 at 17:29