Hermite interpolation

This article is about polynomial interpolation on a segment using derivatives. For the family of orthogonal polynomials on the real line, see Hermite polynomials.

In numerical analysis, Hermite interpolation, named after Charles Hermite, is a method of polynomial interpolation, which generalizes Lagrange interpolation. Lagrange interpolation allows computing a polynomial of degree less than n that takes the same value at n given points as a given function. Instead, Hermite interpolation computes a polynomial of degree fewer than n such that the polynomial and its first few derivatives have the same values at m (fewer than n) given points as the given function and its first few derivatives at those points. The number of pieces of information, function values and derivative values, must add up to $n$.

Hermite's method of interpolation is closely related to the Newton's interpolation method, in that both can be derived from the calculation of divided differences. However, there are other methods for computing a Hermite interpolating polynomial. One can use linear algebra, by taking the coefficients of the interpolating polynomial as unknowns, and writing as linear equations the constraints that the interpolating polynomial must satisfy. For another method, see § Using Chinese remainder theorem. For yet another method, see § Using contour integration.

Hermite interpolating polynomials are also related to matrix generalized eigenspace projection via generalized Frobenius covariants.

Statement of the problem

In a restricted formulation, Hermite interpolation consists, given a differentiable function, of computing a polynomial function of degree as low as possible such that the two functions and their first r derivatives take the same values at m different points. This means that one starts from ⁠$n=m(r+1)$⁠ constraints, each resulting of one of the following given pairs of values: $${\begin{matrix}(x_{1},y_{1}),&(x_{2},y_{2}),&\ldots ,&(x_{m},y_{m}),\\[1ex](x_{1},y_{1}'),&(x_{2},y_{2}'),&\ldots ,&(x_{m},y_{m}'),\\[1ex]\vdots &\vdots &&\vdots \\[1.2ex](x_{1},y_{1}^{(r)}),&(x_{2},y_{2}^{(r)}),&\ldots ,&(x_{m},y_{m}^{(r)})\end{matrix}}$$ As a polynomial of degree less than n is defined by n coefficients, the resulting polynomial must have a degree less than n, since, otherwise, the solution could not be unique. In a more general case, there is no need for r to have a fixed value; that is, some points may have more known derivatives than others. In this case the resulting polynomial has a degree less than the total number of constraints n.

A theorem, commonly known as the Hermite interpolation theorem, asserts the existence and uniqueness of the resulting polynomial.

Let us consider a polynomial ⁠$h(x)$⁠ of degree less than n with indeterminate coefficients; that is, the coefficients of ⁠$h(x)$⁠ are n new variables. Then, by writing the constraints that the interpolating polynomial must satisfy, one gets a system of n linear equations in n unknowns.

Hermite used contour integration to prove that this is effectively the case here, and to find the unique solution, provided that the ⁠$x_{i}$⁠ are pairwise different. Hermite interpolation is a problem of linear algebra that has the coefficients of the interpolation polynomial as unknown variables and a confluent Vandermonde matrix as its matrix. The general methods of linear algebra, and specific methods for confluent Vandermonde matrices are often used for computing the interpolation polynomial, as detailed in § Using linear algebra. Other methods are described below.

Explicit formula

For a given set of m distinct nodes ⁠$x_{1},\dots ,x_{m}$⁠ with associated maximum derivative orders ⁠$r_{1},\dots ,r_{m}$⁠, let ⁠$n_{i}=r_{i}+1$⁠ be the multiplicity of node ⁠$x_{i}$⁠, and let ⁠$\textstyle n=\sum _{i=1}^{m}n_{i}$⁠ be the total number of conditions. Given a target function ⁠$f$⁠, the Hermite interpolating polynomial ⁠$h$⁠ is the unique polynomial of degree at most ⁠$n-1$⁠ that satisfies the interpolation conditions: $$h^{(j)}(x_{i})=f^{(j)}(x_{i})=y_{i,j}\quad {\text{for }}i=1,\dots ,m{\text{ and }}j=0,\dots ,r_{i}.$$

The explicit closed-form formula for the Hermite interpolating polynomial ⁠$h(x)$⁠ is given by: $$h(x)=\sum _{i=1}^{m}\underbrace {\left[\sum _{j=0}^{r_{i}}{\frac {1}{j!}}\left.\left({\frac {f(\xi )}{q_{i}(\xi )}}\right)^{(j)}\right|_{\xi =x_{i}}(x-x_{i})^{j}\right]} _{{\mathcal {T}}_{x_{i}}^{r_{i}}\left[f(x)/q_{i}(x)\right]}q_{i}(x),$$ where ⁠$\textstyle q_{i}(x)=\prod _{j\neq i}(x-x_{j})^{n_{j}}$⁠ and ⁠$\textstyle {\mathcal {T}}_{a}^{k}[f(x)]$⁠ denotes the kth degree Taylor operator of a function ⁠$f(x)$⁠ centered at ⁠$x=a$⁠.

When only the prescribed interpolation conditions ⁠$y_{i,j}$⁠ are given, for each node ⁠$x_{i}$⁠, let ⁠$g_{i}(x)$⁠ be the prescribed local ⁠$r_{i}$⁠th degree Taylor polynomial of ⁠$h$⁠ at ⁠$x=x_{i}$⁠: $$g_{i}(x)={\mathcal {T}}_{x_{i}}^{r_{i}}[h(x)]=\sum _{j=0}^{r_{i}}{\frac {y_{i,j}}{j!}}(x-x_{i})^{j}.$$ Using the Taylor composition theorem, ⁠$f(x)$⁠ can be directly substituted with the prescribed local ⁠$r_{i}$⁠th degree Taylor polynomial ⁠$g_{i}(x)$⁠ inside any ⁠$r_{i}$⁠th degree Taylor operators at ⁠$x=x_{i}$⁠, giving the explicit closed-form formula in terms of the prescribed interpolation conditions: $$h(x)=\sum _{i=1}^{m}{\mathcal {T}}_{x_{i}}^{r_{i}}\left[{\frac {g_{i}(x)}{q_{i}(x)}}\right]q_{i}(x).$$

When ⁠$r_{i}=0$⁠ for all distinct nodes ⁠$x_{i}$⁠, meaning no derivative conditions are specified, this explicit closed-form formula reduces directly to the explicit formula of Lagrange polynomials.

Basis polynomials in closed form

The Hermite interpolating polynomial ⁠$h(x)$⁠ can also be represented as a linear combination of Hermite basis polynomials ⁠$h_{i,j}(x)$⁠: $$h(x)=\sum _{i=1}^{m}\sum _{j=0}^{r_{i}}y_{i,j}h_{i,j}(x).$$

An explicit closed-form expression for each basis polynomial ⁠$h_{i,j}(x)$⁠ can be obtained directly by prescribing ⁠$y_{k,l}=\delta _{ik}\delta _{jl}$⁠ for ⁠$k=1,\dots ,m$⁠ and ⁠$l=0,\dots ,r_{k}$⁠, where ⁠$\delta _{ab}$⁠ denotes the Kronecker delta, and constructing the corresponding prescribed local Taylor polynomials: $${\begin{aligned}g_{k}(x)&={(x-x_{i})}^{j}/j!&&{\text{for }}k=i{\text{ and}}\\g_{k}(x)&=0&&{\text{for }}k\neq i.\end{aligned}}$$ Substituting ⁠$g_{k}(x)$⁠ into the explicit formula yields: $$h_{i,j}(x)={\mathcal {T}}_{x_{i}}^{r_{i}}\left[{\frac {1}{j!}}{\frac {(x-x_{i})^{j}}{q_{i}(x)}}\right]q_{i}(x).$$

Following the § Alternative basis construction via degree-0 basis polynomials, the full set of basis polynomials can be constructed entirely using only the degree-0 basis polynomial ⁠$h_{i,0}(x)$⁠. For ⁠$h_{i,0}(x)$⁠, the local Taylor polynomials simplify to ⁠$g_{k}(x)=\delta _{ik}$⁠, giving the explicit closed-form expression: $$h_{i,0}(x)={\mathcal {T}}_{x_{i}}^{r_{i}}\left[{\frac {1}{q_{i}(x)}}\right]q_{i}(x).$$

An alternative recursively-defined explicit basis polynomial expression is given by: $$h_{i,j}(x)={\frac {(x-x_{i})^{j}}{j!}}q_{i}(x)-\sum _{k=j+1}^{r_{i}}{\binom {k}{j}}q_{i}^{(k-j)}(x_{i})h_{i,k}(x).$$

Using Chinese remainder theorem

See also: Chinese remainder theorem § Hermite interpolation

Let m be a positive integer representing the number of distinct nodes, ⁠$r_{1},\ldots ,r_{m}$⁠ be nonnegative integers specifying the maximum derivative order at each node, ⁠$n_{i}=r_{i}+1$⁠ be the multiplicity of node ⁠$x_{i}$⁠, and the nodes ⁠$x_{1},\ldots ,x_{m}$⁠ be real numbers or elements of any other field of characteristic zero. The Hermite interpolation problem consists of finding a polynomial ⁠$h(x)$⁠ such that $$h(x_{i})=y_{i,0},h'(x_{i})=y_{i,1},\ldots ,h^{(r_{i})}(x_{i})=y_{i,r_{i}}$$ for ⁠$i=1,\ldots ,m$⁠, where the ⁠$y_{i,j}$⁠ are given values in the same field as the ⁠$x_{i}$⁠.

These conditions implies that the Taylor polynomial of ⁠$h(x)$⁠ of degree ⁠$r_{i}$⁠ at ⁠$x_{i}$⁠ is $$g_{i}(x)=\sum _{j=0}^{r_{i}}{\frac {y_{i,j}}{j!}}(x-x_{i})^{j}.$$ In other words, the desired polynomial ⁠$h(x)$⁠ is congruent to this polynomial ⁠$g_{i}(x)$⁠ modulo ⁠$\textstyle p_{i}(x)=(x-x_{i})^{n_{i}}$⁠.

The Chinese remainder theorem for polynomials implies that there is exactly one solution of degree less than ⁠$\textstyle n=\sum _{i=1}^{m}n_{i}$⁠.

Moreover, this solution can be computed with ⁠$\textstyle O(n^{2})$⁠ arithmetic operations, or even faster with fast polynomial multiplication.

This approach does not work in positive characteristic, because of the denominators of the coefficients of the Taylor polynomial. The approach through divided differences, below, works in every characteristic.

Explicit construction via Chinese remainder theorem

For general Hermite interpolation, an explicit solution satisfying the system of simultaneous congruences can be constructed directly by following the constructive proof of the Chinese remainder theorem.

The target polynomial ⁠$h(x)$⁠ must satisfy the system of modular equivalences: $$h(x)\equiv g_{i}(x){\pmod {p_{i}(x)}}\quad {\text{for }}i=1,\dots ,m,$$ where ⁠$\textstyle p_{i}(x)=(x-x_{i})^{n_{i}}$⁠ is the modulus polynomial at node ⁠$x_{i}$⁠.

The solution is constructed iteratively by building a sequence of polynomials ⁠$a_{1}(x),a_{2}(x),\dots ,a_{m}(x)$⁠ such that ⁠$a_{i}(x)$⁠ satisfies the first i congruences.

1. Base Step: Set the initial partial solution for ⁠$i=1$⁠: $$a_{1}(x)=g_{1}(x).$$

2. Inductive Step: Define ⁠$\textstyle b_{i}(x)=\prod _{j=1}^{i}p_{j}(x)$⁠ for ⁠$i=1,\dots ,m$⁠. For ⁠$i=2,\dots ,m$⁠, since ⁠$b_{i-1}(x)$⁠ and ⁠$p_{i}(x)$⁠ share no roots and thus are coprime. By Bézout's identity, there exist polynomials ⁠$u_{i}(x)$⁠ and ⁠$v_{i}(x)$⁠ satisfying: $$u_{i}(x)b_{i-1}(x)+v_{i}(x)p_{i}(x)=1.$$

The polynomials ⁠$u_{i}(x)$⁠ and ⁠$v_{i}(x)$⁠ can be solved together using the polynomial extended Euclidean algorithm.

Define the ith partial solution ⁠$a_{i}(x)$⁠ as: $$a_{i}(x)=g_{i}(x)u_{i}(x)b_{i-1}(x)+a_{i-1}(x)v_{i}(x)p_{i}(x).$$

This construction satisfies all conditions up to node i:

• Modulo ⁠$p_{i}(x)$⁠, the term containing ⁠$p_{i}(x)$⁠ vanishes while ⁠$u_{i}(x)b_{i-1}(x)\equiv 1{\pmod {p_{i}(x)}}$⁠, giving ⁠$a_{i}(x)\equiv g_{i}(x){\pmod {p_{i}(x)}}$⁠.
• Modulo ⁠$p_{j}(x)$⁠ for any ⁠$j<i$⁠, ⁠$b_{i-1}(x)\equiv 0{\pmod {p_{j}(x)}}$⁠ and ⁠$v_{i}(x)p_{i}(x)\equiv 1{\pmod {p_{j}(x)}}$⁠, giving ⁠$a_{i}(x)\equiv a_{i-1}(x)\equiv g_{j}(x){\pmod {p_{j}(x)}}$⁠.

After ⁠$m-1$⁠ inductive steps, the final interpolating polynomial ⁠$h(x)$⁠ is obtained by reducing ⁠$a_{m}(x)$⁠ modulo the total degree product polynomial ⁠$\textstyle b_{m}(x)=\prod _{i=1}^{m}p_{i}(x)$⁠: $$h(x)=a_{m}(x){\bmod {b}}_{m}(x).$$

Basis polynomials via Chinese remainder theorem

Following § Alternative basis construction via degree-0 basis polynomials, a full basis can be constructed by only solving the degree-0 basis polynomials ⁠$h_{i,0}(x)$⁠, which can be obtained in a simplified manner.

For each node ⁠$x_{i}$⁠, the degree-0 basis polynomial ⁠$h_{i,0}(x)$⁠ must satisfy:

• ⁠$h_{i,0}(x)\equiv 0{\pmod {p_{j}(x)}}\quad {\text{for }}j\neq i$⁠,
• ⁠$h_{i,0}(x)\equiv 1{\pmod {p_{i}(x)}}$⁠.

To satisfy the zero-remainder conditions at all other nodes ⁠$j\neq i$⁠, ⁠$h_{i,0}(x)$⁠ must contain the product factor ⁠$\textstyle q_{i}(x)=\prod _{j\neq i}p_{j}(x)$⁠. Thus, ⁠$h_{i,0}(x)=u_{i}(x)q_{i}(x)$⁠ for some polynomial ⁠$u_{i}(x)$⁠.

Since ⁠$p_{i}(x)$⁠ and ⁠$q_{i}(x)$⁠ share no roots and thus are coprime, Bézout's identity guarantees the existence of polynomials ⁠$u_{i}(x)$⁠ and ⁠$v_{i}(x)$⁠ satisfying: $$u_{i}(x)q_{i}(x)+v_{i}(x)p_{i}(x)=1.$$

Taking this relation modulo ⁠$p_{i}(x)$⁠, the term ⁠$v_{i}(x)p_{i}(x)$⁠ vanishes, showing that ⁠$u_{i}(x)$⁠ is a modular multiplicative inverse of ⁠$q_{i}(x)$⁠: $$u_{i}(x)q_{i}(x)\equiv 1{\pmod {p_{i}(x)}}.$$

Since the polynomial ⁠$v_{i}(x)$⁠ is not required for the construction, ⁠$u_{i}(x)$⁠ can be computed directly using half extended Euclidean algorithm that tracks only one set of Bézout coefficients.

Because the degree of the modular inverse ⁠$u_{i}(x)$⁠ obtained via the extended Euclidean algorithm is strictly less than ⁠$\deg(p_{i})=n_{i}$⁠, and ⁠$\deg(q_{i})=n-n_{i}$⁠, the product degree is bounded by: $$\deg(u_{i}q_{i})=\deg(u_{i})+\deg(q_{i})<n_{i}+(n-n_{i})=n.$$ Thus, ⁠$\deg(u_{i}q_{i})\leq n-1$⁠, meaning the product ⁠$u_{i}(x)q_{i}(x)$⁠ already satisfies the degree constraint without requiring any further modular reduction, yielding directly the unique degree-0 basis polynomial ⁠$h_{i,0}(x)=u_{i}(x)q_{i}(x)$⁠.

For any individual simple node ⁠$x_{i}$⁠ with no derivative conditions (⁠$r_{i}=0$⁠), ⁠$p_{i}(x)=x-x_{i}$⁠. By the polynomial remainder theorem, the modular condition ⁠$u_{i}(x)q_{i}(x)\equiv 1{\pmod {x-x_{i}}}$⁠ simplifies to the point evaluation ⁠$u_{i}(x_{i})q_{i}(x_{i})=1$⁠. In this case, ⁠$u_{i}(x)$⁠ reduces to a scalar constant that can be evaluated directly in closed form: $$u_{i}(x)={\frac {1}{q_{i}(x_{i})}}=\prod _{j\neq i}{\frac {1}{p_{j}(x_{i})}}.$$

Using divided differences

See also: Divided differences and Newton polynomial

For a given set of m distinct nodes ⁠$x_{1},\dots ,x_{m}$⁠ with associated maximum derivative orders ⁠$r_{1},\dots ,r_{m}$⁠, let ⁠$n_{i}=r_{i}+1$⁠ be the multiplicity of node ⁠$x_{i}$⁠, and let ⁠$\textstyle n=\sum _{i=1}^{m}n_{i}$⁠ be the total number of conditions. An efficient way to construct the Hermite interpolating polynomial ⁠$h(x)$⁠ is by extending Newton's divided difference method. Since Newton interpolation and Lagrange interpolation construct the exact same unique polynomial for a distinct set of nodes, Newton's method can be generalized to the confluent case where distinct nodes coalesce via confluent divided differences. To interpolate a function ⁠$f(x)$⁠ or a given set of derivative conditions ⁠$y_{i,j}$⁠ (for ⁠$i=1,\dots ,m$⁠ and ⁠$j=0,\dots ,r_{i}$⁠), the original set of nodes is transformed into a dataset ⁠$z_{1},\dots ,z_{n}$⁠ by repeating each node ⁠$x_{i}$⁠ according to its multiplicity ⁠$n_{i}$⁠ times: $$(z_{1},\dots ,z_{n})=(\underbrace {x_{1},\dots ,x_{1}} _{n_{1}{\text{ times}}},\underbrace {x_{2},\dots ,x_{2}} _{n_{2}{\text{ times}}},\dots ,\underbrace {x_{m},\dots ,x_{m}} _{n_{m}{\text{ times}}}).$$ When constructing the divided difference table for this dataset, the standard finite difference quotient for identical confluent nodes would result in an undefined ⁠$0/0$⁠ expression; these entries are directly replaced by the corresponding scaled derivative values ⁠$\textstyle f^{(j)}(x_{i})/j!$⁠ (or ⁠$y_{i,j}/j!$⁠). The interpolating polynomial is then represented in Newton form: $$h(x)=\sum _{i=1}^{n}f[z_{1},\dots ,z_{i}]\prod _{j=1}^{i-1}(x-z_{j}).$$ Once the target function and derivative values are known, this generalized construction requires only basic algebraic operations (addition, subtraction, multiplication, and division) to compute the divided difference table and construct the polynomial.

Through Cauchy's integral formula for confluent divided differences, constructing the Hermite interpolating polynomial via divided differences can also be viewed as a contour integration approach, providing an alternative formulation to the integral representation described in § Using contour integral.

Simple case with only first derivatives

Here we will consider only the simplest case ⁠$r_{i}=1$⁠ (i.e., ⁠$n_{i}=2$⁠) for all distinct nodes. Therefore, given m distinct nodes ⁠$x_{1},\ldots ,x_{m}$⁠, and values ⁠$f(x_{1}),\ldots ,f(x_{m})$⁠ and ⁠$f'(x_{1}),\ldots ,f'(x_{m})$⁠ for a function ⁠$f(x)$⁠ that we want to interpolate, we create a new dataset ⁠$z_{1},\ldots ,z_{2m}$⁠ such that $$z_{2i-1}=z_{2i}=x_{i}.$$

Now, we create a divided differences table for the nodes ⁠$z_{1},\ldots ,z_{2m}$⁠. However, for some divided differences, $$z_{i}=z_{i+1}\implies f[z_{i},z_{i+1}]={\frac {f(z_{i+1})-f(z_{i})}{z_{i+1}-z_{i}}}={\frac {0}{0}}$$ which is undefined. In this case, the divided difference is replaced by ⁠$f'(z_{i})$⁠. All others are calculated normally.

General case for higher-order derivatives

In the general case, suppose a given node ⁠$x_{i}$⁠ has up to ⁠$r_{i}$⁠th order derivatives. Then the dataset ⁠$z_{1},\ldots ,z_{n}$⁠ contains ⁠$n_{i}$⁠ adjacent identical copies of ⁠$x_{i}$⁠. When creating the table, divided differences of ⁠$j+1$⁠ (for ⁠$j=0,\ldots ,r_{i}$⁠) identical values will be calculated as $$f[\underbrace {x_{i},\ldots ,x_{i}} _{j+1\,{\text{times}}}]={\frac {f^{(j)}(x_{i})}{j!}}.$$

For example, $$f[x_{i},x_{i},x_{i}]={\frac {f''(x_{i})}{2}}$$ $$f[x_{i},x_{i},x_{i},x_{i}]={\frac {f^{(3)}(x_{i})}{6}}$$ etc.

A fast algorithm for the fully general case is given by Schneider and Werner. A slower but more numerically stable algorithm is given by Corless and Fillion.

Example

Case of uniform derivative orders

Consider the function ⁠$\textstyle f(x)=x^{8}+1$⁠. Evaluating the function and its first two derivatives at ⁠$x\in \{-1,0,1\}$⁠, we obtain the following data:

⁠$x$⁠	⁠$f(x)$⁠	⁠$f'(x)$⁠	⁠$f''(x)$⁠
−1	2	−8	56
0	1	0	0
1	2	8	56

Since we have two derivatives to work with, we construct the dataset $$(z_{i},\ldots ,z_{n})=(-1,-1,-1,0,0,0,1,1,1).$$ Our divided difference table is then: $${\begin{array}{llcclrrrrr}z_{1}=-1&f[z_{1}]=2&&&&&&&&\\&&{\frac {f'(z_{1})}{1}}=-8&&&&&&&\\z_{2}=-1&f[z_{2}]=2&&{\frac {f''(z_{2})}{2}}=28&&&&&&\\&&{\frac {f'(z_{2})}{1}}=-8&&f[z_{4},z_{3},z_{2},z_{1}]=-21&&&&&\\z_{3}=-1&f[z_{3}]=2&&f[z_{4},z_{3},z_{2}]=7&&15&&&&\\&&f[z_{4},z_{3}]=-1&&f[z_{5},z_{4},z_{3},z_{2}]=-6&&-10&&&\\z_{4}=0&f[z_{4}]=1&&f[z_{5},z_{4},z_{3}]=1&&5&&4&&\\&&{\frac {f'(z_{4})}{1}}=0&&f[z_{6},z_{5},z_{4},z_{3}]=-1&&-2&&-1&\\z_{5}=0&f[z_{5}]=1&&{\frac {f''(z_{5})}{2}}=0&&1&&2&&1\\&&{\frac {f'(z_{5})}{1}}=0&&f[z_{7},z_{6},z_{5},z_{4}]=1&&2&&1&\\z_{6}=0&f[z_{6}]=1&&f[z_{7},z_{6},z_{5}]=1&&5&&4&&\\&&f[z_{7},z_{6}]=1&&f[z_{8},z_{7},z_{6},z_{5}]=6&&10&&&\\z_{7}=1&f[z_{7}]=2&&f[z_{8},z_{7},z_{6}]=7&&15&&&&\\&&{\frac {f'(z_{7})}{1}}=8&&f[z_{9},z_{8},z_{7},z_{6}]=21&&&&&\\z_{8}=1&f[z_{8}]=2&&{\frac {f''(z_{8})}{2}}=28&&&&&&\\&&{\frac {f'(z_{8})}{1}}=8&&&&&&&\\z_{9}=1&f[z_{9}]=2&&&&&&&&\\\end{array}}$$ and by taking the coefficients from the diagonal of the divided difference table, and multiplying the ith coefficient by ⁠$\textstyle \prod _{j=1}^{i-1}(x-z_{j})$⁠, as we would when constructing Newton form to generate a Newton polynomial, the generated polynomial ⁠$h(x)$⁠ is $${\begin{aligned}h(x)&=2-8(x+1)+28(x+1)^{2}-21(x+1)^{3}+15x(x+1)^{3}-10x^{2}(x+1)^{3}\\&\quad {}+4x^{3}(x+1)^{3}-1x^{3}(x+1)^{3}(x-1)+x^{3}(x+1)^{3}(x-1)^{2}\\&=2-8+28-21-8x+56x-63x+15x+28x^{2}-63x^{2}+45x^{2}-10x^{2}-21x^{3}\\&\quad {}+45x^{3}-30x^{3}+4x^{3}+x^{3}+x^{3}+15x^{4}-30x^{4}+12x^{4}+2x^{4}+x^{4}\\&\quad {}-10x^{5}+12x^{5}-2x^{5}+4x^{5}-2x^{5}-2x^{5}-x^{6}+x^{6}-x^{7}+x^{7}+x^{8}\\&=x^{8}+1.\end{aligned}}$$

Case of non-uniform derivative orders

Consider the following set of derivative conditions ⁠$y_{i,j}$⁠ at nodes ⁠$x_{i}\in \{-1,0,1\}$⁠:

⁠$i$⁠	⁠$x_{i}$⁠	⁠$y_{i,0}$⁠	⁠$y_{i,1}$⁠
1	−1	1	2
2	0	1
3	1	-1	-2

For this example with non-uniform derivative orders, we construct the dataset $$(z_{1},\ldots ,z_{n})=(-1,-1,0,1,1).$$ Denoting the divided differences as ⁠$y[z_{a},\dots ,z_{b}]$⁠ when constructed from the derivative conditions ⁠$y_{i,j}$⁠, the constructed divided difference table is given by: $${\begin{array}{llclrr}z_{1}=-1&y[z_{1}]=y_{1,0}=1&&&&\\&&{\frac {y_{1,1}}{1}}=2&&&\\z_{2}=-1&y[z_{2}]=y_{1,0}=1&&y[z_{3},z_{2},z_{1}]=-2&&\\&&y[z_{3},z_{2}]={\phantom {-}}0&&{\frac {1}{2}}&\\z_{3}=0&y[z_{3}]=y_{2,0}=1&&y[z_{4},z_{3},z_{2}]=-1&&0\\&&y[z_{4},z_{3}]=-2&&{\frac {1}{2}}&\\z_{4}=1&y[z_{4}]=y_{3,0}=-1&&y[z_{5},z_{4},z_{3}]=0&&\\&&{\frac {y_{3,1}}{1}}=-2&&&\\z_{5}=1&y[z_{5}]=y_{3,0}=-1&&&&\\\end{array}}$$ and the generated polynomial ⁠$h(x)$⁠ in Newton form is $${\begin{aligned}h(x)&=1+2(x+1)-2(x+1)^{2}+{\frac {1}{2}}x(x+1)^{2}+0x(x+1)^{2}(x-1)\\&={\frac {1}{2}}x^{3}-x^{2}-{\frac {3}{2}}x+1.\end{aligned}}$$

Using linear algebra

See also: Vandermonde matrix, Confluent Vandermonde matrix, and Lagrange polynomial § A perspective from linear algebra

The interpolation problem can also be formulated and solved directly using linear algebra in a similar way to Lagrange interpolation, by replacing the standard Vandermonde matrix with a confluent Vandermonde matrix to account for both function values and derivative conditions at each node.

Let m be the number of distinct nodes ⁠$x_{1},\ldots ,x_{m}$⁠, let ⁠$n_{i}$⁠ be the multiplicity of node ⁠$x_{i}$⁠ representing the number of specified values including derivatives, and let ⁠$\textstyle n=\sum _{i=1}^{m}n_{i}$⁠ be the total number of data conditions. Let ⁠$\mathbf {V}$⁠ be the constructed confluent Vandermonde matrix, let ⁠$\textstyle \mathbf {H} =\mathbf {V} ^{-1}$⁠ be its inverse, and let ⁠$d=n-1$⁠ be the degree bound of the interpolating polynomial. The matrix ⁠$\mathbf {H}$⁠ is given by: $$\mathbf {H} ={\begin{pmatrix}H_{00}&\cdots &H_{0d}\\\vdots &\ddots &\vdots \\H_{d0}&\cdots &H_{dd}\end{pmatrix}}.$$

The sample vector ⁠$\mathbf {y}$⁠ stacks all target function values and derivative constraints into a single column vector of total dimension n. The pair of node index i and derivative index j is mapped to a single flattened index ⁠$k\in \{0,\dots ,d\}$⁠ via the indexing function: $$k=\sigma (i,j)=j+\sum _{l=1}^{i-1}n_{l}.$$

Given the interpolation condition vector ⁠$\mathbf {y}$⁠, the coefficient vector ⁠$\textstyle \mathbf {a} ={\begin{pmatrix}a_{0}&a_{1}&\cdots &a_{d}\end{pmatrix}}^{T}$⁠ of the Hermite interpolating polynomial ⁠$\textstyle h(x)=\sum _{k=0}^{d}a_{k}x^{k}$⁠ can be computed as: $$\mathbf {a} =\mathbf {H} \mathbf {y} .$$

Furthermore, the individual Hermite interpolating basis polynomials ⁠$h_{i,j}(x)$⁠ correspond to the basis polynomials ⁠$h_{k}(x)$⁠ in the flattened indexing under ⁠$k=\sigma (i,j)$⁠, allowing the interpolating polynomial to be expressed equivalently as: $$h(x)=\sum _{k=0}^{d}y_{k}h_{k}(x)=\sum _{i=1}^{m}\sum _{j=0}^{n_{i}-1}y_{\sigma (i,j)}h_{i,j}(x).$$ The basis polynomials ⁠$h_{i,j}(x)=h_{k}(x)$⁠ can be constructed directly from the columns of ⁠$\mathbf {H}$⁠: $$h_{i,j}(x)=h_{k}(x)=H_{0k}+H_{1k}x+\cdots +H_{dk}x^{d}.$$

In numerical practice, directly computing ⁠$\mathbf {H}$⁠ by inverting the confluent Vandermonde matrix is generally avoided because confluent Vandermonde matrices are notoriously ill-conditioned. As the number of nodes or their degree increases, their condition number grows exponentially, making direct matrix inversion numerically unstable and highly sensitive to floating-point errors.

Using contour integration

An elegant representation and derivation of the Hermite interpolating polynomial can be obtained using contour integration in the complex plane, a formulation originally introduced by Charles Hermite in 1878.

Let m be the number of distinct interpolation nodes ⁠$x_{1},\dots ,x_{m}$⁠, with corresponding multiplicities ⁠$n_{1},\dots ,n_{m}$⁠. The total number of interpolation conditions is ⁠$\textstyle n=\sum _{i=1}^{m}n_{i}$⁠. Define the product polynomial of all node factors as $$q(x)=\prod _{i=1}^{m}(x-x_{i})^{n_{i}}.$$

Suppose the target function f is analytic inside and on a positively oriented simple closed contour ⁠$\Gamma$⁠ in the complex plane that encloses the evaluation point ⁠$x$⁠ as well as all nodes ⁠$x_{i}$⁠. The Hermite interpolating polynomial ⁠$h(x)$⁠ of degree at most ⁠$n-1$⁠ can be expressed as $$h(x)=f(x)-{\frac {q(x)}{2\pi i}}\oint _{\Gamma }{\frac {f(\zeta )}{q(\zeta )(\zeta -x)}}\,d\zeta .$$

Alternatively, by expressing ⁠$f(x)$⁠ via Cauchy's integral formula as $$f(x)={\frac {1}{2\pi i}}\oint _{\Gamma }{\frac {f(\zeta )}{\zeta -x}}\,d\zeta$$ and absorbing it into the integral, ⁠$h(x)$⁠ can be written equivalently as $$h(x)={\frac {1}{2\pi i}}\oint _{\Gamma }{\frac {f(\zeta )\left(q(\zeta )-q(x)\right)}{q(\zeta )(\zeta -x)}}\,d\zeta .$$

Derivation using residues

Consider the meromorphic function $$F(\zeta )={\frac {f(\zeta )}{q(\zeta )(\zeta -x)}}.$$

Assuming ⁠$x$⁠ is not one of the nodes ⁠$x_{i}$⁠, ⁠$F(\zeta )$⁠ has the following poles inside ⁠$\Gamma$⁠:

• A simple pole at ⁠$\zeta =x$⁠.
• Poles at ⁠$\zeta =x_{i}$⁠ with order ⁠$n_{i}$⁠ for each ⁠$i=1,\dots ,m$⁠.

By the residue theorem, integrating ⁠$F(\zeta )$⁠ along ⁠$\Gamma$⁠ yields $${\frac {1}{2\pi i}}\oint _{\Gamma }{\frac {f(\zeta )}{q(\zeta )(\zeta -x)}}\,d\zeta =\operatorname {Res} _{\zeta =x}F(\zeta )+\sum _{i=1}^{m}\operatorname {Res} _{\zeta =x_{i}}F(\zeta ).$$

The residue at the simple pole ⁠$\zeta =x$⁠ is easily evaluated as $$\operatorname {Res} _{\zeta =x}F(\zeta )=\lim _{\zeta \to x}\,(\zeta -x){\frac {f(\zeta )}{q(\zeta )(\zeta -x)}}={\frac {f(x)}{q(x)}}.$$

Substituting this residue back into the integral equation gives $${\frac {1}{2\pi i}}\oint _{\Gamma }{\frac {f(\zeta )}{q(\zeta )(\zeta -x)}}\,d\zeta ={\frac {f(x)}{q(x)}}+\sum _{i=1}^{m}\operatorname {Res} _{\zeta =x_{i}}F(\zeta ).$$

Rearranging terms to isolate ⁠$f(x)$⁠ gives $$f(x)={\frac {q(x)}{2\pi i}}\oint _{\Gamma }{\frac {f(\zeta )}{q(\zeta )(\zeta -x)}}\,d\zeta -q(x)\sum _{i=1}^{m}\operatorname {Res} _{\zeta =x_{i}}F(\zeta ).$$

Each residue ⁠$\operatorname {Res} _{\zeta =x_{i}}F(\zeta )$⁠ is a rational function in ⁠$x$⁠ whose denominator factors are cancelled out by ⁠$q(x)$⁠. Summing these terms yields a polynomial in ⁠$x$⁠ of degree at most ⁠$n-1$⁠ that satisfies all interpolation conditions. Thus, it corresponds precisely to the Hermite interpolating polynomial ⁠$h(x)$⁠: $$h(x)=-q(x)\sum _{i=1}^{m}\operatorname {Res} _{\zeta =x_{i}}F(\zeta ).$$ Substituting this definition into the expression for ⁠$f(x)$⁠ and rearrange yields $$h(x)=f(x)-{\frac {q(x)}{2\pi i}}\oint _{\Gamma }{\frac {f(\zeta )}{q(\zeta )(\zeta -x)}}\,d\zeta .$$

Consequently, the interpolation error ⁠$f(x)-h(x)$⁠ (see also § Error) is given directly by the integral term: $$f(x)-h(x)={\frac {q(x)}{2\pi i}}\oint _{\Gamma }{\frac {f(\zeta )}{q(\zeta )(\zeta -x)}}\,d\zeta .$$

Using orthogonal projection

The Hermite interpolating polynomial can also be constructed via an orthogonal projection onto the space of polynomials of degree at most ⁠$n-1$⁠. By defining a suitable inner product based on the nodes and their derivative conditions, an orthogonal polynomial basis can be generated, allowing the interpolant to be expressed as a series of orthogonal projections.

Let ⁠$x_{1},\dots ,x_{m}$⁠ be the distinct interpolation nodes, let ⁠$r_{i}$⁠ denote the maximum derivative order prescribed at node ⁠$x_{i}$⁠, and let ⁠$n_{i}=r_{i}+1$⁠ be the multiplicity of the node ⁠$x_{i}$⁠. The total number of conditions is ⁠$\textstyle n=\sum _{i=1}^{m}n_{i}$⁠, and the Hermite interpolating polynomial ⁠$h$⁠ belongs to the space of polynomials of degree at most ⁠$n-1$⁠.

We define an inner product on this polynomial space by: $$\langle a,b\rangle =\sum _{i=1}^{m}\sum _{j=0}^{r_{i}}a^{(j)}(x_{i})b^{(j)}(x_{i}),$$ with ⁠$\textstyle \lVert a\rVert ={\sqrt {\langle a,a\rangle }}$⁠ denoting the induced norm.

An orthogonal basis ⁠$\textstyle \{h_{k}\}_{k=1}^{n}$⁠ for the polynomial space can be constructed using the Gram-Schmidt process. Starting with ⁠$h_{1}(x)=1$⁠, subsequent basis polynomials are constructed recursively for ⁠$k=1,\dots ,n-1$⁠ by: $$h_{k+1}(x)=xh_{k}(x)-\sum _{l=1}^{k}{\frac {\langle xh_{k}(x),h_{l}\rangle }{{\phantom {^{2}}}\lVert h_{l}\rVert ^{2}}}h_{l}(x)$$

The Hermite interpolating polynomial ⁠$h$⁠ is given by the orthogonal projection of the target function ⁠$f$⁠ onto the basis: $$h(x)=\sum _{k=1}^{n}{\frac {\langle f,h_{k}\rangle }{{\phantom {^{2}}}\lVert h_{k}\rVert ^{2}}}h_{k}(x),$$ or using the normalized orthonormal basis polynomials ⁠$\textstyle {\hat {h}}_{k}=h_{k}/\lVert h_{k}\rVert$⁠, $$h(x)=\sum _{k=1}^{n}\langle f,{\hat {h}}_{k}\rangle {\hat {h}}_{k}(x).$$

If only the prescribed derivative conditions ⁠$y_{i,j}$⁠ are given instead of an explicit target function ⁠$f$⁠, the inner product terms ⁠$\langle f,g\rangle$⁠ can be evaluated directly using ⁠$\textstyle y_{i,j}=f^{(j)}(x_{i})$⁠ via: $$\langle f,g\rangle =\sum _{i=1}^{m}\sum _{j=0}^{r_{i}}y_{i,j}g^{(j)}(x_{i})$$

Common low-degree cases

Cubic Hermite spline

Main article: Cubic Hermite spline

Cubic Hermite splines are widely used in numerical analysis, computer graphics, and motion planning for keyframe animation, path trajectories, and geometric modeling because they ensure smooth ⁠$\textstyle C^{1}$⁠ continuity while providing precise local control over both endpoint positions and tangents.

On the unit interval ⁠$x\in [0,1]$⁠, the two endpoints are defined as ⁠$x_{1}=0$⁠ and ⁠$x_{2}=1$⁠, the cubic Hermite interpolant ⁠$h(x)$⁠ is uniquely determined by boundary values and first derivatives ⁠$f(0),f(1),f'(0),f'(1)$⁠: $$h(x)=f(0)h_{1,0}(x)+f'(0)h_{1,1}(x)+f(1)h_{2,0}(x)+f'(1)h_{2,1}(x)$$

where the four basis polynomials are: $${\begin{aligned}h_{1,0}(x)&=2x^{3}-3x^{2}+1,\\h_{1,1}(x)&=x^{3}-2x^{2}+x,\\h_{2,0}(x)&=-2x^{3}+3x^{2},\\h_{2,1}(x)&=x^{3}-x^{2}.\end{aligned}}$$

For rescaling to arbitrary intervals, piecewise multi-interval construction, or extension to 2D/3D parametric curves, see Cubic Hermite spline.

Quintic Hermite interpolation

The quintic Hermite interpolation based on the function ⁠$f(x)$⁠, its first ⁠$f'(x)$⁠ and second derivatives ⁠$f''(x)$⁠ at two different points (⁠$x_{1}$⁠ and ⁠$x_{2}$⁠) can be used for example to interpolate the position of an object based on its position, velocity and acceleration. The general form is given by $${\begin{aligned}f(x)&=f(x_{1})+f'(x_{1})(x-x_{1})+{\frac {1}{2}}f''(x_{1})(x-x_{1})^{2}+{\frac {f(x_{2})-y(x_{1})-f'(x_{1})(x_{2}-x_{1})-{\frac {1}{2}}f''(x_{1})(x_{2}-x_{1})^{2}}{(x_{2}-x_{1})^{3}}}(x-x_{1})^{3}\\&+{\frac {3f(x_{1})-3f(x_{2})+2\left(f'(x_{1})+{\frac {1}{2}}f'(x_{2})\right)(x_{2}-x_{1})+{\frac {1}{2}}f''(x_{1})(x_{2}-x_{1})^{2}}{(x_{2}-x_{1})^{4}}}(x-x_{1})^{3}(x-x_{2})\\&+{\frac {6f(x_{2})-6f(x_{1})-3\left(f'(x_{1})+f'(x_{2})\right)(x_{2}-x_{1})+{\frac {1}{2}}\left(f''(x_{2})-f''(x_{1})\right)(x_{2}-x_{1})^{2}}{(x_{2}-x_{1})^{5}}}(x-x_{1})^{3}(x-x_{2})^{2}.\end{aligned}}$$

Basis polynomials

See also: § Basis polynomials in closed form, and § Basis polynomials via Chinese remainder theorem

For a given set of m distinct nodes ⁠$x_{1},\dots ,x_{m}$⁠ with associated maximum derivative orders ⁠$r_{1},\dots ,r_{m}$⁠, let ⁠$n_{i}=r_{i}+1$⁠ be the multiplicity of node ⁠$x_{i}$⁠, and let ⁠$\textstyle n=\sum _{i}n_{i}$⁠ be the total number of conditions. The Hermite interpolating polynomial ⁠$h(x)$⁠ can be expressed as a linear combination of basis polynomials ⁠$h_{i,j}(x)$⁠: $$h(x)=\sum _{i=1}^{m}\sum _{j=0}^{r_{i}}y_{i,j}h_{i,j}(x),$$ where each basis polynomial ⁠$h_{i,j}(x)$⁠ is uniquely defined as the interpolating polynomial of degree at most ⁠$n-1$⁠ for the basis conditions given by the Kronecker delta: $$\left.{\frac {d^{l}}{dx^{l}}}h_{i,j}(x)\right|_{x=x_{k}}=\delta _{i,k}\delta _{j,l}={\begin{cases}1&{\text{if }}k=i{\text{ and }}l=j,\\0&{\text{otherwise.}}\end{cases}}$$

To construct a specific basis polynomial ⁠$h_{i,j}(x)$⁠ using general methods for solving Hermite interpolation (such as those described in this article), one specifies a system of target conditions ⁠$\textstyle y_{k,l}^{(i,j)}$⁠ across all nodes ⁠$x_{k}$⁠ and derivative orders ⁠$0\leq l\leq r_{k}$⁠: $$y_{k,l}^{(i,j)}=\delta _{i,k}\delta _{j,l}={\begin{cases}1&{\text{if }}k=i{\text{ and }}l=j,\\0&{\text{otherwise.}}\end{cases}}$$

By passing these target conditions ⁠$\textstyle y_{k,l}^{(i,j)}$⁠ as input values into any general method for solving Hermite interpolation, the resulting unique interpolating polynomial generated by the method is precisely the basis polynomial ⁠$h_{i,j}(x)$⁠. Repeating this process for each pair of indices ⁠$(i,j)$⁠ constructs the full set of n basis polynomials.

Alternative basis construction via degree-0 basis polynomials

Let ⁠$\textstyle q(x)=\prod _{i}(x-x_{i})^{n_{i}}$⁠ be the polynomial of degree n, which vanishes at each node ⁠$x_{i}$⁠ with multiplicity ⁠$n_{i}$⁠. Once the degree-0 basis polynomials ⁠$h_{i,0}(x)$⁠ (corresponding to the function value condition ⁠$y_{i,0}=1$⁠ at node ⁠$x_{i}$⁠) are obtained, the higher-order derivative basis polynomials ⁠$h_{i,j}(x)$⁠ for ⁠$1\leq j\leq r_{i}$⁠ can be constructed directly via polynomial reduction: $$h_{i,j}(x)=\left({\frac {(x-x_{i})^{j}}{j!}}h_{i,0}(x)\right){\bmod {q}}(x).$$

Taking the remainder modulo ⁠$q(x)$⁠ guarantees that ⁠$\deg(h_{i,j})<\deg(q)=n$⁠, projecting the polynomial into the subspace of degree at most ⁠$n-1$⁠. Because ⁠$q(x)$⁠ vanishes at every node ⁠$x_{i}$⁠ to order ⁠$n_{i}$⁠, all target derivative evaluation conditions at the nodes are preserved under this modular reduction. This bounds the degree and ensures that each constructed ⁠$h_{i,j}(x)$⁠ is the unique polynomial of degree at most ⁠$n-1$⁠ satisfying its corresponding interpolation conditions.

This construction is particularly advantageous when direct evaluation of higher-order basis polynomials ⁠$h_{i,j}(x)$⁠ is computationally expensive or numerically unstable. In some cases, the degree-0 basis polynomials ⁠$h_{i,0}(x)$⁠ are simply cleaner and easier to construct and compute, allowing all remaining derivative basis polynomials for that node to be generated efficiently through straightforward polynomial multiplication and modular reduction.

Error

Let m be the number of distinct nodes ⁠$x_{1},\ldots ,x_{m}$⁠, let ⁠$n_{i}$⁠ be the multiplicity of node ⁠$x_{i}$⁠ representing the number of specified values including derivatives, and let ⁠$\textstyle n=\sum _{i=1}^{m}n_{i}$⁠ be the total number of interpolation conditions. Denote the constructed Hermite interpolating polynomial as ⁠$h(x)$⁠ and the original function as ⁠$f(x)$⁠. Consider first the real-valued case. Let ⁠$\textstyle a=\min _{i}x_{i}$⁠ and ⁠$\textstyle b=\max _{i}x_{i}$⁠. Evaluating at a point ⁠$x\in [a,b]$⁠, the error function is $$f(x)-h(x)={\frac {f^{(n)}(c)}{n!}}\prod _{i=1}^{m}(x-x_{i})^{n_{i}},$$ where c is an unknown value within the range ⁠$[a,b]$⁠, derived using a generalization of the mean value theorem for divided differences to confluent nodes. The degree of the product polynomial on the right is n, which is one higher than the degree bound ⁠$n-1$⁠ for ⁠$h(x)$⁠. Furthermore, the error and all its derivatives up to order ⁠$n_{i}-1$⁠ are zero at each node ⁠$x_{i}$⁠, as expected.

In the complex case, as described for example on p. 360 in , $$f(z)-h(z)={\frac {q(z)}{2\pi i}}\oint _{C}{\frac {f(\zeta )}{q(\zeta )(\zeta -z)}}d\zeta$$ where the contour C encloses z and all the nodes ⁠$x_{i}$⁠, and ⁠$\textstyle q(z)=\prod _{i=1}^{m}(z-x_{i})^{n_{i}}$⁠ is the product polynomial of all node factors.

See also

• Newton series, also known as finite differences
• Neville's schema
• Bernstein polynomials
• Generalized Frobenius covariant
• Generalized Sylvester's formula

References

• de Boor, Carl (2005). "Divided Differences". Surveys in Approximation Theory. 1: 46–69. arXiv:math/0502036. ISSN 1555-578X.
• Higham, Nicholas J. (2008). Functions of Matrices: Theory and Computation. Philadelphia, PA, USA: Society for Industrial and Applied Mathematics. doi:10.1137/1.9780898717778. ISBN 978-0-89871-646-7.
• Dubeau, François (2020). "On Hermite interpolation and divided differences" (PDF). Surveys in Mathematics and its Applications. 15: 257–279. ISSN 1842-6298.

External links

• Hermites Interpolating Polynomial at Mathworld