Newton's divided difference interpolation

You fitted a curve through your data points. Then another measurement came in. With most methods you now throw the whole thing away and start again. Newton's version lets you keep everything you already did and add one term.

The curve you end up with is identical to the one Lagrange interpolation gives you, because only one curve of that degree passes through a given set of points. What differs is how you get there.

The quick answer

Build a table of slopes. Then slopes of the slopes. Then slopes of those. Take the top number from each column and you have your polynomial, written as a running sum where each term corrects the one before it.

What a divided difference is

The first level is just your y value:

f[x_i] = y_i

The second level is an ordinary slope between two neighbouring points, rise over run:

f[x_i, x_i+1] = (y_i+1 − y_i) ÷ (x_i+1 − x_i)

Every level after that follows the same shape: take the difference between two entries in the level below, and divide by the total x distance those entries span.

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)

So the second level measures how fast y is changing, the third measures how fast that rate is changing, and so on. The polynomial uses the top entry from each column:

P(x) = f[x₀] + f[x₀,x₁](x − x₀) + f[x₀,x₁,x₂](x − x₀)(x − x₁) + …

Read that as a series of corrections. The first term is your starting value. The second adds a straight line through the first two points. The third bends that line to catch the third point, and each new term fixes up the curve to pass through one more point without disturbing the ones already handled.

Worked example

Three points taken from y = x³: (1, 1), (2, 8) and (4, 64).

Build the table

First level, the ordinary slopes:

f[x₀,x₁] = (8 − 1) ÷ (2 − 1) = 7
f[x₁,x₂] = (64 − 8) ÷ (4 − 2) = 56 ÷ 2 = 28

The curve is getting steeper, 7 then 28, which fits data taken from a cube. Second level, the difference between those slopes divided by the full span from x = 1 to x = 4:

f[x₀,x₁,x₂] = (28 − 7) ÷ (4 − 1) = 21 ÷ 3 = 7

The numbers you need are the ones along the top row: 1, 7 and 7.

Write it out and use it

P(x) = 1 + 7(x − 1) + 7(x − 1)(x − 2)

At x = 3:

P(3) = 1 + 7(2) + 7(2)(1) = 1 + 14 + 14 = 29

That matches what Lagrange gives on the same three points, as it has to. The true value of 3³ is 27, so a quadratic through these points is out by 2.

Multiplying out confirms it is the same polynomial:

1 + 7x − 7 + 7(x² − 3x + 2) = 7x² − 14x + 8

Now add a fourth point

A new measurement arrives: (5, 125). With Lagrange you would rebuild all four helper functions from scratch. Here you just extend the table down one diagonal.

f[x₂,x₃] = (125 − 64) ÷ (5 − 4) = 61
f[x₁,x₂,x₃] = (61 − 28) ÷ (5 − 2) = 33 ÷ 3 = 11
f[x₀,x₁,x₂,x₃] = (11 − 7) ÷ (5 − 1) = 4 ÷ 4 = 1

Every term you already had stays exactly as it was. You append one:

P(x) = 1 + 7(x − 1) + 7(x − 1)(x − 2) + 1(x − 1)(x − 2)(x − 4)

At x = 3 the new term contributes (2)(1)(−1) = −2, so:

P(3) = 29 − 2 = 27

Exactly right, and not by luck. The data came from x³, four points pin down a cubic completely, so the curve is now the original function. Notice the last divided difference came out to 1, which is the coefficient of x³ in the function the data came from. That happens in general: the top divided difference tells you the leading coefficient of the polynomial that fits.

Reading the table tells you things

The size of the numbers in each column is useful information about your data, not just working.

Columns shrinking quickly toward zero means a low degree curve already describes your data well. Adding more terms will not buy you anything.

Columns staying large or growing means a polynomial is a poor model for what you have. Adding degree makes it oscillate rather than improve, which is Runge's phenomenon showing up in the arithmetic before you ever plot anything.

One unusually large entry high up the table usually means a bad data point rather than real curvature. People used divided difference tables to find typos in astronomical data long before computers existed.

Practical notes

The order you list your points does not change the resulting curve at all. It does change how accurate the arithmetic is. Putting the points nearest your target first keeps the running total close to the final answer the whole way, which loses less precision.

Evaluate with nested brackets rather than expanding everything out. It takes fewer multiplications and is more accurate.

If your points are evenly spaced, all the divisions become the same number repeated, which collapses into the forward and backward difference formulas from older textbooks. Same method, fewer keystrokes.

Repeated x values are allowed, which is unusual. Letting two points slide together turns the divided difference into a derivative, and that is how Hermite interpolation matches both the value and the slope at a point.

As with any polynomial method, watch the degree. Past five or six points the fit gets worse between your data rather than better, and a cubic spline is the right tool. For small sets, check your table against the interpolation calculator.