Numerical Methods with Programming

Given \(n\) data points \((x_k, y_k)\) with \(k=1, \dots, n\text{,}\) where all \(x_k\) are different, we are looking for a polynomial \(p\) of degree \(\lt n\) such that \(p(x_k) = y_k\) for \(k=1, \dots, n\text{.}\) Such interpolating polynomial exists and is unique. Existence will be shown below, and uniqueness is a consequence of the fact that a polynomial of degree less than \(n\) cannot have \(n\) distinct zeros. Outside of exceptional cases, the interpolating polynomial will have degree \(n-1\text{.}\) The main question is how to construct it.

A straightforward way is to do this is to write \(p(x) = \sum_{j=0}^{n-1} c_j x^j\) and solve a linear system for \(c_j\text{;}\) it has \(n\) equations with \(n\) unknowns. But this process is both slow and numerically fragile; the matrix of the linear system is likely to be ill-conditioned. And even if we could find the coefficients of monomials \(x_j\) easily, we do not necessarily want to, for the following reason.

We are used to seeing polynomials written in terms of the monomial basis \(x^0, x^1, \dots, x^d\text{,}\) that is, as a linear combination of the monomials. But this basis is ill-suited for numerical computations away from the point \(x=0\text{,}\) as the following example shows.

c = [1, -105, 5145, -156065, 3277365, -50471421, 588833245, -5299499205, ...
    37096494435, -201969803035, 848273172747, -2699051004195, 6297785676455, ...
    -10173346092735, 10173346092735, -4747561509943];
d = 15;
x = linspace(6, 8, 1000);
p = zeros(size(x));
for k=1:d+1
    p = p + c(k)*x.^(d+1-k);
end
subplot(1, 2, 1)
plot(x, p)
subplot(1, 2, 2)
plot(x, (x-7).^15)