13
$\begingroup$

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?

Improve this question
$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – rcollyer
    Commented Oct 9, 2012 at 19:05
  • $\begingroup$ 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. $\endgroup$
    – Searke
    Commented Oct 9, 2012 at 19:11
  • $\begingroup$ To clarify: do you want to do an interpolation, or a regression? PlatoManiac's answer shows both routes. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Oct 9, 2012 at 22:06
  • $\begingroup$ @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. $\endgroup$
    – BOMEz
    Commented Oct 10, 2012 at 12:45
  • 2
    $\begingroup$ @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. $\endgroup$
    – rcollyer
    Commented Oct 10, 2012 at 12:48

3 Answers 3

Reset to default
7
$\begingroup$

First take your data

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]

enter image description here

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.

enter image description here

Improve this answer
$\endgroup$
3
  • 4
    $\begingroup$ 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. $\endgroup$
    – rcollyer
    Commented Oct 9, 2012 at 19:07
  • $\begingroup$ @rcollyer thanks for pointing out! I tried to modify the answer.. $\endgroup$
    – PlatoManiac
    Commented Oct 9, 2012 at 19:33
  • $\begingroup$ @PlatoManiac I suggest deleting the LinearModelFit part, as it could be misleading for future readers $\endgroup$
    – Dr. belisarius
    Commented Oct 10, 2012 at 3:58
8
$\begingroup$

You may want a quick sketch :)

Mathematica graphics

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

Mathematica graphics

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

Improve this answer
$\endgroup$
2
  • 3
    $\begingroup$ ...or, one could use the (undocumented) function LinearAlgebra`VandermondeSolve[]. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Oct 10, 2012 at 11:26
  • 2
    $\begingroup$ @J.M. One could only if one knows it :) . Thanks $\endgroup$
    – Dr. belisarius
    Commented Oct 10, 2012 at 11:28
4
$\begingroup$

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

enter image description here

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

enter image description here

Improve this answer
$\endgroup$
7
  • 1
    $\begingroup$ The output of InterpolatingPolynomial[] is in Newton form, not Horner form. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Oct 10, 2012 at 12:54
  • 1
    $\begingroup$ @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. $\endgroup$
    – rcollyer
    Commented Oct 10, 2012 at 19:10
  • $\begingroup$ @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. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Oct 10, 2012 at 23:55
  • $\begingroup$ 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... $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Oct 11, 2012 at 1:28
  • $\begingroup$ @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. $\endgroup$
    – rcollyer
    Commented Oct 11, 2012 at 1:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.