07 Piecewise Polynomial Interpolation and Splines

Recap: The Limits of Global Polynomials

In our previous lectures, we explored global polynomial interpolation in depth. We established two main approaches:

1. Barycentric Interpolation on Chebyshev Nodes: This is a robust and efficient method. Using the barycentric formula on Chebyshev nodes provides a numerically stable way to interpolate data that converges well for smooth functions. The cost is $O(n^2)$ for the initial weight calculation and then a fast $O(n)$ for each subsequent evaluation.

2. Chebyshev Basis Interpolation via FFT: This is an even faster alternative. By representing the polynomial in the Chebyshev basis, we can find all the coefficients in a remarkable $O(n\log n)$ time using the Fast Fourier Transform (FFT). Evaluation is then done efficiently with Clenshaw’s algorithm.

However, both of these powerful methods share a fundamental characteristic: they are global.

The Global Nature of Polynomial Interpolation

What does “global” mean in this context? If you change just one data point $(x_i,y_i)$, the entire interpolating polynomial changes everywhere.

• Lagrange/Barycentric: Changing one $y_i$ alters one term in the sum, affecting the value of $p(x)$ for all $x$.
• Fourier/Chebyshev: Changing one data point alters all the Fourier/Chebyshev coefficients (since each coefficient is an integral or sum over all data points), again changing the entire function.

This global dependency can be a significant drawback. A single measurement error can corrupt the entire model. Furthermore, high-degree polynomials can introduce unwanted oscillations (“wiggles”) between data points, even with well-chosen Chebyshev nodes, especially if the underlying function is not smooth. This is known as a lack of shape preservation.

To overcome these limitations, we shift our paradigm from using a single high-degree polynomial to stitching together multiple low-degree polynomials. This is the world of piecewise polynomial interpolation.

Piecewise Linear Interpolation: The Simplest Case

The most intuitive way to connect a set of data points is with straight lines. This is piecewise linear interpolation.

Given a set of nodes, or “breakpoints,” $a=x_0<x_1<\cdots <x_N=b$ and corresponding data values $y_j=f(x_j)$, we define a separate linear function $s_j(x)$ on each subinterval $[x_{j-1},x_j]$.

The Reference Interval Approach: A Powerful Construction Principle

While we can easily write the formula for the line segment on $[x_{j-1},x_j]$, it’s more instructive to use a powerful and generalizable technique: the reference interval method.

The "Divide and Conquer" of Function Construction

The core idea is to separate the creative work from the mechanical work.

1. Creative Step (Reference Interval): We solve the simplest possible version of our problem on a standard reference interval, which we’ll take as $t\in [0,1]$.
2. Mechanical Step (Mapping): We use a simple affine map to stretch and shift this “master solution” from the reference interval to any “real” or “current” interval $[x_{j-1},x_j]$.

Step 1: The Solution on the Reference Interval On $[0,1]$, we want a line connecting $(0,y_0)$ and $(1,y_1)$. We can construct this using two simple “hat” basis functions:

• $\chi (1-t)$: This is 1 at $t=0$ and 0 at $t=1$.
• $\chi (t)$: This is 0 at $t=0$ and 1 at $t=1$. The interpolant is then a simple linear combination:

$$f(t)=y_0\cdot \chi (1-t)+y_1\cdot \chi (t)$$

Step 2: Mapping to the Real Interval We map from the reference variable $t$ to the real-world variable $x$ using the linear transformation:

$$t=\frac{x-x_{j-1}}{h_j}⟺x=x_{j-1}+h_{j}t\text{where }{h}_j=x_j-x_{j-1}$$

Substituting this into our reference solution gives the formula for the piecewise linear interpolant $s_j(x)$ on the interval $[x_{j-1},x_j]$:

$$s_j(x)=y_{j-1}\left(\frac{x_j-x}{h_j}\right)+y_j\left(\frac{x-x_{j-1}}{h_j}\right)$$

This method of building a solution on a simple domain and then mapping it is a cornerstone of numerical methods, especially the Finite Element Method.

Properties of Piecewise Linear Interpolation

• Continuity: The resulting global function $s(x)$ is continuous ($C^0$) but is not differentiable at the breakpoints (it has “corners”).
• Locality: This is a crucial advantage. If we change a data value $y_j$, only the two adjacent segments, on $[x_{j-1},x_j]$ and $[x_j,x_{j+1}]$, are affected. The rest of the function remains unchanged. This makes the method robust to local errors.

Cubic Hermite Interpolation: Achieving Smoothness

Piecewise linear interpolation is simple and local, but its lack of smoothness is often a problem. We want a function that is not only continuous but also has a continuous first derivative ($C^1$). This means the slopes of the polynomial pieces must match up at the breakpoints.

To achieve this, we need more flexibility, so we move from degree 1 (linear) to degree 3 (cubic) polynomials on each segment. A cubic polynomial has four degrees of freedom ($c_0+c_{1}x+c_2{x}^2+c_3{x}^3$). On each interval $[x_{j-1},x_j]$, we have four pieces of information we can use to pin it down:

1. The value at the left endpoint: $s(x_{j-1})=y_{j-1}$
2. The value at the right endpoint: $s(x_j)=y_j$
3. The derivative at the left endpoint: $s^{\prime }(x_{j-1})=c_{j-1}$
4. The derivative at the right endpoint: $s^{\prime }(x_j)=c_j$

This is called Cubic Hermite Interpolation. By explicitly matching both the values and the derivatives at each breakpoint, we guarantee that the resulting piecewise function is $C^1$ smooth.

s

The Problem: Where Do the Derivatives Come From?

This construction is beautiful, but it relies on a critical assumption: that we know the desired derivative values $c_j=f^{\prime }(x_j)$ at each breakpoint. In most real-world scenarios, we only have the data values $y_j$. The derivatives are unknown.

The central challenge of creating smooth, piecewise interpolants is therefore: How do we choose or estimate the derivative values $c_j$ in a sensible way?

This choice is a degree of freedom that gives rise to a whole family of methods.

• Naive Approach: We could estimate the derivative using a finite difference, e.g., $c_j\approx (y_{j+1}-y_{j-1})/(x_{j+1}-x_{j-1})$. This often works but can lead to undesirable oscillations.
• Shape-Preserving Methods: More sophisticated methods choose the $c_j$ values to ensure the resulting curve preserves properties of the data, like monotonicity or convexity. If the data is monotonic, the interpolant should also be monotonic.

Shape-Preserving Splines: Akima, PCHIP

Global polynomial interpolation often fails at shape preservation, introducing “wiggles” where none should exist. This is where specialized cubic Hermite methods shine.

As seen on, for a function that should be monotonic, the global barycentric interpolant (blue) exhibits significant over- and undershoots. In contrast, methods designed for shape preservation produce a much more plausible curve.

• Akima & Makima: Popular methods that estimate the derivatives $c_j$ based on a local, weighted average of secant slopes. They are known for producing visually pleasing curves.
• PCHIP (Piecewise Cubic Hermite Interpolating Polynomial): This is the method implemented in MATLAB’s pchip. It guarantees that if the data is monotonic, the interpolant will also be monotonic. It achieves this by carefully constraining the derivative estimates. PCHIP produces a shape-preserving curve, avoiding the oscillations of other methods.

Cubic Splines: Achieving Even Greater Smoothness

Cubic Hermite methods give us a $C^1$ function. What if we need even more smoothness? A cubic spline is a piecewise cubic polynomial 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 degree of smoothness, we must impose an additional constraint. Instead of choosing the first derivatives $c_j$ freely, we enforce the condition that the second derivatives match at each interior breakpoint:

$$s_{j-1}^{\prime \prime }(x_j)=s_j^{\prime \prime }(x_j)\text{for }j=1,\dots ,N-1$$

This condition links all the segments together. It turns the local problem of choosing derivatives into a global problem. Applying this condition at every interior node (using the chain rule and a lot of algebra) results in a tridiagonal system of linear equations for the unknown derivative values $c_j$.

This system is diagonally dominant and can be solved very efficiently in $O(N)$ time.

Boundary Conditions: Tying Up the Loose Ends

The system of equations for the second derivative continuity gives us $N-1$ equations for the $N+1$ unknown derivatives $c_0,\dots ,c_N$. We are two equations short. We need to specify two additional boundary conditions to get a unique solution. Common choices include:

1. “Clamped” Spline: We specify the true derivatives at the endpoints: $s^{\prime }(x_0)=f^{\prime }(x_0)$ and $s^{\prime }(x_N)=f^{\prime }(x_N)$. This is the most accurate if the endpoint derivatives are known.
2. “Natural” Spline: We assume the second derivatives at the endpoints are zero: $s^{\prime \prime }(x_0)=0$ and $s^{\prime \prime }(x_N)=0$. This corresponds to the physical behavior of a flexible ruler being laid across the points, which would be straight (zero curvature) beyond the last points.
3. “Not-a-Knot” Spline: This condition forces the third derivative to be continuous at the first and last interior knots ($x_1$ and $x_{N-1}$). This effectively means that the first two polynomial pieces are the same cubic, and the last two are the same cubic. It’s a popular default choice when no other information is available.
4. Periodic Spline: If the function is known to be periodic, we enforce $s^{\prime }(x_0)=s^{\prime }(x_N)$ and $s^{\prime \prime }(x_0)=s^{\prime \prime }(x_N)$.

The Optimality of Natural Splines

Why are splines so important? The “natural” cubic spline has a remarkable physical and mathematical property.

Theorem: Minimum Curvature Property

Among all twice-differentiable functions $g(x)$ that interpolate the given data points, the natural cubic spline is the one that minimizes the total squared curvature.

$${\int }_a^b|s^{\prime \prime }(x){|}^{2}dx\le {\int }_a^b|g^{\prime \prime }(x){|}^{2}dx$$

Physically, this means the natural spline is the shape a flexible, elastic ruler (a draftsman’s spline) would take if pinned at the data points. It is the “smoothest” possible interpolant in the sense of minimizing bending energy.

The convergence plots show that these various spline methods provide excellent accuracy, especially for smooth functions, demonstrating high orders of convergence ($O(h^4)$ for the function value). They successfully combine the local flexibility of piecewise polynomials with the high accuracy that comes from enforcing global smoothness conditions.

Continue here: 08 A Grand Review of Interpolation Methods and Splines