Hermite interpolation

Short description: Polynomial interpolation using derivative values

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, which uses contour integration.

Hermite interpolating polynomials are also deeply related to matrix diagonalization, Jordan normal form and matrix functions via generalized Frobenius covariants and generalized Sylvester's formula.

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{array}{cccc}(x_1,y_1),& (x_2,y_2),& \dots ,& (x_m,y_m),\\ (x_1,{y_1}^{\prime }),& (x_2,{y_2}^{\prime }),& \dots ,& (x_m,{y_m}^{\prime }),\\ ⋮& ⋮& & ⋮\\ (x_1,y_1^{(r)}),& (x_2,y_2^{(r)}),& \dots ,& (x_m,y_m^{(r)})\end{array}$$ 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.

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 § A perspective from 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 $n={\sum }_{i=1}^m{n}_i$ be the total number of conditions. The Hermite interpolating polynomial $h(x)$ is the unique polynomial of degree at most $n-1$ that satisfies the interpolation conditions: $$h^{(j)}(x_i)=y_{i,j}\text{for }i=1,\dots ,m\text{ and }j=0,\dots ,r_i.$$

Define the complementary factor polynomial $q_i(x)$ at node $x_i$ as: $$q_i(x)=\prod _{j\ne i}(x-x_j{)}^{n_j}.$$

For each node $x_i$, let $g_i(x)$ be the local $r_i$th degree Taylor polynomial matching the target derivative conditions: $$g_i(x)=\sum _{j=0}^{r_i}\frac{y_{i,j}}{j!}(x-x_i{)}^j.$$

The explicit closed-form formula for the Hermite interpolating polynomial $h(x)$ is given by: $$h(x)=\sum _{i=1}^m\underset{{{𝒯}}_{x_i}^{m_i}\left[\frac{g_i(x)}{q_i(x)}\right]}{\underbrace{[\sum _{j=0}^{r_i}\frac{1}{j!}{{\left(\frac{g_i(x)}{q_i(x)}\right)}^{(j)}|}_{x=x_i}(x-x_i{)}^j]}}{q}_i(x),$$ where ${{𝒯}}_a^k[f(x)]$ denotes the kth degree Taylor operator of a function $f(x)$ centered at $x=a$.

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.

Interpolating a function in closed form

When interpolating a function $f(x)$ that is $r_i$ times differentiable at each node $x_i$, the interpolation conditions are determined by the function's derivatives (i.e., $y_{i,j}=f^{(j)}(x_i)$). In this case, the local Taylor polynomial $g_i(x)$ represents the $r_i$th degree Taylor expansion of $f(x)$ at $x_i$.

Since the operator ${{𝒯}}_{x_i}^{r_i}$ isolates terms only up to order $r_i$ at $x_i$, substituting $f(x)$ in place of its local Taylor expansion $g_i(x)$ yields the identical Taylor expansion up to order $r_i$: $${{𝒯}}_{x_i}^{r_i}\left[\frac{g_i(x)}{q_i(x)}\right]={{𝒯}}_{x_i}^{r_i}\left[\frac{f(x)}{q_i(x)}\right].$$

The explicit closed-form Hermite interpolating polynomial for $f(x)$ can thus be expressed concisely as: $$h(x)=\sum _{i=1}^m{{𝒯}}_{x_i}^{r_i}\left[\frac{f(x)}{q_i(x)}\right]q_i(x).$$

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 choosing the local Taylor polynomial $g_i(x)=(x-x_i{)}^j/j!$ and setting $g_k(x)=0$ for all $k\ne i$. Substituting $g_i(x)$ into the explicit formula yields: $$h_{i,j}(x)={{𝒯}}_{x_i}^{r_i}\left[\frac{(x-x_i{)}^j}{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 polynomial simplifies to $g_i(x)=1$, giving the explicit closed-form expression: $$h_{i,0}(x)={{𝒯}}_{x_i}^{r_i}\left[\frac{1}{q_i(x)}\right]q_i(x).$$

Proof of the closed form

The proof proceeds in three steps. First, it shows that evaluating derivatives at any node $x_i$ isolates only the ith term in the summation. Next, it verifies that this term matches the required derivative values via local Taylor expansion. Finally, it establishes that the resulting polynomial satisfies the degree bound and uniqueness.

Vanishing of cross-terms at distinct nodes

Fix an arbitrary target node $x_i$ and derivative order $j\in \{0,\dots ,r_i\}$. Evaluating the jth derivative of $h(x)$ at $x=x_i$ yields:

$$h^{(j)}(x_i)={\frac{d^j}{d{x}^j}(\sum _{k=1}^m{{𝒯}}_{x_k}^{r_k}\left[\frac{g_k(x)}{q_k(x)}\right]q_k(x))|}_{x=x_i}$$

For any node index $k\ne i$, the factor polynomial $q_k(x)$ contains $(x-x_i{)}^{n_i}$ as an explicit factor:

$$q_k(x)=(x-x_i{)}^{n_i}\prod _{l\ne i,k}(x-x_l{)}^{n_l}$$

By the general Leibniz rule, evaluating the jth derivative of ${{𝒯}}_{x_k}^{r_k}[g_k(x)/q_k(x)]q_k(x)$ at $x=x_i$ produces a linear combination of terms involving ${d^l/d{x}^l(x-x_i{)}^{n_i}|}_{x=x_i}$ for $0\le l\le j$. Because $l\le j\le r_i<n_i$, each derivative evaluates to 0 at $x=x_i$:

$${\frac{d^l}{d{x}^l}(x-x_i{)}^{n_i}|}_{x=x_i}={\frac{n_i!}{(n_i-l)!}(x-x_i{)}^{n_i-l}|}_{x=x_i}=0\text{for all }0\le l\le j$$

Consequently, all terms $k\ne i$ in the summation vanish at $x=x_i$.

Verification of target interpolation conditions

Evaluating the remaining term $k=i$ at $x=x_i$:

$$h^{(j)}(x_i)={\frac{d^j}{d{x}^j}({{𝒯}}_{x_i}^{r_i}\left[\frac{g_i(x)}{q_i(x)}\right]q_i(x))|}_{x=x_i}$$

Since the nodes are distinct, $q_i(x_i)=\prod _{k\ne i}(x_i-x_k{)}^{n_k}\ne 0$, so the rational function $g_i(x)/q_i(x)$ is analytic at $x=x_i$. By definition of the Taylor operator ${{𝒯}}_{x_i}^{r_i}$, it expands locally around $x_i$ as:

$${{𝒯}}_{x_i}^{r_i}\left[\frac{g_i(x)}{q_i(x)}\right]=\frac{g_i(x)}{q_i(x)}+O((x-x_i{)}^{n_i})$$

Multiplying by $q_i(x)$ yields:

$${{𝒯}}_{x_i}^{r_i}\left[\frac{g_i(x)}{q_i(x)}\right]q_i(x)=(\frac{g_i(x)}{q_i(x)}+O((x-x_i{)}^{n_i}))q_i(x)=g_i(x)+O((x-x_i{)}^{n_i})$$

Because the big O term represents an analytic function near $x_i$ with a zero of order at least $n_i$, all of its derivatives of orders up to $r_i$ vanish at $x=x_i$. Taking the jth derivative ($0\le j\le r_i$) and evaluating at $x=x_i$ yields:

$$h^{(j)}(x_i)={\frac{d^j}{d{x}^j}(g_i(x)+O((x-x_i{)}^{n_i}))|}_{x=x_i}=g_i^{(j)}(x_i)+0=y_{i,j}$$

Thus, $h(x)$ satisfies all n interpolation conditions $h^{(j)}(x_i)=y_{i,j}$.

Degree bound and uniqueness

For each distinct node $x_i$, the degree of ${{𝒯}}_{x_i}^{r_i}[g_i(x)/q_i(x)]$ is at most $r_i$, and $\mathrm{deg}(q_i)=\sum _{j\ne i}{n}_j=n-n_i$. Therefore, each term in $h(x)$ has degree at most $r_i+n-n_i=n-1$, so $\mathrm{deg}(h)\le n-1$.

To establish uniqueness, assume another polynomial $p(x)$ of degree $\le n-1$ satisfies the same n conditions. Because $h^{(j)}(x_i)=p^{(j)}(x_i)=y_{i,j}$ for all $1\le i\le m$ and $0\le j\le r_i$, the difference $d(x)=f(x)-p(x)$ satisfies $d^{(j)}(x_i)=0$ for all $0\le j\le r_i$, meaning each node $x_i$ is a root of multiplicity at least $n_i$. This gives n roots counted with multiplicity across the m nodes. A nonzero polynomial of degree at most $n-1$ cannot have n roots counted with multiplicity. Therefore, $f(x)-p(x)$ must be the zero polynomial, so $f(x)=p(x)$.

Using Chinese remainder theorem

Let m be a positive integer representing the number of distinct nodes, $r_1,\dots ,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,\dots ,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^{\prime }(x_i)=y_{i,1},\dots ,h^{(r_i)}(x_i)=y_{i,r_i}$$ for $i=1,\dots ,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 $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 $n={\sum }_{i=1}^m{n}_i$.

Moreover, this solution can be computed with $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)(\mathrm{mod}{p}_i(x))\text{for }i=1,\dots ,m,$$ where $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 $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(\mathrm{mod}{p}_i(x))$, giving $a_i(x)\equiv g_i(x)(\mathrm{mod}{p}_i(x))$.
• Modulo $p_j(x)$ for any $j<i$, $b_{i-1}(x)\equiv 0(\mathrm{mod}{p}_j(x))$ and $v_i(x)p_i(x)\equiv 1(\mathrm{mod}{p}_j(x))$, giving $a_i(x)\equiv a_{i-1}(x)\equiv g_j(x)(\mathrm{mod}{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 $b_m(x)={\prod }_{i=1}^m{p}_i(x)$: $$h(x)=a_m(x){modb}_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(\mathrm{mod}{p}_j(x))\text{for }j\ne i$,
• $h_{i,0}(x)\equiv 1(\mathrm{mod}{p}_i(x))$.

To satisfy the zero-remainder conditions at all other nodes $j\ne i$, $h_{i,0}(x)$ must contain the product factor $q_i(x)=\prod _{j\ne 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(\mathrm{mod}{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 $\mathrm{deg}(p_i)=n_i$, and $\mathrm{deg}(q_i)=n-n_i$, the product degree is bounded by: $$\mathrm{deg}(u_i{q}_i)=\mathrm{deg}(u_i)+\mathrm{deg}(q_i)<n_i+(n-n_i)=n.$$ Thus, $\mathrm{deg}(u_i{q}_i)\le 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(\mathrm{mod}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\ne 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 $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. 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)=(\underset{n_1\text{ times}}{\underbrace{x_1,\dots ,x_1}},\underset{n_2\text{ times}}{\underbrace{x_2,\dots ,x_2}},\dots ,\underset{n_m\text{ times}}{\underbrace{x_m,\dots ,x_m}}).$$ 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 $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.

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,\dots ,x_m$, and values $f(x_1),\dots ,f(x_m)$ and $f^{\prime }(x_1),\dots ,f^{\prime }(x_m)$ for a function $f(x)$ that we want to interpolate, we create a new dataset $z_1,\dots ,z_{2m}$ such that $$z_{2i-1}=z_{2i}=x_i.$$

Now, we create a divided differences table for the nodes $z_1,\dots ,z_{2m}$. However, for some divided differences, $$z_i=z_{i+1}⟹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^{\prime }(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,\dots ,z_n$ contains $n_i$ adjacent identical copies of $x_i$. When creating the table, divided differences of $j+1$ (for $j=0,\dots ,r_i$) identical values will be calculated as $$f[\underset{j+1\text{times}}{\underbrace{x_i,\dots ,x_i}}]=\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 $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^{\prime }(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,\dots ,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^{\prime }(z_1)}{1}=-8& & & & & & & \\ z_2=-1& f[z_2]=2& & \frac{f^{″}(z_2)}{2}=28& & & & & & \\ & & \frac{f^{\prime }(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^{\prime }(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^{\prime }(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^{\prime }(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^{\prime }(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 ${\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{array}{rl}h(x)& =2-8(x+1)+28(x+1{)}^2-21(x+1{)}^3+15x(x+1{)}^3-10x^2(x+1{)}^3\\ & +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\\ & +45x^3-30x^3+4x^3+x^3+x^3+15x^4-30x^4+12x^4+2x^4+x^4\\ & -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{array}$$

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,\dots ,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{array}{rl}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{array}$$

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,\dots ,x_m$, let $n_i$ be the multiplicity of node $x_i$ representing the number of specified values including derivatives, and let $n={\sum }_{i=1}^m{n}_i$ be the total number of data conditions. Let V be the constructed confluent Vandermonde matrix, let $H=V^{-1}$ be its inverse, and let $d=n-1$ be the degree bound of the interpolating polynomial. The matrix H is given by: $$H=\left(\begin{array}{ccc}{H}_{00}& \cdots & H_{0d}\\ ⋮& \ddots & ⋮\\ H_{d0}& \cdots & H_{dd}\end{array}\right).$$

The sample vector $\mathbf{𝐲}$ 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{𝐲}$, the coefficient vector $\mathbf{𝐚}={\left(\begin{array}{cccc}{a}_0& a_1& \cdots & a_d\end{array}\right)}^T$ of the Hermite interpolating polynomial $h(x)={\sum }_{k=0}^d{a}_k{x}^k$ can be computed as: $$\mathbf{𝐚}=H\mathbf{𝐲}.$$

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 H: $$h_{i,j}(x)=h_k(x)=H_{0k}+H_{1k}x+\cdots +H_{dk}{x}^d.$$

In numerical practice, directly computing 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.

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 $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^{\prime }(0),f^{\prime }(1)$: $$h(x)=f(0)h_{1,0}(x)+f^{\prime }(0)h_{1,1}(x)+f(1)h_{2,0}(x)+f^{\prime }(1)h_{2,1}(x)$$

where the four basis polynomials are: $$\begin{array}{rl}{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{array}$$

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^{\prime }(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{array}{rl}f(x)& =f(x_1)+f^{\prime }(x_1)(x-x_1)+\frac{1}{2}{f}^{″}(x_1)(x-x_1{)}^2+\frac{f(x_2)-y(x_1)-f^{\prime }(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(f^{\prime }(x_1)+\frac{1}{2}{f}^{\prime }(x_2))(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(f^{\prime }(x_1)+f^{\prime }(x_2))(x_2-x_1)+\frac{1}{2}(f^{″}(x_2)-f^{″}(x_1))(x_2-x_1{)}^2}{(x_2-x_1{)}^5}(x-x_1{)}^3(x-x_2{)}^2.\end{array}$$

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 $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: $${\frac{d^l}{d{x}^l}{h}_{i,j}(x)|}_{x=x_k}={\delta }_{i,k}{\delta }_{j,l}=\{\begin{array}{ll}1& \text{if }k=i\text{ and }l=j,\\ 0& \text{otherwise.}\end{array}$$

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 $y_{k,l}^{(i,j)}$ across all nodes $x_k$ and derivative orders $0\le l\le r_k$: $$y_{k,l}^{(i,j)}={\delta }_{i,k}{\delta }_{j,l}=\{\begin{array}{ll}1& \text{if }k=i\text{ and }l=j,\\ 0& \text{otherwise.}\end{array}$$

By passing these target conditions $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 $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\le j\le r_i$ can be constructed directly via polynomial reduction: $$h_{i,j}(x)=(\frac{(x-x_i{)}^j}{j!}{h}_{i,0}(x))modq(x).$$

Taking the remainder modulo $q(x)$ guarantees that $\mathrm{deg}(h_{i,j})<\mathrm{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,\dots ,x_m$, let $n_i$ be the multiplicity of node $x_i$ representing the number of specified values including derivatives, and let $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 $a=\underset{i}{\min }{x}_i$ and $b=\underset{i}{\max }{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 $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

External links

• Hermites Interpolating Polynomial at Mathworld