Hermite interpolation
Hermite interpolation is a method in numerical analysis for constructing a polynomial that matches both the values and specified derivatives of a given function at a set of distinct points, thereby providing a smoother approximation than standard interpolation techniques when derivative information is available.[1] Named after the French mathematician Charles Hermite, who first posed and solved the problem in 1878, this approach generalizes classical Lagrange interpolation by incorporating osculatory conditions, where the interpolating polynomial and its derivatives up to certain orders agree with those of the original function at the interpolation nodes.[2] For n+1 distinct points with function values and first derivatives specified, the unique interpolating polynomial has degree at most $2n+1, satisfying $2n+2 conditions in total.[3] This method is particularly valuable in applications requiring high accuracy and continuity, such as cubic Hermite splines in computer graphics and animation for smooth curve fitting, numerical integration schemes like Hermite's rule, and simulations of wave propagation in scientific computing.[4][5][6] The error term for Hermite interpolation can be expressed using the next higher derivative of the function, analogous to the Newton form, enabling precise bounds on approximation quality.[7]Problem Formulation
Statement of the Problem
Hermite interpolation is a method in numerical analysis for constructing a polynomial that matches both the values and specified derivatives of a given function at a set of distinct points. Given a smooth function f and n distinct interpolation points x_1, x_2, \dots, x_n \in \mathbb{R}, along with non-negative integers k_1, k_2, \dots, k_n specifying the orders of derivatives to match at each point, the goal is to find a polynomial p of minimal degree such that p^{(j)}(x_i) = f^{(j)}(x_i), \quad j = 0, 1, \dots, k_i, \quad i = 1, 2, \dots, n, where p^{(0)}(x_i) := p(x_i) = f(x_i) and p^{(j)} denotes the j-th derivative of p. This imposes a total of m + 1 := \sum_{i=1}^n (k_i + 1) interpolation conditions, so the polynomial p belongs to the space \Pi_m of polynomials of degree at most m = \sum_{i=1}^n (k_i + 1) - 1.[8] The notation for these conditions highlights the Hermite data: at each x_i, the function value f(x_i) and the derivatives f'(x_i), \dots, f^{(k_i)}(x_i) are prescribed, allowing for varying derivative orders across points (e.g., k_i = 1 for first derivatives only, or higher for greater smoothness). This setup generalizes the classical Lagrange interpolation problem, which corresponds to the special case where all k_i = 0, matching only function values without derivatives. By incorporating derivative information, Hermite interpolation produces smoother approximations, particularly useful for applications requiring continuity in both the interpolant and its derivatives, such as in spline construction or numerical differentiation.[8] In essence, the Hermite interpolant at each point x_i agrees with the local Taylor polynomial of f centered at x_i up to order k_i, but globally combines these local behaviors into a single polynomial of controlled degree.[8]Uniqueness Theorem
The uniqueness of the Hermite interpolating polynomial is a fundamental result that guarantees the existence of a single polynomial satisfying the specified interpolation conditions. Consider distinct points x_1, \dots, x_s \in \mathbb{R} and nonnegative integers k_1, \dots, k_s such that m = \sum_{i=1}^s (k_i + 1) - 1. For a given smooth function f, there exists a unique polynomial p \in \Pi_m (the space of polynomials of degree at most m) satisfying p^{(j)}(x_i) = f^{(j)}(x_i), \quad 0 \leq j \leq k_i, \quad 1 \leq i \leq s. [9][10] This uniqueness follows from a dimension-counting argument combined with the linear independence of the associated evaluation functionals. The vector space \Pi_m has dimension m+1. The interpolation conditions impose exactly m+1 linear constraints on p, given by the total number of specified values and derivatives: \sum_{i=1}^s (k_i + 1) = m + 1. These constraints define a linear map from \Pi_m to \mathbb{R}^{m+1} via the evaluation functionals \{\delta_{x_i}^{(j)} \mid 0 \leq j \leq k_i, \, 1 \leq i \leq s \}, where \delta_x^{(j)}(p) = p^{(j)}(x). Since the domain and codomain have the same finite dimension, it suffices to show that this map is injective (or surjective) to establish bijectivity and thus uniqueness.[11][9] To prove injectivity, suppose q \in \Pi_m satisfies the homogeneous conditions q^{(j)}(x_i) = 0 for all specified i, j. Then q has zeros of multiplicity at least k_i + 1 at each x_i, so q(x) = [\omega(x)](/page/Omega_X) r(x) where \omega(x) = \prod_{i=1}^s (x - x_i)^{k_i + 1} has degree m + 1 and r \in \Pi_m. But \deg q \leq m, which forces r \equiv 0 and hence q \equiv 0. Alternatively, the linear independence of the functionals \{\delta_{x_i}^{(j)}\} on \Pi_m ensures that the kernel is trivial: if \sum c_{i,j} \delta_{x_i}^{(j)}(p) = 0 for all p \in \Pi_m, then the coefficients c_{i,j} = 0. This independence holds because the points x_i are distinct and the derivatives up to order k_i form a basis for the dual space under these conditions.[11][10] From a linear algebra perspective, represent p(x) = \sum_{\ell=0}^m a_\ell x^\ell. The interpolation conditions yield the linear system V \mathbf{a} = \mathbf{b}, where \mathbf{a} = (a_0, \dots, a_m)^T, \mathbf{b} collects the values f^{(j)}(x_i), and V is the (m+1) \times (m+1) confluent Vandermonde matrix with entries involving powers and derivatives evaluated at the x_i. For distinct x_i, this matrix is nonsingular, confirming a unique solution \mathbf{a}.[10][9]Construction Methods
Lagrange Basis Method
The Lagrange basis method constructs the Hermite interpolating polynomial using a collection of basis polynomials l_{i,j}(x) tailored to satisfy both function values and specified derivatives at the interpolation points x_1, \dots, x_n, where up to the k_i-th derivative is given at each x_i. These basis polynomials generalize the classical Lagrange basis by enforcing zero conditions with appropriate multiplicities at all points except x_i, while at x_i, the j-th basis satisfies l_{i,j}^{(s)}(x_i) = j! \, \delta_{s j} for s = 0, \dots, k_i, ensuring the correct scaling for the coefficients involving derivatives. The Hermite interpolant p(x) of degree at most \sum_{i=1}^n (k_i + 1) - 1 is expressed explicitly as p(x) = \sum_{i=1}^n \left[ f(x_i) \, l_{i,0}(x) + \sum_{j=1}^{k_i} \frac{f^{(j)}(x_i)}{j!} \, l_{i,j}(x) \right]. This form leverages the linearity of polynomials: each term f(x_i) \, l_{i,0}(x) contributes the value at x_i while vanishing with derivatives at other points, and the derivative terms \frac{f^{(j)}(x_i)}{j!} \, l_{i,j}(x) contribute the j-th derivative at x_i without affecting lower or higher specified derivatives there or any conditions elsewhere.[10] To build the basis, first define the auxiliary function \omega_i(x) = \prod_{\substack{m=1 \\ m \neq i}}^n \frac{(x - x_m)^{k_m + 1}}{(x_i - x_m)^{k_m + 1}}, which satisfies \omega_i(x_i) = 1 and \omega_i^{(s)}(x_m) = 0 for s = 0, \dots, k_m and all m \neq i. This \omega_i(x) encodes the required zero multiplicities at the other interpolation points. The basis polynomials are then constructed starting from the highest order at each x_i and proceeding recursively downward. For the highest derivative, l_{i,k_i}(x) = (x - x_i)^{k_i} \, \omega_i(x), which yields l_{i,k_i}^{(s)}(x_i) = 0 for s < k_i and l_{i,k_i}^{(k_i)}(x_i) = k_i!, as required. For lower orders j = 0, \dots, k_i - 1, l_{i,j}(x) = (x - x_i)^j \, \omega_i(x) - \sum_{s = j+1}^{k_i} \binom{s}{j} \, \omega_i^{(s - j)}(x_i) \, l_{i,s}(x). This recursion subtracts projections onto higher-order bases to enforce l_{i,j}^{(s)}(x_i) = 0 for s \neq j and l_{i,j}^{(j)}(x_i) = j!, while preserving the zero conditions at other points via the factor \omega_i(x). Unrolling the recursion yields an explicit sum form for each l_{i,j}(x) = \omega_i(x) \sum_{s=0}^{k_i - j} c_{j,s} (x - x_i)^s, where the coefficients c_{j,s} depend on powers of (x - x_i) and derivatives of \omega_i evaluated at x_i.[10] The resulting basis polynomials \{ l_{i,j}(x) \} form a Schauder basis for the space of polynomials of degree at most \sum (k_i + 1) - 1, with each satisfying the isolated interpolation conditions at x_i and complete vanishing (up to the specified multiplicities) at all other x_m. This ensures the uniqueness of the Hermite interpolant, as referenced in the uniqueness theorem, since the basis spans the full solution space and matches all \sum (k_i + 1) conditions. The method, while explicit, can be computationally intensive for high multiplicities due to the need to compute derivatives of \omega_i(x), but it provides a direct product-based formula analogous to the classical Lagrange approach.[10]Newton Divided Difference Method
The Newton divided difference method provides an efficient and numerically stable approach to constructing the Hermite interpolating polynomial by extending the classical divided difference interpolation to handle derivative conditions through confluent divided differences. In this framework, the interpolation points x_i are treated with multiplicity m_i = k_i + 1, where k_i is the highest derivative order specified at x_i, leading to a total of N+1 = \sum m_i conditions for a polynomial of degree at most N. The method leverages a recursive computation of divided differences, where repeated evaluations at the same point incorporate higher-order derivatives, ensuring the polynomial satisfies both function values and derivatives at the nodes.[12] For confluent divided differences at a repeated point x_i, the k-th order divided difference is defined as f[x_i, x_i, \dots, x_i] (k+1 times) = \frac{f^{(k)}(x_i)}{k!}, which arises as the limit of the standard divided difference formula when distinct points converge to x_i. This adaptation allows the divided difference table to be built by listing each x_i repeated m_i times in the sequence of nodes z_0, z_1, \dots, z_N, where the z_j follow the order of the points with their multiplicities. The zeroth-order differences are set to f[z_j] = f(z_j) for the initial occurrence in each group, while higher-order confluent differences are directly assigned using the derivative formula when consecutive z_j are equal; otherwise, the standard recursive formula f[z_j, \dots, z_{j+\ell}] = \frac{f[z_{j+1}, \dots, z_{j+\ell}] - f[z_j, \dots, z_{j+\ell-1}]}{z_{j+\ell} - z_j} is applied, with care taken for the indeterminate form when denominators vanish by using the derivative definitions.[8][12] The divided difference table is constructed iteratively, starting from the zeroth column with function values and derivatives, and filling subsequent diagonal entries (the coefficients) level by level. For instance, in the case of first derivatives at two points (m_0 = m_1 = 2), the node sequence is z_0 = z_1 = x_0, z_2 = z_3 = x_1; the table begins with f[z_0] = f(x_0), f[z_1] = f(x_0), f[z_2] = f(x_1), f[z_3] = f(x_1), then the first-order differences include f[z_0, z_1] = f'(x_0), f[z_2, z_3] = f'(x_1), and mixed differences are computed recursively, yielding coefficients on the main diagonal for the Newton basis. This process scales linearly in the number of conditions for sequential computation, promoting stability over direct methods like Lagrange basis products, especially for higher multiplicities.[13][8] The resulting Hermite interpolant in Newton form is expressed as p(x) = \sum_{k=0}^{N} a_k \prod_{j=0}^{k-1} (x - z_j), where the coefficients a_k = f[z_0, z_1, \dots, z_k] are the entries from the k-th diagonal of the table, and the product incorporates the repeated nodes to account for multiplicities, ensuring the polynomial matches the specified derivatives. To compute the coefficients algorithmically, initialize the table with the repeated nodes and data: set the zeroth column for function values, then iteratively compute higher-order differences using the recursive formula, substituting derivative-based confluent values whenever consecutive nodes coincide (e.g., for the second-order at a double point, f[z_i, z_i, z_{i+1}] = \frac{f[z_i, z_{i+1}] - f'(z_i)}{z_{i+1} - z_i} if z_{i+1} \neq z_i, but directly f[z_i, z_i, z_i] = f''(z_i)/2 for triples). This yields the coefficients a_k efficiently, with the nested product form allowing hierarchical evaluation and addition of terms for adaptive implementations. When all k_i = 0, the method reduces to the standard Newton divided difference interpolation for distinct points.[12][14]Chinese Remainder Theorem Method
The Chinese Remainder Theorem provides an algebraic framework for constructing the Hermite interpolating polynomial by decomposing the problem into local solutions modulo pairwise coprime ideals in the polynomial ring R, where R is the base field of real or complex numbers. The interpolation conditions require the polynomial p(x) to match the function f(x) and its derivatives up to order k_i at distinct points x_1, \dots, x_n, which translates to the congruence p(x) \equiv T_i(x) \pmod{(x - x_i)^{k_i + 1}} for each i = 1, \dots, n, where T_i(x) is the Taylor polynomial of f at x_i of degree at most k_i. The relevant ideal is I = \prod_{i=1}^n (x - x_i)^{k_i + 1}, and since the factors (x - x_i)^{k_i + 1} are pairwise coprime (their generators share no common roots), the Chinese Remainder Theorem implies that R/I \cong \prod_{i=1}^n R / (x - x_i)^{k_i + 1}. This isomorphism ensures a unique solution p(x) of degree less than \deg I = \sum_{i=1}^n (k_i + 1) satisfying all local conditions, aligning with the uniqueness theorem for Hermite interpolation. To construct p(x), first compute the local Taylor polynomials q_i(x) = T_i(x) for each i, each of degree less than k_i + 1, so that q_i(x_i^{(j)}) = f^{(j)}(x_i) for j = 0, \dots, k_i (where x_i^{(j)} denotes the j-th derivative evaluation). These satisfy q_i(x) \equiv T_i(x) \pmod{(x - x_i)^{k_i + 1}}. Next, for each i, define w_i(x) = \prod_{m \neq i} (x - x_m)^{k_m + 1}, which vanishes to the required multiplicity at all points except x_i, ensuring w_i(x) \equiv 0 \pmod{\prod_{m \neq i} (x - x_m)^{k_m + 1}}. Compute the modular inverse s_i(x) of w_i(x) modulo (x - x_i)^{k_i + 1}, a polynomial of degree less than k_i + 1 such that w_i(x) s_i(x) \equiv 1 \pmod{(x - x_i)^{k_i + 1}}; this inverse exists because w_i(x_i) \neq 0. The global interpolant is then formed by linear combination: p(x) = \sum_{i=1}^n q_i(x) \cdot w_i(x) \cdot s_i(x). Each term q_i(x) w_i(x) s_i(x) satisfies the local condition at x_i (since q_i w_i s_i \equiv T_i \cdot 1 \equiv T_i \pmod{(x - x_i)^{k_i + 1}}) and vanishes at all other points to the required orders (since w_i \equiv 0 there, and s_i is low-degree). The inverses s_i(x) can be found efficiently using the extended Euclidean algorithm adapted to power series truncation modulo (x - x_i)^{k_i + 1}. This method excels in theoretical contexts, as the ring decomposition directly proves existence and uniqueness without explicit basis constructions, facilitating analysis in commutative algebra. It is also advantageous computationally when interpolation points are clustered—meaning higher k_i at fewer distinct x_i—since local computations modulo high powers can leverage efficient series expansions or modular lifting techniques, reducing overall complexity compared to global basis methods for ill-conditioned cases.Special Cases
Case with First Derivatives Only
In the case where only function values and first derivatives are specified at each of n distinct interpolation points x_0 < x_1 < \dots < x_{n-1}, the Hermite interpolation problem seeks a unique polynomial p(x) of degree at most $2n-1 satisfying p(x_i) = f(x_i) and p'(x_i) = f'(x_i) for i = 0, \dots, n-1.[8] This setup provides $2n conditions, ensuring existence and uniqueness of the interpolant within the space of polynomials of degree less than $2n.[15] Such interpolation is foundational for constructing smooth approximations, particularly in spline-based methods where continuity of the function and its derivative is required. The Newton divided difference form offers a practical construction for this case, building on the general divided difference method by incorporating repeated points to account for derivatives.[8] Specifically, the divided difference table alternates function values and scaled derivatives: the zeroth-order differences are f[x_i] = f(x_i), the first-order confluent differences at each repeated point are f[x_i, x_i] = f'(x_i), and higher-order differences are computed via the standard recursive formula f[x_j, \dots, x_{j+k}] = \frac{f[x_{j+1}, \dots, x_{j+k}] - f[x_j, \dots, x_{j+k-1}]}{x_{j+k} - x_j} across the points, with confluent values set directly for repeated arguments.[15] This yields the coefficients for the Newton polynomial p(x) = a_0 + a_1 (x - z_0) + \dots + a_{2n-1} \prod_{j=0}^{2n-2} (x - z_j), where the z_j is the sequence of nodes with each x_i repeated twice (i.e., z = [x_0, x_0, x_1, x_1, \dots, x_{n-1}, x_{n-1}]) and the a_k are the divided differences f[z_0, \dots, z_k].[8] For the common subcase of n=2 points x_0 < x_1, the interpolant is a cubic polynomial of degree at most 3, often expressed in terms of four explicit basis functions that isolate contributions from each endpoint value and slope.[15] Let h = x_1 - x_0 and t = \frac{x - x_0}{h}. The basis functions are: \begin{align*} h_{0,0}(t) &= (1 + 2t)(1 - t)^2, \\ h_{1,0}(t) &= t^2 (3 - 2t), \\ h_{0,1}(t) &= t (1 - t)^2, \\ h_{1,1}(t) &= t^2 (t - 1). \end{align*} The interpolating polynomial is then p(x) = h_{0,0}(t) f(x_0) + h_{1,0}(t) f(x_1) + h f_{0,1}(t) f'(x_0) + h h_{1,1}(t) f'(x_1), where the factor h scales the derivative basis to match units.[15] These basis functions satisfy h_{i,0}(t_i) = 1 and h_{i,0}(t_j) = 0 for i \neq j (with t_0 = 0, t_1 = 1), and h_{i,0}'(t_k) = 0 for k = 0,1, and similarly for the derivative bases with zero value at both endpoints, zero derivative at the opposite endpoint, and unit derivative at their respective endpoint.[16] When applied piecewise over consecutive intervals with consistent derivative values at shared knots, this construction yields cubic Hermite splines that achieve C^1 continuity, meaning the function and its first derivative are continuous across the points while higher derivatives may discontinue.[8] This property makes the method ideal for applications requiring smooth curves, such as computer graphics and numerical simulation, without the added complexity of higher derivatives.[15]Cases with Higher Derivatives
In Hermite interpolation, cases involving higher derivatives extend the conditions beyond function values and first derivatives by specifying derivatives up to order k_i \geq 2 at interpolation points x_i. This requires treating each point with multiplicity \mu_i = k_i + 1 in the divided difference table, where the divided differences for coincident arguments incorporate the higher-order derivatives via the relation f[x_i, x_i, \dots, x_i] (with m+1 repetitions) equals f^{(m)}(x_i)/m!. For instance, when k_i = 2, the triple divided difference satisfies f[x_i, x_i, x_i] = f''(x_i)/2!, and higher multiplicities follow analogously from Taylor expansion limits of the divided difference definition.[7][12] The resulting interpolating polynomial has degree at most \sum_i (k_i + 1) - 1, reflecting the total number of interpolation conditions minus one. For a global polynomial, this ensures exact matching of the specified derivatives, but in piecewise constructions—such as higher-order splines—the approach allows for C^{\max k_i} smoothness across knots, provided the pieces are joined with compatible higher derivatives. This elevated smoothness is particularly valuable in applications requiring precise curvature or acceleration modeling, like trajectory planning.[7][1] However, incorporating higher k_i introduces challenges, including increased condition numbers due to the sensitivity of higher derivatives to perturbations in data or points, which can amplify numerical instability in computations. Additionally, the higher polynomial degree exacerbates oscillatory behavior, akin to the Runge phenomenon, potentially leading to poor approximation away from the data points. Such higher multiplicities prove useful near singularities, where they enable correction terms to mitigate Gibbs-like oscillations and preserve jump conditions in the function or its derivatives, achieving improved convergence rates like O(h^4) for adapted interpolants.[17][18] The Newton divided difference form adapts naturally to these cases by constructing the sequence of nodes with repetitions according to each \mu_i, yielding basis polynomials with repeated factors (x - x_i)^{k_i + 1} in the higher-order products \prod_{j=0}^{k-1} (x - z_j), where the z_j include \mu_i copies of x_i. This structure facilitates efficient evaluation and extension from the first-derivative case, though care must be taken in the divided difference table to handle the confluent entries accurately.[7][12]Illustrative Examples
Cubic Hermite Interpolation Example
Consider a simple example of cubic Hermite interpolation using two interpolation points x_0 = 0 and x_1 = 1, with specified function values f(0) = 1, f'(0) = 0, f(1) = 2, and f'(1) = 1. This setup requires a cubic polynomial p(x) of degree at most 3 that satisfies these four conditions. The polynomial can be constructed using the standard cubic Hermite basis functions, which are derived for the interval [x_0, x_1] with h = x_1 - x_0 = 1. Let t = (x - x_0)/h = x. The basis functions are: h_{0,0}(t) = 2t^3 - 3t^2 + 1, \quad h_{0,1}(t) = t^3 - 2t^2 + t, h_{1,0}(t) = -2t^3 + 3t^2, \quad h_{1,1}(t) = t^3 - t^2. These satisfy h_{i,j}(x_k) = \delta_{ik} for function values (j=0) and the derivative conditions scaled by h at the endpoints (for j=1). The interpolant is then p(x) = h_{0,0}(x) \cdot 1 + h_{0,1}(x) \cdot 0 + h_{1,0}(x) \cdot 2 + h_{1,1}(x) \cdot 1. Substituting the basis functions yields p(x) = (2x^3 - 3x^2 + 1) + 2(-2x^3 + 3x^2) + (x^3 - x^2) = -x^3 + 2x^2 + 1. Verification confirms p(0) = 1, p'(0) = 0, p(1) = 2, and p'(1) = 1, where p'(x) = -3x^2 + 4x. Alternatively, the same polynomial can be obtained via the Newton divided-difference form, using confluent divided differences to incorporate the derivatives. For the repeated nodes z_0 = z_1 = 0 and z_2 = z_3 = 1, the divided-difference table is constructed as follows:| Order | z_0 = 0 | z_1 = 0 | z_2 = 1 | z_3 = 1 |
|---|---|---|---|---|
| 0 | f[z_0] = 1 | f[z_2] = 2 | ||
| 1 | f[z_0, z_1] = f'(0) = 0 | f[z_1, z_2] = (2 - 1)/1 = 1 | f[z_2, z_3] = f'(1) = 1 | |
| 2 | f[z_0, z_1, z_2] = (1 - 0)/1 = 1 | f[z_1, z_2, z_3] = (1 - 1)/1 = 0 | ||
| 3 | f[z_0, z_1, z_2, z_3] = (0 - 1)/1 = -1 |