Lagrange Interpolation

The Lagrange form writes the unique interpolating polynomial as a weighted sum of basis polynomials, each one at its own node and zero at the others.

Basis polynomials that select one point

The Lagrange form builds the interpolating polynomial directly from the data. For nodes x_0..x_n it defines basis polynomials L_i(x) that equal 1 at x_i and 0 at every other node. The interpolant is then simply the sum of y_i L_i(x). No linear system needs to be solved.

def lagrange(xs, ys, x):
    n = len(xs); total = 0.0
    for i in range(n):
        term = ys[i]
        for j in range(n):
            if j != i:
                term *= (x - xs[j])/(xs[i] - xs[j])
        total += term
    return total

Elegant but costly to update

The Lagrange form is conceptually clean and excellent for theory and derivations, but it has drawbacks in practice. Adding a new data point requires recomputing every basis polynomial, and naive evaluation costs O(n^2) per point. The Newton divided-difference form is preferred when points are added incrementally.

The barycentric improvement

The barycentric Lagrange formula rewrites the same polynomial using precomputed weights, reducing evaluation to O(n) per point after an O(n^2) setup. It is numerically stable and is the recommended way to evaluate polynomial interpolants in modern practice, especially with Chebyshev nodes.

Error behavior

The interpolation error at x depends on the (n+1)th derivative of the underlying function and on the product of distances to all nodes. This product is what grows near interval ends for equally spaced nodes, producing Runge's phenomenon. Clustering nodes toward the endpoints controls the error.

Lagrange interpolation underpins many quadrature rules and finite-element shape functions used in physics codes, including the element formulations behind breeder Hyperion field solvers.