Hermite interpolation

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$ . The resulting polynomial is commonly known as the Hermite interpolating polynomial or the Hermite interpolant.

There are various methods for computing and formulating the Hermite interpolating polynomial. 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. Alternatively, 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. Other approaches include explicit closed-form formulas and the barycentric form, the latter of which enables efficient updating of the interpolating polynomial when adding new distinct nodes or higher-order derivative conditions. Further advanced or specialized techniques include the Chinese remainder theorem approach, the orthogonal projection method, and Hermite’s original approach using contour integration, which also yields an error function formulation.

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$ of degree less than $n$ with indeterminate coefficients; that is, the coefficients of $h$ 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$ ⁠ 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/q_{i}\right](x)}q_{i}(x),$ where ⁠ $\textstyle q_{i}(x)=\prod _{j\neq i}(x-x_{j})^{n_{j}}$ ⁠ and ⁠ $\textstyle {\mathcal {T}}_{a}^{k}[f]$ ⁠ denotes the $k$ th degree Taylor operator of a function ⁠ $f$ ⁠ centered at ⁠ $a$ ⁠.

For each node ⁠ $x_{i}$ ⁠, let ⁠ $g_{i}$ ⁠ be the local ⁠ $r_{i}$ ⁠th degree Taylor polynomial of both ⁠ $f$ ⁠ and ⁠ $h$ ⁠ at ⁠ $x_{i}$ ⁠, constructed from the interpolation conditions ⁠ $y_{i,j}$ ⁠: $g_{i}(x)=\sum _{j=0}^{r_{i}}{\frac {y_{i,j}}{j!}}(x-x_{i})^{j}={\mathcal {T}}_{x_{i}}^{r_{i}}[f](x)={\mathcal {T}}_{x_{i}}^{r_{i}}[h](x).$

Consider the case when only the prescribed interpolation conditions ⁠ $y_{i,j}$ ⁠ are given. Using the Taylor composition theorem, ⁠ $f$ ⁠ can be directly substituted with the prescribed local ⁠ $r_{i}$ ⁠th degree Taylor polynomial ⁠ $g_{i}$ ⁠ inside any occurrence of the ⁠ $r_{i}$ ⁠th degree Taylor operator at ⁠ $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}}{q_{i}}}\right]\!(x)\,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$ ⁠ can also be represented as a linear combination of Hermite basis polynomials ⁠ $h_{i,j}$ ⁠: $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}$ ⁠ 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}$ ⁠ into the explicit formula yields: $h_{i,j}(x)={\mathcal {T}}_{\xi =x_{i}}^{r_{i}}\!\left[{\frac {1}{j!}}{\frac {(\xi -x_{i})^{j}}{q_{i}(\xi )}}\right]\!(x)\,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}$ ⁠. For ⁠ $h_{i,0}$ ⁠, 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}}}\right]\!(x)\,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

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$ ⁠ such that $h(x_{i})=y_{i,0},\quad h'(x_{i})=y_{i,1},\quad \ldots ,\quad 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$ ⁠ 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$ ⁠ is congruent to this polynomial ⁠ $g_{i}$ ⁠ modulo ⁠ $p_{i}$ ⁠, with ⁠ $\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$ ⁠ must satisfy the system of modular equivalences: $h\equiv g_{i}{\pmod {p_{i}}}\quad {\text{for }}i=1,\dots ,m,$

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

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

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

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

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

This construction satisfies all conditions up to node $i$ :

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

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

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}$ ⁠, which can be obtained in a simplified manner.

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

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

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

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

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

Since the polynomial ⁠ $v_{i}$ ⁠ is not required for the construction, ⁠ $u_{i}$ ⁠ 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}$ ⁠ 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}q_{i}$ ⁠ already satisfies the degree constraint without requiring any further modular reduction, yielding directly the unique degree- $0$ basis polynomial ⁠ $h_{i,0}=u_{i}q_{i}$ ⁠.

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 ⁠ $\textstyle 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}$ ⁠ 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

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$ ⁠ 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$ ⁠ 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}]\,\underbrace {\textstyle \prod _{j=1}^{i-1}(x-z_{j})} _{h_{i}(x)},$ where ⁠ $\textstyle \{h_{i}\}_{i=1}^{n}$ ⁠ is commonly known as the Newton basis. 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.

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})$ ⁠, 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 $i$ th 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$ ⁠ 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 ⁠ $f[z_{a},\dots ,z_{b}]$ ⁠ when constructed from the prescribed derivative conditions ⁠ $y_{i,j}$ ⁠, the constructed divided difference table is given by: ${\begin{array}{llclrr}z_{1}=-1&f[z_{1}]=y_{1,0}=1&&&&\\&&{\frac {y_{1,1}}{1!}}=2&&&\\z_{2}=-1&f[z_{2}]=y_{1,0}=1&&f[z_{3},z_{2},z_{1}]=-2&&\\&&f[z_{3},z_{2}]={\phantom {-}}0&&{\frac {1}{2}}&\\z_{3}=0&f[z_{3}]=y_{2,0}=1&&f[z_{4},z_{3},z_{2}]=-1&&0\\&&f[z_{4},z_{3}]=-2&&{\frac {1}{2}}&\\z_{4}=1&f[z_{4}]=y_{3,0}=-1&&f[z_{5},z_{4},z_{3}]=0&&\\&&{\frac {y_{3,1}}{1!}}=-2&&&\\z_{5}=1&f[z_{5}]=y_{3,0}=-1&&&&\\\end{array}}$ and the generated polynomial ⁠ $h$ ⁠ 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}}$

Confluent divided differences via contour integration

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

When ⁠ $f$ ⁠ is analytic on a simply connected domain ⁠ $\Omega \subset \mathbb {C}$ ⁠ containing the dataset nodes ⁠ $z_{1},\dots ,z_{n}$ ⁠, substituting Cauchy's integral formula into the Newton form yields: $h(x)={\frac {1}{2\pi i}}\sum _{i=1}^{n}\oint _{\Gamma _{i}}{\frac {f(z)}{\prod _{j=1}^{i}(z-z_{j})}}\,dz\prod _{j=1}^{i-1}(x-z_{j}),$ where ⁠ $\Gamma _{i}$ ⁠ is a positively oriented, simple closed contour in ⁠ $\Omega$ ⁠ enclosing the first ⁠ $i$ ⁠ dataset nodes ⁠ $z_{1},\dots ,z_{i}$ ⁠ while excluding the remaining nodes.

Using 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}$ ⁠ correspond to the basis polynomials ⁠ $h_{k}$ ⁠ 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}=h_{k}$ ⁠ 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$ ⁠ 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$ ⁠ 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$ ⁠ 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$ ⁠ has the following poles inside ⁠ $\Gamma$ ⁠:

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

By the residue theorem, integrating ⁠ $F$ ⁠ 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 ⁠ $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$ ⁠: $h(x)=-q(x)\sum _{i=1}^{m}\operatorname {Res} _{\zeta =x_{i}}F(\zeta ).$ Substituting this definition into the expression for ⁠ $f(x)$ ⁠ above 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-h$ ⁠ (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})$

Using barycentric form

Barycentric Hermite interpolation is a formulation allowing for the updating of the resulting Hermite interpolating polynomial, such as adding new interpolation nodes or increasing the derivative degree at existing nodes, without having to reconstruct the entire polynomial from scratch.

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}$ ⁠. Let ⁠ $h$ ⁠ be the resulting unique interpolating polynomial of degree at most ⁠ $n-1$ ⁠, and let ⁠ $f$ ⁠ be the target function, with ⁠ $\textstyle y_{i,j}=f^{(j)}(x_{i})$ ⁠ denoting the prescribed derivatives.

The first barycentric form is given by: $h(x)=\sum _{i=1}^{m}\sum _{j=0}^{r_{i}}w_{i,j}\,q(x)\sum _{k=0}^{r_{i}-j}{\frac {y_{i,k}}{k!}}(x-x_{i})^{(j+k)-n_{i}}$ where ⁠ $\textstyle q(x)=\prod _{i=1}^{m}(x-x_{i})^{n_{i}}$ ⁠ and ⁠ $w_{i,j}$ ⁠ for ⁠ $i=1,\dots ,m$ ⁠ and ⁠ $j=0,\dots ,r_{i}$ ⁠ are the barycentric weights to be determined. The barycentric form can be viewed as a linear combination whose coefficients are given by the barycentric weights. Solving the first barycentric form for the derivative conditions ⁠ $y_{i,j}=\delta _{j,0}$ ⁠, where ⁠ $\delta _{a,b}$ ⁠ is the Kronecker delta, yields the unique polynomial ⁠ $g(x)=1$ ⁠; dividing the first barycentric form by ⁠ $g$ ⁠ yields the second barycentric form: $h(x)=\sum _{i=1}^{m}\sum _{j=0}^{r_{i}}w_{i,j}{\frac {\textstyle \sum _{k=0}^{r_{i}-j}{\frac {y_{i,k}}{k!}}(x-x_{i})^{(j+k)-n_{i}}}{\textstyle \sum _{k=1}^{m}\sum _{l=0}^{r_{k}}w_{k,l}(x-x_{k})^{l-n_{k}}}}.$ Note that the factor ⁠ $q(x)$ ⁠ is cancelled out in the second barycentric form, which improves numerical stability compared to the first barycentric form. For evaluation points ⁠ $x$ ⁠ which are close to one of the nodes ⁠ $x_{i}$ ⁠, catastrophic cancellation would ordinarily be a problem for the value ⁠ $(x-x_{i})$ ⁠; however, this quantity appears in both numerator and denominator and they cancel out, leaving good relative accuracy in the final result. Using this formula to evaluate ⁠ $h(x)$ ⁠ at one of the nodes ⁠ $x_{i}$ ⁠ will result in the indeterminate ⁠ $\textstyle \sum _{j=0}^{r_{i}}\infty y_{i,j}/\infty$ ⁠; computer implementations must replace such results by the prescribed local Taylor polynomial ⁠ $\textstyle h(x_{i})=\sum _{j=0}^{r_{i}}y_{i,j}(x-x_{i})^{j}/j!$ ⁠.

There are several ways to determine the barycentric weights ⁠ $w_{i,j}$ ⁠. One approach is to solve: ${\mathcal {T}}_{x_{i}}^{r_{i}}\!\left[{\frac {1}{q_{i}}}\right]\!(x)=\sum _{j=0}^{r_{i}}w_{i,j}(x-x_{i})^{j}\implies w_{i,j}={\frac {1}{j!}}\left.{\frac {d^{j}}{dx^{j}}}{\frac {1}{q_{i}(x)}}\right|_{x=x_{i}}$ where ⁠ $\textstyle q_{i}(x)=\prod _{j\neq i}(x-x_{j})^{n_{j}}$ ⁠ and ⁠ $\textstyle {\mathcal {T}}_{a}^{k}[f]$ ⁠ denotes the $k$ th degree Taylor operator of a function ⁠ $f$ ⁠ centered at ⁠ $a$ ⁠. Another approach is to use the partial fraction decomposition: ${\frac {1}{q(x)}}=\sum _{i=1}^{m}\sum _{j=0}^{r_{i}}{\frac {w_{i,j}}{(x-x_{i})^{n_{i}-j}}}.$

Incremental updates for new nodes and higher derivative orders

There exists an algorithm with ⁠ $O(n)$ ⁠ space complexity and ⁠ $O(n)$ ⁠ time complexity for adding new nodes or adding higher derivative orders to existing nodes. For initial distinct nodes and derivative conditions, with ⁠ ${\mathcal {I}}_{i,0}=1$ ⁠, the initial intermediate data and barycentric weights are constructed as follows: ${\begin{array}{lll}C_{i}&=\textstyle \prod _{j\neq i}(x_{i}-x_{j})^{-n_{j}}&{\text{for }}i=1,\dots ,m\\{\mathcal {P}}_{i,j}&=\textstyle \sum _{k\neq i}n_{k}(x_{k}-x_{i})^{-j}&{\text{for }}i=1,\dots ,m{\text{ and }}j=1,\dots ,r_{i}\\{\mathcal {I}}_{i,j}&=\textstyle {\frac {1}{j}}\sum _{k=1}^{j}{\mathcal {P}}_{i,k}{\mathcal {I}}_{i,j-k}&{\text{for }}i=1,\dots ,m{\text{ and }}j=1,\dots ,r_{i}\\w_{i,j}&=C_{i}{\mathcal {I}}_{i,j}&{\text{for }}i=1,\dots ,m{\text{ and }}j=0,\dots ,r_{i}\end{array}}$

When either a new distinct node is added at ⁠ $\xi$ ⁠ or a higher derivative order is added at ⁠ $\xi =x_{k}$ ⁠ for some ⁠ $1\leq k\leq m$ ⁠, with ⁠ $\textstyle {\mathcal {I}}'_{i,0}={\mathcal {I}}_{i,0}$ ⁠ and taking ⁠ $w'_{i,-1}=0$ ⁠, the intermediate data and barycentric weights at existing distinct nodes ⁠ $x_{i}\neq \xi$ ⁠, omitting ⁠ $x_{k}$ ⁠ when ⁠ $\xi =x_{k}$ ⁠, are updated according to: ${\begin{array}{lll}C'_{i}&=C_{i}/(x_{i}-\xi )\\{\mathcal {P}}'_{i,j}&={\mathcal {P}}_{i,j}+(\xi -x_{i})^{-j}&{\text{for }}j=1,\dots ,r_{i}\\{\mathcal {I}}'_{i,j}&={\mathcal {I}}_{i,j}+{\mathcal {I}}'_{i,j-1}/(\xi -x_{i})&{\text{for }}j=1,\dots ,r_{i}\\w'_{i,j}&=(w_{i,j}-w'_{i,j-1})/(x_{i}-\xi )&{\text{for }}j=0,\dots ,r_{i}\end{array}}$

In the case of adding a new distinct node at ⁠ $x_{m+1}=\xi$ ⁠, noting that there is no ⁠ ${\mathcal {P}}_{m+1,0}$ ⁠ and that ⁠ ${\mathcal {I}}_{i,0}=1$ ⁠ for ⁠ $i=m+1$ ⁠, the remaining new intermediate data and barycentric weight follow the initial construction rules, which are given by: ${\begin{array}{ll}C_{m+1}&=\textstyle \prod _{i=1}^{m}(\xi -x_{i})^{-n_{i}}\\w_{m+1,0}&=C_{m+1}{\mathcal {I}}_{m+1,0}=C_{m+1}\end{array}}$

In the case of adding a new derivative order at an existing distinct node ⁠ $x_{i}=\xi$ ⁠, existing intermediate data and barycentric weights at node ⁠ $x_{i}$ ⁠ do not require updating, and the new intermediate data and barycentric weight corresponding to the new derivative order ⁠ $r_{i}+1$ ⁠ follow the initial construction rules, which are given by: ${\begin{array}{ll}{\mathcal {P}}_{i,r_{i}+1}&=\textstyle \sum _{k\neq i}n_{k}(x_{k}-x_{i})^{-(r_{i}+1)}\\{\mathcal {I}}_{i,r_{i}+1}&=\textstyle {\frac {1}{r_{i}+1}}\sum _{k=1}^{r_{i}+1}{\mathcal {P}}_{i,k}{\mathcal {I}}_{i,(r_{i}+1)-k}\\w_{i,r_{i}+1}&=C_{i}{\mathcal {I}}_{i,r_{i}+1}\end{array}}$

The updated Hermite interpolating polynomial ⁠ $h$ ⁠ can be constructed using the updated variables by applying either the first barycentric form with an updated ⁠ $q'(x)=q(x)(x-\xi )$ ⁠ or the second barycentric form directly.

Basis polynomials in barycentric form

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

In the first barycentric form: $h_{i,j}(x)={\frac {q(x)}{j!}}\sum _{k=0}^{r_{i}-j}w_{i,k}(x-x_{i})^{(j+k)-n_{i}}$ In the second barycentric form: $h_{i,j}(x)={\frac {1}{j!}}{\frac {\textstyle \sum _{k=0}^{r_{i}-j}w_{i,k}(x-x_{i})^{(j+k)-n_{i}}}{\textstyle \sum _{k=1}^{m}\sum _{l=0}^{r_{k}}w_{k,l}(x-x_{k})^{l-n_{k}}}}.$

Common low-degree cases

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$ ⁠ 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$ ⁠, its first ⁠ $f'$ ⁠ and second derivatives ⁠ $f''$ ⁠ 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 of the resulting quintic polynomial ⁠ $h$ ⁠ is given by ${\begin{aligned}h(x)&=f(x_{1})+f'(x_{1})(x-x_{1})+{\tfrac {1}{2}}f''(x_{1})(x-x_{1})^{2}\\[4pt]&\qquad +{\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}\\[2pt]&\qquad +{\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})\\[2pt]&\qquad +{\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

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$ ⁠ can be expressed as a linear combination of basis polynomials ⁠ $h_{i,j}$ ⁠: $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}$ ⁠ is uniquely defined as the interpolating polynomial of degree at most ⁠ $n-1$ ⁠ satisfying the derivative conditions given by the Kronecker delta: $h_{i,j}^{(l)}(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}$ ⁠, one can apply any general method for solving Hermite interpolation using the corresponding derivative conditions prescribed by the Kronecker delta above. By passing these required derivative values as input into the selected interpolation method, the resulting unique polynomial is precisely ⁠ $h_{i,j}$ ⁠. Repeating this process for each pair of indices ⁠ $(i,j)$ ⁠ constructs the full set of $n$ basis polynomials.

The above basis is commonly known as the generalized Lagrange basis. Other notable common bases are the Newton basis of the confluent divided differences method, the monomial basis of the confluent Vandermonde matrix, and the orthogonal basis of the orthogonal projection method. There is also the barycentric basis that expresses the Hermite interpolating polynomial in terms of barycentric weights.

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}$ ⁠ are obtained, the Hermite interpolation polynomial ⁠ $h$ ⁠ can be constructed directly via polynomial reduction: $h(x)=\sum _{i=1}^{m}\underbrace {\left[\sum _{j=0}^{r_{i}}{\frac {y_{i,j}}{j!}}(x-x_{i})^{j}\right]} _{g_{i}(x)}h_{i,0}(x){\bmod {q}}(x),$ where ⁠ $g_{i}$ ⁠ is the local ⁠ $r_{i}$ ⁠⁠th degree Taylor polynomial at ⁠ $x_{i}$ ⁠⁠, constructed from the interpolation conditions ⁠ $y_{i,j}$ ⁠. Taking the remainder modulo ⁠ $q$ ⁠ guarantees that ⁠ $\deg(h)<\deg(q)=n$ ⁠, projecting the polynomial into the subspace of degree at most ⁠ $n-1$ ⁠. Because ⁠ $q$ ⁠ 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 the constructed ⁠ $h$ ⁠ is the unique polynomial of degree at most ⁠ $n-1$ ⁠ satisfying its corresponding interpolation conditions.

Alternatively, the higher-order derivative basis polynomials ⁠ $h_{i,j}$ ⁠ for ⁠ $1\leq j\leq r_{i}$ ⁠ can be expressed directly via polynomial reduction: $h_{i,j}(x)={\frac {1}{j!}}(x-x_{i})^{j}h_{i,0}(x){\bmod {q}}(x).$

This construction is particularly advantageous when direct evaluation of higher-order basis polynomials ⁠ $h_{i,j}$ ⁠ is computationally expensive or numerically unstable. In some cases, the degree- $0$ basis polynomials ⁠ $h_{i,0}$ ⁠ 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$ ⁠ and the original function as ⁠ $f$ ⁠. 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$ ⁠. 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 _{\Gamma }{\frac {f(\zeta )}{q(\zeta )(\zeta -z)}}d\zeta$ where the contour ⁠ $\Gamma$ ⁠ 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.

References

↑ Traub, J. F. (December 1964). "On Lagrange—Hermite interpolation". Journal of the Society for Industrial and Applied Mathematics. 12 (4): 886–891. doi:10.1137/0112076. JSTOR 2946361.
↑ Dubeau 2020, p. 258.
1 2 Hermite, Charles (1878). "Sur la formule d'interpolation de Lagrange". Journal für die reine und angewandte Mathematik. 84: 70–79.
↑ Spitzbart, A (January 1960). "A Generalization of Hermite Interpolation". American Mathematical Monthly. 67 (1): 42–46. doi:10.1080/00029890.1960.11989446. JSTOR 2308924. Retrieved 2 June 2024.
↑ Higham 2008, Eq. (1.8).
↑ Buchheim, Arthur (August 1886). "An Extension of a Theorem of Professor Sylvester's relating to Matrices". The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science. Fifth Series. 22 (135): 173–174. doi:10.1080/14786448608627914. ISSN 1941-5982.
↑ Malgrange, Bernard (1966). Ideals of Differentiable Functions. Tata Institute of Fundamental Research Studies in Mathematics. Oxford: Oxford University Press. p. 16. ISBN 978-0-19-635241-1.
1 2 3 Dubeau 2020, pp. 273–274.
↑ Rosen, Kenneth H. (1993). Elementary Number Theory and its Applications (3rd ed.). Addison-Wesley. p. 136. ISBN 978-0201-57889-8.
↑ de Boor 2005.
↑ Higham 2008, Appendix B.16.
↑ Dubeau 2020, pp. 274–276.
↑ de Boor 2005, Proposition 7.
↑ Higham 2008, Eq. (1.10).
↑ de Boor 2005, p. 56.
↑ Schneider, C; Werner, W (1991). "Hermite Interpolation: The Barycentric Approach". Computing. 46: 35–51. doi:10.1007/BF02239010.
1 2 Corless, Robert M; Fillion, Nicolas (2013). A Graduate Introduction to Numerical Methods. New York: Springer. ISBN 978-1-4614-8452-3.
↑ de Boor 2005, Eq. (51).
↑ Higham 2008, Eq. (B.29).
↑ Kalman, D. (1984). "The Generalized Vandermonde Matrix". Mathematics Magazine. 57 (1): 15–21. doi:10.1080/0025570X.1984.11977069.
↑ Dubeau 2020, p. 272.
↑ Dubeau 2020, pp. 276–277.
1 2 3 Sadiq, Burhan; Viswanath, Divakar (2013). "Barycentric Hermite Interpolation". SIAM Journal on Scientific Computing. 35: A1254–A1270. arXiv:1105.3466. doi:10.1137/110833221. ISSN 1064-8275.
1 2 Corless, Robert M. (2024). "Barycentric Hermite Interpolation". Maple Transactions. 4 (2): 19001:1–19001:32. doi:10.5206/mt.v4i2.19001. ISSN 2564-3029.
↑ Dubeau 2020, sec. 5.
↑ Bossoney, Simon; Troyanov, Marc (2026). "A Direct Approach to Hermite Interpolation". arXiv:2607.17998 [math.NA].
↑ Dubeau 2020, p. 261.

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

[1] Traub, J. F. (December 1964). "On Lagrange—Hermite interpolation". Journal of the Society for Industrial and Applied Mathematics. 12 (4): 886–891. doi:10.1137/0112076. JSTOR 2946361.

[FOOTNOTEDubeau2020258-2] Dubeau 2020, p. 258.

[Hermite1878-3] 1 2 Hermite, Charles (1878). "Sur la formule d'interpolation de Lagrange". Journal für die reine und angewandte Mathematik. 84: 70–79.

[4] Spitzbart, A (January 1960). "A Generalization of Hermite Interpolation". American Mathematical Monthly. 67 (1): 42–46. doi:10.1080/00029890.1960.11989446. JSTOR 2308924. Retrieved 2 June 2024.

[FOOTNOTEHigham2008Eq._(1.8)-5] Higham 2008, Eq. (1.8).

[6] Buchheim, Arthur (August 1886). "An Extension of a Theorem of Professor Sylvester's relating to Matrices". The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science. Fifth Series. 22 (135): 173–174. doi:10.1080/14786448608627914. ISSN 1941-5982.

[7] Malgrange, Bernard (1966). Ideals of Differentiable Functions. Tata Institute of Fundamental Research Studies in Mathematics. Oxford: Oxford University Press. p. 16. ISBN 978-0-19-635241-1.

[FOOTNOTEDubeau2020273–274-8] 1 2 3 Dubeau 2020, pp. 273–274.

[9] Rosen, Kenneth H. (1993). Elementary Number Theory and its Applications (3rd ed.). Addison-Wesley. p. 136. ISBN 978-0201-57889-8.

[FOOTNOTEde_Boor2005-10] Boor 2005.

[FOOTNOTEHigham2008Appendix_B.16-11] Higham 2008, Appendix B.16.

[FOOTNOTEDubeau2020274–276-12] Dubeau 2020, pp. 274–276.

[FOOTNOTEde_Boor2005Proposition_7-13] Boor 2005, Proposition 7.

[FOOTNOTEHigham2008Eq._(1.10)-14] Higham 2008, Eq. (1.10).

[FOOTNOTEde_Boor200556-15] Boor 2005, p. 56.

[16] Schneider, C; Werner, W (1991). "Hermite Interpolation: The Barycentric Approach". Computing. 46: 35–51. doi:10.1007/BF02239010.

[CF2013-17] 1 2 Corless, Robert M; Fillion, Nicolas (2013). A Graduate Introduction to Numerical Methods. New York: Springer. ISBN 978-1-4614-8452-3.

[FOOTNOTEde_Boor2005Eq._(51)-18] Boor 2005, Eq. (51).

[FOOTNOTEHigham2008Eq._(B.29)-19] Higham 2008, Eq. (B.29).

[20] Kalman, D. (1984). "The Generalized Vandermonde Matrix". Mathematics Magazine. 57 (1): 15–21. doi:10.1080/0025570X.1984.11977069.

[FOOTNOTEDubeau2020272-21] Dubeau 2020, p. 272.

[FOOTNOTEDubeau2020276–277-22] Dubeau 2020, pp. 276–277.

[Sadiq2013-23] 1 2 3 Sadiq, Burhan; Viswanath, Divakar (2013). "Barycentric Hermite Interpolation". SIAM Journal on Scientific Computing. 35: A1254–A1270. arXiv:1105.3466. doi:10.1137/110833221. ISSN 1064-8275.

[Corless2024-24] 1 2 Corless, Robert M. (2024). "Barycentric Hermite Interpolation". Maple Transactions. 4 (2): 19001:1–19001:32. doi:10.5206/mt.v4i2.19001. ISSN 2564-3029.

[FOOTNOTEDubeau2020sec._5-25] Dubeau 2020, sec. 5.

[26] Bossoney, Simon; Troyanov, Marc (2026). "A Direct Approach to Hermite Interpolation". arXiv:2607.17998 [math.NA].

[FOOTNOTEDubeau2020261-27] Dubeau 2020, p. 261.

Authority control databases
International	GND
Other	Yale LUX