Interpolation and Approximation

Interpolation constructs a function passing through given data points. Approximation finds a function that is close to the data (not necessarily passing through it). Both are essential for data fitting, function evaluation, and numerical methods.

Polynomial Interpolation

Given n + 1 points (x₀, y₀), ..., (xₙ, yₙ) with distinct xᵢ, there exists a unique polynomial of degree ≤ n passing through all points.

Lagrange Interpolation

P(x) = Σᵢ₌₀ⁿ yᵢ · Lᵢ(x)

where the Lagrange basis polynomials are:

Lᵢ(x) = Πⱼ≠ᵢ (x - xⱼ) / (xᵢ - xⱼ)

Properties: Lᵢ(xⱼ) = δᵢⱼ (1 if i=j, 0 otherwise).

Advantages: Conceptually simple, explicit formula. Disadvantages: O(n²) to evaluate. Adding a point requires recomputation. Numerically unstable for large n.

Newton Interpolation

P(x) = c₀ + c₁(x-x₀) + c₂(x-x₀)(x-x₁) + ... + cₙ(x-x₀)...(x-xₙ₋₁)

where cₖ = f[x₀, x₁, ..., xₖ] are divided differences.

Divided differences (recursion):

f[xᵢ] = yᵢ
f[xᵢ, xᵢ₊₁, ..., xᵢ₊ₖ] = (f[xᵢ₊₁,...,xᵢ₊ₖ] - f[xᵢ,...,xᵢ₊ₖ₋₁]) / (xᵢ₊ₖ - xᵢ)

Advantages: O(n) to add a new point. O(n) to evaluate (Horner-like nesting). Disadvantages: Still ill-conditioned for large n with bad nodes.

Neville's Algorithm

Recursively build interpolating polynomials of increasing degree. Useful for computing P(x) at a single point without finding the polynomial explicitly.

Interpolation Error

The error of interpolating f(x) with polynomial P(x) through n+1 points:

f(x) - P(x) = f^(n+1)(ξ) / (n+1)! · Πᵢ₌₀ⁿ (x - xᵢ)

for some ξ in the interval containing x and all xᵢ.

The error depends on:

1. The (n+1)-th derivative of f (smoothness)
2. The product Πᵢ(x - xᵢ) (node placement)

Runge's Phenomenon

Problem: High-degree polynomial interpolation on equally spaced points can oscillate wildly near the endpoints.

Classic example: f(x) = 1/(1 + 25x²) on [-1, 1] with equidistant nodes. As n increases, the interpolant diverges near ±1.

Solution: Use Chebyshev nodes instead of equidistant points.

Chebyshev Nodes

The Chebyshev nodes on [-1, 1]:

xₖ = cos((2k+1)π / (2(n+1))),  k = 0, 1, ..., n

These are the roots of the Chebyshev polynomial Tₙ₊₁(x).

Why they work: Chebyshev nodes minimize max|Πᵢ(x - xᵢ)| over [-1, 1]. They cluster near the endpoints, counteracting the Runge phenomenon.

Near-optimal: Among all possible node sets, Chebyshev nodes minimize the maximum interpolation error (within a factor of ~2 of optimal).

For a general interval [a, b]: map Chebyshev nodes via x = (a+b)/2 + (b-a)/2 · cos(...).

Splines

Instead of one high-degree polynomial, use piecewise polynomials joined smoothly at the data points (knots).

Linear Splines

Connect adjacent points with straight lines. Continuous but not smooth (corners at knots).

Cubic Splines

Piecewise cubic polynomials S(x) satisfying:

1. S(xᵢ) = yᵢ at each data point
2. S, S', S'' are continuous at interior knots
3. Boundary conditions (natural, clamped, not-a-knot, or periodic)

On each interval [xᵢ, xᵢ₊₁]:

Sᵢ(x) = aᵢ + bᵢ(x-xᵢ) + cᵢ(x-xᵢ)² + dᵢ(x-xᵢ)³

4n unknowns (4 coefficients per interval). Conditions give 4n equations → unique solution.

Natural spline: S''(x₀) = S''(xₙ) = 0. The spline minimizes ∫(S'')² dx (minimizes total curvature — smoothest interpolant).

Advantages: Avoids Runge's phenomenon. Smooth (C²). Stable. Local control (changing one data point affects nearby intervals only).

B-Splines

Basis splines (B-splines) are a basis for the space of splines. Any spline can be written as:

S(x) = Σᵢ cᵢ Bᵢ,ₖ(x)

where Bᵢ,ₖ are B-spline basis functions of order k.

Properties: Non-negative, compact support, partition of unity, numerically stable.

De Boor's algorithm: Evaluates B-spline curves efficiently (analog of de Casteljau for Bézier curves).

NURBS (Non-Uniform Rational B-Splines)

Rational B-splines with non-uniform knot vectors. Can exactly represent conic sections (circles, ellipses).

Standard in CAD/CAM (computer-aided design). Used in OpenGL, automotive design, font rendering.

Piecewise Polynomial Interpolation

Hermite Interpolation

Matches both function values and derivatives at each node. Piecewise cubic Hermite (two values + two derivatives per interval = 4 conditions per piece).

PCHIP (Piecewise Cubic Hermite Interpolating Polynomial): Chooses slopes to preserve monotonicity. Avoids overshoot. Used in MATLAB's pchip.

Trigonometric Interpolation

For periodic functions, use sines and cosines instead of polynomials.

Given 2n+1 equally spaced points on [0, 2π):

P(x) = a₀/2 + Σₖ₌₁ⁿ (aₖcos(kx) + bₖsin(kx))

This is the discrete Fourier transform (DFT) and is computed by the FFT in O(n log n).

Padé Approximation

Approximate f(x) as a ratio of polynomials:

f(x) ≈ P_m(x) / Q_n(x)

where P_m has degree m and Q_n has degree n. Matches the first m+n+1 Taylor coefficients.

Advantages: Much better than Taylor series for functions with poles or singularities. Can represent asymptotic behavior.

Example: Padé[1,1] for eˣ: (1 + x/2)/(1 - x/2). This is the basis of the Cayley transform and the Crank-Nicolson method.

Least Squares Approximation

When data has noise, interpolation overfits. Instead, fit a lower-degree polynomial (or other basis) minimizing the sum of squared residuals.

Discrete Least Squares

Fit f(x) = Σⱼ cⱼ φⱼ(x) to data points by minimizing Σᵢ(yᵢ - Σⱼ cⱼ φⱼ(xᵢ))².

Normal equations: (ΦᵀΦ)c = Φᵀy, or use QR decomposition for stability.

Continuous Least Squares

Minimize ∫(f(x) - p(x))² w(x) dx. Leads to orthogonal polynomials:

• Legendre: w(x) = 1 on [-1, 1]
• Chebyshev: w(x) = 1/√(1-x²) on [-1, 1]
• Laguerre: w(x) = e⁻ˣ on [0, ∞)
• Hermite: w(x) = e⁻ˣ² on (-∞, ∞)

Minimax Approximation

Minimize the maximum error instead of the sum of squares:

min_p max_x |f(x) - p(x)|

The Remez algorithm (exchange algorithm) computes the best minimax polynomial. The error equioscillates (Chebyshev equioscillation theorem).

Used in libm implementations: hardware math library functions (sin, cos, exp, log) use minimax polynomials tuned to machine precision.

Applications in CS

• Font rendering: Bézier curves and B-splines define glyph outlines (TrueType, OpenType).
• Computer graphics: Spline curves for animation paths, surface patches for modeling.
• Signal processing: Interpolation for sample rate conversion, image resizing (bilinear, bicubic, Lanczos).
• Scientific computing: Interpolation tables for expensive functions. Spectral methods use Chebyshev interpolation.
• Game development: Spline paths for cameras and objects. Terrain heightmap interpolation.
• Machine learning: Basis function expansion (radial basis functions, polynomial features).
• Math libraries: libm functions (sin, exp, log) use polynomial/rational approximations. Minimax or Chebyshev expansions tuned to machine precision.
• Numerical integration: Quadrature rules are derived from polynomial interpolation.
• Data visualization: Smooth curves through scatter plots.