Lagrange interpolation

So, the problem is to find a polynomial p(x) whose values are specified at finitely many points, ie. p(x_i) = f_i (I'll say that there are n such points). To do this, all you need is the following simple trick:

If the polynomial functions p_i(x) are defined in the following manner:

p_i(x) = (product over j, j!=i) (x-x_j)/(x_i-x_j)

then p_i(x_j) = δ_ij (n.b. each of these polynomials are of order n-1).

Given this, simply define the required polynomial as p(x) = (sum over i) f_ip_i(x), and then a simple check reveals that p(x_i) = f_i, as intended. The p(x) is of order n-1, and is the unique polynomial of this order satisfying the condition (since if p(x) and q(x) both do, then (p-q)(x) has n zeros at the x_i, so is the zero polynomial).
As to whether the interpolation is any good elsewhere, it depends on how widely separated you x_i values are - in short, if they're all bunched up together, then further away from the points, p(x) will be either growning or shrinking enormously, so it's not you're best bet for predicting stock market prices.

Using a Lagrange Polynomial approximation, one can generate a single polynomial expression which passes through every data point given. This requires no additional information about the points (such as what the original function's derivative might be there, etc). As a result of this, it's not a great match, but it gets the job done in some cases.

Basics first. Given n points, there will be a unique polynomial of degree n-1 that passes through all of these points. The simplest case is approximating a line through two points. You may be used to doing this with standard form or slope-intercept form. There's another way. Given a function passing through two points, let's call them a and b, at x_a this function should evaluate to y_a, and at x_b it should evaluate to y_b. If this hypothetical function is called f(x), then we already know f(x_a)=y_a and f(x_b) = y_b. Therefore we just need a function L(x) for each point n which is equal to 1 at that point and 0 at the other, and then our approximating polynomial is:

P(x) = L₀(x)f(x₀) + L₁(x₁)f(x₁)

So if L₀(x₀) is 1, then P(x₀) must be f(x₀)! Maybe you've already figured out what the L_n(x) equations are. If so, you're a better (wo)man than I, because it took me a while to really get ahold of this. For the two-point case:

L₀(x) = (x - x₁) / (x₀ - x₁)

and

L₁(x) = (x - x₀) / (x₁ - x₀)

You'll note if you plug in x₀ and x₁ to those equations, they'll be either 1 or 0, depending on the point. This can be extended to multiple points, which starts to make a mess. The general formula for each L(x) equation from a set of n points at point k is:

                 (x-x₀)...(x-x_k-1)(x-x_k+1)...(x-x_n)
L_n,k(x) = -------------------------------------------
                 (x_k - x₀)...(x_k-x_k-1)(x_k-x_k+1)...(x_k-x_n)

Notice that, at the top, (x-x_k) is skipped, and at the bottom, (x_k - x_k) is skipped. That's good, because the latter would make the value of the function undefined.

And you string n of these together via the previous example to get the new P(x). This can be really bad in some cases, as for large numbers of points you get very high-degree polynomials, which tend to oscillate violently, especially if your points are not so close together. Other methods such as osculating polynomial, hermite polynomials, and cubic splines may yield more accurate results.

Additionally, multiplying these things out is a job best left to Mathematica. There's also a more information-obscuring but easier-to-do method for obtaining this same polynomial known as the Divided Difference Method.