Lagrange Interpolation

The Interpolation Problem¶

Problem: Given $n+1$ distinct points $(x_0, y_0), (x_1, y_1), \ldots, (x_n, y_n)$, find a polynomial $p(x)$ such that:

Lagrange Basis Polynomials¶

For nodes $x_0, x_1, \ldots, x_n$, define the Lagrange basis polynomial:

Key property:

The interpolating polynomial is then:

Barycentric Interpolation¶

Can we do better? Yes—by algebraically reorganizing the Lagrange formula, we obtain the barycentric form, which resolves all three issues while computing the same polynomial.

The Node Polynomial¶

Define the node polynomial:

Its derivative at a node $x_j$ is:

This lets us rewrite the Lagrange basis as:

Barycentric Weights¶

Define the barycentric weights:

Then the Lagrange basis becomes:

First Barycentric Formula¶

Substituting into the interpolant gives the first barycentric formula:

Once the weights $\lambda_j$ are known (computed once in $O(n^2)$ time), evaluating $p(x)$ costs only $O(n)$.

Second Barycentric Formula¶

An even more elegant formula comes from the identity $\sum_{j=0}^{n} L_j(x) = 1$ (the Lagrange polynomials form a partition of unity—they interpolate the constant function 1).

Dividing the interpolant by this identity:

Advantages¶

1. $O(n)$ per evaluation after $O(n^2)$ preprocessing for weights

2. Numerically stable even for large $n$

3. Adding a point only requires updating weights

Implementation¶

def bary_weights(x):
    """Compute barycentric weights for nodes x."""
    n = len(x)
    w = np.ones(n)
    for j in range(n):
        for i in range(n):
            if i != j:
                w[j] /= (x[j] - x[i])
    return w

def bary_interp(xeval, x, y, w):
    """Evaluate interpolant at xeval using barycentric formula."""
    # Handle evaluation at nodes
    for j, xj in enumerate(x):
        if np.isclose(xeval, xj):
            return y[j]

    terms = w / (xeval - x)
    return np.dot(terms, y) / np.sum(terms)

Chebyshev Points and Barycentric Weights¶

For Chebyshev points, the barycentric weights have a remarkably simple closed form:

This makes Chebyshev interpolation especially efficient—no $O(n^2)$ weight computation needed!

def cheb_bary_weights(n):
    """Barycentric weights for Chebyshev points."""
    w = np.ones(n+1)
    w[0] = 0.5
    w[-1] = 0.5
    w[1::2] *= -1
    return w

Numerical Stability¶

The two barycentric formulas have different stability properties. The following results are due to Higham (2004).

Backward Stability of the First Formula¶

Forward Stability of the Second Formula¶

Comparison¶

Why Scale Invariance Matters¶

The second barycentric formula is scale invariant: we can multiply all weights by any nonzero constant without changing the result:

This means we can rescale weights to avoid overflow, which is critical for large $n$.

Warning: Equispaced Points¶

For equispaced points, the weights grow like:

This grows exponentially, making polynomial interpolation through equispaced points numerically unstable—even with the barycentric formula. This is separate from (but compounds) Runge’s phenomenon.

Runge’s Phenomenon: A Warning¶

Consider $f(x) = \frac{1}{1 + 25x^2}$ on $[-1, 1]$.

# Equally spaced: error GROWS
x_eq = np.linspace(-1, 1, n)
f_eq = 1/(1 + 25*x_eq**2)
# Interpolant oscillates wildly near boundaries!

# Chebyshev: error DECREASES
x_cheb = np.cos(np.pi * np.arange(n) / (n-1))
f_cheb = 1/(1 + 25*x_cheb**2)
# Smooth, accurate approximation

The lesson: node placement matters. Chebyshev nodes cluster near the endpoints, exactly where equally spaced nodes cause trouble.

Use barycentric interpolation. It’s the standard in modern software (Chebfun, etc.).

References¶