Approximation With Interpolating Polynomials

Lagrange Form

Say you want to find an function approximating a set of data points$(x_0, f_0), (x_1, f_1), \ldots, (x_n, f_n)$. One way to do this is to find a polynomial $p_n \in \mathbb{P}_n$ satisfying the interpolation property

$p_n(x_i) = f_i\qquad i \in [0, n]$

Playing around a bit, we find that this polynomial can be written as

$\begin{align*} p_n(x) = &\frac{(x - x_1) * (x - x_2) * \ldots * (x - x_n)}{(x_0 - x_1) * (x_0 - x_2) * \ldots * (x_0 - x_n)} f_0 \\ &+ \frac{(x - x_0) * (x - x_2) * \ldots * (x - x_n)}{(x_1 - x_0) * (x_1 - x_2) * \ldots * (x_1 - x_n)} f_1 \\ &+ \ldots \\ &+ \frac{(x - x_0) * (x - x_1) * \ldots * (x - x_{n-1})}{(x_n - x_0) * (x_n - x_1) * \ldots * (x_n - x_{n-1})} f_n \\ = &\sum_{i=0}^n \prod_{j=0\\j\ne i}^n \frac{x-x_j}{x_i-x_j} f_i\qquad j \in [0, n] \end{align*}$

This is known as the Lagrange form of an interpolating polynomial, and is more commonly expressed as

$\begin{align*} p_n(x) &= \sum_{i=0}^n \ell_i(x)f_i \\ \ell_i(x) &= \prod_{j=0\\j\ne i}^n \frac{x-x_j}{x_i-x_j}\qquad j \in [0, n] \end{align*}$

where $\ell_i$ is a Lagrange polynomial of degree $n$.⊕

To see that $p_n$ interpolates correctly, observe that

$\ell_i(x) = \begin{cases} 1 & x = x_i\\ 0 & x \in \left\{x_0,\ldots,x_n\right\}\setminus \left\{x_i\right\} \end{cases}$

and so $p_n(x_i) = f_i + \sum_{j=0\j\ne i}^n 0 f_j = f_i$.

It's not difficult to prove that $p_n$ is unique and of at most degree $n$, but we won't do that here.

Cost

The Lagrange form is not suitable for computation. Its construction requires massive multiplication and division, which can be catastrophic for floating point precision. The form is also not memoizable; evaluating $p_n(a)$ does not help in any part of the evaluation of $p_n(b)$ when $b \ne a$.

With some more play, we can get a computationally-suitable equation known as the Barycentric form of an interpolating polynomial, defined as

$p_n(x) = \frac{\sum_{i=0}^n \frac{\omega_i}{x-x_i} * f_i}{\sum_{i=0}^n \frac{\omega_i}{x-x_i}}$

where $w_i$ is a barycentric weight

$\omega_i(x) = \frac{1}{\prod_{j=0\\j\ne i}^n (x_i - x_j)}$

Because $\omega_i$ relies only on the explicitly-defined nodes $x_0, \ldots, x_n$, it does not need to be recomputed for subsequent evaluations of $p_n$.

Attempt 1

Approximating $f(x) = \sin(4x\cos(8x))$ with various degrees of the Barycentric interpolant, using equidistant interpolating nodes, gives

These polynomials definitely interpolate points of $f(x)$, but have the unfortunate tendency to blow up around the interval endpoints. Note that the magnitude of the error at the edges grows with the degree of the approximation.

Runge

This behavior is Runge's phenomenon, and means that increasing the order of an approximating polynomial interpolant may not improve its accuracyWhy this occurs is beyond the scope of this post; Wikipedia gives a good explanation..

One way to minimize Runge's phenomenon is to densify the interpolating points near the endpoints of the approximation intervalThereby more accurately fitting the parts of the approximation affected by Runge's phenomenon.. Chebyshev nodes exhibit exactly this kind of edge-heavy density, making them very useful for polynomial interpolation.

Creating a Barycentric interpolant with Chebyshev nodes yields a nice improvement.

Convergence

Plotting the pointwise error of the Barycentric interpolantgenerated with Chebyshev nodes approximating $f(x) = \sin(4x\cos(8x))$ suggests the polynomial interpolation is convergent.

Inspection of the error $l^2$ normsNote that many other norms can be used for convergence analysis as well. suggests more:

Order   L2 norm
16      8.562880126489002
32      0.4243100460236342
64      5.324902188796007e-05
128     8.698460604825334e-14

Clearly, the rate of convergence is pretty fast.