Hermite Interpolation Polynomial

• Hermite interpolation polynomial is a method that integrates both function values and derivative data to build a unique interpolant across prescribed nodes.
• It employs algebraic formulations, barycentric forms, and Newton divided differences to ensure computational efficiency and numerical stability.
• Its applications span spectral methods, signal processing, and multivariate data fitting, providing high-accuracy error correction and smooth interpolants.

The Hermite interpolation polynomial generalizes Lagrange interpolation by incorporating not only the values of a function at prescribed nodes, but also the values of its derivatives up to specified orders at those nodes. This construction, which lies at the core of approximation theory, is computationally and analytically central to areas including spectral methods, signal processing, multivariate data fitting, and symbolic computation. Its modern expressions leverage algebraic, barycentric, and residue-theoretic constructions as well as stable, efficient numerical schemes.

1. Classical Hermite Interpolation: Foundations and Explicit Formulas

Given distinct nodes $x_0, x_1, \dots, x_n$ and prescribed values $f^{(k)}(x_i)$ for $k=0,1,\dots,m_i-1$ at each $x_i$, the classical Hermite interpolation problem is to construct the unique polynomial $P(x)$ such that

$$P^{(k)}(x_i) = f^{(k)}(x_i), \qquad k=0,\dots,m_i-1,\;\; i=0,\dots,n,$$

with $\deg P < \sum_i m_i$.

The canonical algebraic formula involves Hermite basis polynomials (Durakov et al., 2022):

$$P(x) = \sum_{i=0}^n \sum_{k=0}^{m_i-1} C_{i,k} \, H_{i,k}(x),$$

where

$$H_{i,k}(x) = \frac{1}{k!} \, \frac{d^k}{dx^k}\left[ (x - x_i)^{m_i} \prod_{j\neq i} (x - x_j)^{m_j} \right] \Bigg/ \prod_{j\neq i} (x_i - x_j)^{m_j}.$$

A more computationally advantageous explicit closed form, involving only polynomial derivatives and no rational function derivatives, is given by (Kechriniotis et al., 2011):

$$P(x) = \sum_{j=0}^{n-1} \sum_{k=0}^{m_j-1} c_{k,j} \frac{(x-d_j)^k}{k!} L_j(x),$$

where $L_j(x) = \prod_{k\neq j} (x-d_k)^{m_k}$, and the coefficients $c_{k,j}$ are determined by inversion of a lower-triangular matrix built from derivatives of $L_j$ at $d_j$.

For cases where each $m_i = 1$, this reduces to the standard Lagrange interpolation polynomial.

In practical computation, the Newton divided-difference form is preferred in certain applications, especially when nodes are repeated, as it leads to a triangular system, minimizing computational effort and numerical instability (Kopotun et al., 2020).

2. Barycentric and Fast Evaluation Techniques

Barycentric forms offer superior numerical stability and computational complexity for high-degree or closely-spaced nodes, especially when high-order derivatives are interpolated (Sadiq et al., 2011, Xiang et al., 2014). For nodes $z_1, \ldots, z_K$ with multiplicities $n_1, \ldots, n_K$, $N = n_1+\cdots+n_K$, the Hermite interpolant of degree $N-1$ admits the rational barycentric form:

$$H(z) = \frac{\sum_{k=1}^K \sum_{r=0}^{n_k-1} \omega_{k,r} f_{k,r} B_{k,r}(z)}{\sum_{k=1}^K \sum_{r=0}^{n_k-1} \omega_{k,r} B_{k,r}(z)},$$

where $B_{k,r}(z) = (z-z_k)^{n_k - r}$ and weights $\omega_{k,r}$ are derived efficiently via Newton identities with $O(N^2)$ preprocessing and $O(N)$ updates per additional derivative.

For Jacobi or Chebyshev nodes (Gauss–Jacobi, Gauss–Lobatto–Jacobi), explicit and recursively stable formulas for barycentric weights are obtained by leveraging the Sturm–Liouville structure of the underlying orthogonal polynomials (Xiang et al., 2014). Cancellation of common exponential growth factors yields simplified, numerically robust weights that enable $O(n)$ precomputation and $O(n)$ evaluation (for fixed $m$).

High-order Hermite–Fejér interpolation, as a special case, enjoys similar stability and fast evaluation. For analytic functions, spectral convergence is observed, while algebraic rates pertain for functions with finite regularity.

3. Error Analysis, Pointwise Estimates, and Optimality

The classical error formula for Hermite interpolation at $n+1$ nodes, considering both values and derivatives, is given by (Prentice, 2023): $$E(x) = f(x) - H_{2n+1}(x) = \frac{f^{(2n+2)}(\xi(x))}{(2n+2)!} \prod_{k=0}^{n} (x-x_k)^2$$ for some $\xi(x)\in (x_0, x_n)$, with the location function $\xi(x)$ ("Rolle function") computable via a first-order ODE derived directly from the error representation and solvable numerically (e.g., via explicit Runge–Kutta integrators).

The error can be further polished by fitting the $(2n+2)$-th derivative composed with $\xi(x)$ by a polynomial or spline, effectively constructing a much more accurate corrected Hermite interpolant: $$\widetilde H(x) = H_{2n+1}(x) + \frac{P(x)}{(2n+2)!} \prod_{k=0}^{n}(x-x_k)^2,$$ dramatically reducing global error—by up to 12 orders of magnitude in concrete experiments (Prentice, 2023).

Sharp pointwise error estimates for Hermite interpolation polynomials are given in terms of the local mesh norm and modulus of smoothness: $$|f(x) - P_n(x)| \leq C A \rho_n^r(x) \omega_k(f^{(r)}, \rho_n(x)),$$ with $\rho_n(x) = n^{-1}\sqrt{1-x^2} + n^{-2}$, and are proven to be unimprovable under only Hermite constraints (Kopotun et al., 2020).

4. Multivariate and Symmetric Hermite Interpolation

Multivariate Hermite interpolation extends univariate theory in both algebraic and coordinate directions.

Algebraic-multivariate:

For nodes $\mathbf{x}_i$ and a set of partial derivatives (multi-indices) to be matched at each node, residue-theoretic formulations yield explicit interpolants: $F(\mathbf{z}) = \sum_w \sum_{\boldsymbol{\alpha}} \frac{1}{\boldsymbol{\alpha}!} \partial^\boldsymbol{\alpha}[\det H_w(\mathbf{z})]_{\mathbf{z}=w} (\mathbf{z}-w)^{\boldsymbol{\alpha}} C_{w,\boldsymbol{\alpha}},$ with local conditions ensured via the multidimensional Grothendieck residue (Durakov et al., 2022).

Coordinate-multivariate:

On $n$-dimensional rectilinear grids,

$$f(\mathbf{x}) = \sum_{\mathbf{a} \in A} \big[\Lambda_\mathbf{a}^{-1} \mathbf{T}_\mathbf{a}\big]^T \mathbf{H}_\mathbf{a}(\mathbf{x}),$$

where $\Lambda_\mathbf{a}$ encodes Hermite moments, maintains single-summation closed form and spline-like continuity across block boundaries (Kechriniotis et al., 2023). This yields practical algorithms matching or surpassing state-of-the-art splines in approximation error and smoothness for both regular and scattered data.

Symmetry constraints:

Multivariate symmetric Hermite interpolation leverages generalized (confluent) Vandermonde determinants, resulting in

$$r_A(h)(X) = \sum_{A'} \mathcal{E}_{A'} \frac{D_{A'}(U_n h)(A')}{v(A)} \det V(X \cup (A \setminus A')),$$

guaranteeing symmetry and all Hermite conditions for any partitioned data set with multiplicities (Krick et al., 30 Jan 2025).

5. Alternative and Specialized Forms

Hermite–trigonometric interpolation for periodic functions on uniform grids is efficiently solvable via FFT-based reduction to small frequency-local linear systems, allowing construction of a trigonometric polynomial that interpolates both function and derivative values (Denysiuk, 2019). Precomputation and FFT methods collapse the cost to $O((r+1)M \log M + M(r+1)^3)$ for $M$ nodes and $r$ derivatives.

Multicentric expansions relate analytic function approximation to Hermite interpolation: truncating the series in $p(z)$ at $n$ yields the Hermite interpolant of degree $d(n+1)-1$ with $d(n+1)$ basis polynomials, leading to improved numerical stability compared to direct evaluation in the Hermite basis, particularly when $n$ is large (Nevanlinna et al., 10 Nov 2025).

6. Extensions: Noncommutative and Skew Polynomial Settings

Hermite interpolation extends to free multivariate skew-polynomial rings $$\mathcal{R} = \mathbb{F}[x_1, ..., x_n; \sigma, \delta]$$ over a division ring, with both “right” and “left” $(\sigma, \delta)$-partial derivatives. Existence and uniqueness are characterized by a noncommutative “DP-independence” of nodes and prescribed monomials, and interpolants are constructed via Newton-type recursions with explicit bases for prescribed derivative directions (Donoso et al., 2022).

7. Practical Applications, Numerical Stability, and Advanced Algorithms

In applications, the Hermite interpolant provides the foundation for spectral and pseudospectral discretizations in PDEs, the design of interpolatory splines, post-processing for high-accuracy error correction, and computational quadrature. The use of barycentric forms with optimized weights (Sadiq et al., 2011, Xiang et al., 2014), residue-based closed formulas (Durakov et al., 2022), multicentric expansions (Nevanlinna et al., 10 Nov 2025), and explicit triangular systems (Kechriniotis et al., 2011, Kechriniotis et al., 2023) ensures both algorithmic efficiency and numerical robustness. Practical stability for high-order or near-coalescent nodes is best achieved in barycentric or residue-based settings where spurious growth of coefficients is mitigated.

For multidimensional grids, modern Hermite schemes deliver higher-order smoothness and local continuity, with performance comparable or superior to cubic and generalized splines in data-fitting, simulation, and geometric design tasks (Kechriniotis et al., 2023).

In summary, the Hermite interpolation polynomial and its generalizations constitute a central construct with rigorous algebraic, geometric, and analytic underpinnings. Current research continuously refines fast and stable computational approaches, extends classical theory to new algebraic and analytic contexts, and demonstrates practical superiority over more restrictive Lagrange or spline-based interpolants for high-accuracy, high-continuity approximation.