This page introduces the mathematical foundations of B-splines as used throughout FITPACK. All curve and surface fitting routines in the library represent their results as linear combinations of B-spline basis functions, so understanding these foundations is essential for interpreting knot vectors, coefficients, and the various evaluation, differentiation, and integration operations.

The presentation follows Dierckx (1993), Chapters 1 and 6.

Piecewise Polynomials and Knot Vectors

A spline of degree \( k \) is a piecewise polynomial function whose pieces join with \( k - 1 \) continuous derivatives at designated breakpoints called knots. The knots are collected into a non-decreasing knot vector

\[ t_0 \leq t_1 \leq \cdots \leq t_{n+k+1} \]

which partitions the approximation interval \( [a, b] \) into sub-intervals. The number of distinct interior knots determines the flexibility of the spline: more knots allow closer approximation to complicated shapes, while fewer knots enforce smoothness.

The set of all splines of degree \( k \) on a given knot vector \( \mathbf{t} \) forms a linear space \( \mathcal{S}(k, \mathbf{t}) \) of dimension \( n + 1 \), where \( n + 1 \) is the number of B-spline basis functions (and therefore the number of coefficients needed to represent any element of the space).

In FITPACK the first and last \( k + 1 \) knots are clamped to the boundary:

\[ t_0 = t_1 = \cdots = t_k = a, \qquad t_{n+1} = t_{n+2} = \cdots = t_{n+k+1} = b \]

so that the spline interpolates (rather than merely approximates) its coefficient values at the endpoints.

B-Spline Basis Functions

Every element of \( \mathcal{S}(k, \mathbf{t}) \) can be written uniquely as a linear combination of B-spline basis functions \( N_{i,k}(x) \), which are defined recursively by the Cox–de Boor recurrence.

Degree zero

For \( k = 0 \) the basis functions are simple indicator functions:

\[ N_{i,0}(x) = \begin{cases} 1 & \text{if } t_i \leq x < t_{i+1}, \\ 0 & \text{otherwise}. \end{cases} \]

Recurrence for degree \( k \geq 1 \)

Higher-degree basis functions are built from lower-degree ones:

\[ N_{i,k}(x) = \frac{x - t_i}{t_{i+k} - t_i} \, N_{i,k-1}(x) + \frac{t_{i+k+1} - x}{t_{i+k+1} - t_{i+1}} \, N_{i+1,k-1}(x) \]

where any fraction with a zero denominator is taken to be zero (this arises when consecutive knots coincide). The routine fpbspl implements this recurrence.

Key properties

The B-spline basis functions enjoy several properties that make them ideal for numerical computation:

Non-negativity: \( N_{i,k}(x) \geq 0 \) for all \( x \).
Local support: \( N_{i,k}(x) = 0 \) outside the interval \( [t_i,\, t_{i+k+1}] \). Each basis function is nonzero on at most \( k + 1 \) knot spans, so spline evaluation is a local operation.
Partition of unity: At every point in \( [a, b] \) the basis functions sum to one: \( \sum_{i} N_{i,k}(x) = 1 \). This ensures that any convex combination of coefficients produces a spline value within the convex hull of those coefficients.
Smoothness: At a simple knot (multiplicity 1), \( N_{i,k} \) has \( k - 1 \) continuous derivatives. Each additional coincident knot reduces smoothness by one order.

See also: Dierckx, Ch. 1, §1.1–1.2

Spline as a Linear Combination

Any spline \( s \in \mathcal{S}(k, \mathbf{t}) \) is represented as

\[ s(x) = \sum_{i=0}^{n} c_i \, N_{i,k}(x) \]

where the \( c_i \) are the B-spline coefficients (also called control points or de Boor points). The entire fitting problem therefore reduces to determining the \( n + 1 \) coefficients for a given knot vector.

In FITPACK the coefficient array is stored alongside the knot vector in each derived type (e.g., fitpack_curve, fitpack_surface). Once computed by the fitting routines, the same coefficients are used by all evaluation, derivative, and integration routines.

Evaluation: the de Boor Algorithm

Given a point \( x \in [t_j,\, t_{j+1}) \), the spline value \( s(x) \) can be computed without explicitly forming all basis functions. The de Boor algorithm is a triangular recurrence analogous to Horner's method for ordinary polynomials.

Define the initial values

\[ d_i^{[0]} = c_i, \qquad i = j-k, \ldots, j \]

and compute for \( r = 1, \ldots, k \):

\[ d_i^{[r]} = \frac{x - t_i}{t_{i+k+1-r} - t_i} \, d_i^{[r-1]} + \frac{t_{i+k+1-r} - x}{t_{i+k+1-r} - t_i} \, d_{i-1}^{[r-1]}, \qquad i = j-k+r, \ldots, j \]

After \( k \) stages the result is \( s(x) = d_j^{[k]} \). The cost is \( \mathcal{O}(k^2) \) per evaluation point, independent of the total number of knots. This is the core loop in splev and underlies every eval method in the library.

See also: Dierckx, Ch. 1, Eq. 1.11

Derivatives

Differentiation of a B-spline representation yields another B-spline of one degree lower. The derivative of \( s(x) = \sum_i c_i \, N_{i,k}(x) \) is

\[ s'(x) = \sum_{i=1}^{n} c_i^{(1)} \, N_{i,k-1}(x) \]

where the new coefficients are obtained by differencing:

\[ c_i^{(1)} = \frac{k}{t_{i+k} - t_i} \, (c_i - c_{i-1}) \]

Higher-order derivatives follow by repeated application: the \( \nu \)-th derivative is a spline of degree \( k - \nu \) with coefficients obtained from \( \nu \) successive differencing steps. This structure is exploited by splder, which computes derivatives of any order up to \( k \).

Because each differentiation reduces the degree by one, at most \( k \) derivatives are nonzero. A cubic spline ( \( k = 3 \)) therefore has a piecewise-constant third derivative and a zero fourth derivative, regardless of the coefficients.

See also: Dierckx, Ch. 1, §1.3

Integration

The definite integral of a spline over an interval can be expressed directly in terms of its coefficients, without numerical quadrature. For a spline \( s(x) = \sum_i c_i \, N_{i,k}(x) \) on the interval \( [a, b] \):

\[ \int_a^b s(x) \, dx = \sum_{i=0}^{n} c_i \, \frac{t_{i+k+1} - t_i}{k + 1} \]

This formula follows from the property that the integral of each basis function satisfies

\[ \int_a^b N_{i,k}(x) \, dx = \frac{t_{i+k+1} - t_i}{k + 1} \]

when the knot vector is clamped at both ends. The routine splint implements this formula, giving exact results (up to floating-point arithmetic) regardless of the complexity of the spline.

See also: Dierckx, Ch. 1, §1.4

Knot Insertion (Oslo Algorithm)

A fundamental operation on B-spline representations is knot insertion: adding a new knot \( \bar{t} \) to the knot vector without changing the spline function itself. The resulting representation has one more knot, one more coefficient, and the same function values everywhere.

Suppose a knot \( \bar{t} \in [t_j,\, t_{j+1}) \) is inserted. The new coefficients \( \bar{c}_i \) are obtained from the old ones by a convex combination:

\[ \bar{c}_i = \begin{cases} c_i & \text{if } i \leq j - k, \\[4pt] \alpha_i \, c_i + (1 - \alpha_i) \, c_{i-1} & \text{if } j - k < i \leq j, \\[4pt] c_{i-1} & \text{if } i > j, \end{cases} \]

with weights

\[ \alpha_i = \frac{\bar{t} - t_i}{t_{i+k} - t_i}. \]

Because \( 0 \leq \alpha_i \leq 1 \), the new coefficients lie in the convex hull of their neighbours—a key property that preserves shape.

Knot insertion is the basis of the Oslo algorithm (Cohen, Lyche, and Schumaker, 1986), which generalizes the operation to insert several knots at once. In FITPACK the insert routine performs single knot insertion and is used both for refinement and for the zero-finding algorithm in sproot.

See also: Dierckx, Ch. 6, §6.2 (pp. 114–118)

Tensor Products for Bivariate Splines

The B-spline framework extends naturally to two dimensions via tensor products. Given knot vectors \( \mathbf{t}_x \) and \( \mathbf{t}_y \) with associated basis functions \( N_{i,k_x}(x) \) and \( M_{j,k_y}(y) \), a bivariate spline is

\[ s(x, y) = \sum_{i=0}^{n_x} \sum_{j=0}^{n_y} c_{ij} \, N_{i,k_x}(x) \, M_{j,k_y}(y) \]

The coefficient array \( c_{ij} \) has \( (n_x + 1)(n_y + 1) \) entries and is stored in column-major order in FITPACK.

Evaluation of a tensor-product spline decomposes into two sequences of univariate evaluations: first evaluate in \( x \) for each row of coefficients, then evaluate the resulting intermediate values in \( y \). This separability is the reason tensor-product splines scale well to two dimensions.

The routines bispev (grid evaluation) and bispeu (scattered evaluation) implement this strategy. Partial derivatives (parder, pardeu) and double integrals (dblint) follow the same tensor-product decomposition.

See also: Dierckx, Ch. 1, §1.5