Chapter 4 - Piecewise Polynomial Interpolation

In the last two chapters, we explored the power and peril of using a single, high-degree polynomial to fit a set of points. We saw that with a clever choice of nodes, the Chebyshev points, we could achieve remarkable accuracy. However, this approach has two major weaknesses. First, we are not always free to choose our data points; they are often given to us, perhaps from an experiment at evenly spaced time intervals. Second, a single global polynomial is a rigid tool. A single sharp corner or discontinuity anywhere in the function can send ruinous oscillations across the entire domain.

This calls for a more flexible, robust, and practical philosophy: divide and conquer.

Instead of trying to fit one complex, high-degree curve to all the data, we will fit many simple, low-degree curves to small sections of it. We will break our interval $[a,b]$ into smaller subintervals and use a low-degree polynomial on each piece. This is the world of piecewise polynomial interpolation.

The error formula from Chapter 2 provides the mathematical motivation. For a polynomial of degree $n$ on an interval of length $h$, the error is roughly:

$$|f(x)-s(x)|\le C\frac{h^{n+1}}{(n+1)!}\Vert f^{(n+1)}{\Vert }_{\infty }$$

This formula tells us we can control the error in two ways: by increasing the polynomial degree $n$ or by decreasing the interval size $h$. The global polynomial approach focuses on increasing $n$, which we saw can be dangerous. The piecewise approach keeps $n$ small (typically 1 or 3) and focuses on making $h$ small. This strategy is far more stable and adaptable.

Piecewise Linear Interpolation: Connecting the Dots

The simplest idea is to just connect the dots with straight lines. This creates a function that is continuous, but not smooth.

Given a set of $N+1$ data points $(x_j,y_j)$, which we now call knots or breakpoints, we define a separate linear function on each subinterval $[x_{j-1},x_j]$. Let $h_j=x_j-x_{j-1}$ be the length of the subinterval. The formula for the line segment is a simple weighted average:

$$s_j(x)=y_{j-1}\frac{x_j-x}{h_j}+y_j\frac{x-x_{j-1}}{h_j},\text{for }x\in [x_{j-1},x_j]$$

This can be understood more formally by defining a linear “hat function” on a reference interval $[0,1]$, $\chi (t)=t$. Then, by mapping each subinterval $[x_{j-1},x_j]$ to $[0,1]$ via $t=(x-x_{j-1})/h_j$, our interpolant becomes:

$$s_j(x)=y_{j-1}(1-\chi (t))+y_j\chi (t)$$

The resulting global function, $s(x)$, is continuous (we say $s\in C^0[a,b]$), but its derivative is discontinuous at the knots, the function has sharp “kinks.” This lack of smoothness is its primary drawback. It’s simple and robust, but it’s not elegant. To do better, we need to control the derivatives.

Cubic Hermite Interpolation: Building Smoothness

A function appears smooth if its slope doesn’t change abruptly. To achieve this, we need to build an interpolant whose first derivative is also continuous, making it a $C^1$ function.

The idea of Hermite interpolation is to match not only the function values at the knots, but also the derivative values. Let’s assume that for each knot $x_j$, we are given both a value $y_j=s(x_j)$ and a slope $c_j=s^{\prime }(x_j)$.

On a single interval $[x_{j-1},x_j]$, we now have four pieces of information: $y_{j-1},c_{j-1},y_j,c_j$. This is exactly the right amount of information to uniquely determine a cubic polynomial (which has four coefficients, $a{x}^3+b{x}^2+cx+d$).

How do we construct this cubic? We can define a set of basis functions on the reference interval $[0,1]$. We need four functions that handle our four constraints:

• Two for the values: ${\varphi }_0(t)$ and ${\varphi }_1(t)$
• Two for the slopes: ${\psi }_0(t)$ and ${\psi }_1(t)$

These basis polynomials are chosen to isolate each condition:

• ${\varphi }_0(t)$: Must be 1 at $t=0$, 0 at $t=1$, and have zero slope at both ends. The unique cubic is ${\varphi }_0(t)=2t^3-3t^2+1$.
• ${\varphi }_1(t)$: Must be 0 at $t=0$, 1 at $t=1$, and have zero slope at both ends. The unique cubic is ${\varphi }_1(t)=-2t^3+3t^2$.
• ${\psi }_0(t)$: Must be 0 at both ends, have slope 1 at $t=0$, and slope 0 at $t=1$. The unique cubic is ${\psi }_0(t)=t^3-2t^2+t$.
• ${\psi }_1(t)$: Must be 0 at both ends, have slope 0 at $t=0$, and slope 1 at $t=1$. The unique cubic is ${\psi }_1(t)=t^3-t^2$.

With these building blocks, the cubic interpolant on the interval $[x_{j-1},x_j]$ is a combination that satisfies all four conditions:

$$s_j(x)=y_{j-1}{\varphi }_0(t)+y_j{\varphi }_1(t)+h_j{c}_{j-1}{\psi }_0(t)+h_j{c}_j{\psi }_1(t)\text{where }t=\frac{x-x_{j-1}}{h_j}$$

The factors of $h_j$ are crucial; they come from the chain rule when converting derivatives with respect to $t$ to derivatives with respect to $x$.

By construction, the resulting piecewise cubic function $s(x)$ will be $C^1$ smooth. But this raises a crucial question: where do the slope values $c_j$ come from?

Case 1: Slopes are Known (True Hermite Interpolation)

If the derivatives of the underlying function are known, we can use them directly. In this ideal case, the error is very small.

Theorem 4.2.2: Error of Cubic Hermite Interpolation If $f\in C^4[a,b]$ and $s(x)$ is the piecewise cubic Hermite interpolant using the exact derivatives $c_j=f^{\prime }(x_j)$, then the error is bounded by:

$$\Vert f^{(k)}-s^{(k)}{\Vert }_{L^{\infty }}\le C_k{h}^{4-k}\Vert f^{(4)}{\Vert }_{L^{\infty }}\text{for }k=0,1$$

where $h={\max }_j{h}_j$. The function value converges as $O(h^4)$ and the derivative as $O(h^3)$. This is excellent convergence.

Case 2: Slopes are Estimated (PCHIP)

In most real-world problems, we only have the data values $y_j$. We must estimate the slopes. Different estimation schemes lead to different types of cubic interpolants. One of the most useful is the Piecewise Cubic Hermite Interpolating Polynomial (PCHIP).

PCHIP estimates the slope at a knot $x_j$ by looking at the slopes of the lines connecting its neighbors. It uses a carefully chosen weighted average designed to prevent overshoots and preserve the shape of the data.

The Power of Shape Preservation

The key feature of PCHIP is that it is monotonicity-preserving. If your data points are monotonically increasing, the PCHIP interpolant will also be monotonically increasing. It will not introduce new “bumps” or wiggles where none should exist. This is an invaluable property when interpolating physical data that is known to be non-oscillatory.

The price for this robustness is a lower order of accuracy compared to using the true derivatives, but its shape-preserving nature often makes it the right tool for the job.

from scipy.interpolate import pchip_interpolate
# PchipInterpolator is also available for an object-oriented approach
 
# Data that is monotonic
x_data = np.array([0, 1, 2, 4, 5, 6])
y_data = np.array([0, 0, 1, 1, 2, 2])
 
x_plot = np.linspace(0, 6, 400)
y_pchip = pchip_interpolate(x_data, y_data, x_plot)
 
plt.plot(x_data, y_data, 'ro', label='Data')
plt.plot(x_plot, y_pchip, 'b-', label='PCHIP Interpolant')
plt.title("PCHIP Preserves Monotonicity")
plt.legend(); plt.grid(True); plt.show()

Cubic Splines: The Gold Standard of Smoothness

PCHIP gives us a $C^1$ function, which has a continuous slope. But if you look closely, the curvature (the second derivative) can still jump at the knots. For applications in computer graphics, engineering design, or simply creating a visually perfect curve, we often want an even higher degree of smoothness. We want a function that is twice continuously differentiable, or $C^2$.

This is the domain of the cubic spline.

The idea is beautiful. We start, as before, with the Hermite cubic form on each interval, but now the slopes $c_j=s^{\prime }(x_j)$ are unknown variables. We will determine them all at once by enforcing a global smoothness condition.

At each of the $N-1$ interior knots, we demand that the second derivative of the cubic piece on the left is equal to the second derivative of the cubic piece on the right:

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

After some calculus (taking two derivatives of the Hermite basis form and applying the chain rule), this condition simplifies into a linear relationship between three neighboring slopes: $c_{j-1},c_j,c_{j+1}$.

$$\frac{1}{h_j}{c}_{j-1}+2\left(\frac{1}{h_j}+\frac{1}{h_{j+1}}\right)c_j+\frac{1}{h_{j+1}}{c}_{j+1}=3\left(\frac{y_j-y_{j-1}}{h_j^2}+\frac{y_{j+1}-y_j}{h_{j+1}^2}\right)$$

This gives us $N-1$ equations for our $N+1$ unknown slopes. We are two equations short! To get a unique solution, we must specify two additional boundary conditions.

Types of Spline Boundary Conditions

1. Clamped Spline: We explicitly specify the slopes at the endpoints, $s^{\prime }(x_0)=c_0$ and $s^{\prime }(x_N)=c_N$. This is useful if we know the derivative of the underlying function at the boundaries.
2. Natural Spline: This is the most common choice. It assumes the second derivative is zero at the endpoints: $s^{\prime \prime }(x_0)=s^{\prime \prime }(x_N)=0$. This corresponds physically to a flexible draftsman’s ruler (a “spline”) being held straight at its ends.
3. “Not-a-Knot” Spline: A clever condition that forces the third derivative to be continuous at the first and last interior knots ($x_1$ and $x_{N-1}$). Since our pieces are cubics, this forces the first two cubic pieces (and the last two) to be the same polynomial. This is the default in many software packages, including scipy.interpolate.CubicSpline.
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)$.

Once we’ve chosen our boundary conditions, we have a system of $N+1$ linear equations for the $N+1$ slopes. The matrix for this system is tridiagonal (or nearly so for periodic splines). This is fantastic news, because a tridiagonal system can be solved with a specialized, fast algorithm in just $O(N)$ operations, not the usual $O(N^3)$ for a dense matrix.

from scipy.interpolate import CubicSpline
 
x_data = np.array([0, 1, 3, 4, 6])
y_data = np.array([1, 4, 2, 2, 3])
 
# Create the spline objects with different boundary conditions
cs_notaknot = CubicSpline(x_data, y_data, bc_type='not-a-knot')
cs_natural = CubicSpline(x_data, y_data, bc_type='natural')
cs_clamped = CubicSpline(x_data, y_data, bc_type='clamped', extrapolate=False) # Clamped to zero slope
 
x_plot = np.linspace(0, 6, 400)
 
plt.plot(x_data, y_data, 'ro', label='Data Points')
plt.plot(x_plot, cs_notaknot(x_plot), label='Not-a-Knot (default)')
plt.plot(x_plot, cs_natural(x_plot), '--', label='Natural Spline')
plt.plot(x_plot, cs_clamped(x_plot), ':', label='Clamped Spline')
plt.title("Comparison of Cubic Spline Boundary Conditions")
plt.legend(); plt.grid(True); plt.show()

The Optimality of Splines

The natural cubic spline has a beautiful optimality property. Among all functions that are twice-differentiable and pass through the data points, the natural spline is the one that minimizes the total “bending energy,” defined as ${\int }_a^b(f^{\prime \prime }(x){)}^{2}dx$. It is, in a very precise mathematical sense, the “smoothest” possible interpolant.

Convergence Comparison

The choice of method matters. For smooth functions, the higher-order smoothness of splines pays off in faster convergence.

The plots show that for smooth functions, both clamped splines and true Hermite splines achieve the excellent $O(h^4)$ convergence rate. PCHIP, which prioritizes shape over smoothness, converges more slowly but is more robust against oscillations. The “not-a-knot” spline also achieves $O(h^4)$ convergence.

The overarching lesson is powerful. While a single, global polynomial is a delicate and sometimes dangerous tool, a chain of low-degree piecewise polynomials is robust, flexible, and forms the foundation of modern interpolation, approximation, and computer-aided design.