Mathematical foundations of Lagrange Polynomial Interpolation.
Polynomial interpolation is a bedrock discipline of numerical analysis and scientific computing, forming the mathematical foundation for finite element shape functions, numerical quadrature rules, and predictor-corrector differential equation solvers. Given n + 1 distinct points (x_0, y_0), ..., (x_n, y_n), the Weierstrass approximation theorem guarantees that continuous functions can be approximated arbitrarily well by polynomials, but polynomial interpolation places the stronger demand of exact nodal colocation P_n(x_j) = y_j. In this laboratory, we study the construction of Lagrange interpolating polynomials, analyze the cardinal basis functions, and rigorously track the attenuation of the supremum norm error as the mesh spacing h is systematically halved, validating the theoretical O(h^{n+1}) error asymptotic rate.
Viva voce questions & analytical theory
What is the classical error representation theorem for polynomial interpolation?
For f ∈ C^{n+1}[a, b], the pointwise interpolation error satisfies f(x) - P_n(x) = [f^{(n+1)}(ξ_x) / (n+1)!] · ∏_{j=0}^n (x - x_j), where ξ_x lies in the minimal interval containing x and the interpolation nodes {x_j}.
Why is the naive monomial basis 1, x, x², ... ill-conditioned for interpolation?
The monomial basis generates a Vandermonde coefficient matrix whose condition number grows exponentially with n (κ(V) ~ e^{C·n}), leading to catastrophic loss of numerical precision during matrix inversion.
How does the barycentric formula improve computational efficiency and stability?
The second barycentric interpolation formula computes P_n(x) = [∑ (w_j / (x - x_j)) f_j] / [∑ (w_j / (x - x_j))] in O(n) operations rather than O(n²), with forward backward stability matching best-practice floating point algorithms.
What is the Lebesgue constant and why does it matter for uniform grids?
The Lebesgue constant Λ_n = max_{x} ∑ |L_{n,j}(x)| governs interpolation stability via ||f - P_n||_∞ ≤ (1 + Λ_n) inf_p ||f - p||_∞. For uniform nodes, Λ_n grows exponentially (~ 2^n / (e·n·ln(n))), causing severe edge oscillations.
How does the Faber theorem restrict polynomial interpolation on general grids?
Faber's theorem (1914) proves that for any predefined triangular nodal array on [a, b], there always exists a continuous function f ∈ C[a, b] whose interpolating polynomials fail to converge uniformly to f.
Why does the convergence order α = 4.0 emerge in this fixed-degree mesh refinement experiment?
When the polynomial degree is fixed at n = 3 (cubic interpolation) and the domain sub-interval spacing h is systematically halved, the error bound satisfies ||f - P_3||_∞ ≤ C · h^4 · ||f^{(4)}||_∞, yielding asymptotic convergence order α = 4.0.
What techniques resolve boundary divergence in high-degree polynomial fitting?
Switching from equidistant nodes to Chebyshev or Gauss-Lobatto nodes reduces Lebesgue constant growth to logarithmic (Λ_n ~ (2/π) ln(n)), or using piecewise spline interpolation with low-degree polynomials avoids high-degree instability altogether.
Algorithm precautions and stability rules
- Always verify that all nodal coordinates x_0, x_1, ..., x_n are strictly distinct to prevent division by zero in the basis denominators.
- Utilize the barycentric formulation rather than the naive product formulation to maintain 64-bit numerical stability and avoid catastrophic underflow.
- When evaluating errors, sample across interior non-nodal points; evaluating only at nodal abscissas trivially yields zero error.
- Ensure floating-point evaluation handles the case where the evaluation target x coincides exactly with a node x_j without generating a NaN.
- Maintain uniform test grid density across all refinement steps to prevent undersampling of localized high-frequency error peaks.
- Do not extrapolate outside the convex hull [x_0, x_n] of the interpolation nodes, as polynomial divergence accelerates dramatically outside the domain.
Theoretical limitations
- Equidistant high-degree polynomial interpolation suffers from Runge's phenomenon, restricting uniform mesh interpolation to low polynomial degrees.
- Finite precision IEEE 754 arithmetic introduces machine roundoff ε_mach ≈ 2.22 × 10⁻¹⁶, capping the minimum observable interpolation error.
- Global polynomial interpolation requires full domain recomputation if an additional data point is introduced, motivating Newton's divided difference formulation.