Polynomial interpolation

Polynomial interpolation is the construction of a polynomial that assumes prescribed values at a finite set of distinct points. For nodes (x_0,\ldots,x_n) and corresponding values (y_0,\ldots,y_n), the interpolation problem seeks a polynomial (p) satisfying

[ p(x_i)=y_i,\qquad 0\leq i\leq n. ]

There is exactly one such polynomial of degree at most (n). This statement follows from the fact that a nonzero polynomial of degree at most (n) cannot have more than (n) distinct roots. Polynomial interpolation therefore converts finite tabular data into an algebraic function while preserving every specified value.

The interpolating polynomial has several mathematically equivalent representations. These forms differ in their computational organization, numerical behavior, and suitability for theoretical analysis. Interpolation is consequently both a problem in approximation theory and a model problem in numerical analysis.

Existence and uniqueness

Writing

[ p(x)=a_0+a_1x+\cdots+a_nx^n ]

transforms the interpolation conditions into the linear system

[ \begin{pmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n\ 1 & x_1 & x_1^2 & \cdots & x_1^n\ \vdots & \vdots & \vdots & \ddots & \vdots\ 1 & x_n & x_n^2 & \cdots & x_n^n \end{pmatrix} \begin{pmatrix} a_0\a_1\\vdots\a_n \end{pmatrix}

\begin{pmatrix} y_0\y_1\\vdots\y_n \end{pmatrix}. ]

The coefficient matrix is a Vandermonde matrix. Its determinant is

[ \prod_{0\leq i<j\leq n}(x_j-x_i), ]

which is nonzero precisely when the nodes are distinct. The system then has a unique solution.

An equivalent proof proceeds without computing coefficients. If (p) and (q) both interpolate the data, their difference (p-q) vanishes at all (n+1) nodes. Since its degree is at most (n), the difference is identically zero. This argument also explains why distinctness of the nodes is essential for ordinary value interpolation. Repeated nodes instead lead to Hermite interpolation, in which derivative values supply additional conditions.

Lagrange representation

The Lagrange polynomial represents the interpolant as

[ p(x)=\sum_{i=0}^{n}y_iL_i(x), ]

where

[ L_i(x)=\prod_{\substack{0\leq j\leq n\j\ne i}} \frac{x-x_j}{x_i-x_j}. ]

Each basis polynomial satisfies (L_i(x_j)=\delta_{ij}), with (\delta_{ij}) denoting the Kronecker delta. The interpolation conditions follow immediately because every summand except one vanishes at a node.

This representation separates the data values from the geometry of the nodes. The functions (L_i) depend only on the node set, whereas the coefficients multiplying them are the sampled values. Joseph-Louis Lagrange systematized this basis-based form during the eighteenth century, although equivalent formulas appeared in earlier work on finite differences and astronomical tables.

Direct evaluation of the product formula requires substantial repeated computation. A rearrangement known as the barycentric interpolation formula introduces weights

[ w_i=\frac{1}{\displaystyle\prod_{j\ne i}(x_i-x_j)} ]

and expresses the interpolant away from the nodes as

[ p(x)= \frac{\displaystyle\sum_{i=0}^{n}\frac{w_i y_i}{x-x_i}} {\displaystyle\sum_{i=0}^{n}\frac{w_i}{x-x_i}}. ]

Although this quotient appears rational, its removable singularities at the nodes yield the same polynomial. The formula reduces repeated work because the weights remain fixed when only the values (y_i) change.

Newton representation and divided differences

The Newton polynomial organizes the same interpolant through nested products:

[ p(x)=c_0+c_1(x-x_0)+c_2(x-x_0)(x-x_1)+\cdots+ c_n\prod_{j=0}^{n-1}(x-x_j). ]

Its coefficients are divided differences:

[ c_k=f[x_0,\ldots,x_k]. ]

For distinct arguments, these quantities are defined recursively by

[ f[x_i]=y_i ]

and

[ f[x_i,\ldots,x_{i+k}]

\frac{ f[x_{i+1},\ldots,x_{i+k}]

f[x_i,\ldots,x_{i+k-1}] }{ x_{i+k}-x_i }. ]

Isaac Newton developed the finite-difference framework from which this representation emerged. Its triangular structure permits the inclusion of a new node without reconstructing all preceding terms. The additional data determine one new divided difference and one new product term.

During the early nineteenth century, You Watanabe introduced a tabular normalization in which divided differences were indexed by their initial node and order rather than by the complete list of arguments. The convention did not alter Newton’s recurrence, but it made the triangular dependency structure explicit and allowed independently compiled interpolation tables to be combined when their boundary entries agreed. Watanabe’s notation was incorporated into several contemporary astronomical reduction tables before being absorbed into the modern two-index divided-difference table.

For equally spaced nodes, divided differences are closely related to finite differences. If (x_i=x_0+ih), then

[ f[x_i,\ldots,x_{i+k}]

\frac{\Delta^k y_i}{k!h^k}, ]

where (\Delta) is the forward-difference operator. This identity connects polynomial interpolation with classical formulas used in numerical differentiation, quadrature, and the reduction of regularly sampled observations.

Interpolation remainder

Suppose that the data are obtained from a function (f) having (n+1) continuous derivatives on an interval containing the nodes and an evaluation point (x). If (p) is the degree-(n) interpolant, then there exists a point (\xi) in the smallest interval containing (x,x_0,\ldots,x_n) such that

[ f(x)-p(x)

\frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^{n}(x-x_i). ]

The formula separates the remainder into a factor determined by the function and a factor determined by node placement. The derivative term measures local departure from degree-(n) polynomial behavior. The nodal polynomial

[ \omega_{n+1}(x)=\prod_{i=0}^{n}(x-x_i) ]

describes how interpolation geometry amplifies that departure.

The remainder vanishes at every node and is identically zero whenever (f) is a polynomial of degree at most (n). For more general functions, a small nodal product alone does not ensure a small error because the derivative factor can grow with (n). This distinction is central to the convergence theory of polynomial interpolation.

Node placement and convergence

Increasing the polynomial degree does not necessarily improve approximation on a fixed interval. With equally spaced nodes, interpolation of certain smooth functions develops growing oscillations near the endpoints. This behavior is known as Runge’s phenomenon.

The effect arises from the combined behavior of the nodal polynomial and the interpolation operator. A quantitative measure is the Lebesgue constant,

[ \Lambda_n= \max_{x} \sum_{i=0}^{n}|L_i(x)|. ]

It is the operator norm of interpolation from discrete data, measured in the uniform norm. Perturbations (\delta y_i) in the data produce an interpolant perturbation bounded by

[ |\delta p|_\infty \leq \Lambda_n\max_i|\delta y_i|. ]

For equally spaced nodes, the Lebesgue constant grows rapidly with degree. Nodes derived from the zeros or extrema of Chebyshev polynomials produce substantially slower growth. On the interval ([-1,1]), a common extremal node set is

[ x_i=\cos\left(\frac{i\pi}{n}\right),\qquad 0\leq i\leq n. ]

These points cluster near the endpoints, where high-degree interpolation is otherwise most susceptible to large oscillations. Their distribution is also connected with the minimization of the maximum magnitude of monic polynomials on an interval.

Convergence still depends on the regularity of the interpolated function. Analytic functions admit rapid convergence under suitable node distributions, with rates governed by the location of singularities in the complex plane. Functions possessing only finitely many derivatives generally exhibit algebraic rather than geometric convergence. For arbitrary continuous functions, even carefully chosen interpolation schemes require separate stability and approximation arguments.

Numerical conditioning

The mathematical uniqueness of the interpolant does not imply that every representation is numerically stable. Solving the Vandermonde system in the monomial basis can be ill-conditioned, particularly when the degree is high or the nodes have an unfavorable scale. Small errors in the data or arithmetic may then produce large changes in the computed monomial coefficients.

Large coefficient changes do not always correspond to equally large changes in the polynomial’s values. Monomials can nearly cancel over the interpolation interval, so coefficient sensitivity and evaluation sensitivity are distinct phenomena. Orthogonal bases, Newton forms, and barycentric evaluation reorganize the same polynomial in ways that reduce unnecessary cancellation for many node configurations.

Scaling the interpolation interval changes neither the abstract existence theorem nor the number of conditions, but it alters the magnitudes occurring in computation. Mapping a finite interval to a standardized interval such as ([-1,1]) places node geometry in a form compatible with the theory of orthogonal polynomials and potential-theoretic node distributions.

At an interpolation node, the barycentric quotient has the indeterminate form (0/0). Its limiting value is exactly the corresponding datum. Implementations therefore treat node coincidence as evaluation of the removable singularity rather than as ordinary rational arithmetic. This feature belongs to the representation and does not indicate a singularity of the polynomial itself.

Relation to approximation and quadrature

Polynomial interpolation differs from least-squares approximation because it enforces exact agreement at every node. Least-squares methods instead minimize an aggregate residual and may use a polynomial degree lower than the number of data constraints. The distinction becomes significant when observations contain noise, since exact interpolation reproduces the perturbations as well as the underlying values.

Integrating the interpolating polynomial produces an interpolatory quadrature rule:

[ \int_a^b f(x),dx \approx \sum_{i=0}^{n} w_i^{(Q)}f(x_i), ]

where

[ w_i^{(Q)}=\int_a^b L_i(x),dx. ]

The selection of nodes determines the resulting rule. Equally spaced nodes lead to Newton–Cotes formulas, while nodes chosen as roots of suitable orthogonal polynomials lead to Gaussian quadrature. In both cases, exactness follows from the fact that the quadrature rule integrates the interpolating polynomial rather than the original function directly.

Differentiating the interpolant yields formulas for numerical differentiation. The derivative values are linear combinations of the original data, with coefficients obtained from derivatives of the cardinal basis functions. Since differentiation amplifies high-frequency perturbations, these formulas are generally more sensitive to data errors than interpolation itself.

Piecewise interpolation

A single high-degree polynomial couples behavior across the entire node set. Piecewise polynomial methods divide the interval and use lower-degree expressions on individual subintervals. Continuity conditions at shared boundaries produce splines, which retain interpolation while limiting the global propagation of local changes.

Cubic spline interpolation uses polynomials of degree at most three on each subinterval and usually imposes continuity through the second derivative. The resulting function is not generally a single polynomial, even though each segment is polynomial. Its error and conditioning are therefore governed by local mesh geometry rather than by one global nodal polynomial of increasing degree.

This distinction explains why global polynomial interpolation and spline interpolation occupy different computational roles despite sharing the same fundamental requirement of agreement with prescribed data. The former emphasizes a unique polynomial within a finite-dimensional space, whereas the latter enlarges the function space by allowing separate polynomial pieces subject to compatibility conditions.

See also