JOVANA
Explore Library Glossary Getting Started Three Levels Fields How it works Mission
Join the mission
All guides

Polynomial Interpolation: Lagrange and Newton

You have a handful of data points and you need a curve that passes exactly through every one of them. The cleanest answer is a single polynomial — and there is exactly one. This guide shows two beautiful, equivalent ways to write it down: Lagrange's direct recipe and Newton's grow-it-one-point-at-a-time form.

The problem: a curve through every dot

Suppose a weather station logs the temperature once an hour, and you need a sensible estimate for 2:37 in the afternoon — a moment nobody measured. Or a table of physical constants gives a material's density at 0, 100, and 200 degrees, and you want the value at 150. The recurring task is the same: you are handed a few points (x_0, y_0), (x_1, y_1), ..., (x_n, y_n) and asked for a smooth curve that passes EXACTLY through all of them, so you can read off values in between. That is interpolation, and the simplest curve to use is a polynomial — a sum of powers a_0 + a_1 x + a_2 x^2 + .... Polynomials are wonderfully convenient: trivial to evaluate, differentiate, and integrate, which is why polynomial interpolation sits at the foundation of so much numerical work — quadrature and ODE solvers are quietly built on it.

Before we hunt for the curve, one fact makes the whole subject solid: it is UNIQUE. Given n + 1 points with DISTINCT x-values, there is one and only one polynomial of degree at most n that threads all of them. Two points fix a straight line (degree 1), three fix a parabola (degree 2), four fix a cubic, and so on — n + 1 conditions pin down the n + 1 coefficients of a degree-n polynomial. This is not several competing answers to choose between; it is a single object we are about to write down in three different handwritings. The catch in the fine print — 'degree AT MOST n' — matters: three collinear points give back a straight line, not a forced parabola.

The honest-but-treacherous way: solve for coefficients

The most literal idea is to find the coefficients directly. Write p(x) = a_0 + a_1 x + ... + a_n x^n and demand p(x_i) = y_i at each of the n + 1 data points. That gives n + 1 linear equations in the n + 1 unknown coefficients — a system A a = y you can solve with the Gaussian elimination from earlier in this ladder. The matrix here, whose row i is (1, x_i, x_i^2, ..., x_i^n), is the famous Vandermonde matrix. Its determinant is the product of (x_j - x_i) over all i < j, which is nonzero exactly when the nodes are distinct — so this viewpoint gives a clean linear-algebra PROOF of the uniqueness we just claimed.

Conceptually lovely, computationally dreadful. The Vandermonde matrix becomes severely ILL-CONDITIONED as the degree grows: for equispaced nodes its condition number explodes roughly like 2^n. Recall the rule of thumb from the conditioning rung — a condition number near 10^8 already costs you about 8 of your roughly 16 double-precision digits. By degree 20 or so on an equispaced grid the recovered coefficients are numerical garbage, even though the true interpolating polynomial is perfectly fine. This is the difference between the PROBLEM (interpolation, which is well behaved) and a bad ALGORITHM for it (forming and solving A a = y). So we keep the Vandermonde matrix as an explanation and reach for better recipes to actually compute — exactly the spirit of preferring a backward-stable method over a naive one.

Lagrange's recipe: blend the y-values directly

Here is a way to write the answer with no system to solve at all. Lagrange interpolation builds the polynomial as a weighted blend of the y-values, where each weight is a little custom-made polynomial. For the i-th node you cook up a basis polynomial L_i(x) designed to equal 1 at x_i and 0 at every OTHER node. The construction is direct: multiply together the factors (x - x_j) for all j not equal to i (this kills it at every other node), then divide by that same product evaluated at x_i (this normalizes it to 1 at x_i). Then the whole interpolant is simply p(x) = y_0 L_0(x) + y_1 L_1(x) + ... + y_n L_n(x).

Why does this thread every point? Plug in any node x_k. Every basis polynomial L_i with i not equal to k was built to vanish at x_k, so all those terms disappear; only L_k survives and it equals 1, leaving p(x_k) = y_k. It hits every data point automatically, by construction — no equations, no elimination, no Vandermonde matrix. With two points this just reproduces the straight line through them; with three, the parabola. It is the same unique polynomial as before, written in a 'partition of the work among the y-values' handwriting.

Through (0,1), (1,3), (2,2):

        (x-1)(x-2)              (x-0)(x-2)             (x-0)(x-1)
L_0 = --------------    L_1 = --------------   L_2 = --------------
      (0-1)(0-2)              (1-0)(1-2)             (2-0)(2-1)

p(x) = 1*L_0(x) + 3*L_1(x) + 2*L_2(x)

Check x = 1:  L_0(1)=0,  L_1(1)=1,  L_2(1)=0  ->  p(1) = 3
Each basis polynomial is 1 at its own node and 0 at all the others, so the y-values pass through untouched.

Lagrange's form is beautiful for theory and proofs — the answer is written out explicitly with no linear algebra in sight. Its honest drawbacks are practical. In this naive form, evaluating p(x) at one query point costs O(n^2) work, and worse, ADDING a new data point forces you to rebuild every single basis polynomial from scratch, because the (x - x_j) products all change. Newton's form fixes the second complaint, and a clever rewrite called the barycentric formula fixes both at once — but that polished, fast, stable version is the topic of the Chebyshev guide later in this rung.

Newton's form: grow the curve one point at a time

Now imagine you want to build the curve incrementally, so that each new data point just ADDS a correction without disturbing the fit you already have — like appending a finer term to a series. The Newton divided-difference form does exactly this. It writes the SAME unique polynomial as a nested sum: p(x) = c_0 + c_1 (x - x_0) + c_2 (x - x_0)(x - x_1) + c_3 (x - x_0)(x - x_1)(x - x_2) + .... Look at the structure: each new factor (x - x_j) is zero at all the EARLIER nodes, so tacking on a fresh term cannot break the points the earlier terms already nailed. That is exactly why adding a data point only appends one more coefficient and one more factor — cheap, where Lagrange demanded a full rebuild.

The coefficients c_k are the divided differences, computed in a small triangular table. The zeroth differences are just the y-values: f[x_i] = y_i. The first differences are slopes between neighbours: f[x_i, x_{i+1}] = (f[x_{i+1}] - f[x_i]) / (x_{i+1} - x_i). Higher ones recurse — each is a difference of two lower divided differences, divided by the span of x-values they cover: 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 edge of the table — f[x_0], f[x_0,x_1], f[x_0,x_1,x_2], ... — reads off as c_0, c_1, c_2, ... directly.

  1. Write the y-values down the first column of the table: f[x_0] = y_0, f[x_1] = y_1, ..., f[x_n] = y_n.
  2. Form the first divided differences from each adjacent pair: f[x_i, x_{i+1}] = (f[x_{i+1}] - f[x_i]) / (x_{i+1} - x_i).
  3. Recurse up the orders: each higher difference is (the one below-right minus the one below-left) divided by (x_{i+k} - x_i), filling shorter and shorter columns.
  4. Read the coefficients c_0, c_1, ..., c_n straight off the top diagonal of the table.
  5. Evaluate p(x) by nested (Horner-style) multiplication, working from the deepest term outward — O(n) work per query point.

A tiny worked example with the same three points (0, 1), (1, 3), (2, 2): the y-values give f[0] = 1, f[1] = 3, f[2] = 2. First differences: f[0,1] = (3-1)/1 = 2 and f[1,2] = (2-3)/1 = -1. Second difference: f[0,1,2] = (-1 - 2)/(2 - 0) = -1.5. So p(x) = 1 + 2(x - 0) - 1.5(x - 0)(x - 1). Expand it and you get exactly x^2 + ... — the very same parabola the Lagrange form produced, and the same one the Vandermonde solve would have given. One curve, three handwritings. Building the table is a one-time O(n^2) cost; every evaluation afterward is O(n).

Same curve, three handwritings — and one honest warning

It bears repeating because it surprises people: Vandermonde, Lagrange, and Newton do NOT give three different interpolants. By the uniqueness theorem they cannot — there is only one polynomial of degree at most n through n + 1 distinct points, so all three must produce it. What differs is only the bookkeeping: the cost to set up, the cost to evaluate, the numerical stability, and how gracefully you can add a new point. Vandermonde is the clearest explanation but the worst computation; Lagrange is the cleanest for theory; Newton is the practical choice when data arrives incrementally or when you want the error estimate that the next guide builds on.

There is a deeper bridge worth noticing. The leading divided difference is tied to a derivative: f[x_0, ..., x_n] = f^(n)(xi) / n! for some point xi inside the data span. That is no accident — it is the seed of the interpolation error formula, where the next neglected divided difference behaves exactly like the gap between the polynomial and the true function. Newton's form, in other words, carries its own error estimate in its very next coefficient, which is part of why it remains the workhorse for adaptive schemes.