Skip to main content

Section 21.1 Interpolation in the monomial basis

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.

Consider the polynomial \(p(x) = (x-7)^{15}\text{.}\) Expand this formula in the monomial basis and try using it to plot the polynomial on the interval \([6, 8]\text{.}\) Also plot the original formula for comparison.

Answer
In the code below, c is a vector of coefficients; the ellipsis notation allows it to by split between lines for readability (without …, a linebreak would be understood as the beginning of a new row of a matrix).
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)

There is a catastrophic loss of significance here. If one uses p = polyval(c, x), the result is not much better, even though this command tries to avoid the loss of significance. Expanding this polynomial in the monomial basis is numerically unwise.