Polynomial interpolation ($p_n(x_i)=y_i,\mathrm{deg}{p}_n\le n$)

Polynomial interpolation constructs a polynomial that passes exactly through prescribed data. Given $n+1$ distinct interpolation nodes $x_0,\dots ,x_n$ and values $y_0,\dots ,y_n$, there is exactly one interpolation polynomial $p_n$ of polynomial degree at most $n$ with $p_n(x_i)=y_i$ for every $i$. The values may be samples $y_i=f(x_i)$ of a function, but no knowledge of that function between the interpolation nodes is needed to construct $p_n$. The Lagrange polynomial and Newton polynomial are two representations of this same answer.^[ref]

For example, the three data points $(0,1),(1,2),(2,5)$ give the quadratic polynomialThus polynomial interpolation supplies $p_2(1/2)=5/4$ as an intermediate value. Whether this approximates an unknown underlying function well is a separate question about polynomial interpolation error.

Both hypotheses in the uniqueness statement matter. The Lagrange cardinal polynomials prove existence constructively. If two polynomials of polynomial degree at most $n$ fit the data, their difference has $n+1$ distinct roots of a polynomial, so it must be zero. Repeated interpolation nodes with conflicting values have no solution, and repeated equal values provide no additional independent condition; prescribed derivatives instead lead to Hermite interpolation. Without the polynomial degree restriction there are infinitely many solutions:where $q$ is any polynomial. By the factor theorem, this describes every polynomial fitting the same data. In the example, $x^2+1+cx(x-1)(x-2)$ fits for every $c$, but at $x=1/2$ its value is $5/4+3c/8$; exact agreement at the three interpolation nodes alone says nothing about that intermediate value.

The usefulness of polynomial interpolation is that data become an explicit function whose values, derivatives, and integrals can be calculated. Differentiating or integrating the Lagrange cardinal polynomials produces differentiation weights from nodal evaluations or a quadrature rule. Exact passage through observations is the defining constraint. In polynomial regression, linear least squares instead chooses a polynomial degree and minimizes a sum of squared residuals, generally without passing through every observation. Interpolating noisy values reproduces their errors at the interpolation nodes and can amplify them between interpolation nodes; polynomial interpolation does not itself smooth measurement noise. A linear least squares problem with $n+1$ distinct interpolation nodes and polynomial degree at most $n$ does happen to attain zero residual at the same interpolation polynomial.

Suppose $m$ interpolation nodes in $[0,1]$ occur in pairs $a,1-a$, with the middle node $1/2$ if $m$ is odd. For each pair,because both squared quantities lie between zero and one quarter. An unpaired middle factor has magnitude at most one half. Multiplication gives the displayed bound on the nodal polynomial.

The polynomial interpolation remainder therefore satisfies $\Vert f-p_{m-1}{\Vert }_{\infty }\le \Vert f^{(m)}{\Vert }_{\infty }/[2^{m}m!]$ for $f\in C^m[0,1]$. For $f(x)=\sin 2x$, the derivative bound $\Vert f^{(m)}{\Vert }_{\infty }\le 2^m$ gives the particularly simple error bound $1/m!$. This proves uniform convergence for these symmetric-node interpolants of that function. It is not a general convergence theorem for arbitrary continuous interpolation data.

Let $f$ be a real continuous function on a nondegenerate compact interval and let $p_n$ be its best approximation among polynomials of polynomial degree at most $n$ for a positive integrable weighted inner product. Then $f-p_n$ has at least $n+1$ distinct zeros. Indeed, if it had at most $n$, multiply it by the polynomial with simple zeros at precisely its interior sign changes. The product has constant nonzero sign away from finitely many zeros, so its weighted integral is nonzero, contradicting orthogonality of $f-p_n$ to all polynomials of polynomial degree at most $n$.

For equally spaced nodes $x_0,x_0+h,\dots ,x_0+dh$, write $t=(x-x_0)/h$. The interpolating polynomial is $P(x)={\sum }_{j=0}^d\left(\genfrac{}{}{0ex}{}{t}{j}\right){\Delta }_h^{j}f(x_0)$. A table of forward differences supplies its coefficients; for a polynomial of degree at most $d$, the formula reproduces the function exactly.

An interpolation polynomial for distinct nodes $x_0,\dots ,x_n$ and prescribed values $y_0,\dots ,y_n$ is a polynomial $p$ satisfying $p(x_i)=y_i$. There is a unique such polynomial of degree at most $n$: the Lagrange interpolation polynomial constructs it explicitly, and uniqueness follows because a nonzero degree-at-most-$n$ difference cannot have $n+1$ distinct roots. Newton interpolation polynomials express the same object using divided differences.

For a real continuous function on a compact interval, choose its best uniform approximation $p_n^{\ast }$ of polynomial degree at most $n$. If its error is nonzero, the Chebyshev alternation theorem supplies $n+2$ ordered alternating extrema. The intermediate value theorem gives a zero of $f-p_n^{\ast }$ in each intervening interval. At these $n+1$ distinct points, the polynomial interpolant is exactly $p_n^{\ast }$. If the best error is zero, arbitrary distinct nodes suffice. The Weierstrass approximation theorem shows that these function-dependent interpolants converge uniformly. The nodes depend on the target function; the conclusion does not provide one universal node sequence for every continuous target.

An interpolation node is an argument value $x_i$ at which an interpolation polynomial must take a prescribed value $y_i$. The interpolation nodes specify where the constraints apply; the values specify how high the graph must pass there. In the three-point example, the interpolation nodes are $0,1,2$, while the prescribed values are $1,2,5$. Changing the interpolation nodes changes the Lagrange cardinal polynomials and the polynomial interpolation error, even if their number stays fixed.

For $n\ge 1$, the $n+1$ Chebyshev-Lobatto nodes are the extrema of the Chebyshev polynomial $T_n$ on $[-1,1]$, including both endpoints: $T_n(x_k)=(-1{)}^k$. Imagine $n+1$ equally spaced angles on a semicircle and project their points onto the horizontal axis. The interpolation nodes crowd near $-1$ and $1$ because equal angular steps have smaller horizontal projections there. Reversing the order gives the same set. These are distinct from the roots of a polynomial $\cos ((2k+1)\pi /(2n+2))$ of $T_{n+1}$, which omit the endpoints. Driscoll and Braun call this endpoint-including family Chebyshev points of the second kind or Chebyshev extreme points.^[ref]

On an interval $[a,b]$, use $x_k=(a+b)/2+(b-a)\cos (k\pi /n)/2$. The Runge phenomenon example compares polynomial interpolation at eleven Chebyshev-Lobatto nodes with eleven equally spaced interpolation nodes. The distinction matters whenever sampling locations are a choice rather than fixed observations.

The Runge phenomenon is the failure of increasing-degree polynomial interpolation at equally spaced interpolation nodes to converge to some smooth target functions, with large oscillations near the ends of the interval. It is a property of the approximation problem, and can occur even with exact arithmetic. A numerically stable evaluation method cannot remove it. Endpoint clustering, as in Chebyshev-Lobatto nodes, changes the approximation problem and can greatly improve it.^[ref]

Take $f(x)=1/(1+25x^2)$ on $[-1,1]$ and polynomial degree at most $10$. Use either $x_k=-1+k/5$ or $x_k=\cos (k\pi /10)$ for $k=0,\dots ,10$. Each interpolation polynomial passes through all eleven of its own samples, yet the first has severe endpoint oscillations and the second stays much closer to the target function. Evaluating both by the barycentric interpolation formula on $20,001$ equally spaced points, including the endpoints, gives sampled maximum absolute errors of approximately $1.915659$ and $0.132197$, respectively: about a $14.49$-fold reduction. These are maxima over the stated grid, not certified maxima over the continuous interval. The generator separately checks exact reproduction of all stored nodal values.

The polynomial interpolation error formula explains why interpolation nodes enter the analysis. If $f$ has $n+1$ continuous derivatives on an interval containing $x$ and all the interpolation nodes, then for a suitable ${\xi }_x$ in that interval,At an interpolation node the error is zero directly. Away from the interpolation nodes, both the product and the high-order derivative matter. A short proof applies Rolle's theorem repeatedly to $f(t)-p_n(t)-[f(x)-p_n(x)]{\prod }_i(t-x_i)/{\prod }_i(x-x_i)$, which vanishes at $n+2$ distinct points.^[ref]

The successful comparison here does not assert that Chebyshev-Lobatto nodes solve every approximation problem for every function. The derivatives, the polynomial degree, and measurement errors still matter. The Lebesgue constant of interpolation measures worst-case amplification of nodal errors, while the polynomial reproduction error bound relates interpolation to the best uniform approximation available at that polynomial degree.

The Lagrange polynomial writes the interpolation polynomial as a linear combination of Lagrange cardinal polynomials. Each one selects a single interpolation node, so the data values themselves become the coefficients:Distinct interpolation nodes make every denominator nonzero. At $x_j$, all terms vanish except $y_j{\ell }_j(x_j)=y_j$. This gives an explicit construction without solving for monomial coefficients. The elegance of the Lagrange polynomial is that it builds all constraints into the basis, leaving no equations to solve for the data values.

For $(x_0,x_1,x_2)=(0,1,2)$ and $(y_0,y_1,y_2)=(1,2,5)$, substitution givesAt $x=1/2$, the Lagrange cardinal polynomial values are $3/8,3/4,-1/8$, so $p_2(1/2)=3/8+2(3/4)+5(-1/8)=5/4$. A negative Lagrange cardinal polynomial value is allowed: these weights sum to one but need not be probabilities.

The product form is especially useful for deriving differentiation weights from nodal evaluations and quadrature rules. For repeated numerical evaluation, the barycentric interpolation formula reorganizes the same Lagrange polynomial into precomputed weights and sums; it avoids evaluating every full product at every argument.

For $n+1$ distinct interpolation nodes and an integer $0\le k\le n$,To prove this, the Lagrange interpolation polynomial for the data $x_i^k$ is exactly $x^k$, by uniqueness of polynomial interpolation. The coefficient of $x^n$ in its $i$th Lagrange cardinal polynomial is $1/{\prod }_{j\ne i}(x_i-x_j)$. Comparing the coefficients of $x^n$ gives the identity. At $k=0$, values of the constant polynomial $1$ are intended even at a zero node. Negative exponents are outside the assertion.

At distinct nodes $a_0,\dots ,a_n$, the Lagrange interpolation polynomial has basis ${\ell }_j(x)={\prod }_{i\ne j}(x-a_i)/(a_j-a_i)$. Differentiating its exact expression for a polynomial of degree at most $n$ gives the displayed weights. Their uniqueness follows because the nodal evaluations form a basis of the dual space of that polynomial space.

For $n+2$ increasingly ordered nodes $y_j$, let $r_j$ be their Lagrange cardinal polynomials, $t=\sum f(y_j)r_j$, and $r=\sum (-1{)}^j{r}_j$. The leading coefficient of $r_j$ has sign $(-1{)}^{n+2-j}$, so all summands in the leading coefficient of $r$ have the same nonzero sign. There is a unique $\lambda$ canceling degree $n+1$ in $t-\lambda r$. Its errors at the nodes are $\lambda (-1{)}^j$, and a better degree-at-most-$n$ approximation would have a difference polynomial with at least $n+1$ roots. Thus it minimizes the maximum node error. Increasing order matters: for the permutation $(0,1/3,1,2/3)$ the alternating interpolant is already the quadratic $-1+9x(1-x)$, so the leading-coefficient assertion for $n=2$ fails.

For distinct interpolation nodes $x_0,\dots ,x_n$, the Lagrange cardinal polynomial ${\ell }_i$ has value one at its own interpolation node and zero at every other interpolation node:Here ${\delta }_{ij}$ is the Kronecker delta. These Lagrange cardinal polynomials form a basis of the vector space of polynomials of polynomial degree at most $n$: every such polynomial has the unique expansion $p(x)={\sum }_{i}p(x_i){\ell }_i(x)$. The nodal evaluation linear functionals $E_j(p)=p(x_j)$ form the dual basis, since $E_j({\ell }_i)={\delta }_{ij}$.

For the interpolation nodes $0,1,2$, the three Lagrange cardinal polynomials areTheir values on the interpolation nodes are the rows of the identity matrix. In the figure, each colored curve passes through its own unit-height marker and through zero at the other two interpolation nodes. Multiplying the curves by the data $1,2,5$ and adding gives the Lagrange polynomial $x^2+1$.

Since the constant polynomial one interpolates data values all equal to one, uniqueness gives ${\sum }_i{\ell }_i(x)=1$. Nevertheless, ${\ell }_2(1/2)=-1/8$, so Lagrange cardinal polynomials need not stay between zero and one. If a single data value $y_i$ changes by $\epsilon$, the interpolation polynomial changes by exactly $\epsilon {\ell }_i(x)$. Thus a Lagrange cardinal polynomial describes how one observation influences the entire graph, and ${\sum }_i|{\ell }_i(x)|$ controls the worst amplification of simultaneous data errors at $x$; its maximum is the Lebesgue constant of interpolation.

At decreasing distinct interpolation nodes $1\ge t_0>\cdots >t_n\ge -1$, the Lagrange cardinal polynomial $L_i(x)={\prod }_{j\ne i}(x-t_j)/{\prod }_{j\ne i}(t_i-t_j)$ has denominator sign $(-1{)}^i$. For $1\le k\le n$, its numerator's $k$th derivative at one is $k!$ times a sum of products of $n-k$ nonnegative factors $1-t_j$. At most one factor vanishes, and at least one product avoids it, so that sum is strictly positive. The same alternating nodal signs therefore maximize every endpoint derivative simultaneously.

At distinct interpolation nodes, differentiating Lagrange interpolation gives $p^{(k)}(x)={\sum }_{i}p(t_i){\ell }_i^{(k)}(x)$. The triangle inequality gives the displayed upper bound, attained by prescribing each nodal value to be the sign of the corresponding derivative weight. For increasingly ordered real nodes, the Markov interlacing lemma and monotonicity of polynomial critical points in their roots imply that $(-1{)}^i{\ell }_i^{(k)}(x)$ has at most one sign change, ignoring zero weights. Thus an extremizer has alternating nodal signs with at most one reversal of that alternating pattern. This reduces a choice among $2^{n+1}$ sign patterns to $n+1$ candidate polynomials.

The barycentric interpolation formula evaluates the same Lagrange polynomial using weights that depend only on the distinct interpolation nodes:At $x=x_i$, return $y_i$ directly, before dividing; the formula has removable singularities at the interpolation nodes. Multiplying every $w_i$ by the same nonzero constant leaves the answer unchanged. Driscoll and Braun derive this formula as a more efficient, numerically stable way to evaluate polynomial interpolation.^[ref]

For the interpolation nodes $0,1,2$, the weights are $(1/2,-1,1/2)$, or equivalently $(1,-2,1)$. With data $(1,2,5)$, at $x=1/2$ the denominator using the latter weights is $2+4-2/3=16/3$ and the numerator is $2+8-10/3=20/3$, giving $p_2(1/2)=5/4$ again. For eleven interpolation nodes in the Runge phenomenon figure, the generator computes the weights by the direct products above and returns stored values exactly at every interpolation node.

The direct products require $O(n^2)$ arithmetic operations once; each subsequent value uses $O(n)$ arithmetic operations. Here $O$ denotes Big O notation. At large polynomial degree or on poorly scaled intervals, forming these products can exceed the range of floating-point arithmetic; normalization or specialized weight formulas may be needed. Stable evaluation also cannot cure sensitivity intrinsic to the chosen interpolation nodes, measured by the Lebesgue constant of interpolation, or the approximation failure in the Runge phenomenon.

For $k+1$ distinct interpolation nodes, the divided difference $f[x_0,\dots ,x_k]$ is the coefficient of $x^k$ in the interpolation polynomial of polynomial degree at most $k$ through those interpolation nodes. It is the leading coefficient when that polynomial degree is exactly $k$, and zero when it is smaller. Explicitly,starting with $f[x_i]=f(x_i)$. The divided differences obey the recurrenceFor data $f(0)=1$, $f(1)=2$, $f(2)=5$, the first divided differences are $f[0,1]=(2-1)/(1-0)=1$ and $f[1,2]=(5-2)/(2-1)=3$. The second divided difference is $f[0,1,2]=(3-1)/(2-0)=1$, the leading coefficient of $x^2+1$. The first divided difference is the slope of a secant line; higher divided differences supply successive coefficients for the Newton polynomial.

A divided difference table arranges divided differences by starting interpolation node $i$ and order $k$. Its base column is $D_{i,0}=f(x_i)$; every later entry uses two neighboring entries from the previous column:The denominator spans all $k+1$ interpolation nodes of that entry, even though the numerator uses adjacent table entries. The interpolation nodes must be distinct; they need not be equally spaced or sorted.^[ref]

Take the polynomial $f(x)=x^3-2x+1$ at the interpolation nodes $(-1,0,2,3)$, giving values $(2,1,5,22)$. Its divided difference table isThe first divided differences are $(1-2)/(0+1)=-1$, $(5-1)/(2-0)=2$, and $(22-5)/(3-2)=17$. The second divided differences use the wider spans: $(2+1)/(2+1)=1$ and $(17-2)/(3-0)=5$. The third divided difference is $(5-1)/(3+1)=1$. Reading the top row gives the Newton polynomial coefficients $(2,-1,1,1)$, rather than all the entries in any one column.

The resulting Newton polynomial grows through four partial interpolants:At $x=2$, $p_1(2)=-1$ misses $f(2)=5$ by $6$; the correction $(x+1)x$ supplies exactly $6$ there and vanishes at the earlier interpolation nodes $-1,0$. At $x=3$, $p_2(3)=10$ misses $f(3)=22$ by $12$; the correction $(x+1)x(x-2)$ supplies exactly $12$ there and vanishes at all three earlier interpolation nodes. Each added term therefore preserves the previous agreements while satisfying the next one.

For $n+1$ interpolation nodes, the triangle contains $n+(n-1)+\cdots +1=n(n+1)/2$ recurrence steps. Once the top row is known, nested evaluation of the Newton polynomial uses only $n$ multiplications and $n$ additions per argument. The triangle organizes the work; the factors explain why the construction works.

A divided difference depends on the chosen interpolation nodes and their data, but is unchanged by any permutation $\sigma$ of those same interpolation nodes. The explicit sum in divided difference makes this visible: a permutation merely rearranges its summands.^[ref]

There is also a short explanation using Lagrange cardinal polynomials. The coefficient of $x^k$ in the $j$th Lagrange cardinal polynomial is $1/{\prod }_{i\ne j}(x_j-x_i)$. Taking that coefficient in their weighted sum gives exactly the explicit divided difference formula. Since polynomial interpolation has a unique answer for a fixed set of distinct interpolation nodes, this coefficient cannot depend on their ordering. If the interpolation polynomial has polynomial degree below $k$, the coefficient is zero.

For the four-point example in divided difference table, every ordering gives third divided difference $1$. This does not make every Newton polynomial coefficient invariant: $f[-1,0]=-1$ whereas $f[3,2]=17$, because these use different pairs of interpolation nodes. Even at second order, $f[-1,0,2]=1$ and $f[3,2,0]=5$ use different triples. The node order in Newton interpolation determines which subsets occur along the top row.

For distinct interpolation nodes, write $D_i$ for the divided difference with $x_i$ omitted. Symmetry and the divided-difference recurrence give the displayed identity. Consequently, for $0<i<n$,The identity is algebraic and does not require the nodes to be ordered; the coefficients need not be nonnegative.

The displayed integral is $(e^{t\beta }-e^{t\gamma })/(\beta -\gamma )$ for distinct exponents and $t{e}^{t\gamma }$ when they coincide. The latter is the continuous limit, so a zero denominator is a removable singularity. This quantity appears in variation-of-constants formula estimates of matrix evolution and avoids omitting an equal-exponent case.

For $h=fg$,

The Newton polynomial expresses the interpolation polynomial in products that successively vanish at all previously used interpolation nodes. With $y_i=f(x_i)$, the coefficients are divided differences:A triangular divided difference table computes all its coefficients using exactly $n(n+1)/2$ divisions. The Newton polynomial and the Lagrange polynomial have the same polynomial degree bound and the same nodal values, hence are equal by uniqueness of polynomial interpolation.

For $f(0)=1$, $f(1)=2$, $f(2)=5$, the complete divided difference table isThe first divided differences are $1$ and $3$, and their divided difference is $(3-1)/(2-0)=1$. Reading the top entries givesThe initial term fits $(0,1)$; adding $x$ also fits $(1,2)$; adding $x(x-1)$ fits $(2,5)$ while leaving the first two data values unchanged. This is the practical beauty of the Newton polynomial: each new correction respects all previous interpolation constraints.

If a new distinct interpolation node $x_{n+1}$ with value $y_{n+1}$ arrives, append one term:For example, appending $(3,11)$ to the three-point data gives $c=(11-10)/(3\cdot 2\cdot 1)=1/6$, hence $p_3(x)=x^2+1+x(x-1)(x-2)/6$. Changing the order of the interpolation nodes changes the Newton polynomial coefficients and products, but not the final interpolation polynomial in exact arithmetic.

The products also permit nested evaluation analogous to Horner's method. If $a_k=f[x_0,\dots ,x_k]$, start with $b_n=a_n$ and work backwards by $b_k=a_k+(x-x_k)b_{k+1}$; the answer is $b_0$. The three-point example becomes $1+x[1+(x-1)]$, needing two multiplications for each evaluation. This organization is useful when adding observations incrementally or evaluating a stored Newton polynomial many times; the Runge phenomenon still depends on the interpolation nodes rather than on which representation is used.

To justify the update formula, its added product vanishes at every previous interpolation node. Choosing $c$ by the new residual forces agreement at $x_{n+1}$, so uniqueness of polynomial interpolation identifies the result. For $n\ge 0$, $p_n$ has polynomial degree at most $n$, and the added product has leading coefficient one. Thus the new coefficient $c$ is exactly $f[x_0,\dots ,x_{n+1}]$. The initial state is $p_0(x)=f(x_0)$, with the empty product interpreted as one.^[ref]

Appending $(4,57)$ to the four-point example in divided difference table leaves the Newton polynomial unchanged: $p_3(4)=57$, so the new fourth divided difference is zero. If the appended value is instead $58$, the new coefficient isThis polynomial still takes values $2,1,5,22$ at the original interpolation nodes, and now takes value $58$ at $4$. A zero update is a valid completed step: adding another exact sample of the original polynomial need not increase the polynomial degree.

Reordering the interpolation nodes changes the successive subsets used for the divided differences and the factors of the Newton polynomial. The final interpolation polynomial remains unchanged in exact arithmetic. For $f(x)=x^3-2x+1$, reversing $(-1,0,2,3)$ gives $(3,2,0,-1)$ and values $(22,5,1,2)$. The successive divided difference columns are nowFor example, the first second divided difference is $(2-17)/(0-3)=5$. The top row becomes $(22,17,5,1)$, givingThe earlier partial interpolants differ because they fit different subsets: $p_2$ fits $-1,0,2$, whereas $q_2$ fits $3,2,0$. The highest divided difference is invariant because both complete top rows finish with the same set of four interpolation nodes; this is symmetry of divided differences.

This distinction is useful when checking a divided difference table: permuting all the interpolation nodes must preserve its final divided difference, even though the intermediate entries change. It also explains why data can be appended in their arrival order. In floating-point arithmetic, the subtractions and divisions may accumulate different rounding errors in different orders; exact invariance describes the mathematical polynomial, rather than guaranteeing bit-for-bit equality of numerical results.