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.
Example 21.1.1. Loss of significance in the monomial basis.
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.
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.