2
$\begingroup$

I've read the answers of

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

enter image description here

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

$\endgroup$
2
$\begingroup$

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

enter image description here

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – LLlAMnYP Sep 25 '15 at 16:17
2
$\begingroup$

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]

plot of the spline

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – anderstood Oct 8 '15 at 15:06
  • $\begingroup$ Actually FullSimplify[] does a better job, but it takes longer. $\endgroup$ – J. M. will be back soon Oct 8 '15 at 17:29
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – anderstood Sep 25 '15 at 16:24

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.