Prev
Unit 7 · Interpolation and Approximation
Polynomial Interpolation
60m100 XP
Next
📈Interpolation and Approximation

Polynomial Interpolation

Existence and uniqueness, the Lagrange and Newton forms, the barycentric formula and the error term

Every 15 minutes, to a few centimetres
The International GNSS Service combines tracking data from hundreds of ground stations into precise orbits for every GPS, Galileo, GLONASS and BeiDou satellite. Its final product gives each satellite's position to about 2–3 cm — but only at tabulated epochs, every 15 minutes. A receiver doing precise positioning needs the satellite's position at the exact instant each signal left it — a time that never coincides with a table entry. In 15 minutes a GPS satellite travels about 3,500 km along a curved orbit. A straight line between two table entries misses the true position by tens of kilometres. Yet processing software routinely recovers the position between epochs to millimetres, from nothing but the tabulated numbers. The method is older than the satellites by two centuries: pass a polynomial through nearby table entries, and evaluate it in between.
💡Existence and uniqueness
Given n+1n+1 points (x0,y0),…,(xn,yn)(x_0, y_0), \dots, (x_n, y_n) with distinct xjx_j, there is exactly one polynomial pp of degree at most nn with p(xj)=yjp(x_j) = y_j for every jj. Existence: the construction below. Uniqueness: if two such polynomials existed, their difference would have degree at most nn and n+1n+1 roots — so it is zero. Whatever form we write pp in — monomials, Lagrange, Newton — it is the same polynomial; the forms differ only in cost and numerical stability.
p(x)=∑j=0nyj ℓj(x),ℓj(x)=∏m≠jx−xmxj−xmp(x) = \sum_{j=0}^{n} y_j\,\ell_j(x), \qquad \ell_j(x) = \prod_{m \ne j}\frac{x - x_m}{x_j - x_m}

The Lagrange form. Each basis polynomial ℓⱼ equals 1 at its own node and 0 at all the others, so p takes the value yⱼ at xⱼ by construction.

Newton's form builds the same polynomial one node at a time:
p(x)=c0+c1(x−x0)+c2(x−x0)(x−x1)+…p(x) = c_0 + c_1(x - x_0) + c_2(x - x_0)(x - x_1) + \dots
The coefficients are divided differences, ck=f[x0,…,xk]c_k = f[x_0, \dots, x_k], computed by a triangular table: f[xi,xi+1]=(yi+1−yi)/(xi+1−xi)f[x_i, x_{i+1}] = (y_{i+1} - y_i)/(x_{i+1} - x_i), and so on. Adding a new data point adds one coefficient and leaves the others unchanged, and evaluation by nested multiplication costs O(n)O(n). Interpolating ISA air density from a table at 0, 2, …, 10 km gives 0.73645 kg/m³ at 5 km against the tabulated 0.73612 — and adding the 12 km entry keeps the first six coefficients exactly.
Derivation
Why Newton's coefficients are divided differences
Impose the interpolation conditions one node at a time; each new node brings in exactly one new coefficient, because every later term of Newton's form vanishes at the earlier nodes.
Node x0x_0: every term after c0c_0 contains (x−x0)(x - x_0), so c0=y0c_0 = y_0.
Node x1x_1: y0+c1(x1−x0)=y1y_0 + c_1(x_1 - x_0) = y_1, so c1=y1−y0x1−x0=f[x0,x1]c_1 = \dfrac{y_1 - y_0}{x_1 - x_0} = f[x_0, x_1] — the slope of the chord.
Node x2x_2: y2=y0+f[x0,x1](x2−x0)+c2(x2−x0)(x2−x1)y_2 = y_0 + f[x_0, x_1](x_2 - x_0) + c_2(x_2 - x_0)(x_2 - x_1). Write y2−y0=f[x1,x2](x2−x1)+f[x0,x1](x1−x0)y_2 - y_0 = f[x_1, x_2](x_2 - x_1) + f[x_0, x_1](x_1 - x_0); then most terms cancel and
c2=f[x1,x2]−f[x0,x1]x2−x0=f[x0,x1,x2],c_2 = \frac{f[x_1, x_2] - f[x_0, x_1]}{x_2 - x_0} = f[x_0, x_1, x_2],
a difference of slopes divided by the width they span — a discrete second derivative (about f′′/2f''/2). The pattern continues:
f[xi,…,xi+k]=f[xi+1,…,xi+k]−f[xi,…,xi+k−1]xi+k−xi,f[x_i, \dots, x_{i+k}] = \frac{f[x_{i+1}, \dots, x_{i+k}] - f[x_i, \dots, x_{i+k-1}]}{x_{i+k} - x_i},
which is the triangular table. A new node xn+1x_{n+1} adds one new term and leaves c0,…,cnc_0, \dots, c_n alone, because they were fixed by the earlier nodes only.
💡The interpolation error
If ff has n+1n+1 continuous derivatives, then for each xx there is a ξ\xi in the span of the nodes and xx with f(x)−p(x)=f(n+1)(ξ)(n+1)!∏j=0n(x−xj).f(x) - p(x) = \dfrac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{j=0}^{n}(x - x_j). Two factors: how wiggly ff is (its high derivatives) and where the nodes are (the node polynomial). For a smooth orbit the derivatives are modest and (n+1)!(n+1)! wins, so the error falls fast with nn. Outside the span of the nodes the node polynomial grows like ∣x∣n+1|x|^{n+1}: extrapolation is dangerous, which is why ephemeris software always centres its window of epochs on the requested time.
Derivation
The error formula, by Rolle's theorem
Fix a point xx that is not a node, and write w(t)=∏j=0n(t−xj)w(t) = \prod_{j=0}^{n}(t - x_j) for the node polynomial. Choose the constant KK so that the function
g(t)=f(t)−p(t)−K w(t)g(t) = f(t) - p(t) - K\,w(t)
is zero at t=xt = x; this is possible because w(x)≠0w(x) \ne 0. Now count the zeros of gg: it vanishes at all n+1n + 1 nodes (where f=pf = p and w=0w = 0) and at xx — n+2n + 2 zeros.
Rolle's theorem: between two zeros of a smooth function its derivative has a zero. So g′g' has at least n+1n + 1 zeros, g′′g'' at least nn, …, and g(n+1)g^{(n+1)} at least one, at some point ξ\xi. Differentiate gg n+1n + 1 times: pp has degree nn, so p(n+1)=0p^{(n+1)} = 0; ww is tn+1+…t^{n+1} + \dots, so w(n+1)=(n+1)!w^{(n+1)} = (n+1)!. Hence
0=g(n+1)(ξ)=f(n+1)(ξ)−K (n+1)!⟹K=f(n+1)(ξ)(n+1)!.0 = g^{(n+1)}(\xi) = f^{(n+1)}(\xi) - K\,(n+1)! \qquad\Longrightarrow\qquad K = \frac{f^{(n+1)}(\xi)}{(n+1)!}.
Finally g(x)=0g(x) = 0 says f(x)−p(x)=K w(x)f(x) - p(x) = K\,w(x), which is the formula above. The same argument with a single node is Taylor's theorem with the Lagrange remainder (lesson 1-3).
📐Interpolating an orbit
A GPS-like orbit (aa = 26,560 km, ee = 0.01, ii = 55°) tabulated every 15 minutes, interpolated at 06:07:30 with a window of epochs centred on that time:
  • 2 epochs (linear): 56 km error;
  • 4 epochs (cubic): 170 m;
  • 6 epochs: 0.47 m;
  • 8 epochs: 0.19 mm;
  • 10 epochs (degree 9): 0.02 mm.
The pure Keplerian orbit is smoother than a real one, which carries perturbations; real processing uses degree 8–10 windows and reaches millimetre consistency, well below the 2–3 cm accuracy of the tabulated positions themselves.
Science is built up with facts, as a house is with stones. But a collection of facts is no more a science than a heap of stones is a house.
– Henri Poincaré, Science and Hypothesis (1902)
How to evaluate it. Do not solve the Vandermonde system Vc=yV\mathbf c = \mathbf y for monomial coefficients: VV becomes hopelessly ill-conditioned — κ≈1019\kappa \approx 10^{19} for 81 points on [−1,1][-1, 1] — and the result loses most of its digits. The barycentric formula
p(x)=∑jwjx−xjyj∑jwjx−xj,wj=1∏m≠j(xj−xm)p(x) = \dfrac{\sum_j \dfrac{w_j}{x - x_j}y_j}{\sum_j \dfrac{w_j}{x - x_j}}, \qquad w_j = \dfrac{1}{\prod_{m\ne j}(x_j - x_m)}
is the Lagrange form rearranged: O(n)O(n) per evaluation once the weights are known, and numerically stable. On 81 Chebyshev points it is accurate to 2.4×10−92.4 \times 10^{-9} while the Vandermonde route is off by 4×10−44 \times 10^{-4}. scipy.interpolate.BarycentricInterpolator implements it.
Derivation
The barycentric formula is Lagrange's, rearranged
Every Lagrange basis polynomial contains all the node factors except one. With ℓ(x)=∏m(x−xm)\ell(x) = \prod_m(x - x_m) and the weights wjw_j above,
ℓj(x)=∏m≠jx−xmxj−xm=ℓ(x) wjx−xj,sop(x)=ℓ(x)∑jwjx−xj yj.\ell_j(x) = \prod_{m \ne j}\frac{x - x_m}{x_j - x_m} = \ell(x)\,\frac{w_j}{x - x_j}, \qquad\text{so}\qquad p(x) = \ell(x)\sum_j\frac{w_j}{x - x_j}\,y_j.
Now interpolate the constant function 1. The interpolant of a constant is that constant, so with every yj=1y_j = 1:
1=ℓ(x)∑jwjx−xj.1 = \ell(x)\sum_j\frac{w_j}{x - x_j}.
Divide the first line by the second: ℓ(x)\ell(x) cancels and leaves the barycentric formula. The weights are computed once in O(n2)O(n^2); each evaluation then costs O(n)O(n) and never forms a large monomial coefficient. For the Chebyshev extreme points of lesson 7-2 the weights are even known in closed form — up to a common factor they are ±1\pm1 with alternating signs, halved at the two ends.
💡Where you will use this
GNSS Navigation: Lagrange interpolation of precise ephemerides and clock tables, exactly as in the lab. Rocket Propulsion and Aircraft Engines: thermochemical tables (enthalpy, specific heat against temperature). Flight software: aerodynamic coefficient tables interpolated at every control cycle. Numerical integration and ODE solvers (Units 8–9): Simpson's rule, Gauss quadrature and Adams methods are all built by integrating interpolating polynomials.
📚Further reading
  • Jean-Paul Berrut and Lloyd N. Trefethen, 'Barycentric Lagrange interpolation', SIAM Review 46(3), 501–517 (2004) — why the barycentric formula should replace every other way of evaluating an interpolant.
  • Michael Heath, Scientific Computing: An Introductory Survey (2nd ed.; SIAM Classics 2018), chapter 7, 'Interpolation'.
  • Endre Süli and David Mayers, An Introduction to Numerical Analysis (Cambridge University Press, 2003), chapter 6, 'Polynomial interpolation' — the error formula and divided differences with full proofs.
↓But there is a problem
For the orbit, more points gave a better answer every time. It seems obvious that this always holds: more data, more accuracy. In 1901 the German mathematician Carl Runge tried the most natural experiment imaginable — interpolate the smooth, bell-shaped function 1/(1+25x2)1/(1 + 25x^2) at more and more equally spaced points on [−1,1][-1, 1]. In the middle the fit improved. Near the ends the polynomial began to oscillate, and the oscillations grew without bound as points were added.