doc/euclidean/natural-uniform.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "This notebook is part of https://github.com/AudioSceneDescriptionFormat/splines, see also https://splines.readthedocs.io/.\n",
    "\n",
    "[back to overview](natural.ipynb) -\n",
    "[properties](natural-properties.ipynb) -\n",
    "[derivation (non-uniform)](natural-non-uniform.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Uniform Natural Splines\n",
    "\n",
    "For deriving natural splines,\n",
    "we first look at the *uniform* case,\n",
    "which means that the parameter interval in each segment\n",
    "is chosen to be $1$.\n",
    "\n",
    "The more general case with arbitrary parameter intervals\n",
    "is derived in a separate\n",
    "[notebook about non-uniform natural splines](natural-non-uniform.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing(order='grevlex')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We import some helpers from [utility.py](utility.py):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utility import NamedExpression, dotproduct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = sp.symbols('t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get started, let's look at two neighboring segments:\n",
    "Let's say the fourth segment,\n",
    "from $\\boldsymbol{x}_3$ to $\\boldsymbol{x}_4$,\n",
    "defined by the polynomial $\\boldsymbol{p}_3$,\n",
    "and the fifth segment, from $\\boldsymbol{x}_4$ to $\\boldsymbol{x}_5$,\n",
    "defined by the polynomial $\\boldsymbol{p}_4$.\n",
    "In both cases, we use $0 \\le t \\le 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coefficients3 = sp.symbols('a:dbm3')[::-1]\n",
    "coefficients4 = sp.symbols('a:dbm4')[::-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We apply these coefficients to the [monomial basis](polynomials.ipynb) ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b_monomial = t**3, t**2, t, 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... to define the two polynomials ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p3 = NamedExpression('pbm3', dotproduct(b_monomial, coefficients3))\n",
    "p4 = NamedExpression('pbm4', dotproduct(b_monomial, coefficients4))\n",
    "display(p3, p4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and we calculate their first derivatives:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd3 = p3.diff(t)\n",
    "pd4 = p4.diff(t)\n",
    "display(pd3, pd4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this,\n",
    "we obtain 8 equations\n",
    "containing the 8 yet unknown coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "equations = [\n",
    "    p3.evaluated_at(t, 0).with_name('xbm3'),\n",
    "    p3.evaluated_at(t, 1).with_name('xbm4'),\n",
    "    p4.evaluated_at(t, 0).with_name('xbm4'),\n",
    "    p4.evaluated_at(t, 1).with_name('xbm5'),\n",
    "    pd3.evaluated_at(t, 0).with_name('xbmdot3'),\n",
    "    pd3.evaluated_at(t, 1).with_name('xbmdot4'),\n",
    "    pd4.evaluated_at(t, 0).with_name('xbmdot4'),\n",
    "    pd4.evaluated_at(t, 1).with_name('xbmdot5'),\n",
    "]\n",
    "display(*equations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can solve the system of equations\n",
    "to get an expression for each coefficient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coefficients = sp.solve(equations, coefficients3 + coefficients4)\n",
    "for c, e in coefficients.items():\n",
    "    display(NamedExpression(c, e))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, this is the same as we have done in\n",
    "[the notebook about uniform Hermite splines](hermite-uniform.ipynb).\n",
    "In fact, the above constants are the same as in $M_H$!\n",
    "\n",
    "An additional constraint for natural splines\n",
    "is that the second derivatives are continuous,\n",
    "so let's calculate those derivatives ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdd3 = pd3.diff(t)\n",
    "pdd4 = pd4.diff(t)\n",
    "display(pdd3, pdd4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and set them to be equal at the segment border:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.Eq(pdd3.expr.subs(t, 1), pdd4.expr.subs(t, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inserting the equations from above\n",
    "leads to this equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_.subs(coefficients).simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can generalize this expression by renaming index $4$ to $i$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "\\dot{\\boldsymbol{x}}_{i-1}\n",
    "+\n",
    "4 \\dot{\\boldsymbol{x}}_{i}\n",
    "+\n",
    "\\dot{\\boldsymbol{x}}_{i+1}\n",
    "=\n",
    "3 (\\boldsymbol{x}_{i+1} - \\boldsymbol{x}_{i-1})\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can be used for each segment\n",
    "-- except for the very first and last one --\n",
    "yielding a matrix with $N$ columns and $N-2$ rows:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "1 & 4 & 1 && \\cdots & 0 \\\\\n",
    "& 1 & 4 & 1 && \\vdots \\\\\n",
    "&& \\ddots & \\ddots && \\\\\n",
    "\\vdots && 1 & 4 & 1 & \\\\\n",
    "0 & \\cdots && 1 & 4 & 1\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "\\dot{\\boldsymbol{x}}_0\\\\\n",
    "\\dot{\\boldsymbol{x}}_1\\\\\n",
    "\\vdots\\\\\n",
    "\\dot{\\boldsymbol{x}}_{N-2}\\\\\n",
    "\\dot{\\boldsymbol{x}}_{N-1}\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "3 (\\boldsymbol{x}_2 - \\boldsymbol{x}_0)\\\\\n",
    "3 (\\boldsymbol{x}_3 - \\boldsymbol{x}_1)\\\\\n",
    "\\vdots\\\\\n",
    "3 (\\boldsymbol{x}_{N-2} - \\boldsymbol{x}_{N-4})\\\\\n",
    "3 (\\boldsymbol{x}_{N-1} - \\boldsymbol{x}_{N-3})\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## End Conditions\n",
    "\n",
    "We need a first and last row for this matrix\n",
    "to be able to fully define a natural spline.\n",
    "The following subsections show a selection of a few end conditions\n",
    "which can be used to obtain the missing rows of the matrix.\n",
    "End conditions (except \"closed\") can be mixed,\n",
    "e.g. \"clamped\" at the beginning and \"natural\" at the end.\n",
    "The Python class [splines.Natural](../python-module/splines.rst#splines.Natural)\n",
    "uses \"natural\" end conditions by default."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Natural\n",
    "\n",
    "Natural end conditions are commonly used for natural splines,\n",
    "which is probably why they are named that way.\n",
    "\n",
    "There is a\n",
    "[separate notebook about \"natural\" end conditions](end-conditions-natural.ipynb),\n",
    "from which we can get the uniform case\n",
    "by setting $\\Delta_i = 1$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "2 \\dot{\\boldsymbol{x}}_0\n",
    "+\n",
    "\\dot{\\boldsymbol{x}}_1\n",
    "&=\n",
    "3 (\\boldsymbol{x}_1 - \\boldsymbol{x}_0)\n",
    "\\\\\n",
    "\\dot{\\boldsymbol{x}}_{N-2}\n",
    "+\n",
    "2 \\dot{\\boldsymbol{x}}_{N-1}\n",
    "&=\n",
    "3 (\\boldsymbol{x}_{N-1} - \\boldsymbol{x}_{N-2})\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding this to the matrix from above leads to a full $N \\times N$ matrix:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "2 & 1 &&& \\cdots & 0\\\\\n",
    "1 & 4 & 1 &&& \\vdots \\\\\n",
    "& 1 & 4 & 1 && \\\\\n",
    "&& \\ddots & \\ddots && \\\\\n",
    "&& 1 & 4 & 1 & \\\\\n",
    "\\vdots &&& 1 & 4 & 1\\\\\n",
    "0 & \\cdots &&& 1 & 2\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "\\dot{\\boldsymbol{x}}_0\\\\\n",
    "\\dot{\\boldsymbol{x}}_1\\\\\n",
    "\\vdots\\\\\n",
    "\\dot{\\boldsymbol{x}}_{N-2}\\\\\n",
    "\\dot{\\boldsymbol{x}}_{N-1}\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "3 (\\boldsymbol{x}_1 - \\boldsymbol{x}_0)\\\\\n",
    "3 (\\boldsymbol{x}_2 - \\boldsymbol{x}_0)\\\\\n",
    "3 (\\boldsymbol{x}_3 - \\boldsymbol{x}_1)\\\\\n",
    "\\vdots\\\\\n",
    "3 (\\boldsymbol{x}_{N-2} - \\boldsymbol{x}_{N-4})\\\\\n",
    "3 (\\boldsymbol{x}_{N-1} - \\boldsymbol{x}_{N-3})\\\\\n",
    "3 (\\boldsymbol{x}_{N-1} - \\boldsymbol{x}_{N-2})\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clamped\n",
    "\n",
    "We can simply provide arbitrarily chosen values\n",
    "$D_\\text{begin}$ and $D_\\text{end}$\n",
    "for the end tangents.\n",
    "This is called *clamped* end conditions.\n",
    "\n",
    "\\begin{align*}\n",
    "\\dot{\\boldsymbol{x}}_0 &= D_\\text{begin}\\\\\n",
    "\\dot{\\boldsymbol{x}}_{N-1} &= D_\\text{end}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This leads to a very simple first and last line:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "1 &&&& \\cdots & 0\\\\\n",
    "1 & 4 & 1 &&& \\vdots \\\\\n",
    "& 1 & 4 & 1 && \\\\\n",
    "&& \\ddots & \\ddots && \\\\\n",
    "&& 1 & 4 & 1 & \\\\\n",
    "\\vdots &&& 1 & 4 & 1\\\\\n",
    "0 & \\cdots &&&& 1\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "\\dot{\\boldsymbol{x}}_0\\\\\n",
    "\\dot{\\boldsymbol{x}}_1\\\\\n",
    "\\vdots\\\\\n",
    "\\dot{\\boldsymbol{x}}_{N-2}\\\\\n",
    "\\dot{\\boldsymbol{x}}_{N-1}\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "D_\\text{begin}\\\\\n",
    "3 (\\boldsymbol{x}_2 - \\boldsymbol{x}_0)\\\\\n",
    "3 (\\boldsymbol{x}_3 - \\boldsymbol{x}_1)\\\\\n",
    "\\vdots\\\\\n",
    "3 (\\boldsymbol{x}_{N-2} - \\boldsymbol{x}_{N-4})\\\\\n",
    "3 (\\boldsymbol{x}_{N-1} - \\boldsymbol{x}_{N-3})\\\\\n",
    "D_\\text{end}\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Closed\n",
    "\n",
    "We can close the spline by connecting $\\boldsymbol{x}_{N-1}$\n",
    "with $\\boldsymbol{x}_0$.\n",
    "This can be realized by cyclically extending the matrix\n",
    "in both directions:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "4 & 1 && \\cdots & 0 & 1\\\\\n",
    "1 & 4 & 1 && 0 & 0 \\\\\n",
    "& 1 & 4 & 1 && \\vdots \\\\\n",
    "&& \\ddots & \\ddots && \\\\\n",
    "\\vdots && 1 & 4 & 1 & \\\\\n",
    "0 & 0 && 1 & 4 & 1\\\\\n",
    "1 & 0 & \\cdots && 1 & 4\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "\\dot{\\boldsymbol{x}}_0\\\\\n",
    "\\dot{\\boldsymbol{x}}_1\\\\\n",
    "\\vdots\\\\\n",
    "\\dot{\\boldsymbol{x}}_{N-2}\\\\\n",
    "\\dot{\\boldsymbol{x}}_{N-1}\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "3 (\\boldsymbol{x}_1 - \\boldsymbol{x}_{N-1})\\\\\n",
    "3 (\\boldsymbol{x}_2 - \\boldsymbol{x}_0)\\\\\n",
    "3 (\\boldsymbol{x}_3 - \\boldsymbol{x}_1)\\\\\n",
    "\\vdots\\\\\n",
    "3 (\\boldsymbol{x}_{N-2} - \\boldsymbol{x}_{N-4})\\\\\n",
    "3 (\\boldsymbol{x}_{N-1} - \\boldsymbol{x}_{N-3})\\\\\n",
    "3 (\\boldsymbol{x}_{0} - \\boldsymbol{x}_{N-2})\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving the System of Equations\n",
    "\n",
    "The matrices above are *tridiagonal* and can therefore\n",
    "be solved efficiently with a\n",
    "[tridiagonal matrix algorithm](https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm).\n",
    "The class [splines.Natural](../python-module/splines.rst#splines.Natural),\n",
    "however,\n",
    "is not very concerned about efficiency and simply uses NumPy's\n",
    "[linalg.solve()](https://numpy.org/doc/stable/reference/generated/numpy.linalg.solve.html)\n",
    "function to solve the system of equations."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}