Barycentric Interpolation

A numerically stable and fast reformulation of polynomial interpolation through a set of points.

The trouble with Lagrange interpolation

The classical Lagrange form of the interpolating polynomial is elegant on paper but poor in practice: evaluating it directly costs order N-squared per point, and adding a new data point requires redoing everything. The barycentric formula rewrites the same polynomial in a form that is both faster to evaluate and numerically stable.

The barycentric formula

Each interpolation node is assigned a fixed weight that depends only on the node locations, not on the data values. The interpolant is then a ratio: the sum over nodes of the weight divided by (x minus node) times the data value, all divided by the sum of the weight divided by (x minus node). Once the weights are precomputed (a one-time order N-squared step), each evaluation costs only order N, and changing the data values requires no recomputation of weights.

import numpy as np

def barycentric(xnodes, w, f, x):
    num = np.zeros_like(x, dtype=float)
    den = np.zeros_like(x, dtype=float)
    for xj, wj, fj in zip(xnodes, w, f):
        d = x - xj
        # exact-node hits handled separately in production code
        num += wj / d * fj
        den += wj / d
    return num / den

Choosing the nodes and weights

For Chebyshev points the barycentric weights have a simple closed form (alternating plus-minus one, with halved endpoints), making Chebyshev-barycentric interpolation both trivial to set up and provably stable. This is the backbone of the Chebfun approach, which represents smooth functions to machine precision by such interpolants. Equally spaced nodes have weights that grow exponentially, reflecting the underlying ill-conditioning of that node set.

Why it matters

Barycentric interpolation is the recommended way to evaluate polynomial interpolants in scientific computing: it is fast, stable, and the foundation of Chebyshev spectral methods and high-order function approximation. It also generalizes to rational interpolation and to constructing quadrature and differentiation rules from the same nodes and weights.