Where the last four guides leave us
By now you can build the unique degree-n polynomial through n+1 points in either the Lagrange or the Newton form, and you can read the error formula: f(x) - p(x) = f^{(n+1)}(xi) / (n+1)! times the node polynomial (x - x_0)(x - x_1)...(x - x_n). Guide 3 then delivered the gut-punch — Runge's phenomenon. On equispaced points, raising the degree on a perfectly smooth function like 1/(1 + 25 x^2) makes the interpolant oscillate ever more violently near the ends of the interval, and the error diverges as n grows. More points, more wrong.
Guide 4 gave one escape: stop using one high-degree polynomial and stitch together many low-degree pieces — the cubic spline and its cousins. Splines are the right tool when you are stuck with data given to you. But there is a second, very different escape that this guide is about: keep the single global polynomial, keep the high degree — just stop putting the points in the wrong places. The villain in the error formula is not the degree; it is the node polynomial (x - x_0)...(x - x_n), and you get to choose the x_i. That freedom is the whole game.
Chebyshev nodes: crowd the edges
Look hard at the node polynomial w(x) = (x - x_0)(x - x_1)...(x - x_n). With equispaced points it is small in the middle but grows huge near the two ends — and the error inherits that shape, which is exactly why Runge's wiggles live at the edges. So the fix writes itself: choose the x_i to keep the peak of |w(x)| as small as possible across the interval. The points that do this are not spread evenly. They cluster toward the endpoints and thin out in the middle. These are the Chebyshev nodes.
The cleanest picture: take n+1 equally spaced points around the upper half of a circle, then drop each straight down onto the diameter. The projections bunch up near the two ends — that is the Chebyshev distribution. Formally, on the interval [-1, 1] the Chebyshev nodes are x_k = cos((2k + 1) / (2n + 2) * pi), the cosines of equally spaced angles. With these points the peak of |w(x)| is as small as it can possibly be, smaller than for any other choice of n+1 nodes — that is a theorem, not a heuristic.
The payoff against Runge is dramatic. For 1/(1 + 25 x^2) the equispaced interpolant's error grows without bound as you add points; the Chebyshev interpolant's error shrinks toward zero. For a function that is genuinely smooth (analytic) the Chebyshev interpolation error decays geometrically — like c^n for some c < 1 — so a few dozen well-placed points can reach machine precision. The honest caveat: this rescue assumes the function is smooth. If f has a kink or a jump, no node placement makes a single polynomial converge fast — that is when you reach back for the splines of guide 4.
The barycentric formula: fast and stable
Choosing good nodes is only half the job; you still have to evaluate the interpolant without wrecking it in floating point. The naive Lagrange form re-multiplies n+1 basis polynomials from scratch at every x — O(n^2) work per point and numerically clumsy. The fix is the barycentric formula, which rewrites the very same polynomial as a weighted ratio. You precompute one weight w_j per node once, in O(n^2) total, and after that every evaluation costs only O(n) and is beautifully stable.
barycentric form of the interpolant at a point x:
sum_j [ w_j / (x - x_j) ] * y_j
p(x) = ---------------------------------------
sum_j [ w_j / (x - x_j) ]
w_j = 1 / prod_{k != j} (x_k - x_j) (precompute once)
if x equals some node x_j exactly -> return y_j directlyTwo things make this more than a tidy rewrite. First, Chebyshev nodes have known barycentric weights — the simple alternating +1, -1 pattern — so you skip the O(n^2) setup entirely and the method becomes O(n) start to finish. Second, the barycentric form is backward stable on a sensible node set: small relative changes in the y_j produce only small relative changes in p(x). Adding a new data point is also cheap — you update the weights in O(n) rather than rebuilding everything. This is the formula real software uses; it is what makes high-degree Chebyshev interpolation practical instead of merely theoretical.
Interpolation versus best approximation
Now the deeper distinction. An interpolant is pinned to touch its data points exactly — but who said the best approximation has to touch anything? Best approximation in the minimax (or uniform) sense asks for the degree-n polynomial p that makes the worst-case error max over x of |f(x) - p(x)| as small as possible. This minimax polynomial generally does not pass through any prescribed points; it floats freely, balancing its error above and below the function so that no single peak is taller than it has to be. It is the genuinely optimal answer, the bar everything else is measured against.
How would you even recognize the best one? The equioscillation theorem of Chebyshev gives a startlingly crisp answer. The degree-n minimax polynomial is the unique one whose error curve f(x) - p(x) attains its maximum magnitude at least n+2 times, with the sign alternating +, -, +, -, ... at those extrema. Picture the error as a wave that touches the same height up and down, n+2 times, like a perfectly level row of equal peaks and troughs. If your candidate's error wave is lopsided — one peak taller than the rest — you can always lower that peak by nudging the polynomial, so you are not yet optimal.
This is exactly why Chebyshev nodes are so good and not a coincidence of names. The node polynomial built on Chebyshev points is itself an equioscillating Chebyshev polynomial — its peaks are all the same height. So Chebyshev interpolation is, in effect, mimicking the equioscillation that defines the true optimum. The result is a near-best theorem: the Chebyshev interpolant's worst error is larger than the unbeatable minimax error by only a slowly growing factor (it creeps up like log n), so for practical degrees it is within a few percent of optimal. You get almost all of best approximation, just by sampling at the right points.
If you truly need the exact minimax polynomial, the classic tool is the Remez exchange algorithm: guess a set of n+2 points, fit a polynomial that equioscillates on them, find where the actual error peaks really are, swap your points toward those peaks, and repeat until the peaks level out. It converges fast, but it is iterative and fiddly — which is exactly why, in practice, people so often just take the Chebyshev interpolant and call it close enough. Optimal is nice; near-optimal-for-free is usually better engineering.
Orthogonal polynomials, and the bigger picture
There is a second, equally important way to measure error: not the worst-case peak, but the integrated squared error — the least-squares sense from the previous rung. Minimizing that over polynomials leads straight to orthogonal polynomials. Just as orthogonal vectors made least-squares solving clean, a family of polynomials that are mutually orthogonal (under an integral inner product) lets you find the best least-squares polynomial approximation coefficient by coefficient, with no ill-conditioned system to solve. The Chebyshev polynomials are one such family; the Legendre polynomials are another.
These families behave well for the same structural reason Gram-Schmidt orthogonalization tamed least squares: orthogonality decouples the coefficients, so adding one more term never disturbs the ones you already found, and the linear system is diagonal instead of a near-singular Vandermonde. They are also the engine of two later rungs — Gauss-Legendre quadrature places its sample points at the roots of the Legendre polynomial to integrate as accurately as possible, and spectral methods expand the solution of a differential equation in exactly such a basis. Orthogonal polynomials are one of the great connecting threads of numerical analysis.
Underwriting all of this is the Weierstrass approximation theorem: any continuous function on a closed interval can be approximated as closely as you like by some polynomial. That is the license — it guarantees a good polynomial always exists. But it is only existence; it does not tell you how to find it or how high a degree you need, and it says nothing about equispaced interpolation, which (Runge) can fail spectacularly. The art of this rung is turning Weierstrass's promise into a real, computable, stable approximation — and Chebyshev nodes with the barycentric formula are the cleanest way to cash that promise in.
Carry three things up the ladder. First, the degree is rarely the problem — node placement is; equispaced points cause Runge, Chebyshev points cure it by shrinking the node polynomial, and the barycentric formula evaluates the result in O(n) with backward stability. Second, interpolation and best approximation are different questions, joined by the equioscillation theorem, and Chebyshev interpolation is near-best for free — within a log-n whisker of the unbeatable minimax answer. Third, every one of these is an approximation: smoothness sets the convergence rate, a kink or jump caps it (use splines then), and in finite-precision arithmetic the right basis and the right points are what stand between you and a number you can trust.