Interpolation: Piecewise & Shape-Preserving

Overview

Interpolation is the construction of a function that passes exactly through given nodes $(x_i, y_i)$ ($i = 0, \ldots, n$), used to estimate values between the nodes. The most naive approach is to build a single degree-$n$ polynomial interpolant through all $n+1$ points, but this breaks down as the degree grows.

The Runge phenomenon

When high-degree polynomial interpolation is performed on equally spaced nodes, the interpolating polynomial oscillates violently near the ends of the interval, and the oscillation worsens as more nodes are added. The classic example is $f(x) = 1/(1 + 25 x^2)$ interpolated with equally spaced nodes on $[-1, 1]$, where the error diverges exponentially at the ends. This is the Runge phenomenon.

The cause is that the amplitude of the Lagrange basis $\ell_i(x)$ at equally spaced nodes blows up explosively near the ends. There are broadly two remedies:

• Cluster the nodes toward the ends (e.g. Chebyshev nodes) — suppress the oscillation while keeping a single polynomial.
• Make it piecewise and local — connect each interval with a low-degree (1st to 3rd) polynomial and localize the error. This is piecewise interpolation.

Piecewise interpolation has a locality property: the value is determined solely by the coefficients of the single interval to which the evaluation point belongs. Moving one node affects only its neighbourhood, which makes it easy to work with in data-driven applications. This page covers the piecewise and shape-preserving methods that sangi implements — cubic splines, PCHIP, modified Akima, Catmull-Rom, and quintic Hermite — together with the underlying Lagrange and Newton interpolation on which they rest.

Related API: Interpolation.

Lagrange Interpolation and Newton Divided Differences

Before turning to the piecewise methods, we review the two classical forms of global polynomial interpolation on which they are founded. Both represent the same interpolating polynomial $p(x)$, but they differ in the efficiency of construction and evaluation.

The Lagrange basis

For nodes $x_0, \ldots, x_n$, define the Lagrange basis as

$$\ell_i(x) = \prod_{\substack{j=0 \\ j \neq i}}^{n} \frac{x - x_j}{x_i - x_j}$$

so that $\ell_i(x_k) = \delta_{ik}$ (1 when $i = k$, 0 otherwise). The interpolating polynomial is then

$$p(x) = \sum_{i=0}^{n} y_i\, \ell_i(x)$$

sangi's lagrange_interpolation follows this definition directly, requiring a double loop of $O(n^2)$ work per evaluation point $x_i$.

Newton divided differences

Writing the same polynomial in a form that adds nodes one at a time gives the Newton form:

$$p(x) = \sum_{i=0}^{n} c_i \prod_{j=0}^{i-1} (x - x_j)$$

The coefficients $c_i = f[x_0, \ldots, x_i]$ are the divided differences, computed by the recurrence:

$$f[x_i] = y_i, \qquad f[x_i, \ldots, x_{i+k}] = \frac{f[x_{i+1}, \ldots, x_{i+k}] - f[x_i, \ldots, x_{i+k-1}]}{x_{i+k} - x_i}$$

sangi's newton_interpolation_coefficients initializes the coefficient array with $y$ and constructs the divided-difference table by overwriting in place with a double loop ($O(n^2)$). Evaluation is handled by newton_interpolation_evaluate using Horner's method, at $O(n)$ per evaluation:

$$p(x) = c_0 + (x - x_0)\bigl(c_1 + (x - x_1)\bigl(c_2 + \cdots\bigr)\bigr)$$

The advantage of the Newton form is that when one node is added, the existing coefficients can be reused (one only appends a column to the difference table).

The barycentric Lagrange form

The barycentric form speeds Lagrange interpolation up to $O(n)$ per evaluation while also stabilizing it numerically. First define the weight of each node as

$$w_j = \frac{1}{\displaystyle\prod_{i \neq j} (x_j - x_i)}$$

(sangi's barycentricLagrangeWeights, computed once in $O(n^2)$). The interpolating polynomial can then be written as

$$p(x) = \frac{\displaystyle\sum_{j=0}^{n} \frac{w_j}{x - x_j}\, y_j}{\displaystyle\sum_{j=0}^{n} \frac{w_j}{x - x_j}}$$

(barycentricLagrangeEvaluate). The advantages of this second (barycentric) form:

• $O(n)$ evaluation: because the weights are reused, it is faster than the ordinary form when many points are evaluated on the same nodes.
• Numerical stability: the same $w_j/(x - x_j)$ appears in both the numerator and the denominator, so the common factor cancels and cancellation error is unlikely. When $x$ coincides with a node $x_j$, $1/(x - x_j)$ diverges, so the implementation detects $x - x_j = 0$ and returns $y_j$ directly.

That said, because the barycentric form is still a global polynomial, it cannot escape the Runge phenomenon on equally spaced nodes. It is the piecewise methods of the following sections that break this oscillation.

Cubic Splines

Connecting each interval $[x_i, x_{i+1}]$ with a cubic polynomial, continuous up to the first and second derivatives at the nodes, gives a cubic spline. Write the polynomial of interval $i$ in the local coordinate $t = x - x_i$ as

$$S_i(x) = a_i + b_i (x - x_i) + c_i (x - x_i)^2 + d_i (x - x_i)^3$$

With $n$ intervals there are $4n$ coefficients. They are determined by the following conditions.

Continuity conditions

• Interpolation condition: each interval passes through both endpoint nodes ($S_i(x_i) = y_i,\ S_i(x_{i+1}) = y_{i+1}$).
• $C^1$ continuity: the first derivative matches at interior nodes ($S_{i-1}'(x_i) = S_i'(x_i)$).
• $C^2$ continuity: the second derivative matches at interior nodes ($S_{i-1}''(x_i) = S_i''(x_i)$).

These alone fall two conditions short, so a boundary condition is imposed at each end. sangi's SplineBoundaryCondition offers two kinds:

• Natural spline (Natural): $S''(x_0) = S''(x_n) = 0$ at both ends. The most standard choice, setting the curvature to 0 at the ends.
• Clamped (Clamped): specify the values of the first derivative $S'(x_0), S'(x_n)$ at both ends. Used when you want to fix the end slopes to known values.

Deriving the tridiagonal system

Taking the coefficients $c_i$ ($i = 0, \ldots, n$) associated with the second derivative as unknowns, the $C^2$ continuity condition yields one equation at each interior node. With interval width $h_i = x_{i+1} - x_i$, at interior node $i$ we have

$$h_{i-1}\, c_{i-1} + 2(h_{i-1} + h_i)\, c_i + h_i\, c_{i+1} = 3\!\left(\frac{y_{i+1} - y_i}{h_i} - \frac{y_i - y_{i-1}}{h_{i-1}}\right)$$

Since each equation involves only the three adjacent unknowns $c_{i-1}, c_i, c_{i+1}$, the coefficient matrix is tridiagonal. Adding the boundary conditions (for the natural case, $c_0 = c_n = 0$) gives the linear system for the coefficient vector $c = (c_0, \ldots, c_n)^\top$

$$M\,c = d$$

Here $M$ is a diagonally dominant tridiagonal matrix and is nonsingular.

$O(n)$ solution via the Thomas algorithm

A tridiagonal system can be solved in $O(n)$ with the Thomas algorithm (forward elimination plus back substitution specialized for tridiagonal systems), without resorting to general $O(n^3)$ Gaussian elimination. sangi's cubic_spline_coefficients consists of these two stages:

• Forward elimination: update the diagonal entries as $l_i = 2(h_{i-1} + h_i) - h_{i-1}\,\mu_{i-1}$ and the superdiagonal as $\mu_i = h_i / l_i$, while computing the intermediate solution $z_i$ of the right-hand side in order.
• Back substitution: solve in reverse from $c_{n} = z_n$ via $c_i = z_i - \mu_i\, c_{i+1}$. Together with this, fix the interval coefficients $b_i, d_i$ by $b_i = \dfrac{y_{i+1} - y_i}{h_i} - \dfrac{h_i (c_{i+1} + 2 c_i)}{3}$, $d_i = \dfrac{c_{i+1} - c_i}{3 h_i}$.

Evaluation is handled by cubic_spline_evaluate, which finds the interval $i$ to which the evaluation point $x$ belongs and then computes $S_i(x)$ in a Horner-like fashion. Because the coefficient notation $\{a_i, b_i, c_i, d_i\}$ is shared with the other piecewise methods (PCHIP, Akima, Catmull-Rom), they can all obtain values from the same evaluation function.

Shape-Preserving Interpolation

The cubic spline is the smoothest ($C^2$), but on steeply rising data or step-like data it creates peaks or valleys (overshoot) that are not in the data. When you do not want to disturb the monotonicity or sign of measured values, you use a method that preserves the shape, even at the cost of a little smoothness. All of the following are based on the per-interval cubic Hermite polynomial

$$p_i(t) = y_i\,(1-t)^2(1+2t) + y_{i+1}\,t^2(3-2t) + d_i\,h_i\,t(1-t)^2 + d_{i+1}\,h_i\,t^2(t-1), \quad t = \frac{x - x_i}{h_i}$$

and differ only in how the slope $d_i$ at each node is chosen. If the slopes are converted into the per-interval coefficients $\{a,b,c,d\}$, they can be handled by the same evaluation function as the cubic spline.

PCHIP (Fritsch-Carlson)

PCHIP (Piecewise Cubic Hermite Interpolating Polynomial) is a method that preserves monotonicity using the Fritsch-Carlson algorithm (sangi's pchipCoefficients). With interval slope $\delta_i = (y_{i+1} - y_i)/h_i$, the slope at an interior node is determined as follows:

• The neighbouring interval slopes have the same sign ($\delta_{i-1}\,\delta_i > 0$): the harmonic mean weighted by interval width, $$d_i = \frac{w_1 + w_2}{\dfrac{w_1}{\delta_{i-1}} + \dfrac{w_2}{\delta_i}}, \quad w_1 = 2h_i + h_{i-1},\ w_2 = h_i + 2h_{i-1}$$ Because it is a harmonic mean, if either slope is 0 (flat) then $d_i = 0$.
• The slopes have opposite signs, or one of them is 0 (an extremum node): set $d_i = 0$, so that no peak or valley is created at that point.

For the endpoints, an initial estimate is made with the one-sided three-point Bessel formula; if $d_0$ and $\delta_0$ have opposite signs it is set to 0, and if $|d_0|$ exceeds $3|\delta_0|$ it is limited to $3\delta_0$.

Finally, the Fritsch-Carlson monotonicity condition is imposed on every interval. Setting $\alpha = d_i/\delta_i$ and $\beta = d_{i+1}/\delta_i$, a sufficient condition for the interval to be monotone is

$$\alpha^2 + \beta^2 \le 9$$

On any interval that violates this, multiply $d_i, d_{i+1}$ by $\tau = 3/\sqrt{\alpha^2 + \beta^2}$ to pull them back inside the circle of radius 3. Thanks to this shrinking, monotone data produce a monotone interpolant and no overshoot occurs. The result is $C^1$ continuous.

Modified Akima

Akima's method (1970) estimates the slope at each node from the local slopes of the four neighbouring intervals, weighted by the differences of the interval slopes. Its characteristic is that it does not react too sensitively to a single outlier-like point, giving a natural curve with little oscillation. sangi's modifiedAkimaCoefficients constructs the slope weights as

$$w_1 = |\delta_{i+1} - \delta_i| + \varepsilon, \qquad t_i = \frac{w_1\, \delta_{i-1} + w_2\, \delta_i}{w_1 + w_2}$$

In the original Akima method the denominator can become 0 (when successive interval slopes are equal); this modified version avoids and stabilizes that by adding a small quantity $\varepsilon$ ($\varepsilon$ is a scale formed by multiplying the maximum slope difference by machine precision). At both ends, virtual interval slopes are supplied by linear extrapolation and handled with the same formula as the interior. The result is $C^1$ continuous.

Catmull-Rom

The Catmull-Rom spline is a $C^1$ continuous interpolant whose tangent at each node is determined from its two immediate neighbours, widely used for curve interpolation in computer graphics (sangi's catmullRomCoefficients). The slope at an interior node is given by the secant connecting the surrounding nodes:

$$t_i = \frac{y_{i+1} - y_{i-1}}{x_{i+1} - x_{i-1}}$$

The endpoints use one-sided differences. Since it has no slope limiter, it looks smoother than PCHIP but does not guarantee monotonicity. Its advantage is the convenience of passing exactly through the control points, which suits applications where the shape is designed directly.

Quintic Hermite

Specifying at each node the value $f$, the first derivative $f'$, and the second derivative $f''$, and connecting the intervals with a degree-5 polynomial, gives quintic Hermite interpolation (sangi's quinticHermiteCoefficients). Each interval has six coefficients

$$S_i(t) = a + b\,t + c\,t^2 + d\,t^3 + e\,t^4 + f\,t^5, \quad t = x - x_i$$

and achieves $C^2$ continuity. From the boundary conditions $S(0)=y_0,\ S(h)=y_1,\ S'(0)=f'_0,\ \ldots,\ S''(h)=f''_1$ the three low-order coefficients $a, b, c$ are determined directly, and the three high-order coefficients $a_3, a_4, a_5$ are obtained by solving a $3 \times 3$ linear system with Gaussian elimination. When derivative information is at hand, it can deliver higher smoothness and accuracy than cubic Hermite.

Why no overshoot on monotone data

A cubic spline overshoots because, to satisfy $C^2$ continuity, the node slopes $d_i$ are determined globally as a coupled system and may take large slopes that deviate from the local data shape. Shape-preserving methods, by contrast, determine the node slopes locally and with an upper bound:

• Set the slope to 0 at an extremum node ($\delta_{i-1}$ and $\delta_i$ have opposite signs) → no peak or valley absent from the data.
• Limit the magnitude of the slope to within a constant multiple of the interval slope $\delta_i$ (the Fritsch-Carlson $\alpha^2 + \beta^2 \le 9$) → the cubic polynomial stays monotone within the interval.

This combination of "locality + slope limiting" is the theoretical basis for monotonicity (= no overshoot) on monotone data. $C^2$ continuity is sacrificed, however, and shape-preserving methods generally guarantee only up to $C^1$. Smoothness and shape preservation are in a trade-off relationship, and the method is chosen according to the application.

Comparison Table

A rough guide: for maximum smoothness use a cubic spline; to preserve monotonicity and sign use PCHIP; for a natural curve robust to outliers use modified Akima; to handle control points directly in CG use Catmull-Rom; and when derivatives are known and high accuracy is required use quintic Hermite. If you only look up values many times on the same nodes, barycentric Lagrange is fast.

References

• de Boor, C. (2001). A Practical Guide to Splines. Revised ed. Springer.
• Fritsch, F. N., & Carlson, R. E. (1980). "Monotone Piecewise Cubic Interpolation". SIAM Journal on Numerical Analysis, 17(2), 238–246.
• Akima, H. (1970). "A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures". Journal of the ACM, 17(4), 589–602.
• Berrut, J.-P., & Trefethen, L. N. (2004). "Barycentric Lagrange Interpolation". SIAM Review, 46(3), 501–517.
• Catmull, E., & Rom, R. (1974). "A Class of Local Interpolating Splines". In Computer Aided Geometric Design (pp. 317–326). Academic Press.