Fact-checked by Grok 2 weeks ago

Lagrange polynomial

In numerical analysis, the Lagrange interpolating polynomial is the unique polynomial of degree at most n-1 that passes through a given set of n distinct points (x_k, y_k) for k = 1, \dots, n, where the y_k values are typically samples of some underlying function f(x) evaluated at the x_k. This construction provides an explicit way to approximate f(x) at arbitrary points by ensuring exact agreement at the interpolation nodes, making it a cornerstone method for polynomial interpolation. The formula for the Lagrange interpolating polynomial P(x) is given by P(x) = \sum_{i=1}^n y_i \ell_i(x), where each basis polynomial \ell_i(x) is defined as \ell_i(x) = \prod_{\substack{1 \leq j \leq n \\ j \neq i}} \frac{x - x_j}{x_i - x_j}. Each \ell_i(x) satisfies \ell_i(x_k) = \delta_{ik} (the ), ensuring P(x_i) = y_i for all i. Although named after the Italian-French mathematician , who formalized and published the method in his 1795 work Théorie des fonctions analytiques, the formula was originally discovered by Edward Waring in 1779 and independently rediscovered by Leonhard Euler in 1783. Key properties of the Lagrange polynomial include its explicit computability without solving systems of equations, unlike Newton interpolation, though it can suffer from —oscillations near the endpoints—for high-degree approximations on equispaced nodes. It also admits efficient barycentric forms for and evaluation. Applications span numerical methods, where it facilitates and derivative estimation from discrete data; for ; and physics simulations for modeling continuous phenomena from sampled measurements. In engineering and scientific computing, it underpins techniques like and finite element methods.

Fundamentals

Definition

Lagrange interpolation constructs a P(x) of degree at most n that passes through n+1 given distinct points (x_0, y_0), \dots, (x_n, y_n), where the x_i are distinct real or numbers and the y_i are the corresponding values, often y_i = f(x_i) for some f. The interpolating is explicitly given by the P(x) = \sum_{i=0}^n y_i \ell_i(x), where the \ell_i(x) are the Lagrange basis polynomials defined as \ell_i(x) = \prod_{\substack{j=0 \\ j \neq i}}^n \frac{x - x_j}{x_i - x_j}. Each basis polynomial \ell_i(x) is of n and satisfies \ell_i(x_k) = \delta_{ik}, the , ensuring P(x_k) = y_k for each k. This form assumes the interpolation points x_i are distinct, as the denominators in the basis polynomials would otherwise be zero. The Lagrange formulation offers an explicit, direct expression for the without requiring the computation of , unlike the form.

Illustrative Example

To illustrate the construction of a Lagrange , consider the task of finding a polynomial P(x) of degree at most 2 that interpolates the points (0, 1), (1, 2), and (2, 0). The basis polynomials are computed as follows: \ell_0(x) = \frac{(x - 1)(x - 2)}{(0 - 1)(0 - 2)} = \frac{(x - 1)(x - 2)}{2} \ell_1(x) = \frac{(x - 0)(x - 2)}{(1 - 0)(1 - 2)} = \frac{x(x - 2)}{-1} = -x(x - 2) \ell_2(x) = \frac{(x - 0)(x - 1)}{(2 - 0)(2 - 1)} = \frac{x(x - 1)}{2} The interpolating polynomial is then P(x) = 1 \cdot \ell_0(x) + 2 \cdot \ell_1(x) + 0 \cdot \ell_2(x) = \frac{(x - 1)(x - 2)}{2} + 2 \left[ -x(x - 2) \right]. Expanding this yields P(x) = \frac{x^2 - 3x + 2}{2} - 2x(x - 2) = \frac{x^2 - 3x + 2}{2} - 2x^2 + 4x = -\frac{3}{2}x^2 + \frac{5}{2}x + 1. To verify, evaluate P(x) at the interpolation points:
  • P(0) = -\frac{3}{2}(0)^2 + \frac{5}{2}(0) + 1 = 1
  • P(1) = -\frac{3}{2}(1)^2 + \frac{5}{2}(1) + 1 = -\frac{3}{2} + \frac{5}{2} + 1 = 2
  • P(2) = -\frac{3}{2}(4) + \frac{5}{2}(2) + 1 = -6 + 5 + 1 = 0
The values match the given points, confirming the polynomial passes through them. For further intuition, the following table shows P(x) evaluated at additional points between 0 and 2:
xP(x)
0.51.875
1.51.375
This quadratic curve connects the points smoothly, demonstrating how Lagrange interpolation produces an exact fit at the specified locations.

Formulation

Basis Polynomials

The Lagrange basis polynomials, denoted \ell_i(x) for i = 0, 1, \dots, n, serve as the building blocks for the Lagrange scheme given n+1 distinct nodes x_0, x_1, \dots, x_n. Each \ell_i(x) is the unique of at most n such that \ell_i(x_i) = 1 and \ell_i(x_j) = 0 for all j \neq i. This uniqueness follows from the fact that the of of at most n has n+1, and the n+1 interpolation conditions imposed on \ell_i(x) form a linearly set. The explicit construction of \ell_i(x) takes the product form \ell_i(x) = \prod_{\substack{0 \leq j \leq n \\ j \neq i}} \frac{x - x_j}{x_i - x_j}. This formula directly enforces the required properties: when x = x_i, each factor in the numerator equals the corresponding factor in the denominator, yielding \ell_i(x_i) = 1; when x = x_k for k \neq i, the numerator includes a factor of zero while the denominator does not, yielding \ell_i(x_k) = 0. Each basis polynomial \ell_i(x) has exact degree n, arising from the product of n linear terms in the numerator. A key characteristic of the Lagrange basis polynomials is their from the interpolated values y_i = f(x_i); they depend only on the locations x_i. Consequently, the \ell_i(x) can be precomputed for a fixed set of nodes and reused to construct interpolants for any values at those nodes, enhancing computational efficiency in repeated tasks on the same grid.

Interpolation Formula

The Lagrange interpolating polynomial P(x) for a f(x) at n+1 distinct points (x_0, y_0), (x_1, y_1), \dots, (x_n, y_n), where y_i = f(x_i), is constructed by linearly combining the Lagrange basis polynomials \ell_i(x) weighted by the corresponding data values y_i. Specifically, the formula is given by P(x) = \sum_{i=0}^n y_i \ell_i(x), where each basis polynomial is \ell_i(x) = \prod_{\substack{j=0 \\ j \neq i}}^n \frac{x - x_j}{x_i - x_j}. This form arises directly from the properties of the basis polynomials, which are designed such that \ell_i(x_k) = \delta_{ik} (the Kronecker delta, equal to 1 if i = k and 0 otherwise) for k = 0, \dots, n. Substituting x = x_k into the formula yields P(x_k) = \sum_{i=0}^n y_i \delta_{ik} = y_k, ensuring that P(x) exactly matches f(x) at each interpolation node x_k. The formula effectively blends the data values y_i according to their proximity to the evaluation point x: the weight \ell_i(x) is close to 1 when x is near x_i (since the product emphasizes small denominators for nearby points) and close to 0 when x is far from x_i, providing a transition between the given points. The resulting P(x) is a of degree at most n, as each \ell_i(x) is of degree n and the preserves this bound. Direct evaluation of P(x) at a single point requires computing each of the n+1 basis polynomials, each involving n multiplications and divisions, leading to O(n^2) arithmetic operations overall without any preprocessing.

Properties

Uniqueness

The uniqueness theorem states that for any set of n+1 distinct points (x_0, y_0), \dots, (x_n, y_n) in the plane, there exists at most one polynomial P of degree at most n such that P(x_i) = y_i for each i = 0, \dots, n. To prove this, suppose there are two polynomials P and Q, both of at most n, satisfying the interpolation conditions: P(x_i) = y_i and Q(x_i) = y_i for i = 0, \dots, n. Consider the difference R(x) = P(x) - Q(x). Then \deg R \leq n, and R(x_i) = 0 for each of the n+1 distinct points x_i. A nonzero polynomial of at most n can have at most n roots, unless it is the zero polynomial. Since R has n+1 roots, it must be identically zero, so P(x) = Q(x) for all x. This establishes . The proof relies on the fundamental property that a polynomial of degree at most n is uniquely determined by its values at n+1 distinct points, a direct consequence of the and the fact that polynomials are analytic functions with isolated roots unless identically zero. As the Lagrange interpolation formula explicitly constructs a polynomial of degree at most n that satisfies the given conditions, it follows that this construction yields the unique such interpolating polynomial. An alternative perspective views the space \mathcal{P}_n of all polynomials of degree at most n over the real numbers as of n+1, with basis \{1, x, x^2, \dots, x^n\}. The map sending P \in \mathcal{P}_n to the vector (P(x_0), \dots, P(x_n)) \in \mathbb{R}^{n+1} is a linear transformation. For distinct x_i, this map is represented by the Vandermonde matrix, which has full rank n+1 and is thus invertible. The invertibility implies that the map is bijective, so for any target vector (y_0, \dots, y_n), there is exactly one P \in \mathcal{P}_n mapping to it, confirming both existence and uniqueness of the interpolant.

Derivatives

The k-th derivative of the Lagrange interpolating P(x) = \sum_{i=0}^n y_i \ell_i(x) is given by term-by-term : P^{(k)}(x) = \sum_{i=0}^n y_i \ell_i^{(k)}(x), where \ell_i^{(k)}(x) is the k-th of the i-th basis \ell_i(x). This formula holds because is linear and commutes with . Each basis \ell_i(x) has n, so its k-th \ell_i^{(k)}(x) has at most n - k, implying that P^{(k)}(x) has at most n - k. The explicit form of \ell_i^{(k)}(x) arises from repeated application of the to \ell_i(x) = \prod_{j \neq i} \frac{x - x_j}{x_i - x_j}, resulting in expressions that are sums of polynomials of n - k. These terms resemble Lagrange basis polynomials of lower constructed over reduced sets of the original points excluding certain x_j. The Lagrange framework extends naturally to Hermite interpolation, which incorporates derivative values at the interpolation points as additional conditions. For the case of n+1 distinct points x_0, \dots, x_n where both function values f(x_i) and first derivatives f'(x_i) are specified, the unique Hermite interpolating polynomial H(x) of degree at most 2n + 1 takes the form H(x) = \sum_{i=0}^n \left[ f(x_i) \left( 1 - 2 \ell_i'(x_i) (x - x_i) \right) \ell_i^2(x) + f'(x_i) (x - x_i) \ell_i^2(x) \right]. This construction ensures H(x_i) = f(x_i) and H'(x_i) = f'(x_i) for each i, as the squared basis \ell_i^2(x) vanishes to second order at all other points x_j (j ≠ i), while the linear factors adjust the derivative conditions. For higher-order derivatives up to order m_i - 1 at each point (with total degree at most \sum (m_i) - 1), the generalized Hermite form uses higher powers of \ell_i(x) multiplied by polynomials that satisfy the moment conditions at x_i. A practical application of these derivatives is , where P^{(k)}(x) or H^{(k)}(x) approximates the k-th derivative of an underlying f within the , leveraging the known values and the of the basis.

Error Analysis

Remainder Term

When interpolating a f that is n+1 times continuously differentiable on an containing the distinct nodes x_0, x_1, \dots, x_n and the evaluation point x, the Lagrange polynomial P_n(x) of degree at most n satisfies the conditions P_n(x_i) = f(x_i) for i = 0, \dots, n. The , or term, is given by f(x) - P_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^n (x - x_i), where \xi lies in the smallest containing x_0, \dots, x_n, x. The product \prod_{i=0}^n (x - x_i) is known as the nodal polynomial, which quantifies the geometric from x to the interpolation nodes; it vanishes at each x_i, ensuring the error is zero at the nodes as required. The factor \frac{f^{(n+1)}(\xi)}{(n+1)!} captures the local smoothness of f, with the remainder depending on the of the (n+1)-th and the node spacing. This form implies that the approximation error can be bounded by \frac{M}{(n+1)!} \max_{x \in I} \left| \prod_{i=0}^n (x - x_i) \right|, where M = \max |f^{(n+1)}(\xi)| over the relevant interval I, highlighting how denser or better-distributed nodes reduce the error for smooth functions. However, for equispaced nodes over large intervals, the nodal polynomial grows rapidly near the endpoints, leading to the Runge phenomenon—large oscillations in the interpolant that fail to converge uniformly to f, even for analytic functions like f(x) = \frac{1}{1 + 25x^2} on [-1, 1].

Derivation of Remainder

To derive the remainder term in the Lagrange interpolation error, assume that the interpolation points x_0, x_1, \dots, x_n are distinct and lie in an interval containing the evaluation point x, and that the function f is (n+1)-times continuously differentiable on this interval, i.e., f \in C^{n+1}. Let P(x) denote the unique Lagrange interpolating polynomial of degree at most n such that P(x_i) = f(x_i) for i = 0, 1, \dots, n. The error function is defined as e(x) = f(x) - P(x), which satisfies e(x_i) = 0 for each i = 0, 1, \dots, n. Consider the auxiliary function \phi(t) = e(t) - \frac{e(x)}{\omega(x)} \omega(t), where \omega(t) = \prod_{i=0}^n (t - x_i) is the of degree n+1 with roots at the interpolation points. By construction, \phi(x_i) = 0 for i = 0, 1, \dots, n, and also \phi(x) = 0, so \phi(t) has n+2 distinct zeros (assuming x \neq x_i for all i). Since \phi(t) is (n+1)-times differentiable and has n+2 zeros, repeated application of implies that there exists some point \xi in the interval spanned by x_0, \dots, x_n, x such that \phi^{(n+1)}(\xi) = 0. Differentiating \phi(t) n+1 times yields \phi^{(n+1)}(t) = f^{(n+1)}(t) - \frac{e(x)}{\omega(x)} \cdot (n+1)!, because the (n+1)-th derivative of \omega(t) is the constant (n+1)! (as \omega(t) is monic of degree n+1) and P(t) is a of degree at most n, so its (n+1)-th derivative vanishes. Setting \phi^{(n+1)}(\xi) = 0 gives $0 = f^{(n+1)}(\xi) - \frac{e(x)}{\omega(x)} \cdot (n+1)!, so e(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!} \omega(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^n (x - x_i) for some \xi between the minimum and maximum of \{x_0, \dots, x_n, x\}.

Representations

Barycentric Form

The barycentric form provides a computationally efficient and numerically stable alternative to the standard Lagrange interpolation formula for evaluating the interpolating polynomial P(x) at a point x. It reformulates the interpolation by introducing precomputed barycentric weights w_i = \frac{1}{\prod_{j=0, j \neq i}^n (x_i - x_j)} for i = 0, \dots, n, which capture the denominators of the Lagrange basis polynomials. The first barycentric form expresses the interpolant as P(x) = \frac{\sum_{i=0}^n \frac{w_i}{x - x_i} y_i}{\sum_{i=0}^n \frac{w_i}{x - x_i}}, where the sums are taken over the distinct interpolation nodes x_0, \dots, x_n and corresponding values y_0, \dots, y_n. This form requires O(n^2) operations to compute the weights w_i once, followed by O(n) operations per evaluation of P(x), making it suitable for repeated evaluations. Compared to the direct Lagrange form, the barycentric approach is superior in , particularly for equispaced nodes or when evaluating near the interpolation points, as it sidesteps the explicit formation of large intermediate products that can lead to subtractive cancellation. The preprocessing cost of O(n^2) is offset by the linear-time evaluation, rendering it the preferred method for practical implementations of .

Linear Algebra View

The interpolating P(x) of degree at most n that passes through the distinct points (x_i, y_i) for i = 0, \dots, n can be expressed in the as P(x) = \sum_{k=0}^n c_k x^k, where the coefficient vector \mathbf{c} = (c_0, \dots, c_n)^T satisfies the V \mathbf{c} = \mathbf{y}. Here, V is the (n+1) \times (n+1) with entries V_{i,j} = x_i^j for i,j = 0, \dots, n, and \mathbf{y} = (y_0, \dots, y_n)^T is the vector of data values. This system has a unique solution because V is invertible whenever the x_i are distinct, ensuring the existence and uniqueness of the interpolating of degree at most n. The Lagrange form provides an explicit expression for P(x) without directly solving for \mathbf{c}: P(x) = \sum_{i=0}^n y_i \ell_i(x), where each Lagrange basis polynomial \ell_i(x) satisfies \ell_i(x_j) = \delta_{i j}. In the , \ell_i(x) = \sum_{k=0}^n a_{k i} x^k, and the matrix A = (a_{k i}) is precisely the V^{-1}, with the coefficients of \ell_i(x) given by the i-th column of V^{-1}. Consequently, the coefficients of P(x) are \mathbf{c} = V^{-1} \mathbf{y}, highlighting how the Lagrange basis implicitly applies the of the to the in the . The condition number \kappa(V) of the Vandermonde matrix grows exponentially with n for equispaced interpolation points, often as O(2^n / \sqrt{n}) or worse, which explains the severe numerical instability when solving V \mathbf{c} = \mathbf{y} directly for the coefficients, as small perturbations in \mathbf{y} or rounding errors amplify dramatically in \mathbf{c}. This ill-conditioning is closely tied to the Lebesgue constant \Lambda_n, which measures the maximum amplification of errors in the operator; specifically, \kappa(V_n) = \|V_n\| \Lambda_n, where \|V_n\| is the matrix , and both quantities become prohibitively large for equispaced points, contributing to phenomena like the Runge oscillation. More generally, expressing the interpolating polynomial in alternative bases, such as the basis via , corresponds to a from the form, typically involving a triangular factorization of the transposed (e.g., V^T = [LU](/page/Lu)), where L maps Newton polynomials to monomials and U performs the reverse transformation. This basis change often yields more stable algorithms, as the resulting matrices are better conditioned for hierarchical evaluation.

Extensions

Finite Fields

Lagrange interpolation extends naturally to finite fields \mathbb{F}_q, where q is a , allowing the construction of the unique polynomial of degree at most n-1 that passes through n distinct points (x_i, y_i) with x_i \in \mathbb{F}_q and n \leq q. The formula holds identically to the classical case over the reals or complexes, provided the characteristic of the field exceeds n or, more precisely, the points are distinct to ensure all denominators x_i - x_j (for i \neq j) are nonzero and invertible in \mathbb{F}_q. The Lagrange basis polynomials are defined as \ell_i(x) = \prod_{\substack{0 \leq j \leq n-1 \\ j \neq i}} \frac{x - x_j}{x_i - x_j}, with the interpolating polynomial given by p(x) = \sum_{i=0}^{n-1} y_i \ell_i(x). This construction relies on the field properties, where division is multiplication by the modular inverse, and uniqueness follows from the fact that a nonzero polynomial of degree at most n-1 can have at most n-1 roots in \mathbb{F}_q. In coding theory, Lagrange interpolation plays a key role in the decoding of Reed-Solomon codes, which are evaluation codes over \mathbb{F}_q where codewords correspond to evaluations of polynomials of degree less than k at n distinct points. During decoding, if the received word agrees with the codeword on at least k positions, the original message polynomial can be recovered via Lagrange interpolation on those reliable points, enabling error correction up to (n-k)/2 errors. Extensions to list decoding, such as the Guruswami-Sudan algorithm, further employ interpolation techniques—albeit bivariate—to generate a list of possible polynomials consistent with more erroneous positions, achieving rates beyond the unique decoding radius. A prominent cryptographic application arises in scheme, where a secret is encoded as the constant term of a random of t over \mathbb{F}_q, and shares are distributed as evaluations at distinct points; reconstruction requires t+1 shares and uses to recover the and thus the secret. Challenges in finite fields include the non-uniqueness of high- representations as functions, since in \mathbb{F}_q, any function can be represented by a of at most q-1 due to the identity x^q = x; thus, for n > q, yields the minimal- representative, but computations must respect the field size to avoid redundant high- terms. Unlike over the reals, there is no Runge phenomenon, as the finite domain prevents oscillation issues associated with equidistant points on unbounded intervals.

Applications

Lagrange polynomials play a central role in for function approximation, where they construct a unique polynomial of degree at most n-1 that interpolates n given data points, providing a smooth estimate of the underlying function between those points. This approach is particularly useful for data fitting and extrapolation in scientific , as it ensures exact at the interpolation nodes while allowing evaluation at arbitrary points. In numerical , Lagrange interpolants form the basis for methods like Newton-Cotes formulas, where the of the interpolating over an interval approximates the definite of the original , enabling efficient computation for irregular . For solving ordinary differential equations (ODEs), Lagrange interpolation is employed in methods, such as those deriving implicit Runge-Kutta schemes, by assuming the solution follows a that satisfies the ODE at specific points within each step. In , Bézier curves represent a specialized application of principles through barycentric coordinates, where control points are weighted to generate curves used in rendering and surface modeling. Barycentric coordinates, closely tied to the Lagrange basis functions, facilitate by interpolating texture coordinates across triangle vertices during rasterization, ensuring perspective-correct sampling of image data on 3D surfaces. Lagrange interpolation finds extensive use in for resampling discrete-time signals to non-integer rates, where it reconstructs continuous signals from samples and evaluates them at new points to achieve bandlimited without . In , Lagrange-based interpolators serve as (FIR) filters for fractional delay and multirate processing, optimizing flatness and attenuation in applications like audio and communications. Despite these strengths, high-degree Lagrange polynomials exhibit numerical instability, particularly the Runge phenomenon, which causes large oscillations near the endpoints when interpolating equispaced points on functions like $1/(1 + 25x^2). Selecting , clustered near the interval endpoints, significantly reduces these errors and improves convergence for smooth functions. For applications requiring global smoothness, such as over large datasets, piecewise polynomial methods like cubic splines are preferred over single high-degree Lagrange interpolants, as they maintain lower curvature and avoid while ensuring C^2 continuity. In recent contexts, Lagrange polynomials approximate polynomial kernels in support vector machines and Gaussian processes, enabling efficient computation of nonlinear feature mappings without explicit high-dimensional expansions. Post-2020 developments in have integrated Lagrange polynomial encoding into variational quantum algorithms, facilitating hybrid quantum-classical models for tasks like on noisy intermediate-scale quantum hardware.

Historical Context

The development of polynomial interpolation methods traces back to the , when introduced the concept of as part of his work on finite differences for approximating continuous functions, particularly in the context of and astronomical computations. Newton's approach to finite differences for interpolation, developed in the late and detailed in his unpublished manuscripts and later works such as the 1711 Methodus differentialis, provided a foundational framework for constructing interpolating polynomials using successive differences, enabling efficient computation for tabulated data without explicit polynomial forms. In the late , the specific form now known as the Lagrange interpolating emerged. Edward Waring first published the formula in 1779 within his treatise on the resolution of equations, presenting it as a direct method to construct a passing through given points. Leonhard Euler independently rediscovered the method in 1783, applying it to problems in analysis and numerical computation. formalized and popularized the approach in 1795, in his Théorie des fonctions analytiques (published 1797), where he used it to address the numerical solution of algebraic equations through based on finite differences. The saw significant evolutions in the Lagrange method for enhanced and broader applications. In the 1970s, Carl de Boor advanced the barycentric representation of Lagrange interpolation, introducing weights that improve computational efficiency and reduce rounding errors, particularly in the context of spline approximations as detailed in his 1978 book A Practical Guide to Splines. Concurrently, extensions to finite fields gained prominence; Irving S. Reed and Gustave Solomon applied Lagrange interpolation over finite fields in 1960 to develop Reed-Solomon error-correcting codes, revolutionizing for reliable data transmission. Lagrange polynomials have profoundly influenced modern , serving as the basis for techniques refined by de Boor during the 1970s, which enable smooth approximations in and data fitting. Today, the method underpins implementations in leading numerical libraries, such as SciPy's lagrange function in , facilitating widespread use in scientific .