# Get polynomial interpolation formula

I'm attempting to get a polynomial interpolation formula out of Mathematica but I am absolutely lost. I stared out using Wolfram|Alpha, but it seems as if my input had become too large.

I tried using this on WolframAlpha.com :

Interpolate {1,33},{2,80},{5,286},{10,771},{15,1382},{20,2087},{25,2867},{30,3707},{40,5526},
{50,7470},{60,9482},{70,11507},{80,13495},{90,15391},{100,17313},{110,18631},{120,19752},
{125,20064}


but I get an error.

Thinking that Mathematica was the solution I signed up for a trial but I can't seem to get an interpolation formula out of it. Searching Google and this exchange has given me very few results.

Any idea how I can just get a polynomial interpolation formula from my dataset?

-
Have you tried using the Documentation Center under the Help menu? I found the answer in there, but I'd like to know you've done a little searching yourself, first. –  rcollyer Oct 9 '12 at 19:05
Once you write out the equations for coefficients of the polynomial, it becomes a linear algebra problem with a vandermonde matrix en.wikipedia.org/wiki/Vandermonde_matrix reference.wolfram.com/mathematica/ref/… Or you can use the documentation center and search for Interpolation. –  Searke Oct 9 '12 at 19:11
To clarify: do you want to do an interpolation, or a regression? PlatoManiac's answer shows both routes. –  Guess who it is. Oct 9 '12 at 22:06
@rcollyer I did do some searching on my own using Google and ended up on the Mathematica online help. The results I got to though didn't give me what I was looking for.I was confused because I couldn't just get a straight forward formula like I could on Wolframalpha. From the responses below it seems more involved and when I felt I was in over my head I decided to ask. –  BOMEz Oct 10 '12 at 12:45
@BOMEz a useful and necessary skill is to be able to search within the documentation itself, and not via google. Obviously, you know the answer now, but try searching for "interpolating polynomial" in the documentation center and see what you get. I got to it by going through InterpolationFunction and seeing the related functions. –  rcollyer Oct 10 '12 at 12:48

data = {{1, 33}, {2, 80}, {5, 286}, {10, 771}, {15, 1382}, {20,
2087}, {25, 2867}, {30, 3707}, {40, 5526}, {50, 7470}, {60,
9482}, {70, 11507}, {80, 13495}, {90, 15391}, {100, 17313}, {110,
18631}, {120, 19752}, {125, 20064}};


Then we call LinearModelFit and fit a cubic polynomial to your data.

lm = LinearModelFit[data, {x^3, x^2, x}, x];
Show[ListPlot[data, PlotStyle -> Red,Filling->Bottom],Plot[lm[x],{x, 0, 125}],Frame -> True]


And to get the polynomial that best fits your data.

Normal[lm]


-83.6419 + 69.7325 x + 2.19787 x^2 - 0.0116981 x^3

Now you must realize that above polynomial is not an interpolation of your data but a best continuous approximation of it with respect to the Euclidean norm. Forming a interpolating polynomial for a data of $n$ points require at least a $n$-th degree polynomial. This is not practical as higher degree polynomials come with higher and unwanted oscillations. Hence people use polynomials for peace-wise interpolation. HermitePolynomil can be used for this purpose.

But if you want you can get an interpolating polynomial for your data as follows.

InterpolatingPolynomial[data, x] // Expand // N


3.1289 + 18.9719 x + 12.41 x^2 - 1.71825 x^3 + 0.228052 x^4 - 0.0219672 x^5 + 0.00150955 x^6 - 0.0000751269 x^7 + 2.75112*10^-6 x^8 - 7.49912*10^-8 x^9 + 1.53098*10^-9 x^10 - 2.34119*10^-11 x^11 + 2.66291*10^-13 x^12 - 2.216*10^-15 x^13 + 1.30841*10^-17 x^14 - 5.18423*10^-20 x^15 + 1.2348*10^-22 x^16 - 1.33508*10^-25 x^17

You can see the above polynomial has degree 17 and you have 18 data entries. This interpolating polynomial also fails to keep up with the data trend. You can see this phenomena in the boundaries of the following plot. For $x$ outside your data range the mammoth polynomial (dashed one) starts to disagree with given data trend almost exponentially! The approximating cubic polynomial (thick red one) captures the data behavior very well even beyond the given data range.

-
A linear fit is not guaranteed to go through every point, but an interpolating polynomial is. So, while this answer a question, it is not the same one as asked. –  rcollyer Oct 9 '12 at 19:07
@rcollyer thanks for pointing out! I tried to modify the answer.. –  PlatoManiac Oct 9 '12 at 19:33
@PlatoManiac I suggest deleting the LinearModelFit part, as it could be misleading for future readers –  belisarius Oct 10 '12 at 3:58

Evaluating what the OP tried with free-form input (shortcut =) yields a suggestion cell with this input : Expand[ InterpolatingPolynomial[{286, 771}, x]], similarly a Wolfram|Alpha query yields input interpretation interpolating polynomial, so it is not especially difficult to surmise that the desired function is InterpolatingPolynomial. Defining the list l :

l = {{1, 33}, {2, 80}, {5, 286}, {10, 771}, {15, 1382}, {20, 2087}, {25, 2867}, {30, 3707},
{40, 5526}, {50, 7470}, {60, 9482}, {70, 11507}, {80, 13495}, {90, 15391}, {100, 17313},
{110, 18631}, {120, 19752}, {125, 20064}};


the resulting polynomial (written in a standard expanded form ) is :

p[x_] := Expand @ InterpolatingPolynomial[ l, x]


## How could we define an interpolating polynomial having no InterpolatingPolynomial built-in function ?

A natural way is to define recursively a family of polynomials. In order to make it possibly fast we can take advantage of memoization techiques (see e.g. answers to What is “x := x =” trickery? or Why does Expand not work within a function? )

poly[x_, 0]  := l[[1, 2]]
poly[x_, k_] := poly[ x, k] = poly[ x, k-1] +
(l[[k, 2]] - poly[l[[k, 1]], k-1]) Product[(x-l[[j, 1]])/(l[[k, 1]]-l[[j, 1]]),
{j, k - 1}]


let's check if it works properly :

Expand @ poly[ x, Length @ l] === p[x]

True


Well indeed, now we can imagine what is behind InterpolatingPolynomial so let's restrict to p[x] being a seventeenth order polynomial with rational coefficients (in general when the list is of length n then the polynomial will be of order n-1, and if l is a list of rational pairs, p will be a polynomial with rational coefficients). One can reproduce the list l with p[x] :

p[#] & /@ l[[All, 1]] == l[[All, 2]]

True


Since its coefficients are a bit involved we prefer to write them this way :

p[x] // TraditionalForm // N


Plot[ p[x], {x, -12, 128}, PlotStyle -> Thick, Epilog -> {Red, PointSize[0.01], Point[l]}]


-
The output of InterpolatingPolynomial[] is in Newton form, not Horner form. –  Guess who it is. Oct 10 '12 at 12:54
@J.M. according to the docs (4th under more info) the polynomial is output in Horner form so that it is suitable for numerical evaluation. So, however it generates the polynomial, and I suspect that it is Newton-Cotes, it rearranges it into Horner's form. –  rcollyer Oct 10 '12 at 19:10
@rcollyer, the docs are wrong, then. Compare InterpolatingPolynomial[{{a, b}, {c, d}, {p, q}, {r, s}}, x] and HornerForm[InterpolatingPolynomial[{{a, b}, {c, d}, {p, q}, {r, s}}, x], x]. Horner form is nothing but the special case of the Newton form, when all the interpolation points coalesce. It looks like Horner, but it's not really Horner. –  Guess who it is. Oct 10 '12 at 23:55
BTW @rcollyer: "Netwon-Cotes" is the name used for the series of integration rules based on interpolating polynomials on equispaced points; Newton (divided difference) interpolation (which is what InterpolatingPolynomial[] does) is a whole 'nother thing altogether... –  Guess who it is. Oct 11 '12 at 1:28
@J.M. I must disagree. Newton is written as a sum being a linear combination of polynomials. While debating the content of wikipedia is often pointless, your link does not show anything near Horner's form. While InterpolatingPolynomial doesn't return what we'd precisely call Horner's form, it is a generalization that I have not seen associated with Newton, so I'd call it Horner's form, or generalized Horner's, if I'm being pedantic. And, yes, I mis-remembered Newton-Cotes. –  rcollyer Oct 11 '12 at 1:56

You may want a quick sketch :)

Or, you could also solve the Vandermonde matrix:

n = 10;(* points number*)
{xi, yi} = RandomReal[{0, 1}, {2, n}];
a = LinearSolve[Transpose@Table[xi^k, {k, 0, n - 1}], yi]

(* System solved, now plot it *)
p[x_, m_] := Sum[a[[i + 1]] x^i, {i, 0, m}];
Show[Plot[p[x, n - 1], {x, 0, 1},  PlotRange -> {Automatic, {-20, 20}}],
ListPlot[Transpose[{xi, yi}], PlotStyle -> {PointSize[Medium], Red}]]


Just be aware that it may, depending on your points, results in an ill-conditioned system

-
...or, one could use the (undocumented) function LinearAlgebraVandermondeSolve[]`. –  Guess who it is. Oct 10 '12 at 11:26
@J.M. One could only if one knows it :) . Thanks –  belisarius Oct 10 '12 at 11:28