Section 7.2 Polynomial interpolation
The tipical problem of interpolation is the following: given \(n+1\) points in the plane
\begin{equation*}
(x_0,y_0), (x_1,y_1), \dots, (x_n,y_n),
\end{equation*}
with all distinct \(x\) coordinates, find some "simple" function that passes through all of them. This interpolating function can then be used for several purposes: approximating integrals, approximating derivatives, finding minima and maxima, approximating values between points and so on. Of course the simplest type of functions, the only kind we discuss in this chapter, are polynomial ones. Recall that a polynomial of degree \(n\) has \(n+1\) coefficients, so \(n\) is the lowest degree for which we can find a solution for a generic choice of points. Moreover, for each choice of distinct \(n+1\) points, there is a unique polynomial of degree not larger than \(n\) passing through them. Denote by \(P_{n}\) the linear space of all polynomials of degree at most \(n\text{.}\) Recall that this is a linear space of dimension \(n+1\text{.}\) Then, to find the interpolating polynomial relative to the points \((x_k,y_k)\text{,}\) it is enough choosing a basis
\begin{equation*}
p_0(x),\dots,p_n(x)
\end{equation*}
of \(P_{n}\text{,}\) write
\begin{equation*}
p(x)=a_0p_0(x)+\dots+a_np_n(x)
\end{equation*}
and solve the system
\begin{equation*}
\left\{\begin{array}{@{}l@{}}
a_0p_0(x_0)+\dots+a_np_n(x_0)=y_0\cr
\dots\cr
a_0p_0(x_{n})+\dots+a_np_n(x_{n})=y_{n}\cr
\end{array}\right.
\end{equation*}
Subsection 7.2.1 Case 1: \(p_k(x)=x^k\)
In this case the system above writes
\begin{equation*}
\begin{pmatrix}
1&x_0&x_0^2&\dots&x_0^n\\
\vdots&\vdots&\vdots&\vdots&\vdots\\
1&x_{n}&x_{n}^2&\dots&x_{n}^n\\
\end{pmatrix}
\begin{pmatrix}
a_0\\
\vdots\\
a_n\\
\end{pmatrix}
=
\begin{pmatrix}
y_0\\
\vdots\\
y_{n}\\
\end{pmatrix}
\end{equation*}
This matrix is known as Vandermonde matrix and it is ill-conditioned:
The reason behind this behavior is that monomials of high degree look more and more just like each other: Subsection 7.2.2 Case 2: Lagrange basis
The Lagrange basis is at the opposite end of the spectrum: the corresponding matrix is the identity matrix! Lagrange polynomials are defined as follows:
\begin{equation*}
L_k(x) = \prod_{j\neq k}\frac{x-x_j}{x_k-x_j},
\end{equation*}
so that
\begin{equation*}
L_k(x_j)=\delta_{kj}=\left\{\begin{array}{@{}l@{}}1,\hbox{ if }k=j;\cr 0,\hbox{ otherwise}.\cr\end{array}\right.
\end{equation*}
Hence, the system above now becomes simply
\begin{equation*}
\left\{\begin{array}{@{}l@{}}
a_0=y_0\cr
\dots\cr
a_n=y_{n}\cr
\end{array}\right.
\end{equation*}
and so
\begin{equation*}
p(x) = y_0 L_0(x) + \dots + y_n L_n(x).
\end{equation*}
In the code below we implement a MATLAB function defining Lagrange polynomials and use it to plot them in a concrete case. The plot clearly shows that each Lagrange polyomial is quite distinct from all other ones, quite the opposite behavior of the monomials \(x^k\text{.}\)