Chebyshev Polynomial Interpolation

• Polynomial (Chebyshev) interpolation is the construction of an algebraic polynomial in the Chebyshev basis that approximates functions at optimally chosen nodes to achieve near-minimax error control.
• It employs stable barycentric formulations and extends to multivariate cases using tensor grids and non-tensorial node sets such as Padua and Lissajous points.
• Advanced techniques like filtering, adaptive partitioning, and mock-Chebyshev methods further enhance approximation quality and mitigate issues like the Runge phenomenon.

Polynomial (Chebyshev) interpolation refers to the construction of an algebraic polynomial, typically expressed in a Chebyshev polynomial basis, that interpolates a function at a prescribed set of nodes. When the nodes are chosen as the zeros or extrema of Chebyshev polynomials, or their multispace generalizations such as Lissajous–Chebyshev sets, the resulting interpolation schemes possess optimal stability, minimax error properties, and efficient numerical implementations. The theory covers the classical univariate context, tensorized multivariate grids, non-tensorial nodesets (Padua, Lissajous), filtered variants for improved uniform approximation, and adaptive domain decompositions. It also underpins sparse interpolation, error-resilient algorithms, and applications in scientific computing and computational finance.

1. Mathematical Foundations and Node Systems

Chebyshev polynomials of the first kind, defined by $T_n(x) = \cos(n \arccos x)$ on $[-1,1]$, form an orthogonal system with respect to the weight $w(x) = (1-x^2)^{-1/2}$. Their zeros and extrema,

$$x^{(z)}_k = \cos \Big( \frac{(2k+1)\pi}{2n+2} \Big),\quad x^{(e)}_k = \cos \Big( \frac{k\pi}{n} \Big),\quad k=0,\ldots,n,$$

are used as interpolation nodes due to their optimal distribution: the associated Lebesgue constant grows only logarithmically in $n$ and the nodes minimize the maximum of the interpolation error kernel, achieving near-minimax properties (Amartey, 2024).

Multivariate extensions adopt tensor products (Chebyshev grids) or more sophisticated node sets such as Padua points, degenerate Lissajous curves, and general Lissajous–Chebyshev lattices, parameterized via number-theoretic properties (e.g., coprimality conditions) and parity constraints on multi-indices. Lissajous–Chebyshev nodes uniformly subsume classical Padua and Morrow-Patterson-Xu points and allow both dense (tensor-product) and sparse, highly efficient interpolatory schemes (Dencker et al., 2015, Dencker et al., 2017, Erb, 2015).

2. Interpolation Schemes and Stability

2.1 Classical Lagrange–Chebyshev Interpolation

The Lagrange interpolant using Chebyshev nodes—either nodes (roots) or extrema (Gauss-Lobatto)—is given by: $$P_n[f](x) = \sum_{k=0}^n f(x_k) \ell_k(x),\qquad \ell_k(x) = \prod_{j\ne k} \frac{x-x_j}{x_k-x_j}.$$ Its barycentric formulation is numerically stable and allows $O(n)$ evaluation per point (Diehl et al., 2021, Amartey, 2024, Breuß et al., 2016).

For any smooth $f$, the interpolation error at Chebyshev nodes is given by

$$|f(x) - P_n[f](x)| = \frac{|f^{(n+1)}(\xi)|}{(n+1)!} \prod_{k=0}^n |x - x_k|,\quad \xi \in [-1,1].$$

Use of Chebyshev nodes minimizes $$\max_{x\in[-1,1]} \left| \prod_{k=0}^n (x-x_k) \right|$$, thereby almost minimizing the interpolation error for any given degree.

2.2 Lebesgue Constants and Minimax Properties

The Lebesgue constant,

$[-1,1]$0

governs the uniform stability of interpolation. Chebyshev nodes yield $[-1,1]$1. For multivariate tensor grids, $[-1,1]$2 (Dencker et al., 2016, Bialas-Ciez et al., 2023). Lissajous–Chebyshev nodes achieve the same order, $[-1,1]$3, as full tensor grids, even with "half" the node count (Dencker et al., 2016, Dencker et al., 2017).

Non-tensorial node sets (Padua, rank-1 Lissajous) have Lebesgue constants growing polylogarithmically, e.g., Padua points with $[-1,1]$4 in two dimensions (Erb, 2015).

3. Error Bounds and Convergence Analysis

For analytic $[-1,1]$5 extending to a Bernstein ellipse $[-1,1]$6, Chebyshev interpolation achieves geometric convergence: $[-1,1]$7 In the multivariate tensor setting, the error bound improves to (Glau et al., 2016, Gaß et al., 2015): $[-1,1]$8 If all $[-1,1]$9 are equal and $w(x) = (1-x^2)^{-1/2}$0, the decay remains essentially spectral, while alternative bounds (e.g., anisotropic orderings) can sharpen error guarantees (Glau et al., 2016).

If $w(x) = (1-x^2)^{-1/2}$1 has bounded variation, the error decays algebraically at the optimal rate $w(x) = (1-x^2)^{-1/2}$2 (Xiang, 2014).

4. Advanced Schemes: Filtering, Partitioning, and Least-Squares

4.1 Filtered (de la Vallée–Poussin) Interpolation

Filtered polynomial interpolation via de la Vallée–Poussin means employs a low-pass filter $w(x) = (1-x^2)^{-1/2}$3 applied to the Chebyshev expansion coefficients: $w(x) = (1-x^2)^{-1/2}$4 where $w(x) = (1-x^2)^{-1/2}$5 are the (discrete) Chebyshev coefficients (Occorsio et al., 2020, Occorsio et al., 2021). This scheme preserves values at the nodes, avoids the logarithmic ("Lebesgue penalty") blow-up, and achieves uniform operator-norm boundedness under necessary and sufficient conditions on the weighted norm exponents.

Filtered interpolants admit near-best uniform-norm approximation error, with rates governed by the smoothness of $w(x) = (1-x^2)^{-1/2}$6 and the width $w(x) = (1-x^2)^{-1/2}$7, and suppress Gibbs-type phenomena for non-smooth data. For valid (Chebyshev or Jacobi) weights, the Lebesgue constant of the filtered operator remains uniformly bounded in $w(x) = (1-x^2)^{-1/2}$8, unlike the diverging constants of classical Lagrange interpolation (Occorsio et al., 2020, Occorsio et al., 2021).

4.2 Partition of Unity and Adaptive Splitting

When $w(x) = (1-x^2)^{-1/2}$9 has singularities near the interval or localized features, adaptive domain partitioning accelerates convergence. The partition-of-unity approach covers the domain with overlapping subintervals, constructs high-degree Chebyshev interpolants locally, and blends them via $$x^{(z)}_k = \cos \Big( \frac{(2k+1)\pi}{2n+2} \Big),\quad x^{(e)}_k = \cos \Big( \frac{k\pi}{n} \Big),\quad k=0,\ldots,n,$$0 weights: $$x^{(z)}_k = \cos \Big( \frac{(2k+1)\pi}{2n+2} \Big),\quad x^{(e)}_k = \cos \Big( \frac{k\pi}{n} \Big),\quad k=0,\ldots,n,$$1 The global error is no worse than the worst local error, and by judicious splitting, exponential convergence is restored even when the global analytic domain is narrow (Aiton et al., 2017).

4.3 Mock-Chebyshev Interpolation

When only uniform samples are available, mock-Chebyshev strategies interpolate on subsets of equispaced nodes that mimic Chebyshev–Lobatto points, followed by least-squares correction on the residual. The approach yields substantial suppression of the Runge phenomenon, achieving error rates of $$x^{(z)}_k = \cos \Big( \frac{(2k+1)\pi}{2n+2} \Big),\quad x^{(e)}_k = \cos \Big( \frac{k\pi}{n} \Big),\quad k=0,\ldots,n,$$2 and uniform-norm improvement under analyticity (Marchi et al., 2014).

5. Multivariate, Sparse, and Error-Resilient Interpolation

Multivariate Chebyshev interpolation extends via tensor products or non-tensor Lissajous–Chebyshev sets (Dencker et al., 2015, Dencker et al., 2017). For high-dimensional or sparse regimes, Prony-type sparse interpolation enables the exact recovery of Chebyshev expansions with at most $$x^{(z)}_k = \cos \Big( \frac{(2k+1)\pi}{2n+2} \Big),\quad x^{(e)}_k = \cos \Big( \frac{k\pi}{n} \Big),\quad k=0,\ldots,n,$$3 terms,