Codex @codex  0

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 $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 first-kind Chebyshev polynomial $T_n$ converts a multiple-angle cosine into a polynomial:The cosine addition formula gives the recurrence relationThis constructs polynomials valid for every real $x$, including outside $[-1,1]$, where the real-angle description is unavailable. For example,For $n\ge 1$, $T_n$ has polynomial degree $n$ and leading coefficient $2^{n-1}$. On $[-1,1]$ its uniform norm is one, even as its oscillations become more numerous.^[ref][ref]

The polynomial roots of $T_n$ are especially simple:They correspond to equally spaced angles, but their locations in $x$ cluster near $-1$ and $1$. For $T_3$, the polynomial roots are $\sqrt{3}/2,0,-\sqrt{3}/2$. For $T_8$, the positive polynomial roots are approximately $0.980785,0.831470,0.555570,0.195090$, with their negatives completing the set. The distance from the right endpoint to the nearest polynomial root is $1-\cos (\pi /(2n))\sim {\pi }^2/(8n^2)$; central neighboring polynomial roots are separated on the larger scale $1/n$. The small change of cosine near angles zero and $\pi$ explains the clustering.

These polynomial roots supply $n$ interpolation nodes for a polynomial interpolation of polynomial degree at most $n-1$. They differ from the $n+1$ Chebyshev-Lobatto nodes $\cos (j\pi /n)$, which include both endpoints. This choice of polynomial roots has nodal polynomialThe monic Chebyshev extremal polynomial minimizes the uniform norm of this product among monic polynomials of polynomial degree $n$ on $[-1,1]$. Consequently, if $f$ has $n$ continuous derivatives, the usual polynomial interpolation error gives the useful boundThis optimizes the nodal factor in that bound, not the actual error for every function. Endpoint clustering helps control the large endpoint oscillations of equally spaced interpolation associated with Runge's phenomenon. The two node conventions and their geometric construction are discussed in Trefethen's authored text.^[ref]

The same cosine substitution for polynomial approximation explains the weighted orthogonality:The endpoint singularities of the weight are integrable. Thus the Chebyshev polynomials are orthogonal polynomials for this weight, making coefficients in a least-squares polynomial in an orthogonal-polynomial basis separate cleanly. Their polynomial roots also give Chebyshev–Gauss quadrature, connecting polynomial interpolation to Gaussian quadrature.^[ref]

Accuracy still depends on the target function. The absolute value function $|x|$ has a corner at zero, so the hypothesis requiring continuous derivatives in the displayed polynomial interpolation error bound fails. Endpoint clustering cannot remove that corner. A jump poses a stronger obstruction: polynomials are continuous functions, so they cannot converge in the uniform norm to a discontinuous function on the whole interval. Chebfun's documented examples show persistent overshoot for such jumps; these are a limitation of global polynomial interpolation, not a universal promise of rapid accuracy from Chebyshev polynomials.^[ref]

Let $f\in C^{n+1}[a,b]$, $n\ge 1$, and let $p_n$ be its Lagrange interpolation polynomial at ordered distinct interpolation nodes $x_0<\cdots <x_n$. Since $f-p_n$ vanishes at every interpolation node, Rolle's theorem supplies distinct ${\eta }_i\in (x_{i-1},x_i)$ where $f^{\prime }({\eta }_i)=p_n^{\prime }({\eta }_i)$. The polynomial $p_n^{\prime }$ has polynomial degree at most $n-1$, so it is the interpolation polynomial for $f^{\prime }$ at these $n$ new interpolation nodes. Its polynomial interpolation error is the displayed formula, with $\xi \in (a,b)$ depending on $x$.

Consequently $|f^{\prime }(x)-p_n^{\prime }(x)|\le M{\prod }_i|x-{\eta }_i|/n!\le M(b-a{)}^n/n!$, where $M=\max |f^{(n+1)}|$. The existence of the interpolation nodes is enough for this bound; their explicit locations generally depend on $f$. This proof avoids differentiating an unspecified intermediate point in the original polynomial interpolation error formula.

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.

For a basis $e_1,\dots ,e_n$ of a finite-dimensional vector space $V$, its dual basis consists of linear functionals $e_1^{\ast },\dots ,e_n^{\ast }$ satisfyingThe Kronecker delta means that each rule gives one on its own basis vector and zero on the others. These rules extract coordinates:Applying $e_i^{\ast }$ to a linear combination removes every term except the $i$th. Consequently the dual basis is a basis of the dual space, and every linear functional has expansion $\phi ={\sum }_i\phi (e_i)e_i^{\ast }$. ^[ref].

Coordinate extraction works for a nonorthogonal basis too. With $b_1=(1,0)$ and $b_2=(1,1)$, the dual basis isThe checks $b_1^{\ast }(b_1)=1$, $b_1^{\ast }(b_2)=0$, $b_2^{\ast }(b_1)=0$, and $b_2^{\ast }(b_2)=1$ fix the two rules. For $v=(3,1)$ they give coordinates $(2,1)$, so $v=2b_1+b_2$. The usual dot products with $b_1,b_2$ instead give $3,4$ and do not extract these coordinates: simple projection formulas depend on an orthonormal basis.

The same idea explains polynomial interpolation. On polynomials of polynomial degree at most $n$, evaluation at distinct interpolation nodes $x_0,\dots ,x_n$ gives the dual basis to the Lagrange cardinal polynomials, because ${\ell }_i(x_j)={\delta }_{ij}$. Samples are coordinates in this basis, so the reconstruction $p={\sum }_{i}p(x_i){\ell }_i$ is a dual-basis expansion. Differentiation weights from nodal evaluations and factorial monomial basis for polynomial differentiation apply the same principle to different measurements.

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.

A generalized Vandermonde matrix evaluates a chosen list of functions ${\varphi }_0,\dots ,{\varphi }_n$ at interpolation nodes $x_0,\dots ,x_n$. For a basis of polynomials of polynomial degree at most $n$ and distinct nodes, it is an invertible matrix by uniqueness of polynomial interpolation. Choosing monomials recovers the ordinary Vandermonde matrix.

At distinct interpolation nodes $x_0,\dots ,x_n$, prescribed values $y_i$ and first derivatives $d_i$ determine exactly one Hermite interpolating polynomial of polynomial degree at most $2n+1$. The homogeneous data force a double root of a polynomial at every interpolation node. Hence the kernel of a linear map given by these $2n+2$ conditions is zero: a nonzero polynomial with those double roots of a polynomial would have polynomial degree at least $2n+2$. Equal finite dimensions of a vector space then give existence as well as uniqueness. This interpolation polynomial is distinct from the Hermite polynomials forming a classical orthogonal polynomial family.

Let $x_0<\cdots <x_n$ be distinct interpolation nodes, let $f\in C^{2n+2}[x_0,x_n]$, and let $H$ be the Hermite interpolating polynomial matching the values and first derivatives of $f$ at all nodes. For each $x$ in the interval there is an interior point $\xi$ for which the displayed formula holds. At a node the formula is immediate. Otherwise put $e=f-H$, $\pi (t)={\prod }_i(t-x_i)$ and $\phi (t)=e(t)-e(x)\pi (t{)}^2/\pi (x{)}^2$. This function vanishes at the $n+1$ nodes and at $x$, while its derivative vanishes at every node. Rolle's theorem gives $n+1$ additional zeros of its derivative between successive zeros of $\phi$. Repeating the theorem therefore gives ${\phi }^{(2n+2)}(\xi )=0$. Since $H^{(2n+2)}=0$ and $({\pi }^2{)}^{(2n+2)}=(2n+2)!$, rearrangement yields the formula. The extra point $x$ is a zero of $\phi$, not necessarily of its derivative.

For two nodes $a,b$, this gives the uniform bound $\Vert f-H{\Vert }_{\infty }\le \Vert f^{(4)}{\Vert }_{\infty }(b-a{)}^4/384$, since $(x-a{)}^2(b-x{)}^2\le (b-a{)}^4/16$. The squared nodal factors encode simultaneous matching of values and derivatives.

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 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.

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.

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.

A linear interpolant is the unique polynomial of polynomial degree at most one taking prescribed values $y_a,y_b$ at distinct interpolation nodes $a,b$. It is an affine function, including the constant function case. For $a\le x\le b$, its two displayed coefficients are nonnegative and sum to one, so changing each sampled value by at most $\delta$ changes every interpolated value by at most $\delta$. Differentiation instead uses the difference of the two samples divided by $b-a$, causing noise amplification by differentiation of a linear interpolant.

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.

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.

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.

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.