22
$\begingroup$

I've used Interpolation[] to generate an InterpolatingFunction object from a list of integers.

f = Interpolation[{2, 5, 9, 15, 22, 33, 50, 70, 100, 145, 200, 280, 375, 495, 635, 800,
                   1000, 1300, 1600, 2000, 2450, 3050, 3750, 4600, 5650, 6950}]

I'm using that to generate values like f[27], f[28], ...

Is there any way to print or show the function used by Mathematica that produced the result of f[27]?

$\endgroup$
  • 1
    $\begingroup$ You can view the internal parts of an InterpolatingFunction using "Methods" as I outlined here: (19042) however that only shows the data that is used and the kind of interpolation, not the actual function itself. You might find value in InterpolatingPolynomial, though it is not the same as Interpolation. $\endgroup$ – Mr.Wizard Sep 17 '14 at 19:33
  • 5
    $\begingroup$ See Properties and Relations for an example; see also Some Notes on Internal Implementation. $\endgroup$ – Michael E2 Sep 17 '14 at 19:34
  • $\begingroup$ I am pretty sure that Mathematica uses interpolating polynomials, using Neville's algorithm. It is easy to implement this yourself and see if you get the same value. $\endgroup$ – Igor Rivin Sep 17 '14 at 23:57
  • 2
    $\begingroup$ Maybe you want InterpolatingPolynomial instead of Interpolation. $\endgroup$ – Szabolcs Oct 12 '17 at 13:07
  • 1
    $\begingroup$ @Szabolcs if I'm reading this right, certainly not. That one gives a polynomial of arbitrarily high degree which is probably not that useful. $\endgroup$ – LLlAMnYP Oct 12 '17 at 13:23
22
$\begingroup$

Here is the example from the documentation adapted for the OP's data:

data = MapIndexed[
   Flatten[{#2, #1}] &,
   {2, 5, 9, 15, 22, 33, 50, 70, 100, 145, 200, 280, 375, 495, 635, 800,
    1000, 1300, 1600, 2000, 2450, 3050, 3750, 4600, 5650, 6950}];
f = Interpolation@data
(* InterpolatingFunction[{{1, 26}}, <>] *)
pwf = Piecewise[
     Map[{InterpolatingPolynomial[#, x], x < #[[3, 1]]} &, Most[#]], 
     InterpolatingPolynomial[Last@#, x]] &@Partition[data, 4, 1];

Here is a comparison of the piecewise interpolating polynomials and the interpolating function:

Plot[f[x] - pwf, {x, 1, 28}, PlotRange -> All]

Mathematica graphics

The values of f[27] and f[28] are beyond the domain, which is 1 <= x <= 26, and extrapolation is used. The formula for extrapolation is given by the last InterpolatingPolynomial in pwf:

Last@pwf
(* 3750 + (850 + (100 + 25/3 (-25 + x)) (-24 + x)) (-23 + x) *)

In response to a comment: The error in the plot has to do with round-off error. Apparently the calculation done by InterpolatingFunction, while algebraically equivalent, is not numerically identical. The error was greatest above in the domain 26 < x < 28 where extrapolation is performed. With arbitrary precision, the error is zero, as shown below.

Plot[f[x] - pwf, {x, 1, 28}, PlotRange -> All, 
 WorkingPrecision -> $MachinePrecision, Exclusions -> None, PlotStyle -> Red]

Mathematica graphics

$\endgroup$
  • $\begingroup$ (+1) The data in your answer should be defined as simple 1D list in order your code to work properly. How do you explain relatively large differences (1*10^-12) between your implementation and the built-in? $\endgroup$ – Alexey Popkov Sep 18 '14 at 6:36
  • $\begingroup$ @AlexeyPopkov Thanks. See the update for the error explanation. (I copied the wrong code by mistake. It should work now.) $\endgroup$ – Michael E2 Sep 18 '14 at 11:16
  • 1
    $\begingroup$ Thanks, it is clear now. I would add that plotting relative error would be more correct and it clearly shows that the relative error is of expected magnitude of lesser than 10^-15: Plot[(f[x] - pwf)/pwf, {x, 1, 28}, PlotRange -> All]. $\endgroup$ – Alexey Popkov Sep 18 '14 at 11:32
  • 1
    $\begingroup$ @LLlAMnYP There is no example in the linked dupe, but, yes, the old code assumed the abscissae were successive integers. Fixed now, I think. Let me know if it doesn't work for you. It won't work on all InterpolatingFunction methods, I think, just this type. $\endgroup$ – Michael E2 Oct 12 '17 at 15:39
  • 1
    $\begingroup$ @LLlAMnYP The method above does work on unequally spaced points, but not if the spacing varies wildly. It appears that InterpolatingFunction automatically splits data at points that are too close (I can get it to match if there aren't too many too-close points in a row). Or maybe it adapts the order locally. I don't know the criterion used. It is perhaps if the slope is too great in some relative sense. The result of InterpolatingFunction seems worse than the piecewise one in the example I got from your (random) code. Interesting -- didn't know it did such complicated things. $\endgroup$ – Michael E2 Oct 13 '17 at 11:49
17
$\begingroup$

There is no documented built-in way to convert the InterpolatingFunction object into explicit Piecewise form (thanks to @MichaelE2 for the link!). So the only possibility to get an explicit interpolating function is to re-implement the built-in Interpolation in the high-level Mathematica language. I have already done this for the built-in "Spline" method with InterpolationOrder -> 2 (quadratic spline interpolation with splicing points in the middle of adjacent interpolation points). Spline interpolation in general gives much better results than the default "Hermite" method.

You can use my implementation of quadric spline interpolation in Mathematica to produce an explicit Piecewise function interpolating your data (as opposed to the built-in, it supports arbitrary precision!):

data = Transpose[{Range[Length[#]], #}] &@{2, 5, 9, 15, 22, 33, 50, 70, 100, 145};
spline[\[FormalX]_] = makeSpline[toSplineData[data], \[FormalX]]

screenshot

Here is a comparison with the original data and with the built-in Interpolation:

Table[data[[x, 2]] - spline[x], {x, 10}]
f = Interpolation[data, Method -> "Spline", InterpolationOrder -> 2];
Table[f[x] - spline[x], {x, 10}]
Plot[(f[x] - spline[x])/spline[x], {x, 1, 10}, PlotRange -> All]
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
{0., 0., 0., 1.77636*10^-15, 0., 0., 0., 0., 0., 0.}

plot

$\endgroup$
  • $\begingroup$ Is it equivalent to what Mma implements for Interpolation[pts]? $\endgroup$ – Dr. belisarius Sep 18 '14 at 1:24
  • $\begingroup$ @belisarius It is equivalent to Interpolation[pts, Method -> "Spline", InterpolationOrder -> 2] with the only difference that my implementation supports arbitrary precision while the current Mathematica's built-in for "Spline" allows only MachinePrecision. $\endgroup$ – Alexey Popkov Sep 18 '14 at 1:28
  • $\begingroup$ I answered the same question where you wrote that answer, but I think the OP is asking another specific thing here about the InterpolatingFunctionoject $\endgroup$ – Dr. belisarius Sep 18 '14 at 1:31
  • $\begingroup$ You know that there is no documented built-in way to convert InterpolatingFunction into explicit Piecewise form. I am sure that there is no undocumented way too because (at least "Hermite" method) is written in C, not in the high-level Mathematica language. So the only possibility to answer the question is to provide a high-level (re)implementation of the Interpolation. $\endgroup$ – Alexey Popkov Sep 18 '14 at 1:38
7
$\begingroup$

You could use Series. What's necessary is to know which abscissa values were used for the interpolation. Let's generate some fake data.

xVals = RandomReal[{0, 100}, 35] // Sort;
yVals = Sin[xVals/30] + RandomReal[.1, 35];
ListPlot[{xVals, yVals}\[Transpose]]

enter image description here

Create the interpolation:

iPol = Interpolation[{xVals, yVals}\[Transpose]]

Get the Series expansion at some point:

Series[iPol[x], {x, xVals[[5]], 4}]
0.469053+0.0487866 (x-14.3715)+0.011059 (x-14.3715)^2+0.00100578 (x-14.3715)^3+O[x-14.3715]^5

Note that there is no 4th order term (as expected, since default interpolation order is 3).

Let's generate a piecewise function:

pW[x_] = Piecewise[
  Table[{Normal@Series[iPol[x], {x, x0, 3}], x < x0}, {x0, xVals}]]

enter image description here

Test it out:

Plot[pW[x] - iPol[x], {x, xVals[[1]], xVals[[-1]]} PlotRange->Full]

enter image description here

Looks like (almost) machine precision to me

$\endgroup$
  • $\begingroup$ Thanks for your answer. That's indeed a uniform way to find the explicit representation. $\endgroup$ – user52830 Oct 12 '17 at 13:22
  • $\begingroup$ @user52830 You're welcome. Barring edge cases, if you happen to have just the InterpolatingFunction object, you can modify my Table statement as Table[..., {x0, iPol[[3,1]]}] since iPol[[3,1]] === xVals $\endgroup$ – LLlAMnYP Oct 12 '17 at 13:26
5
$\begingroup$

Here is a (mostly) general routine that (tries to) convert a one-dimensional InterpolatingFunction[] into an equivalent Piecewise[] function:

convertToPiecewise::umet = "Unknown interpolation method `1`.";
SetAttributes[convertToPiecewise, Listable];

convertToPiecewise[iF_InterpolatingFunction, x_,
                   OptionsPattern[{"Extrapolation" -> False,
                                   InterpolationOrder -> Automatic}]] := 
       Module[{bf, extQ, imet, kp, makePP, met, nodes, perQ, pieces, pts, vand, xt},
              Switch[met = iF["InterpolationMethod"],
                     "Hermite" | "Chebyshev",
                     pts = Transpose[{Flatten[iF["Grid"]], iF["ValuesOnGrid"]}],
                     "BSpline",
                     bf = First[Cases[iF, _BSplineFunction, ∞]];
                     pts = {#, bf[#]} & /@ Union[First[bf["Knots"]]],
                     _,
                     Message[convertToPiecewise::umet, met]; Return[$Failed, Module]];
              kp = OptionValue[InterpolationOrder];
              If[kp === Automatic,
                 (* repeated differentiation to determine maximal order *)
                 kp = First[NestWhile[Derivative[1], 
                            Derivative[First[iF["InterpolationOrder"]] + 1][iF],
                            (Norm[#["ValuesOnGrid"], ∞] > 0) &]["DerivativeOrder"]] - 1];
              If[kp > 0, (* normal case *)
                 (* use equispaced nodes in the exact case, and Chebyshev otherwise *)
                 nodes = Range[kp - 1]/kp;
                 If[MatrixQ[pts, InexactNumberQ],
                    nodes = N[Haversine[π nodes], Precision[pts]]];
                 vand = LinearAlgebra`Private`VandermondeSolve[##, Transpose -> True] &;
                 makePP[{{x1_, y1_}, {x2_, y2_}}] := Module[{h = x2 - x1, ip},
                     ip = Transpose[Join[{{x1 - x1, y1}},
                                         {#, iF[x1 + #]} & /@ (h nodes), {{h, y2}}]];
                     (* solve for interpolating polynomial coefficients *)
                     {Fold[(#1 (xt - x1) + #2) &, Reverse[vand @@ ip]], x1 <= xt <= x2}];
                 pieces = makePP /@ Partition[pts, 2, 1],
                 (* zero-order interpolation *)
                 pieces = Transpose[{Rest[pts[[All, -1]]],
                                     #1 <= xt <= #2 & @@@
                                     Partition[pts[[All, 1]], 2, 1]}]];
              perQ = TrueQ[First[iF["Periodicity"]]];
              If[! perQ, extQ = OptionValue["Extrapolation"];
                 If[! ListQ[extQ], extQ = {extQ, extQ}]; extQ = TrueQ /@ extQ;
                 If[extQ[[1]], pieces[[1, 2]] = Delete[pieces[[1, 2]], 1]];
                 If[extQ[[2]], pieces[[-1, 2]] = Delete[pieces[[-1, 2]], 3]]];
                 Piecewise[pieces /. xt ->
                           If[! perQ, x, Mod[x, #2 - #1, #1] & @@ First[iF["Domain"]]]]]

(N.B. replace LinearAlgebra`Private`VandermondeSolve[] with LinearAlgebra`VandermondeSolve[] when using the function in versions before 11.2.)

It should work for InterpolatingFunction[] objects that come from Interpolation[], ListInterpolation[], or FunctionInterpolation[]. It mostly works for InterpolatingFunction[] objects from NDSolve[], but may fail in some cases. (If you find an example, please tell me!)

Some examples:

if1 = Interpolation[{1, 3, 5, 2, 1}, InterpolationOrder -> 1];

Here, we tell convertToPiecewise[] to use the rightmost piece for extrapolation to the right:

convertToPiecewise[if1, x, "Extrapolation" -> {False, True}]

$$\begin{cases} 2 (x-1)+1 & 1\leq x\leq 2 \\ 2 (x-2)+3 & 2\leq x\leq 3 \\ 5-3 (x-3) & 3\leq x\leq 4 \\ 6-x & 4\leq x \\ 0 & \mathtt{True} \end{cases}$$

Convert an InterpolatingFunction[] with irregular spacing:

if2 = Interpolation[{{0, 0}, {0.1, .3}, {0.5, .6}, {1, -.2}, {2, 3}}, Method -> "Spline"];

pw2[x_] = convertToPiecewise[if2, x];

Plot[if4[x] - pw4[x], {x, 0, 2}, PlotRange -> All]

spline interpolant error

Convert the result of NDSolve[]:

if3 = NDSolveValue[{g'[x] == Sin[2 x] - g[x], g[0] == 1}, g, {x, 0, 6}];

pw3[x_] = convertToPiecewise[if3, x];

Plot[if3[x] - pw3[x], {x, 0, 6}, PlotRange -> All]

error of interpolant from NDSolve

Another NDSolve[] example. The previous version of the routine was unable to handle this.

if4 = NDSolveValue[{y''[t] == 10 (1 - y[t]^2) y[t] - y[t], y[0] == 2, y'[0] == 0},
                   y, {t, 0, 6}, Method -> "StiffnessSwitching"];

pw4[x_] = convertToPiecewise[if4, x];

Plot[if4[x] - pw4[x], {x, 0, 6}, PlotRange -> All]

error of interpolant from NDSolve, second order

$\endgroup$
  • $\begingroup$ I seem to have to use LinearAlgebra`Private`VandermondeSolve; I can't get LinearAlgebra`VandermondeSolve to evaluate on anything. (V11.2, Mac). $\endgroup$ – Michael E2 Mar 24 '18 at 22:28
  • $\begingroup$ Ah, you're right, I was doing this in 11.1. I'll add a note... $\endgroup$ – J. M. will be back soon Mar 25 '18 at 1:02
4
$\begingroup$

Instead of trying to come up with a function that replicates the output of the InterpolatingFunction, one could instead "inactivate" it, and then see what kind of arithmetic it does. For instance:

inactiveIF = MapAt[HoldForm, f, {-2, All, 1}];

For example:

r = inactiveIF[10.5]

200-0.5 (-1. (145-200)+0.5 (1. (1. (145-200)+0.5 (-145+280))+0.75 (-1. (-0.333333 (100-280)-0.5 (-145+280))-1. (1. (145-200)+0.5 (-145+280)))))

There are hidden HoldForm wrappers in the above code. If we use ReleaseHold on r we get the same value as using f directly:

ReleaseHold[r]
f[10.5]

170.313

170.313

$\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.