# Spline (mathematics)

Single knots at 1/3 and 2/3 establish a spline of three cubic polynomials meeting with C2 continuity. Triple knots at both ends of the interval ensure that the curve interpolates the end points. The shown curve is a parametric spline and not a polynomial function.

In mathematics, a spline is a function defined piecewise by polynomials. Spline functions are used to interpolate intervals between data points and are used extensively in the computer science subfields of computer-aided design and computer graphics to represent arbitrary curves and surfaces.

Splines are popular curves in these subfields because of the simplicity of their construction, their ease and accuracy of evaluation, and their capacity to approximate complex shapes through curve fitting and interactive curve design.[citation needed] In interpolation problems, spline interpolation is often preferred to polynomial interpolation because polynomial interpolation over a number of data points can lead to polynomials of high degrees and oscillations due to Runge's phenomenon, while a spline is made up of lower degree polynomials which avoid this problem.

The term spline comes from the flexible wooden strips called splines historically used by shipbuilders and draftsmen to draw smooth shapes.[1]

## Introduction

The term "spline" is used to refer to a wide class of functions that are used for two main purposes:

• Interpolation: Given a set of data points, a spline function is fitted to pass through each point while, typically, minimizing some measure of curvature to produce a smooth curve.
• Smoothing: Given a set of data points, often including noise, a spline function is fitted to give a 'line of best fit', in a process similar to least squares approximation. The objectives of fitting the spline are to minimize a combination of some measure of distance between the spline and data points, and some measure of curvature. By changing the weights of each objective, the resulting spline can provide a better fit to the data points, or a smoother curve, respectively.

Splines are described by the type of polynomial functions they are composed of. Common examples are:

• Polynomial spline, which are often described by the degree of the polynomial used, such as linear, quadratic, cubic, etc.
• B-spline or basis spline.
• Bézier spline as a spline constructed of piece-wise Bézier curves.
• Periodic spline, used to approximate periodic functions. The spline is constructed such that the polynomial in the last interval is the same as in the first.[2]
• Cubic Hermite spline or cubic Hermite interpolator is a spline where each piece is a third-degree polynomial specified in Hermite form.
• NURBS, or Non-Uniform Rational B-Spline. In these the piecewise polynomial is a ratio, or rational function, of B-spline functions. NURBS are able to model exactly lines, planes, conic sections and quadric surfaces.[3]

Splines may be univariate (e.g., 2-D curves), bivariate (e.g., 3-D surfaces) or be multivariate to describe functions of higher dimension.

## Definition

The spline ${\displaystyle S}$ maps values from an interval ${\displaystyle [a,b]}$ to the set of real numbers ${\displaystyle \mathbb {R} :}$

${\displaystyle S:[a,b]\to \mathbb {R} .}$

For the piecewise definition of ${\displaystyle S}$ the interval ${\displaystyle [a,b]}$ is partitioned in ${\displaystyle n}$ subintervals that together make up the whole interval. This is done by selecting additional ${\displaystyle n-1}$ numbers ${\displaystyle x_{k}(1\leq k\leq n-1)}$, representing points within this interval, such that

${\displaystyle a=x_{0}\leq x_{1}\leq \cdots \leq x_{n-1}\leq x_{n}=b,}$

yielding ${\displaystyle n}$ subintervals

${\displaystyle [x_{i},x_{i+1}]\quad {\text{for}}\quad 0\leq i\leq n-1,}$

such that

${\displaystyle [a,b]=[x_{0},x_{1}]\cup [x_{1},x_{2}]\cup \cdots \cup [x_{n-2},x_{n-1}]\cup [x_{n-1},x_{n}].}$

For each of these subintervals a polynomial ${\displaystyle P_{i}}$ will be constructed

${\displaystyle P_{i}:[x_{i},x_{i+1}]\to \mathbb {R} \quad {\text{for}}\quad 0\leq i\leq n-1.}$

Excluding the right endpoints of each subinterval (with exception of the rightmost) allows the piecewise definition of ${\displaystyle S}$ on ${\displaystyle [a,b]}$ by the ${\displaystyle n}$ polynomials ${\displaystyle P_{i}}$ on the ${\displaystyle i^{t}h}$ subinterval

${\displaystyle S(x)={\begin{cases}P_{0}(x)\quad {\text{for}}\quad x\in [x_{0},x_{1})\\P_{1}(x)\quad {\text{for}}\quad x\in [x_{1},x_{2})\\\quad \vdots \\P_{n-2}(x)\quad {\text{for}}\quad x\in [x_{n-2},x_{n-1})\\P_{n-1}(x)\quad {\text{for}}\quad x\in [x_{n-1},x_{n}].\end{cases}}}$

If the ${\displaystyle n}$ polynomials ${\displaystyle P_{i}}$ each have degree at most ${\displaystyle m}$, then the spline is said to be of degree ${\displaystyle m}$ (or is of order ${\displaystyle m+1}$).

### Parametric Representation

The definition given above described the spline ${\displaystyle S(x)}$ explicitly in terms of the independent variable ${\displaystyle x}$. Splines may also be represented parametrically in terms of an independent parameter called, for example, ${\displaystyle t}$.

A parametric curve in ${\displaystyle \mathbb {R} ^{n}}$ can be defined as a continuous map from an interval, the domain of the parameter, to the space ${\displaystyle \mathbb {R} ^{n}.}$ For example, a plane curve ${\displaystyle G}$ is a mapping from the interval ${\displaystyle [a,b]}$ to the Euclidean plane ${\displaystyle \mathbb {R} ^{2}}$

{\displaystyle {\begin{aligned}G:\;[a,b]&\to \mathbb {R} ^{2}\\t&\mapsto (X(t),Y(t)),\end{aligned}}}

with ${\displaystyle X}$ and ${\displaystyle Y}$ being real functions on ${\displaystyle [a,b].}$

If both functions ${\displaystyle X}$ and ${\displaystyle Y}$ are simple spline functions of the same degree with the same extended knot vectors on that interval, then

${\displaystyle G(t)=(X(t),Y(t))\quad {\text{for}}\quad t\in [a,b]}$

specifies a spline curve in ${\displaystyle \mathbb {R} ^{2}.}$.

### Continuity

The polynomial functions are continuous within each subinterval of the spline, and by definition meet exactly at the boundaries. This requirement means splines have at least class ${\displaystyle C^{0}}$ parametric continuity (i.e., continuity of position at the sub-interval boundaries). If the first derivatives of the polynomials are equal at subinterval boundaries, the spline has class ${\displaystyle C^{1}}$ continuity (continuity of slope). Equal second derivatives gives class ${\displaystyle C^{2}}$ continuity (continuity of curvature), and so on. As the degree of continuity increases, the spline becomes increasingly smooth.[4]

If the first ${\displaystyle r_{k}}$ derivatives of the two polynomials match at ${\displaystyle x_{k},}$ then ${\displaystyle S}$ has ${\displaystyle C^{r}}$ continuity at ${\displaystyle x_{k}.}$ Assuming a spline of degree ${\displaystyle m}$, the loss of smoothness (with respect to the maximal nontrivial smoothness) at the point ${\displaystyle x_{k}}$ is given by ${\displaystyle m-r_{k}.}$ A loss of ${\displaystyle m}$ at ${\displaystyle x_{k}}$, i.e. ${\displaystyle r_{k}=0,}$ establishes continuity there by imposing the condition ${\displaystyle P_{k-1}(x_{k})=P_{k}(x_{k});}$ a smaller loss, ${\displaystyle (1\leq r_{k}\leq m),}$ yields additional equations ${\displaystyle P_{k-1}^{(d)}(x_{k})=P_{k}^{(d)}(x_{k})}$ for ${\displaystyle 1\leq d\leq r_{k}}$ involving the derivatives up to the desired smoothness. Collecting the ${\displaystyle n-1}$ values ${\displaystyle r_{k}}$ yields the smoothness vector ${\displaystyle {\mathbf {r} }=(r_{1},\dots ,r_{n-1})}$ of the spline ${\displaystyle S}$.

Natural splines are those where second derivatives at ${\displaystyle a}$ and ${\displaystyle b}$ are set to zero.

#### Locally negative continuity

It might be asked what meaning more than n multiple knots in a knot vector have, since this would lead to continuities like

${\displaystyle S(t)\in C^{-m}{\mbox{ , }}m>0}$

at the location of this high multiplicity. By convention, any such situation indicates a simple discontinuity between the two adjacent polynomial pieces. This means that if a knot ti appears more than n + 1 times in an extended knot vector, all instances of it in excess of the (n + 1)th can be removed without changing the character of the spline, since all multiplicities n + 1, n + 2, n + 3, etc. have the same meaning. It is commonly assumed that any knot vector defining any type of spline has been culled in this fashion.

### Knots

The ${\displaystyle n+1}$ points ${\displaystyle x_{j}}$ ${\displaystyle (0\leq j\leq n)}$ are called knots. The vector ${\displaystyle {\mathbf {x} }=(x_{0},\dots ,x_{n})}$ is called a knot vector for the spline. If the knots are equidistantly distributed in the interval ${\displaystyle [a,b]}$ the spline is called uniform, otherwise it is non-uniform.

An inner knot ${\displaystyle x_{k}}$ can be "deleted" by making it equal to a neighboring knot ${\displaystyle x_{k\pm 1}}$. By identifying ${\displaystyle x_{k}}$ with ${\displaystyle x_{k+1}}$ the polynomial piece ${\displaystyle P_{k}}$ disappears, and the pieces ${\displaystyle P_{k-1}}$ and ${\displaystyle P_{k+1}}$ join with the sum of the continuity losses for ${\displaystyle x_{k}}$ and ${\displaystyle x_{k+1}}$. That is,

${\displaystyle S\in C^{m-j_{k}-j_{k+1}}[x_{k}=x_{k+1}],}$ where ${\displaystyle j_{k}=m-r_{k}.}$

This leads to a more general understanding of a knot vector. The continuity loss at any point can be considered to be the result of multiple knots located at that point, and a spline type can be completely characterized by its degree ${\displaystyle m}$ and its extended knot vector

${\displaystyle (x_{0},x_{1},\cdots ,x_{1},x_{2},\cdots ,x_{2},x_{3},\cdots ,x_{n-2},x_{n-1},\cdots ,x_{n-1},x_{n})}$

where ${\displaystyle x_{k}}$ is repeated ${\displaystyle j_{k}}$ times for ${\displaystyle k=1,\dots ,n-1}$. Obviously, ${\displaystyle n}$ does not refer anymore to the number of components in this vector.

### Spline Space

Given a knot vector ${\displaystyle {\mathbf {x} }}$, a degree ${\displaystyle m}$, and a smoothness vector ${\displaystyle {\mathbf {r} }}$ for this knot vector, one can consider the set of all splines of degree ${\displaystyle \leq m}$ having this knot vector ${\displaystyle {\mathbf {x} }}$ and this smoothness vector ${\displaystyle {\mathbf {r} }}$. Equipped with the operation of adding two piecewise defined polynomial functions (pointwise addition) and taking their real multiples, together with the zero function, this set becomes a real vector space. This spline space is commonly denoted by ${\displaystyle S_{m}^{\mathbf {r} }({\mathbf {x} })}$. Briefly this means that adding any two splines of a given type produces spline of that given type, and multiplying a spline of a given type by any constant produces a spline of that given type.

The dimension of the space containing all splines of a certain type can be counted from the extended knot vector:

${\displaystyle a=x_{0}<\underbrace {x_{1}=\cdots =x_{1}} _{j_{1}}<\cdots <\underbrace {x_{k-2}=\cdots =x_{k-2}} _{j_{k-2}}
${\displaystyle j_{i}\leq n+1~,~~i=1,\ldots ,k-2.}$

The dimension is equal to the sum of the degree plus the multiplicities

${\displaystyle d=n+\sum _{i=1}^{k-2}j_{i}.}$

If a type of spline has additional linear conditions imposed upon it, then the resulting spline will lie in a subspace. The space of all natural cubic splines, for instance, is a subspace of the space of all cubic ${\displaystyle C^{2}}$ splines.

## Polynomial Interpolation Splines

Suppose it is desired to provide a continuous interpolation for a data set given by the ${\displaystyle n+1}$ points

${\displaystyle x=(x_{0},x_{1},\cdots ,x_{n})}$
${\displaystyle y=(f_{0},f_{1},\cdots ,f_{n})}$

where

${\displaystyle a=x_{0}\leq x_{1}\leq \cdots \leq x_{n-1}\leq x_{n}=b}$

There will be ${\displaystyle n}$ intervals, each with a polynomial:

${\displaystyle P_{0}(x),x\in [x_{0},x_{1});P_{1}(x),x\in [x_{1},x_{2});\cdots ,P_{n-1}(x),x\in [x_{n-1},x_{n}]}$

If these polynomials have a degree ${\displaystyle \leq k}$, then spline ${\displaystyle S}$ will have a continuous ${\displaystyle (k-1)^{st}}$ derivative on ${\displaystyle [x_{0},x_{n}]}$.

### Linear Spline

A linear spline has degree 1 (i.e., ${\displaystyle k=1}$) and is of the form:

${\displaystyle S(x)={\begin{cases}P_{0}(x)=a_{0}+b_{0}x\quad {\text{for}}\quad x\in [x_{0},x_{1})\\P_{1}(x)=a_{1}+b_{1}x\quad {\text{for}}\quad x\in [x_{1},x_{2})\\\quad \vdots \\P_{n-1}=a_{n-1}+b_{n-1}x\quad {\text{for}}\quad x\in [x_{n-1},x_{n}].\end{cases}}}$

Each polynomial has two coefficients ${\displaystyle (a_{i},b_{i})}$ and two boundary conditions, ${\displaystyle p_{i}(x_{i})=f_{i}{\text{ and }}p_{i}(x_{i+1})=f_{i+1}}$. Solving for the coefficients gives the final form of the spline:

${\displaystyle S(x)={\begin{cases}P_{0}(x)=f_{0}{\tfrac {x-x_{1}}{x_{0}-x_{1}}}+f_{1}{\tfrac {x-x_{0}}{x_{1}-x_{0}}}\quad {\text{for}}\quad x\in [x_{0},x_{1})\\P_{1}(x)=f_{1}{\tfrac {x-x_{2}}{x_{1}-x_{2}}}+f_{2}{\tfrac {x-x_{1}}{x_{2}-x_{1}}}\quad {\text{for}}\quad x\in [x_{1},x_{2})\\\quad \vdots \\P_{n-1}=f_{n-1}{\tfrac {x-x_{n}}{x_{n-1}-x_{n}}}+f_{n}{\tfrac {x-x_{n-1}}{x_{n}-x_{n-1}}}\quad {\text{for}}\quad x\in [x_{n-1},x_{n}].\end{cases}}}$

The quadratic spline has degree 2 and so has continuous first derivatives across subintervals. Its form is

${\displaystyle S(x)={\begin{cases}P_{0}(x)=a_{0}+b_{0}x+c_{0}x^{2}\quad {\text{for}}\quad x\in [x_{0},x_{1})\\P_{1}(x)=a_{1}+b_{1}x+c_{1}x^{2}\quad {\text{for}}\quad x\in [x_{1},x_{2})\\\quad \vdots \\P_{n-1}=a_{n-1}+b_{n-1}x+c_{n-1}x^{2}\quad {\text{for}}\quad x\in [x_{n-1},x_{n}].\end{cases}}}$

Thus there are ${\displaystyle 3n}$ coefficients to solve for. In addition to the ${\displaystyle 2n}$ positional boundary conditions (i.e., ${\displaystyle p_{i}(x_{i})=f_{i}{\text{ and }}p_{i}(x_{i+1})=f_{i+1}}$), the first derivatives of the polynomials must be equal at the interior knots, i.e.,

${\displaystyle p_{i-1}'(x_{i})=p_{i}'(x_{i}),\quad i=1,2,\dots ,n-1.}$

This provides a further ${\displaystyle n-1}$ boundary conditions for a total of ${\displaystyle 3n-1}$ equations with which to solve for the ${\displaystyle 3n}$ coefficients. To derive a unique solution, the analyst must supply one more constraint, a value for the first derivative at either end point of the spline. The values of the coefficients may then be determined to give the final equations for the spline.

### Cubic spline

The cubic spline polynomials are:

${\displaystyle S(x)={\begin{cases}P_{0}(x)=a_{0}+b_{0}x+c_{0}x^{2}+d_{0}x^{3}\quad {\text{for}}\quad x\in [x_{0},x_{1})\\P_{1}(x)=a_{1}+b_{1}x+c_{1}x^{2}+d_{1}x^{3}\quad {\text{for}}\quad x\in [x_{1},x_{2})\\\quad \vdots \\P_{n-1}=a_{n-1}+b_{n-1}x+c_{n-1}x^{2}+d_{n-1}x^{3}\quad {\text{for}}\quad x\in [x_{n-1},x_{n}].\end{cases}}}$

The continuity requirements become

${\displaystyle p_{i}(x_{i})=f_{i}{\text{ and }}p_{i}(x_{i+1})=f_{i+1},\quad i=0,2,\dots ,n-1.}$
${\displaystyle p_{i-1}^{\prime }(x_{i})=p_{i}^{\prime }(x_{i}),\quad i=1,2,\dots ,n-1.}$
${\displaystyle p_{i-1}^{\prime \prime }(x_{i})=p_{i}^{\prime \prime }(x_{i}),\quad i=1,2,\dots ,n-1.}$

Continuity of the resulting spline is of class ${\displaystyle C^{2}}$.

There are a total of ${\displaystyle 3n-2}$ equations here to solve for the ${\displaystyle 3n}$ polynomial coefficients, so two more constraints are needed. Choices available to the analyst include:

• so-called natural cubic spline conditions :${\displaystyle p_{0}^{\prime \prime }(x_{0})=0,\quad p_{n}^{\prime \prime }(x_{n})=0}$
• use estimates of the actual second derivatives at the interval end points :${\displaystyle p_{0}^{\prime \prime }(x_{0})=f^{\prime \prime }(x_{0}),\quad p_{n}^{\prime \prime }(x_{n})=f^{\prime \prime }(x_{n})}$
• so-called not-a-knot condition on the first interior points, i.e., eliminate the knots at ${\displaystyle x_{1}{\text{ and }}x_{n-1}}$ which has the effect of making the polynomials on either side of these knots identical, i.e., ${\displaystyle P_{0}(x)\equiv P_{1}(x){\text{ and }}P_{n-1}(x)\equiv P_{n}(x)}$.
• so-called complete cubic spline condition where slopes of the end points are specified :${\displaystyle p_{0}^{\prime }(x_{0})=f^{\prime }(x_{0}),\quad p_{n}^{\prime }(x_{n})=f^{\prime }(x_{n})}$

Using one of these options, a linear system of equations for the ${\displaystyle 3n}$ unknown polynomial coefficients may be set up and solved.

A somewhat shorter solution can be obtained by applying some of the boundary conditions to the polynomials and providing an initial simplification to reduce the number of coefficients needing to be found. The polynomial form becomes

${\displaystyle P_{i}(x)={\frac {z_{i+1}(x-x_{i})^{3}}{6h_{i}}}+{\frac {z_{i}(x_{i+1}-x)^{3}}{6h_{i}}}+\left[{\frac {f(x_{i+1})}{h_{i}}}-{\frac {z_{i+1}h_{i}}{6}}\right](x-x_{i})+\left[{\frac {f(x_{i})}{h_{i}}}-{\frac {z_{i}h_{i}}{6}}\right](x_{i+1}-x)}$

where

• ${\displaystyle z_{i}=f^{\prime \prime }(x_{i})}$ are the values of the second derivative at the ${\displaystyle i^{th}}$ knot.
• ${\displaystyle h_{i}^{}=x_{i+1}-x_{i}}$
• ${\displaystyle f(x_{i}^{})}$ are the values of the function at the ${\displaystyle i^{th}}$ knot.

The second derivatives at the knots, ${\displaystyle z_{i}}$, are unknown yet, except for values possibly assigned to the end points ${\displaystyle z_{0},z_{n}}$. A system of linear equations in these unknown ${\displaystyle z_{i}}$ values may be set up as

${\displaystyle \mathbf {A} \mathbf {z} =\mathbf {d} }$

where

${\displaystyle \mathbf {z} ={\begin{bmatrix}z_{0},z_{1},\dots ,z_{n}\end{bmatrix}}^{T}}$
${\displaystyle \mathbf {d} ={\begin{bmatrix}d_{0},d_{1},\dots ,d_{n}\end{bmatrix}}^{T}}$
${\displaystyle d_{i}={\begin{cases}{\frac {6}{h_{1}}}(f_{1}-f_{0})-6f_{0}^{\prime }\quad {\text{for}}\quad i=0)\\{\frac {6}{h_{i+1}}}(f_{i+1}-f_{1})-{\frac {6}{h_{1}}}(f_{i}-f_{i-1})\quad {\text{for}}\quad 1\leq i\leq n-12\\6f_{n}^{\prime }-{\frac {6}{h_{n-1}}}(f_{n}-f_{n-1})\quad {\text{for}}\quad i=n\end{cases}}}$

The matrix ${\displaystyle \mathbf {A} }$ has a tridiagonal form

${\displaystyle \mathbf {A} ={\begin{bmatrix}2h_{0}&h_{0}\\h_{0}&2(h_{0}+h_{1})&h_{1}\\&h_{1}&2(h_{1}+h_{2})&h_{2}\\&&\ddots &\ddots &\ddots \\&&&h_{n-2}&2(h_{n-2}+h_{n-1})&h_{n-1}\\&&&&h_{n-1}&2h_{n-1}\end{bmatrix}}}$

which may be solved easily using a tridiagonal matrix algorithm.

## Examples

A quadratic spline is to be created to interpolate the points

${\displaystyle (-1,0),\quad (0,1),\quad (1,3).}$

The polynomial equations are:

${\displaystyle S(x)={\begin{cases}P_{0}(x)=a_{0}+b_{0}x+c_{0}x^{2}\quad {\text{for}}\quad x\in [-1,0)\\P_{1}(x)=a_{1}+b_{1}x+c_{1}x^{2}\quad {\text{for}}\quad x\in [0,1]\\\end{cases}}}$

Substituting the subinterval end points into each polynomial gives

${\displaystyle P_{0}(-1)=a_{0}-b_{0}+c_{0}=0}$
${\displaystyle P_{0}(0)=a_{0}=1}$
${\displaystyle P_{1}(0)=a_{1}=1}$
${\displaystyle P_{1}(1)=a_{1}+b_{1}+c_{1}=3}$

Equating the derivatives of the polynomials ${\displaystyle P_{0}'(0)=P_{1}'(0)}$ gives

${\displaystyle b_{0}=b_{1}}$

For the remaining constraint, we select the derivative at ${\displaystyle x=-1{\text{ to be }}P_{0}'(-1)=0}$ giving

${\displaystyle b_{0}-2c_{0}=0}$

Solving for the unknown coefficients finds their values as

${\displaystyle a_{0}=1,b_{0}=2,c_{0}=1,a_{1}=1,b_{1}=2,c_{1}=0}$

The spline function then becomes

${\displaystyle S(x)={\begin{cases}P_{0}(x)=1+2x+x^{2}\quad {\text{for}}\quad x\in [-1,0)\\P_{1}(x)=1+2x\quad {\text{for}}\quad x\in [0,1]\\\end{cases}}}$

A common spline is the natural cubic spline of degree 3 with continuity C2. The word "natural" means that the second derivatives of the spline polynomials are set equal to zero at the endpoints of the interval of interpolation

${\displaystyle S''(a)\,=S''(b)=0.}$

### Algorithm for computing natural cubic splines

Cubic splines have polynomial pieces of the form ${\displaystyle P_{i}(t)=a_{i}+b_{i}(t-x_{i})+c_{i}(t-x_{i})^{2}+d_{i}(t-x_{i})^{3}.}$ Given ${\displaystyle n+1}$ coordinates ${\displaystyle (x_{0},y_{0}),(x_{1},y_{1}),\dots ,(x_{n},y_{n}),}$ we find ${\displaystyle n}$ polynomials, ${\displaystyle P_{0}}$ and ${\displaystyle P_{k},}$ which satisfy for ${\displaystyle 1\leq k\leq n-1}$:

${\displaystyle P_{0}(x_{0})=y_{0}\quad }$ and ${\displaystyle \quad P_{k-1}(x_{k})=y_{k}=P_{k}(x_{k}),}$
${\displaystyle P'_{k-1}(x_{k})=P'_{k}(x_{k}),}$
${\displaystyle P''_{k-1}(x_{k})=P''_{k}(x_{k}),}$
${\displaystyle P''_{0}(x_{0})=P''_{n-1}(x_{n})=0.}$

One such polynomial ${\displaystyle P}$ is given by a 4-tuple ${\displaystyle (a,b,c,d)}$ where ${\displaystyle a,b,c\,}$ and ${\displaystyle d\,}$ correspond to the coefficients as used above.

Computation of Natural Cubic Splines:
Input: a set of ${\displaystyle k+1}$ coordinates
Output: a spline as a set of polynomial pieces, each represented by a 5-tuple.

1. Create a new array a of size n + 1, and for ${\displaystyle j=0,\ldots ,n}$ set ${\displaystyle a_{j}=y_{j}}$
2. Create new arrays b, d and μ each of size n
3. Create a new array h of size k and for ${\displaystyle i=0,\ldots ,n-1}$ set ${\displaystyle h_{i}=x_{i+1}-x_{i}}$
4. Create a new array β of size n-1 and for ${\displaystyle k=1,\ldots ,n-1}$ set ${\displaystyle \beta _{k}={\tfrac {3}{h_{k}}}(a_{k+1}-a_{k})-{\tfrac {3}{h_{k-1}}}(a_{k}-a_{k-1})}$
5. Create new arrays c, l, and z each of size ${\displaystyle n+1}$.
6. Set ${\displaystyle l_{0}=1,\;\mu _{0}=z_{0}=0}$
7. For ${\displaystyle i=1,\ldots ,n-1\,}$
1. Set ${\displaystyle l_{i}=2(x_{i+1}-x_{i-1})-h_{i-1}\mu _{i-1}.}$
2. Set ${\displaystyle \mu _{i}={\tfrac {h_{i}}{l_{i}}}.}$
3. Set ${\displaystyle z_{i}={\tfrac {\beta _{i}-h_{i-1}z_{i-1}}{l_{i}}}.}$
8. Set ${\displaystyle l_{n}=1;z_{n}=c_{n}=0.}$
9. For ${\displaystyle i=n-1,n-2,\ldots ,0}$
1. Set ${\displaystyle c_{i}=z_{i}-\mu _{i}c_{i+1}}$
2. Set ${\displaystyle b_{i}={\tfrac {a_{i+1}-a_{i}}{h_{i}}}-{\tfrac {h_{i}(c_{i+1}+2c_{i})}{3}}}$
3. Set ${\displaystyle d_{i}={\tfrac {c_{i+1}-c_{i}}{3h_{i}}}.}$
10. Create the spline as a new set of polynomials and call it output_set. Populate it with n 4-tuples for the polynomials Pi.
11. For ${\displaystyle i=0,\ldots ,k-1}$
1. Set Pi,a = ai
2. Set Pi,b = bi
3. Set Pi,c = ci
4. Set Pi,d = di
12. Output output_set

## History

Traditionally, splines were long flexible strips of wood used for lofting the complex curves of boat hulls. These drafting techniques were adopted as well by the aircraft and automobile body styling industries. The strips would be held in place at discrete points by weights (called "ducks", "dogs" or "rats") and between these points would assume shapes of minimum strain energy. It was discovered in the mid-1700s by Euler and the brothers Bernoulli that the shape of the spline forms a piecewise cubic polynomial with no curvature at the ends (i.e., the "natural spline" condition).[5]

A wooden spline

It is commonly accepted that the first use of the term "spline function" is in the 1946 papers by Schoenberg[6][7], which described it in terms of a smooth, piecewise polynomial approximation.[8]

In the foreword to (Bartels et al., 1987)[9], Robin Forrest describes lofting used in the British aircraft industry during World War II to construct templates, or patterns, that contained the full-size geometric data needed for airplane manufacturing. Aircraft factories stored large numbers of these templates, and Forrest suggests one possible impetus for moving from physical to mathematical models for aircraft geometry was the potential catastrophic loss of the critical design data should the loft be hit by an enemy bomb. This desire gave rise to "conic lofting", which used conic sections to model the position of the curve between the ducks. Conic lofting was replaced by what we would call splines in the early 1960s based on work by J. C. Ferguson[10] at Boeing and (somewhat later) by M. A. Sabin at British Aircraft Corporation.

The use of splines for modeling automobile bodies seems to have several independent beginnings. Credit is claimed on behalf of de Casteljau at Citroën, Pierre Bézier at Renault, and Birkhoff[11], Garabedian, and de Boor at General Motors (see Birkhoff and de Boor, 1965), all for work occurring in the very early 1960s or late 1950s. At least one of de Casteljau's papers was published, but not widely, in 1959. De Boor's work at General Motors resulted in a number of papers being published in the early 1960s, including some of the fundamental work on B-splines.

Work was also done at Pratt & Whitney Aircraft, where two of the authors of (Ahlberg et al., 1967)[12]—the first book-length treatment of splines—were employed, and the David Taylor Model Basin, by Feodor Theilheimer. The work at General Motors is detailed nicely in (Birkhoff, 1990)[13] and (Young, 1997)[14]. Davis (1997)[15] summarizes some of this material.

## References

1. ^ Wegman, Edward J., and Ian W. Wright (1983). "Splines in statistics". Journal of the American Statistical Association. 78 (382): 351–365. JSTOR 2288640.CS1 maint: Multiple names: authors list (link)
2. ^ Schumaker, Larry (2007). Spline functions : basic theory (3rd ed.). Cambridge University Press. p. 297. ISBN 978-0521705127.
3. ^ Rogers, David F. (2000). An introduction to NURBS : with historical perspective. Morgan Kaufmann Publishers. p. 129. ISBN 978-1558606692.
4. ^ Schumaker, Larry (2007). Spline functions : basic theory (3rd ed.). Cambridge University Press. p. 13. ISBN 978-0521705127.
5. ^ Schumaker, Larry (2007). Spline functions : basic theory (3rd ed.). Cambridge University Press. p. 6. ISBN 978-0521705127.
6. ^ Schoenberg, Isaac J. (1946). "Contributions to the Problem of Approximation of Equidistant Data by Analytic Functions: Part A.—On the Problem of Smoothing or Graduation. A First Class of Analytic Approximation Formulae" (PDF). Quarterly of Applied Mathematics. 4 (2): 45–99.
7. ^ Schoenberg, Isaac J. (1946). "Contributions to the Problem of Approximation of Equidistant Data by Analytic Functions: Part B.—On the Problem of Osculatory Interpolation. A Second Class of Analytic Approximation Formulae" (PDF). Quarterly of Applied Mathematics. 4 (2): 112–141.
8. ^ Schumaker, Larry (2007). Spline functions : basic theory (3rd ed.). Cambridge University Press. p. 6. ISBN 978-0521705127.
9. ^ Bartels, Richard H.; Beatty, Brian A.; Barsky, John C. (1987). An Introduction to Splines for Use in Computer Graphics and Geometric Modeling. San Mateo: Morgan Kaufmann. ISBN 0-934613-27-3.
10. ^ Ferguson, James C. (1964). "Multivariable Curve Interpolation". Journal of the ACM. 11 (2): 221–228. doi:10.1145/321217.321225.
11. ^ Birkhoff; de Boor (1965). "Piecewise polynomial interpolation and approximation". In Garabedian, H. L. (ed.). Proc. General Motors Symposium of 1964. New York and Amsterdam: Elsevier. pp. 164–190.
12. ^ Ahlberg, J. Harold; Nielson, Edwin N.; Walsh, Joseph L. (1967). The Theory of Splines and Their Applications. New York: Academic Press. ISBN 0-12-044750-9.
13. ^ Birkhoff (1990). "Fluid dynamics, reactor computations, and surface representation". In Nash, Steve (ed.). A History of Scientific Computation. Reading: Addison-Wesley. ISBN 0-201-50814-1.
14. ^ Young (1997). "Garrett Birkhoff and applied mathematics". Notices of the AMS. 44 (11): 1446–1449.
15. ^ Davis (1997). "B-splines and Geometric design". SIAM News. 29 (5).