08 A Grand Review of Interpolation Methods and Splines

A Look Back and a Look Ahead

We have now concluded the most challenging part of this course: interpolation. The difficulty lies in the fact that our unknown is not a single number, but an entire function, a complex, infinite-dimensional object that we are trying to approximate from a finite amount of data.

The upcoming topics, such as numerical integration (quadrature), will be conceptually simpler because the unknown will be just a single number. We will find that many of the powerful ideas we developed for interpolation, like choosing optimal nodes and constructing basis functions, will be directly repurposed for quadrature, making our future work much easier.

Today, we will synthesize everything we’ve learned about interpolation, comparing the various methods to understand their trade-offs. We will then formally introduce the powerful alternative of piecewise polynomial interpolation and splines, which address the fundamental limitations of global methods.

The Landscape of Interpolation Methods: A Showdown

We’ve seen a variety of methods, and it’s crucial to understand that there is no single “best” one. The right choice depends entirely on the problem at hand: what do we know about our function, what are our constraints, and what is our goal?

Let’s compare the main contenders.

Global Polynomial Interpolation

These methods use a single polynomial of degree $n$ to fit $n+1$ data points.

1. The Monomial Basis (Naive Approach)

   • Method: Solve the Vandermonde system for coefficients of $\{1,x,x^2,\dots \}$.
   • Verdict: A numerical disaster. While mathematically sound, the Vandermonde matrix is horribly ill-conditioned. This method should never be used in practice.

2. Barycentric Interpolation on Equidistant Nodes

   • Method: Use the stable barycentric formula on equally spaced points.
   • Verdict: A dangerous trap. While the formula is stable, the choice of equidistant nodes for high-degree polynomials leads to Runge’s phenomenon, wild oscillations near the interval boundaries. The error diverges as the number of points increases.

3. Barycentric Interpolation on Chebyshev Nodes

   • Method: Use the stable barycentric formula, but sample the function at the optimally chosen Chebyshev nodes.
   • Verdict: An excellent, robust choice. This method avoids Runge’s phenomenon and converges rapidly (exponentially) for smooth functions. It is a go-to method when you have the freedom to choose your sample points.

4. Chebyshev Basis Interpolation (via FFT)

   • Method: Represent the polynomial in the Chebyshev basis $\{T_0(x),T_1(x),\dots \}$ and find the coefficients using the FFT.
   • Verdict: The fastest method for smooth functions. This approach leverages the power of the FFT to find the polynomial coefficients in $O(n\log n)$ time. It shares the excellent convergence properties of using Chebyshev nodes.

The Magic of Chebyshev Interpolation

Why does Chebyshev interpolation work so well for non-periodic functions, even though it uses a cosine-based (trigonometric) construction? The key is the variable change $x=\cos (\theta )$. This mapping “tricks” the non-periodic function on $[-1,1]$ into looking like a smooth, even (symmetric) periodic function in the $\theta$ domain. This even symmetry means it can be represented efficiently by a cosine series, which is exactly what a Chebyshev series is! This elegant trick avoids the boundary discontinuities that plague standard trigonometric interpolation for non-periodic functions.

5. Trigonometric (Fourier) Interpolation
   • Method: Approximate the function using a basis of sines and cosines (complex exponentials) on equidistant points.
   • Verdict: Perfect for periodic functions, terrible for non-periodic ones. If the function is naturally periodic and smooth, this method converges exponentially fast. However, if the function is not periodic, its “periodic extension” will have a jump discontinuity at the boundary. This jump introduces high frequencies and destroys the fast convergence, leading to the Gibbs phenomenon and slow, algebraic convergence.

The Big Picture: Global Methods

All these methods, Monomial, Barycentric, Chebyshev, Fourier, are global. Changing a single data point affects the entire interpolating function. This can be a disadvantage if your data is noisy or if the function has localized, non-smooth features.

Piecewise Polynomial Interpolation: The “Divide and Conquer” Philosophy

To overcome the limitations of global methods, we can adopt a “divide and conquer” strategy. Instead of trying to fit one high-degree polynomial to the whole interval $[a,b]$, we break the interval into smaller pieces and fit a low-degree polynomial to each piece.

The Motivation:

1. Locality: Changes to data in one subinterval only affect the interpolant locally. This makes the method robust to noise and measurement errors.
2. Handling Non-Smoothness: If a function has a “kink” or discontinuity in a derivative, we can place a breakpoint at that exact location. This isolates the “bad behavior” to a single point, allowing the interpolant to be highly accurate in the smooth regions on either side.
3. Guaranteed Convergence: The error on each small piece of length $h$ is bounded by a term proportional to $h^{n+1}$, where $n$ is the (low) degree of the polynomial. As we add more points (making $h$ smaller), the error is guaranteed to decrease, avoiding the divergence of Runge’s phenomenon.

Piecewise Linear Interpolation

This is the simplest approach: connect the dots with straight lines.

• Construction: We use “hat functions” (or ReLU functions in machine learning) as a local basis. The interpolant is a sum of these hat functions, weighted by the data values $y_j$.
• Properties: The resulting function is continuous ($C^0$) but has “corners” at the breakpoints, so it is not differentiable. The approximation order is low.

Cubic Hermite Interpolation

To get a smoother curve, we use piecewise cubic polynomials. On each interval, a cubic has four degrees of freedom. We use these to match not only the function values ($y_{j-1},y_j$) but also the first derivatives ($c_{j-1},c_j$) at the endpoints.

• Properties: This construction guarantees a $C^1$ smooth function (continuous function and continuous first derivative).
• The Challenge: We usually don’t know the true derivatives $c_j$. The core problem becomes: how do we choose or estimate these derivative values?

A Family of Methods

The freedom to choose the derivatives $c_j$ gives rise to a whole family of methods, each with different goals:

• High Accuracy: If we could supply the exact derivatives $f^{\prime }(x_j)$, the method would have a very high convergence order ($O(h^4)$).
• Shape Preservation: Methods like PCHIP (Piecewise Cubic Hermite Interpolating Polynomial) choose the $c_j$ values in a way that guarantees the interpolant preserves the monotonicity of the data. If the data is increasing, the curve will not have any “wiggles.” This comes at the cost of a lower convergence order (as seen on Slide 8, the order drops from 4 to ~2.5).

Splines: The Pursuit of Higher Smoothness

What if we need even more smoothness? A spline is a piecewise polynomial that satisfies global smoothness conditions. The most common type is the cubic spline.

Definition: Cubic Spline

A cubic spline is a piecewise cubic function that is twice continuously differentiable ($C^2$). This means the function, its first derivative, and its second derivative are all continuous across the breakpoints.

To achieve this extra level of smoothness, we enforce the condition that the second derivatives of the cubic pieces must match at each interior breakpoint:

$$s_{j-1}^{\prime \prime }(x_j)=s_j^{\prime \prime }(x_j)$$

This constraint links all the pieces together. Instead of choosing the first derivatives $c_j$ freely, they now become the unknowns in a global system of linear equations. This system is tridiagonal, meaning it can be solved very efficiently in $O(N)$ time.

Boundary Conditions

This system gives us $N-1$ equations for the $N+1$ unknown derivatives $c_0,\dots ,c_N$. We need two more conditions to get a unique solution. These are the boundary conditions:

1. Clamped Spline: Specify the derivatives at the endpoints, $c_0$ and $c_N$. (Requires extra information).
2. Natural Spline: Assume the second derivatives at the endpoints are zero, $s^{\prime \prime }(x_0)=s^{\prime \prime }(x_N)=0$.
3. “Not-a-Knot” Spline: A common default that forces the first two pieces (and last two pieces) to be the same cubic polynomial.
4. Periodic Spline: For periodic data, enforce $c_0=c_N$ and $s^{\prime \prime }(x_0)=s^{\prime \prime }(x_N)$.

The Optimality of the Natural Spline

The natural spline is special because it is the “smoothest” possible interpolant.

The Draftsman's Spline

The term “spline” comes from the flexible wooden or metal strips used by draftsmen to draw smooth curves. They would pin the strip at several points, and the natural elastic forces would cause it to bend in a shape that minimizes its total bending energy. The natural cubic spline is the mathematical formalization of this physical object. It is the interpolating function that minimizes the total squared curvature $\int |s^{\prime \prime }(x){|}^{2}dx$.

Final Showdown: Which Method to Use?

Let’s run a final comparison on a few test functions to see these trade-offs in action.

• For the smooth, periodic sin(pi*x) function: Chebyshev interpolation is the clear winner, achieving exponential convergence. Trigonometric interpolation also performs exceptionally well, as expected. Spline methods are very good but show algebraic convergence ($O(h^4)$), which is slower than exponential.
• For the non-smooth |sin(pi*x)|^3 function (which has a $C^2$ “kink”): The performance of the global methods suffers due to the lack of smoothness. Piecewise methods like splines, which can handle such local features gracefully, become much more competitive.
• For the non-periodic Runge function:
  • Barycentric on equidistant nodes fails catastrophically.
  • Trigonometric interpolation is hopelessly bad. Why? Because the function is not periodic. Its periodic extension has a large jump discontinuity at the boundary, which completely destroys the convergence.
  • Chebyshev interpolation is the clear winner. Its variable-change trick creates an effectively smooth, periodic function, leading to spectral convergence.
  • Spline methods (clamped, natural, not-a-knot) all perform very well and are nearly indistinguishable from the Chebyshev solution on a coarse grid. They provide a robust, high-quality approximation without needing special node points.

Conclusion:

• If your function is smooth and you can choose the sample points, use Chebyshev interpolation (either barycentric or FFT-based) for the best performance.
• If your function is smooth and periodic, use trigonometric (Fourier) interpolation.
• If your sample points are fixed (and possibly equidistant), or if your function has kinks, corners, or other non-smooth features, cubic splines are an excellent and robust choice. They provide local control, guaranteed convergence, and high-order accuracy