Given {n+1} distinct real numbers {x_j}, and {n+1} real numbers {y_j}, we wish to find a polynomial {p_n} (where {n} is the degree), such that {p_n(x_j)=y_j, \;j=0,\dots,n}. The polynomial must be of at most degree {n} since we have {n+1} constraints. The real numbers {y_j} could be function values, or table values. Given a function {f}, we can write the constraints as, {p_n(x_j)=f(x_j),\; j=0,\dots,n}.

The most intuitive approach is to start with a general {n}th degree polynomial of the form,


\displaystyle p_n(x) = a_0 + a_1x + a_2x^2 + \cdots + a_nx^n.

which is a linear combination of the basis functions {1,x,x^2,\dots,x_n}. They form a basis for the space of all polynomials of degree {n} or less, which we denote as {\mathcal{P}_n}. Because we are using a monomial basis of {\mathcal{P}_n} to construct the interpolating polynomial, we say that the polynomial is constructed in the monomial basis.

Explicitly writing out all constraints gives the following system of equations,


\displaystyle \begin{array}{rcl} p_n(x_0) &=& a_0 + a_1x_0 + a_2x^2_0 + \cdots + a_nx^n_0 = f(x_0) \\ p_n(x_1) &=& a_0 + a_1x_1 + a_2x^2_1 + \cdots + a_nx^n_1 = f(x_0) \\ &\vdots& \\ p_n(x_n) &=& a_0 + a_1x_n + a_2x^2_n + \cdots + a_nx^n_n = f(x_n), \end{array}

which we can write in matrix notation as,


\displaystyle \left( \begin{matrix} 1 & x_0 & x^2_0 & \cdots & x^n_0 \\ 1 & x_1 & x^2_1 & \cdots & x^n_1 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x^2_n & \cdots & x^n_n \end{matrix} \right) \left(\begin{matrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{matrix}\right) = \left(\begin{matrix} f(x_0) \\ f(x_1) \\ \vdots \\ f(x_n) \end{matrix}\right).

We can solve this system using the QR-decomposition, Gaussian elimination or even specialised algorithms that exploit the Vandermonde structure.

In one of his lecture notes, Mark Embree and the University of Rice poses the following six questions that numerical analysts seek answers too;

1. Does such a polynomial {p_n\in\mathcal{P}_n} exist?
2. If so, is it unique?
3. Does {p_n\in\mathcal{P}_n} behave like {f\in C[a,b]} at points {x\in[a,b]} when {x\neq x_j} for {j=0,\dots,n}?
4. How can we compute {p_n\in\mathcal{P}_n} efficiently on a computer?
5. How can we compute {p_n\in\mathcal{P}_n} accurately on a computer (with floating point arithmetic)?
6. How should the interpolation points {\{x_j\}} be chosen?

We could answer questions 1 and 2 using properties of Linear Algebra, but to avoid that we will leave them open for now. The other questions deserve their own posts, and we will deal with them eventually.