As there is only one interpolation polynomial for a given set of data points it is a bit misleading to call the polynomial the Lagrange interpolation polynomial. The more precise name is interpolation polynomial in the Lagrange form.

Definition

Given a set of k + 1 data points

$(x\_0,\; y\_0),ldots,(x\_k,\; y\_k)$

where no two x_j are the same, the interpolation polynomial in the Lagrange form is a linear combination

$L(x)\; :=\; sum\_\{j=0\}^\{k\}\; y\_j\; ell\_j(x)$

of Lagrange basis polynomials

$ell\_j(x)\; :=\; prod\_\{i=0,,\; ineq\; j\}^\{k\}\; frac\{x-x\_i\}\{x\_j-x\_i\}\; =\; frac\{(x-x\_0)\}\{(x\_j-x\_0)\}\; cdots\; frac\{(x-x\_\{j-1\})\}\{(x\_j-x\_\{j-1\})\}\; frac\{(x-x\_\{j+1\})\}\{(x\_j-x\_\{j+1\})\}\; cdots\; frac\{(x-x\_\{k\})\}\{(x\_j-x\_\{k\})\}.$

Proof

The function we are looking for has to be a polynomial function L(x) of degree less than or equal to k with

$L(x\_j)\; =\; y\_j\; qquad\; j=0,ldots,k$

The Lagrange polynomial is a solution to the interpolation problem.

As can be seen

1. $ell\_j(x)$ is a polynomial and has degree k.
2. $ell\_i(x\_j)\; =\; delta\_\{ij\},quad\; 0\; leq\; i,j\; leq\; k.,$

Thus the function L(x) is a polynomial with degree at most k and

$L(x\_i)\; =\; sum\_\{j=0\}^\{k\}\; y\_j\; ell\_j(x\_i)\; =\; y\_i.$

There can be only one solution to the interpolation problem since the difference of two such solutions would be a polynomial with degree at most k and k+1 zeros. This is only possible if the difference is identically zero, so L(x) is the unique polynomial interpolating the given data.

Main idea

Solving an interpolation problem leads to a problem in linear algebra where we have to solve a matrix. Using a standard monomial basis for our interpolation polynomial we get the very complicated Vandermonde matrix. By choosing another basis, the Lagrange basis, we get the much simpler identity matrix = δ_i,j which we can solve instantly.

Implementation in C

Note : "pos" and "val" arrays are of size "degree".

float lagrangeInterpolatingPolynomial (float pos[], float val[], int degree, float desiredPos) {

  float weight;

  float retVal = 0;

  for (int i=0; i <= degree; i++) {

     weight = 1;

     for (int j=0; j
        weight *= (desiredPos - pos[j]) / (pos[i] - pos[j]);
      }
     for(int j=i+1; j<=degree; j++) {
        weight *= (desiredPos - pos[j]) / (pos[i] - pos[j]);
      }
     retVal += weight * val[i];
   } 
  return retVal;
} 

Usage
 Example 1 
We wish to interpolate $f(x)=tan(x)$ at the points 


 $x\_0=-1.5,$  
  $f(x\_0)=-14.1014,$ 

 $x\_1=-0.75,$ 
  $f(x\_1)=-0.931596,$  

 $x\_2=0,$     
  $f(x\_2)=0,$ 

$x\_3=0.75,$   
  $f(x\_3)=0.931596,$ 

$x\_4=1.5,$    
  $f(x\_4)=14.1014,$ 
The basis polynomials are:
$ell\_0(x)=\{x\; -\; x\_1\; over\; x\_0\; -\; x\_1\}cdot\{x\; -\; x\_2\; over\; x\_0\; -\; x\_2\}cdot\{x\; -\; x\_3\; over\; x\_0\; -\; x\_3\}cdot\{x\; -\; x\_4\; over\; x\_0\; -\; x\_4\}$
             ={1over 243} x (2x-3)(4x-3)(4x+3)
$ell\_1(x)=\{x\; -\; x\_0\; over\; x\_1\; -\; x\_0\}cdot\{x\; -\; x\_2\; over\; x\_1\; -\; x\_2\}cdot\{x\; -\; x\_3\; over\; x\_1\; -\; x\_3\}cdot\{x\; -\; x\_4\; over\; x\_1\; -\; x\_4\}$
             =-{8over 243} x (2x-3)(2x+3)(4x-3)
$ell\_2(x)=\{x\; -\; x\_0\; over\; x\_2\; -\; x\_0\}cdot\{x\; -\; x\_1\; over\; x\_2\; -\; x\_1\}cdot\{x\; -\; x\_3\; over\; x\_2\; -\; x\_3\}cdot\{x\; -\; x\_4\; over\; x\_2\; -\; x\_4\}$
             ={3over 243} (2x+3)(4x+3)(4x-3)(2x-3) 
$ell\_3(x)=\{x\; -\; x\_0\; over\; x\_3\; -\; x\_0\}cdot\{x\; -\; x\_1\; over\; x\_3\; -\; x\_1\}cdot\{x\; -\; x\_2\; over\; x\_3\; -\; x\_2\}cdot\{x\; -\; x\_4\; over\; x\_3\; -\; x\_4\}$
             =-{8over 243} x (2x-3)(2x+3)(4x+3)
$ell\_4(x)=\{x\; -\; x\_0\; over\; x\_4\; -\; x\_0\}cdot\{x\; -\; x\_1\; over\; x\_4\; -\; x\_1\}cdot\{x\; -\; x\_2\; over\; x\_4\; -\; x\_2\}cdot\{x\; -\; x\_3\; over\; x\_4\; -\; x\_3\}$
             ={1over 243} x (2x+3)(4x-3)(4x+3)
Thus the interpolating polynomial then is

&-8f(x_1)x (2x-3)(2x+3)(4x-3) 
&+3f(x_2)(2x+3)(4x+3)(4x-3)(2x-3) 
&-8f(x_3)x (2x-3)(2x+3)(4x+3) 
&+f(x_4)x (2x+3)(4x-3)(4x+3)Big)
&= 4.834848x^3 - 1.477474x.
end{align} 
 Example 2 
We wish to interpolate $f(x)=\; x^2$ over the range $1\; leq\; x\; leq\; 3$:


 $x\_0=1,$  
  $f(x\_0)=1,$ 

 $x\_1=2,$ 
  $f(x\_1)=4,$  

 $x\_2=3,$     
  $f(x\_2)=9,$ 
The interpolating polynomial is:
$begin\{align\}$
L(x) &= {1}cdot{x - 2 over 1 - 2}cdot{x - 3 over 1 - 3}+{4}cdot{x - 1 over 2 - 1}cdot{x - 3 over 2 - 3}+{9}cdot{x - 1 over 3 - 1}cdot{x - 2 over 3 - 2} 
&= x^2. 
end{align} 
 Example 3 
We wish to interpolate $f(x)=\; x^3$ over the range $1\; leq\; x\; leq\; 3$:


 $x\_0=1,$  
  $f(x\_0)=1,$ 

 $x\_1=2,$ 
  $f(x\_1)=8,$  

 $x\_2=3,$     
  $f(x\_2)=27,$ 
The interpolating polynomial is:
$begin\{align\}$
L(x) &= {1}cdot{x - 2 over 1 - 2}cdot{x - 3 over 1 - 3}+{8}cdot{x - 1 over 2 - 1}cdot{x - 3 over 2 - 3}+{27}cdot{x - 1 over 3 - 1}cdot{x - 2 over 3 - 2} 
&=  6x^2 - 11x + 6. 
end{align} 
 Notes 
The Lagrange form of the interpolation polynomial shows the linear character of polynomial interpolation and the uniqueness of the interpolation polynomial. Therefore, it is preferred in proofs and theoretical arguments. But, as can be seen from the construction, each time a node xk changes, all Lagrange basis polynomials have to be recalculated. A better form of the interpolation polynomial for practical (or computational) purposes is the the barycentric form of the Lagrange interpolation (see below) or Newton polynomials. 
Lagrange and other interpolation at equally spaced points, as in the example above, yield a polynomial oscillating above and below the true function. This behaviour is known as the Gibbs phenomenon and can be lessened by choosing interpolation points at Chebyshev nodes. 
The Lagrange basis polynomials can be used in numerical integration to derive the Newton–Cotes formulas.

Barycentric Interpolation
Using the quantity
$ell(x)\; =\; (x\; -\; x\_0)(x\; -\; x\_1)\; dots\; (x\; -\; x\_k)$
we can re-write the Lagrange basis polynomials as
$ell\_j(x)\; =\; frac\{ell(x)\}\{x-x\_j\}\; frac\{1\}\{prod\_\{i=0,i\; neq\; j\}^k(x\_j-x\_i)\}$
or, by defining the barycentric weights
$w\_j\; =\; frac\{1\}\{prod\_\{i=0,i\; neq\; j\}^k(x\_j-x\_i)\}$
we can simply write
$ell\_j(x)\; =\; ell(x)frac\{w\_j\}\{x-x\_j\}$
which is commonly referred to as the first form of the barycentric interpolation formula.
The advantage of this representation is that the interpolation polynomial may now be evaluated as
$L(x)\; =\; ell(x)\; sum\_\{j=0\}^k\; frac\{w\_j\}\{x-x\_j\}f\_j$
which, if the weights $w\_j$ have been pre-computed, requires only $mathcal\; O(n)$ operations (evaluating $ell(x)$ and the weights $w\_j/(x-x\_j)$) as opposed to $mathcal\; O(n^2)$ for evaluating the Lagrange basis polynomials $ell\_j(x)$ individually.
The barycentric interpolation formula can also easily be updated to incorporate a new node $x\_\{k+1\}$ by dividing each of the $w\_j$, $j=0\; dots\; k$ by $(x\_j\; -\; x\_\{k+1\})$ and constructing the new $w\_\{k+1\}$ as above.
We can further simplify the first form by first considering the barycentric interpolation of the constant function $g(x)equiv\; 1$:
$g(x)\; =\; ell(x)\; sum\_\{j=0\}^k\; frac\{w\_j\}\{x-x\_j\}$.
Dividing $L(x)$ by $g(x)$ does not modify the interpolation, yet yields
$L(x)\; =\; frac\{sum\_\{j=0\}^k\; frac\{w\_j\}\{x-x\_j\}f\_j\}\{sum\_\{j=0\}^k\; frac\{w\_j\}\{x-x\_j\}\}$
which is referred to as the second form or true form of the barycentric interpolation formula. This second form has the advantage, that $ell(x)$ need not be evaluated for each evaluation of $L(x)$.
See also
Polynomial interpolation
Newton form of the interpolation polynomial
Bernstein form of the interpolation polynomial
Newton–Cotes formulas
External links

  Lagrange Method of Interpolation - Notes, PPT, Mathcad, Mathematica, Matlab, Maple at  Holistic Numerical Methods Institute
 Lagrange interpolation polynomial on www.math-linux.com
 Module for Lagrange Polynomials by John H. Mathews
  The chebfun Project at Oxford University
References