# 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. – LLlAMnYP Sep 25 '15 at 16:17

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. – anderstood Oct 8 '15 at 15:06
• Actually FullSimplify[] does a better job, but it takes longer. – J. M. is away Oct 8 '15 at 17:29

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. – anderstood Sep 25 '15 at 16:24