# 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. Oct 9, 2012 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. Oct 9, 2012 at 19:11
• To clarify: do you want to do an interpolation, or a regression? PlatoManiac's answer shows both routes. Oct 9, 2012 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. Oct 10, 2012 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. Oct 10, 2012 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. Oct 9, 2012 at 19:07
• @rcollyer thanks for pointing out! I tried to modify the answer.. Oct 9, 2012 at 19:33
• @PlatoManiac I suggest deleting the LinearModelFit part, as it could be misleading for future readers Oct 10, 2012 at 3:58

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[]. Oct 10, 2012 at 11:26
• @J.M. One could only if one knows it :) . Thanks Oct 10, 2012 at 11:28

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. Oct 10, 2012 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. Oct 10, 2012 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. Oct 10, 2012 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... Oct 11, 2012 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. Oct 11, 2012 at 1:56