This page was generated from doc/euclidean/natural-uniform.ipynb. Interactive online version: .

Uniform Natural Splines

For deriving natural splines, we first look at the uniform case, which means that the parameter interval in each segment is chosen to be \(1\).

The more general case with arbitrary parameter intervals is derived in a separate notebook about non-uniform natural splines.

import sympy as sp
sp.init_printing(order='grevlex')

We import some helpers from utility.py:

from utility import NamedExpression, dotproduct

t = sp.symbols('t')

To get started, let’s look at two neighboring segments: Let’s say the fourth segment, from \(\boldsymbol{x}_3\) to \(\boldsymbol{x}_4\), defined by the polynomial \(\boldsymbol{p}_3\), and the fifth segment, from \(\boldsymbol{x}_4\) to \(\boldsymbol{x}_5\), defined by the polynomial \(\boldsymbol{p}_4\). In both cases, we use \(0 \le t \le 1\).

coefficients3 = sp.symbols('a:dbm3')[::-1]
coefficients4 = sp.symbols('a:dbm4')[::-1]

We apply these coefficients to the monomial basis …

b_monomial = t**3, t**2, t, 1

… to define the two polynomials …

p3 = NamedExpression('pbm3', dotproduct(b_monomial, coefficients3))
p4 = NamedExpression('pbm4', dotproduct(b_monomial, coefficients4))
display(p3, p4)

$\displaystyle \boldsymbol{p}_{3} = \boldsymbol{d}_{3} t^{3} + \boldsymbol{c}_{3} t^{2} + \boldsymbol{b}_{3} t + \boldsymbol{a}_{3}$

$\displaystyle \boldsymbol{p}_{4} = \boldsymbol{d}_{4} t^{3} + \boldsymbol{c}_{4} t^{2} + \boldsymbol{b}_{4} t + \boldsymbol{a}_{4}$

… and we calculate their first derivatives:

pd3 = p3.diff(t)
pd4 = p4.diff(t)
display(pd3, pd4)

$\displaystyle \frac{d}{d t} \boldsymbol{p}_{3} = 3 \boldsymbol{d}_{3} t^{2} + 2 \boldsymbol{c}_{3} t + \boldsymbol{b}_{3}$

$\displaystyle \frac{d}{d t} \boldsymbol{p}_{4} = 3 \boldsymbol{d}_{4} t^{2} + 2 \boldsymbol{c}_{4} t + \boldsymbol{b}_{4}$

From this, we obtain 8 equations containing the 8 yet unknown coefficients.

equations = [
    p3.evaluated_at(t, 0).with_name('xbm3'),
    p3.evaluated_at(t, 1).with_name('xbm4'),
    p4.evaluated_at(t, 0).with_name('xbm4'),
    p4.evaluated_at(t, 1).with_name('xbm5'),
    pd3.evaluated_at(t, 0).with_name('xbmdot3'),
    pd3.evaluated_at(t, 1).with_name('xbmdot4'),
    pd4.evaluated_at(t, 0).with_name('xbmdot4'),
    pd4.evaluated_at(t, 1).with_name('xbmdot5'),
]
display(*equations)

$\displaystyle \boldsymbol{x}_{3} = \boldsymbol{a}_{3}$

$\displaystyle \boldsymbol{x}_{4} = \boldsymbol{a}_{3} + \boldsymbol{b}_{3} + \boldsymbol{c}_{3} + \boldsymbol{d}_{3}$

$\displaystyle \boldsymbol{x}_{4} = \boldsymbol{a}_{4}$

$\displaystyle \boldsymbol{x}_{5} = \boldsymbol{a}_{4} + \boldsymbol{b}_{4} + \boldsymbol{c}_{4} + \boldsymbol{d}_{4}$

$\displaystyle \dot{\boldsymbol{x}}_{3} = \boldsymbol{b}_{3}$

$\displaystyle \dot{\boldsymbol{x}}_{4} = \boldsymbol{b}_{3} + 2 \boldsymbol{c}_{3} + 3 \boldsymbol{d}_{3}$

$\displaystyle \dot{\boldsymbol{x}}_{4} = \boldsymbol{b}_{4}$

$\displaystyle \dot{\boldsymbol{x}}_{5} = \boldsymbol{b}_{4} + 2 \boldsymbol{c}_{4} + 3 \boldsymbol{d}_{4}$

We can solve the system of equations to get an expression for each coefficient:

coefficients = sp.solve(equations, coefficients3 + coefficients4)
for c, e in coefficients.items():
    display(NamedExpression(c, e))

$\displaystyle \boldsymbol{a}_{3} = \boldsymbol{x}_{3}$

$\displaystyle \boldsymbol{a}_{4} = \boldsymbol{x}_{4}$

$\displaystyle \boldsymbol{b}_{3} = \dot{\boldsymbol{x}}_{3}$

$\displaystyle \boldsymbol{b}_{4} = \dot{\boldsymbol{x}}_{4}$

$\displaystyle \boldsymbol{c}_{3} = - 3 \boldsymbol{x}_{3} + 3 \boldsymbol{x}_{4} - 2 \dot{\boldsymbol{x}}_{3} - \dot{\boldsymbol{x}}_{4}$

$\displaystyle \boldsymbol{c}_{4} = - 3 \boldsymbol{x}_{4} + 3 \boldsymbol{x}_{5} - 2 \dot{\boldsymbol{x}}_{4} - \dot{\boldsymbol{x}}_{5}$

$\displaystyle \boldsymbol{d}_{3} = 2 \boldsymbol{x}_{3} - 2 \boldsymbol{x}_{4} + \dot{\boldsymbol{x}}_{3} + \dot{\boldsymbol{x}}_{4}$

$\displaystyle \boldsymbol{d}_{4} = 2 \boldsymbol{x}_{4} - 2 \boldsymbol{x}_{5} + \dot{\boldsymbol{x}}_{4} + \dot{\boldsymbol{x}}_{5}$

So far, this is the same as we have done in the notebook about uniform Hermite splines. In fact, the above constants are the same as in \(M_H\)!

An additional constraint for natural splines is that the second derivatives are continuous, so let’s calculate those derivatives …

pdd3 = pd3.diff(t)
pdd4 = pd4.diff(t)
display(pdd3, pdd4)

$\displaystyle \frac{d^{2}}{d t^{2}} \boldsymbol{p}_{3} = 6 \boldsymbol{d}_{3} t + 2 \boldsymbol{c}_{3}$

$\displaystyle \frac{d^{2}}{d t^{2}} \boldsymbol{p}_{4} = 6 \boldsymbol{d}_{4} t + 2 \boldsymbol{c}_{4}$

… and set them to be equal at the segment border:

sp.Eq(pdd3.expr.subs(t, 1), pdd4.expr.subs(t, 0))

$\displaystyle 2 \boldsymbol{c}_{3} + 6 \boldsymbol{d}_{3} = 2 \boldsymbol{c}_{4}$

Inserting the equations from above leads to this equation:

_.subs(coefficients).simplify()

$\displaystyle 3 \boldsymbol{x}_{3} = 3 \boldsymbol{x}_{5} - \dot{\boldsymbol{x}}_{3} - 4 \dot{\boldsymbol{x}}_{4} - \dot{\boldsymbol{x}}_{5}$

We can generalize this expression by renaming index \(4\) to \(i\):

\begin{equation*} \dot{\boldsymbol{x}}_{i-1} + 4 \dot{\boldsymbol{x}}_{i} + \dot{\boldsymbol{x}}_{i+1} = 3 (\boldsymbol{x}_{i+1} - \boldsymbol{x}_{i-1}) \end{equation*}

This can be used for each segment – except for the very first and last one – yielding a matrix with \(N\) columns and \(N-2\) rows:

\begin{equation*} \left[ \begin{matrix} 1 & 4 & 1 && \cdots & 0 \\ & 1 & 4 & 1 && \vdots \\ && \ddots & \ddots && \\ \vdots && 1 & 4 & 1 & \\ 0 & \cdots && 1 & 4 & 1 \end{matrix} \right] \left[ \begin{matrix} \dot{\boldsymbol{x}}_0\\ \dot{\boldsymbol{x}}_1\\ \vdots\\ \dot{\boldsymbol{x}}_{N-2}\\ \dot{\boldsymbol{x}}_{N-1} \end{matrix} \right] = \left[ \begin{matrix} 3 (\boldsymbol{x}_2 - \boldsymbol{x}_0)\\ 3 (\boldsymbol{x}_3 - \boldsymbol{x}_1)\\ \vdots\\ 3 (\boldsymbol{x}_{N-2} - \boldsymbol{x}_{N-4})\\ 3 (\boldsymbol{x}_{N-1} - \boldsymbol{x}_{N-3}) \end{matrix} \right] \end{equation*}

End Conditions

We need a first and last row for this matrix to be able to fully define a natural spline. The following subsections show a selection of a few end conditions which can be used to obtain the missing rows of the matrix. End conditions (except “closed”) can be mixed, e.g. “clamped” at the beginning and “natural” at the end. The Python class splines.Natural uses “natural” end conditions by default.

Natural

Natural end conditions are commonly used for natural splines, which is probably why they are named that way.

There is a separate notebook about “natural” end conditions, from which we can get the uniform case by setting \(\Delta_i = 1\):

\begin{align*} 2 \dot{\boldsymbol{x}}_0 + \dot{\boldsymbol{x}}_1 &= 3 (\boldsymbol{x}_1 - \boldsymbol{x}_0) \\ \dot{\boldsymbol{x}}_{N-2} + 2 \dot{\boldsymbol{x}}_{N-1} &= 3 (\boldsymbol{x}_{N-1} - \boldsymbol{x}_{N-2}) \end{align*}

Adding this to the matrix from above leads to a full \(N \times N\) matrix:

\begin{equation*} \left[ \begin{matrix} 2 & 1 &&& \cdots & 0\\ 1 & 4 & 1 &&& \vdots \\ & 1 & 4 & 1 && \\ && \ddots & \ddots && \\ && 1 & 4 & 1 & \\ \vdots &&& 1 & 4 & 1\\ 0 & \cdots &&& 1 & 2 \end{matrix} \right] \left[ \begin{matrix} \dot{\boldsymbol{x}}_0\\ \dot{\boldsymbol{x}}_1\\ \vdots\\ \dot{\boldsymbol{x}}_{N-2}\\ \dot{\boldsymbol{x}}_{N-1} \end{matrix} \right] = \left[ \begin{matrix} 3 (\boldsymbol{x}_1 - \boldsymbol{x}_0)\\ 3 (\boldsymbol{x}_2 - \boldsymbol{x}_0)\\ 3 (\boldsymbol{x}_3 - \boldsymbol{x}_1)\\ \vdots\\ 3 (\boldsymbol{x}_{N-2} - \boldsymbol{x}_{N-4})\\ 3 (\boldsymbol{x}_{N-1} - \boldsymbol{x}_{N-3})\\ 3 (\boldsymbol{x}_{N-1} - \boldsymbol{x}_{N-2}) \end{matrix} \right] \end{equation*}

Clamped

We can simply provide arbitrarily chosen values \(D_\text{begin}\) and \(D_\text{end}\) for the end tangents. This is called clamped end conditions.

\begin{align*} \dot{\boldsymbol{x}}_0 &= D_\text{begin}\\ \dot{\boldsymbol{x}}_{N-1} &= D_\text{end} \end{align*}

This leads to a very simple first and last line:

\begin{equation*} \left[ \begin{matrix} 1 &&&& \cdots & 0\\ 1 & 4 & 1 &&& \vdots \\ & 1 & 4 & 1 && \\ && \ddots & \ddots && \\ && 1 & 4 & 1 & \\ \vdots &&& 1 & 4 & 1\\ 0 & \cdots &&&& 1 \end{matrix} \right] \left[ \begin{matrix} \dot{\boldsymbol{x}}_0\\ \dot{\boldsymbol{x}}_1\\ \vdots\\ \dot{\boldsymbol{x}}_{N-2}\\ \dot{\boldsymbol{x}}_{N-1} \end{matrix} \right] = \left[ \begin{matrix} D_\text{begin}\\ 3 (\boldsymbol{x}_2 - \boldsymbol{x}_0)\\ 3 (\boldsymbol{x}_3 - \boldsymbol{x}_1)\\ \vdots\\ 3 (\boldsymbol{x}_{N-2} - \boldsymbol{x}_{N-4})\\ 3 (\boldsymbol{x}_{N-1} - \boldsymbol{x}_{N-3})\\ D_\text{end} \end{matrix} \right] \end{equation*}

Closed

We can close the spline by connecting \(\boldsymbol{x}_{N-1}\) with \(\boldsymbol{x}_0\). This can be realized by cyclically extending the matrix in both directions:

\begin{equation*} \left[ \begin{matrix} 4 & 1 && \cdots & 0 & 1\\ 1 & 4 & 1 && 0 & 0 \\ & 1 & 4 & 1 && \vdots \\ && \ddots & \ddots && \\ \vdots && 1 & 4 & 1 & \\ 0 & 0 && 1 & 4 & 1\\ 1 & 0 & \cdots && 1 & 4 \end{matrix} \right] \left[ \begin{matrix} \dot{\boldsymbol{x}}_0\\ \dot{\boldsymbol{x}}_1\\ \vdots\\ \dot{\boldsymbol{x}}_{N-2}\\ \dot{\boldsymbol{x}}_{N-1} \end{matrix} \right] = \left[ \begin{matrix} 3 (\boldsymbol{x}_1 - \boldsymbol{x}_{N-1})\\ 3 (\boldsymbol{x}_2 - \boldsymbol{x}_0)\\ 3 (\boldsymbol{x}_3 - \boldsymbol{x}_1)\\ \vdots\\ 3 (\boldsymbol{x}_{N-2} - \boldsymbol{x}_{N-4})\\ 3 (\boldsymbol{x}_{N-1} - \boldsymbol{x}_{N-3})\\ 3 (\boldsymbol{x}_{0} - \boldsymbol{x}_{N-2}) \end{matrix} \right] \end{equation*}

Solving the System of Equations

The matrices above are tridiagonal and can therefore be solved efficiently with a tridiagonal matrix algorithm. The class splines.Natural, however, is not very concerned about efficiency and simply uses NumPy’s linalg.solve() function to solve the system of equations.