Cubic Spline Interpolation

A cubic spline interpolates data with a sequence of cubic polynomials rather than one high-degree polynomial across the entire interval.

For ordered nodes

$$ x_0<x_1<\cdots<x_n, $$

one cubic is used on each interval

$$ [x_i,x_{i+1}]. $$

The pieces are joined so that the function, first derivative, and second derivative are continuous at every interior knot. The result is a smooth curve with local low-degree behavior.

Why use piecewise cubics?

A straight-line interpolant is simple and local, but it has slope discontinuities at the nodes. A single high-degree polynomial is smooth, but it can oscillate strongly and becomes increasingly sensitive as the number of nodes grows.

Cubic splines occupy a useful middle ground:

• each segment has low degree;
• the full interpolant is smooth;
• the linear system has a sparse tridiagonal structure;
• after setup, evaluating one query requires only one cubic segment.

Segment definition

On interval $[x_i,x_{i+1}]$, write

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

For $n$ intervals there are $4n$ segment coefficients.

They are determined by four types of conditions.

1. Interpolation

Each segment matches its endpoint values:

$$ S_i(x_i) = y_i, \qquad S_i(x_{i+1}) = y_{i+1}. $$

2. First-derivative continuity

At an interior knot,

$$ S_{i-1}'(x_i) = S_i'(x_i). $$

3. Second-derivative continuity

Also require

$$ S_{i-1}''(x_i) = S_i''(x_i). $$

4. Two boundary conditions

The interior conditions alone do not uniquely determine the spline. Two endpoint conditions are required.

This note focuses on the natural spline:

$$ S''(x_0) = 0, \qquad S''(x_n) = 0. $$

Other common choices are clamped, periodic, and not-a-knot boundary conditions.

A natural boundary condition sets the endpoint curvature to zero; it does not set the endpoint slope to zero.

Solving with knot second derivatives

A convenient derivation introduces

$$ M_i = S''(x_i). $$

Let

$$ h_i = x_{i+1} - x_i. $$

For a natural spline,

$$ M_0 = 0, \qquad M_n = 0. $$

The interior second derivatives satisfy the tridiagonal equations

$$ h_{i-1}M_{i-1} + 2(h_{i-1} + h_i)M_i + h_iM_{i+1} = 6 \left[\frac{y_{i+1}-y_i}{h_i} - \frac{y_i-y_{i-1}}{h_{i-1}} \right] $$

for

$$ i = 1,\ldots,n - 1. $$

The right-hand side measures the change in neighboring secant slopes. If the data are locally close to a straight line, that difference is small; if the slope changes strongly, the spline needs more curvature.

Deriving the tridiagonal system

A cubic has a linear second derivative. Interpolating that derivative between its endpoint values gives

$$ S_i''(x)=M_i\frac{x_{i+1}-x}{h_i}+M_{i+1}\frac{x-x_i}{h_i}. $$

Integrate twice and choose the integration constants so that $S_i(x_i)=y_i$ and $S_i(x_{i+1})=y_{i+1}$. The endpoint first derivatives are then

$$ S_i'(x_i)=\frac{y_{i+1}-y_i}{h_i}-\frac{h_i}{6}(2M_i+M_{i+1}), $$

$$ S_i'(x_{i+1})=\frac{y_{i+1}-y_i}{h_i}+\frac{h_i}{6}(M_i+2M_{i+1}). $$

Equating $S_{i-1}'(x_i)$ and $S_i'(x_i)$ and multiplying by six produces the displayed tridiagonal equation. Only three neighboring curvatures appear, explaining the sparse structure.

Counting conditions gives the same result: the $2n$ endpoint-value conditions, $n-1$ first-derivative matches, and $n-1$ second-derivative matches total $4n-2$. The two boundary conditions supply the remaining equations for the $4n$ coefficients. With positive interval widths, the natural-spline interior matrix is strictly diagonally dominant, which supports a pivot-free tridiagonal solve.

Recovering the cubic coefficients

Once the $M_i$ values are known,

$$ a_i = y_i, $$

$$ b_i = \frac{y_{i+1}-y_i}{h_i} - \frac{h_i}{6}(2M_i + M_{i+1}), $$

$$ c_i = \frac{M_i}{2}, $$

and

$$ d_i = \frac{M_{i+1}-M_i}{6h_i}. $$

Thus

$$ S_i(x) = a_i + b_i\Delta x + c_i\Delta x^2 + d_i\Delta x^3, \qquad \Delta x = x - x_i. $$

Why the joins are smooth

The second derivative on a segment is

$$ S_i''(x) = 2c_i + 6d_i(x - x_i). $$

At the left endpoint,

$$ S_i''(x_i) = 2c_i = M_i. $$

At the right endpoint,

$$ S_i''(x_{i+1}) = 2c_i + 6d_ih_i = M_{i+1}. $$

Adjacent segments therefore share the same second derivative at each knot. The tridiagonal system was derived from the matching of the first derivatives, so those are continuous as well.

Worked example

Use

$$ (0,0), \qquad(1,0.5), \qquad(2,0). $$

Both intervals have width

$$ h_0 = h_1 = 1. $$

For a natural spline,

$$ M_0 = M_2 = 0. $$

The only interior equation is

$$ 1\cdot M_0 + 2(1 + 1)M_1 + 1\cdot M_2 = 6 \left[\frac{0-0.5}{1} - \frac{0.5-0}{1} \right]. $$

Therefore,

$$ 4M_1 = -6, $$

so

$$ M_1 = -1.5. $$

For the first interval,

$$ a_0 = 0, \qquad b_0 = 0.75, \qquad c_0 = 0, \qquad d_0 = -0.25. $$

Hence

$$ S_0(x) = 0.75x - 0.25x^3. $$

At $x=0.5$,

$$ S_0(0.5) = 0.75(0.5) - 0.25(0.5)^3 = 0.34375. $$

Thus

$$ \boxed{S(0.5)=0.34375}. $$

The second piece is $S_1(x)=0.5-0.75(x-1)^2+0.25(x-1)^3$. At the shared knot $x=1$, both pieces have value $0.5$, slope $0$, and second derivative $-1.5$. At $x=0$ and $x=2$, the second derivatives are zero. These checks verify both the interpolation and the natural boundary conditions.

Algorithm

For a natural cubic spline:

Step 1. Sort the points by increasing $x_i$;

Step 2. Compute interval widths

$$ h_i = x_{i+1} - x_i; $$

Step 3. Assemble the tridiagonal system for $M_1,\ldots,M_{n-1}$;

Step 4. Solve that system;

Step 5. Set $M_0=M_n=0$;

Step 6. Compute $(a_i,b_i,c_i,d_i)$ for each interval;

VII. for each query, locate its interval and evaluate the corresponding cubic.

For this strictly diagonally dominant natural-spline system, the Thomas algorithm gives a pivot-free solve in $O(n)$ time and $O(n)$ storage. A general nonsingular tridiagonal matrix does not automatically have safe nonzero pivots, so the same claim requires care outside this spline setting.

Evaluation cost

After preprocessing, a query requires:

• interval search: $O(\log n)$ with binary search;
• polynomial evaluation: constant work.

The cubic can be evaluated in nested form:

$$ S_i(x) = a_i + \Delta x \left[b_i + \Delta x(c_i + \Delta x d_i) \right]. $$

Boundary conditions matter

Different endpoint conditions produce different interpolants even when all interior data are identical.

Common options include:

• natural:

$$ S''(x_0) = S''(x_n) = 0; $$

• clamped: prescribed endpoint slopes;
• periodic: endpoint values and derivatives match periodically;
• not-a-knot: the first two and last two cubic pieces are constrained to behave as if the first and last interior knots were not true breaks.

When comparing with a library, make sure the boundary condition matches. A library default may not be natural.

Polynomial degree does not override the boundary conditions

A cubic piece can represent a quadratic, but a natural spline does not necessarily reproduce a quadratic exactly: the imposed endpoint second derivatives may disagree with those of the quadratic. For $f(x)=x^2$ sampled at $0,1,2$, the natural conditions set $M_0=M_2=0$, while the interior equation gives $M_1=3$. The first segment is then $S_0(x)=x/2+x^3/2$, so $S(0.5)=0.3125$, compared with $f(0.5)=0.25$.

The discrepancy comes from the boundary model, not a failure to interpolate the supplied values. Clamping the endpoint slopes to the true values $0$ and $4$ instead reproduces this quadratic. Boundary conditions should therefore reflect available endpoint information, especially when endpoint derivatives or near-boundary accuracy matter.

Numerical behavior and limitations

Cubic splines usually avoid the dramatic oscillations of high-degree global polynomials, but they are not automatically shape preserving.

A cubic spline can:

• overshoot between monotone data points;
• become negative between positive samples;
• change shape globally when a data value changes because the coefficient system is coupled.

If monotonicity is required, use a shape-preserving interpolation method designed for that constraint.