doc/euclidean/end-conditions-natural.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](end-conditions.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Natural End Conditions\n",
    "\n",
    "For the first and last segment, we assume that the inner tangent is known.\n",
    "To find the outer tangent\n",
    "according to *natural* end conditions,\n",
    "the second derivative is set to $0$\n",
    "at the beginning and end of the curve.\n",
    "\n",
    "We are looking only at the non-uniform case here,\n",
    "it's easy to get to the uniform case by setting $\\Delta_i = 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Natural end conditions are naturally a good fit for\n",
    "[natural splines](natural-uniform.ipynb#End-Conditions).\n",
    "And in case you were wondering,\n",
    "natural end conditions are sometimes also called \"relaxed\" end conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing(order='grevlex')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, we are getting some help from [utility.py](utility.py):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utility import NamedExpression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = sp.symbols('t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Begin\n",
    "\n",
    "We are starting with the first polynomial segment $\\boldsymbol{p}_0(t)$,\n",
    "with $t_0 \\le t \\le t_1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t0, t1 = sp.symbols('t:2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coefficients ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a0, b0, c0, d0 = sp.symbols('a:dbm0')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... multiplied with the [monomial basis](polynomials.ipynb)\n",
    "give us the uniform polynomial ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d0 * t**3 + c0 * t**2 + b0 * t + a0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... which we re-scale to the desired parameter range:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p0 = NamedExpression('pbm0', _.subs(t, (t - t0) / (t1 - t0)))\n",
    "p0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need the first derivative\n",
    "(a.k.a. velocity, a.k.a. tangent vector):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd0 = p0.diff(t)\n",
    "pd0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to the\n",
    "[notebook about non-uniform Hermite splines](hermite-non-uniform.ipynb),\n",
    "we are interested in the function values and first derivatives\n",
    "at the control points:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "\\boldsymbol{x}_0 &= \\boldsymbol{p}_0(t_0)\\\\\n",
    "\\boldsymbol{x}_1 &= \\boldsymbol{p}_0(t_1)\\\\\n",
    "\\boldsymbol{\\dot{x}}_0 &= \\boldsymbol{p}_0'(t_0)\\\\\n",
    "\\boldsymbol{\\dot{x}}_1 &= \\boldsymbol{p}_0'(t_1)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "equations_begin = [\n",
    "    p0.evaluated_at(t, t0).with_name('xbm0'),\n",
    "    p0.evaluated_at(t, t1).with_name('xbm1'),\n",
    "    pd0.evaluated_at(t, t0).with_name('xdotbm0'),\n",
    "    pd0.evaluated_at(t, t1).with_name('xdotbm1'),\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get simpler equations,\n",
    "we are substituting $\\Delta_0 = t_1 - t_0$.\n",
    "Note that this is only for display purposes,\n",
    "the calculations are still done with $t_i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_begin = [\n",
    "    (t0, 0),\n",
    "    (t1, sp.Symbol('Delta0')),\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for e in equations_begin:\n",
    "    display(e.subs(delta_begin))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coefficients_begin = sp.solve(equations_begin, [a0, b0, c0, d0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for c, e in coefficients_begin.items():\n",
    "    display(NamedExpression(c, e.subs(delta_begin)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second derivative (a.k.a. acceleration) ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdd0 = pd0.diff(t)\n",
    "pdd0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... at the beginning of the curve ($t = t_0$) ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdd0.evaluated_at(t, t0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... is set to zero ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.Eq(_.expr, 0).subs(coefficients_begin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... leading to an expression for the initial tangent vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xd0 = NamedExpression.solve(_, 'xdotbm0')\n",
    "xd0.subs(delta_begin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can also be written as\n",
    "\n",
    "\\begin{equation*}\n",
    "\\boldsymbol{\\dot{x}}_0 =\n",
    "\\frac{3 \\left(\\boldsymbol{x}_1 - \\boldsymbol{x}_0\\right)}{2 \\Delta_0} -\n",
    "\\frac{\\boldsymbol{\\dot{x}}_1}{2}.\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## End\n",
    "\n",
    "If a spline has $N$ vertices,\n",
    "it has $N-1$ polynomial segments\n",
    "and the last polynomial segment is\n",
    "$\\boldsymbol{p}_{N-2}(t)$, with $t_{N-2} \\le t \\le t_{N-1}$.\n",
    "To simplify the notation a bit,\n",
    "let's assume we have $N = 10$ vertices,\n",
    "which makes $\\boldsymbol{p}_8$ the last polynomial segment.\n",
    "The following steps are very similar\n",
    "to the above derivation of the start conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a8, b8, c8, d8 = sp.symbols('a:dbm8')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t8, t9 = sp.symbols('t8:10')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d8 * t**3 + c8 * t**2 + b8 * t + a8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p8 = NamedExpression('pbm8', _.subs(t, (t - t8) / (t9 - t8)))\n",
    "p8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd8 = p8.diff(t)\n",
    "pd8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "\\boldsymbol{x}_{N-2} &= \\boldsymbol{p}_{N-2}(t_{N-2})\\\\\n",
    "\\boldsymbol{x}_{N-1} &= \\boldsymbol{p}_{N-2}(t_{N-1})\\\\\n",
    "\\boldsymbol{\\dot{x}}_{N-2} &= \\boldsymbol{p}_{N-2}'(t_{N-2})\\\\\n",
    "\\boldsymbol{\\dot{x}}_{N-1} &= \\boldsymbol{p}_{N-2}'(t_{N-1})\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "equations_end = [\n",
    "    p8.evaluated_at(t, t8).with_name('xbm8'),\n",
    "    p8.evaluated_at(t, t9).with_name('xbm9'),\n",
    "    pd8.evaluated_at(t, t8).with_name('xdotbm8'),\n",
    "    pd8.evaluated_at(t, t9).with_name('xdotbm9'),\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define $\\Delta_8 = t_9 - t_8$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_end = [\n",
    "    (t8, 0),\n",
    "    (t9, sp.Symbol('Delta8')),\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for e in equations_end:\n",
    "    display(e.subs(delta_end))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coefficients_end = sp.solve(equations_end, [a8, b8, c8, d8])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for c, e in coefficients_end.items():\n",
    "    display(NamedExpression(c, e.subs(delta_end)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time,\n",
    "the second derivative ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdd8 = pd8.diff(t)\n",
    "pdd8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... *at the end* of the last segment ($t = t_9$) ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdd8.evaluated_at(t, t9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... is set to zero ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.Eq(_.expr, 0).subs(coefficients_end)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... leading to an expression for the final tangent vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xd9 = NamedExpression.solve(_, 'xdotbm9')\n",
    "xd9.subs(delta_end)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Luckily, that's symmetric to the result we got above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equation can be generalized to\n",
    "\n",
    "\\begin{equation*}\n",
    "\\boldsymbol{\\dot{x}}_{N-1} =\n",
    "\\frac{3 \\left(\\boldsymbol{x}_{N-1} - \\boldsymbol{x}_{N-2}\\right)}{2 \\Delta_{N-2}} -\n",
    "\\frac{\\boldsymbol{\\dot{x}}_{N-2}}{2}.\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "\n",
    "We are showing a one-dimensional example where 3 time/value pairs are given.\n",
    "The slope for the middle value is given, the begin and end slopes are calculated using the \"natural\" end conditions as calculated above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "values = 2, 2, 1\n",
    "times = 0, 4, 5\n",
    "slope = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are using a few helper functions from [helper.py](helper.py) for plotting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helper import plot_sympy, grid_lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0, x1 = sp.symbols('xbm0:2')\n",
    "x8, x9 = sp.symbols('xbm8:10')\n",
    "xd1 = sp.symbols('xdotbm1')\n",
    "xd8 = sp.symbols('xdotbm8')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "begin = p0.subs(coefficients_begin).subs_symbols(xd0).subs({\n",
    "    t0: times[0],\n",
    "    t1: times[1],\n",
    "    x0: values[0],\n",
    "    x1: values[1],\n",
    "    xd1: slope,\n",
    "}).with_name(r'p_\\text{begin}')\n",
    "end = p8.subs(coefficients_end).subs_symbols(xd9).subs({\n",
    "    t8: times[1],\n",
    "    t9: times[2],\n",
    "    x8: values[1],\n",
    "    x9: values[2],\n",
    "    xd8: slope,\n",
    "}).with_name(r'p_\\text{end}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_sympy(\n",
    "    (begin.expr, (t, times[0], times[1])),\n",
    "    (end.expr, (t, times[1], times[2])))\n",
    "grid_lines(times, [1, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "begin.diff(t).evaluated_at(t, times[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "end.diff(t).evaluated_at(t, times[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bézier Control Points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Up to now we have assumed\n",
    "that we know one of the tangent vectors\n",
    "and want to find the other tangent vector\n",
    "in order to construct a [Hermite spline](hermite.ipynb).\n",
    "What if we want to construct a [Bézier spline](bezier.ipynb) instead?\n",
    "\n",
    "If the inner Bézier control points\n",
    "$\\boldsymbol{\\tilde{x}}_1^{(-)}$ and\n",
    "$\\boldsymbol{\\tilde{x}}_{N-2}^{(+)}$\n",
    "are given,\n",
    "we can insert the equations for the tangent vectors from the\n",
    "[notebook about non-uniform Bézier splines](bezier-non-uniform.ipynb#Tangent-Vectors)\n",
    "into our tangent vector equations from above\n",
    "and solve them for the outer control points\n",
    "$\\boldsymbol{\\tilde{x}}_0^{(+)}$ and\n",
    "$\\boldsymbol{\\tilde{x}}_{N-1}^{(-)}$, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xtilde0, xtilde1 = sp.symbols('xtildebm0^(+) xtildebm1^(-)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NamedExpression.solve(xd0.subs({\n",
    "    xd0.name: 3 * (xtilde0 - x0) / (t1 - t0),\n",
    "    xd1: 3 * (x1 - xtilde1) / (t1 - t0),\n",
    "}), xtilde0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xtilde8, xtilde9 = sp.symbols('xtildebm8^(+) xtildebm9^(-)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NamedExpression.solve(xd9.subs({\n",
    "    xd8: 3 * (xtilde8 - x8) / (t9 - t8),\n",
    "    xd9.name: 3 * (x9 - xtilde9) / (t9 - t8),\n",
    "}), xtilde9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that all $\\Delta_i$ cancel each other out\n",
    "(as well as the inner vertices\n",
    "$\\boldsymbol{x}_1$ and\n",
    "$\\boldsymbol{x}_{N-2}$)\n",
    "and we get very simple equations for the \"natural\" end conditions:\n",
    "\n",
    "\\begin{align*}\n",
    "\\boldsymbol{\\tilde{x}}_0^{(+)} &=\n",
    "\\frac{\\boldsymbol{x}_0 + \\boldsymbol{\\tilde{x}}_1^{(-)}}{2} \\\\\n",
    "\\boldsymbol{\\tilde{x}}_{N-1}^{(-)} &=\n",
    "\\frac{\\boldsymbol{x}_{N-1} + \\boldsymbol{\\tilde{x}}_{N-2}^{(+)}}{2}\n",
    "\\end{align*}"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}