Newton's Divided Differences

Newton's form builds the interpolating polynomial incrementally from divided-difference coefficients, so new data points extend the polynomial cheaply.

An incremental construction

Newton's form writes the interpolating polynomial as a_0 + a_1(x-x_0) + a_2(x-x_0)(x-x_1) + ..., where the coefficients a_k are divided differences of the data. Each coefficient depends only on points already used, so adding a new point appends one term without recomputing the rest.

The divided-difference table

Divided differences are defined recursively: the zeroth order is f[x_i] = y_i, and higher orders are f[x_i..x_{i+k}] = (f[x_{i+1}..x_{i+k}] - f[x_i..x_{i+k-1}])/(x_{i+k} - x_i). The leading diagonal of the table gives the Newton coefficients.

def divided_diff(xs, ys):
    n = len(xs); c = list(ys)
    for k in range(1, n):
        for i in range(n-1, k-1, -1):
            c[i] = (c[i]-c[i-1])/(xs[i]-xs[i-k])
    return c   # Newton coefficients

def newton_eval(xs, c, x):
    p = c[-1]
    for k in range(len(xs)-2, -1, -1):
        p = p*(x-xs[k]) + c[k]
    return p

Why it is preferred for updates

When data arrives one point at a time, or when the degree is chosen adaptively until accuracy is met, Newton's form is ideal. Its Horner-like evaluation is efficient and reasonably stable, and it produces the same unique polynomial as the Lagrange form, differing only in representation.

Relation to derivatives

Divided differences approximate derivatives: f[x_0,x_1] is a first-difference slope, and the kth divided difference times k! approaches the kth derivative as the nodes coalesce. This link makes Newton's form the natural bridge between interpolation and numerical differentiation.

Incremental interpolation is convenient when refining tabulated physics data adaptively, as done for material and reaction-rate tables in breeder Hyperion simulations.