{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# $\\text{sinc}(x)$ and Its Derivatives\n",
    "\\begin{align}\n",
    "\\text{sinc}(x) &= \\begin{cases} \n",
    "      \\frac{x^4}{120}-\\frac{x^2}{6}+1 & |x| \\leq \\epsilon \\\\\n",
    "      \\frac{\\sin(x)}{x} & \\text{otherwise}\n",
    "   \\end{cases}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "\n",
    "import numpy\n",
    "from sympy import *\n",
    "from sympy.plotting import plot\n",
    "\n",
    "x = Symbol(\"x\")\n",
    "x0 = Symbol(\"x_0\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def taylor_series(f, n=5):\n",
    "    dfs = [f]\n",
    "    for i in range(n):\n",
    "        dfs.append(dfs[i].diff())\n",
    "    return sum([\n",
    "        dfs[i].subs({\"x\": x0}) / factorial(i) * (x - x0)**i\n",
    "        for i in range(len(dfs))\n",
    "    ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\sin{\\left(x \\right)}}{x}$"
      ],
      "text/plain": [
       "sin(x)/x"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sinc = sin(x) / x\n",
    "display(sinc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{x^{4}}{120} - \\frac{x^{2}}{6} + 1$"
      ],
      "text/plain": [
       "x**4/120 - x**2/6 + 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sinc_ts = taylor_series(sinc)\n",
    "display(limit(sinc_ts, x0, 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9999998333333333\n",
      "0.999999833333325\n",
      "2.2204460492503130808472633361816406250000000000000e-16\n",
      "1.4901161193847656250000000000000000e-8\n",
      "1.220703125000000000000000e-4\n"
     ]
    }
   ],
   "source": [
    "x1 = 1e-3\n",
    "print(1.0 - x1**2 / 6)\n",
    "print(1.0 - x1**2 * (1 / 6 + x1**2 / 120))\n",
    "x1 = sys.float_info.epsilon\n",
    "print(N(x1, 50))\n",
    "print(f\"{N(sqrt(x1), 35):e}\")\n",
    "print(f\"{N(sqrt(sqrt(x1)), 25):e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\frac{d}{dx} \\text{sinc}(x)$\n",
    "\\begin{align}\n",
    "\\text{sinc}'(x) &= \\begin{cases} \n",
    "      0 & x = 0 \\\\\n",
    "      \\frac{x\\cos(x) - \\sin(x)}{x^2} & \\text{otherwise}\n",
    "   \\end{cases}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{x \\cos{\\left(x \\right)} - \\sin{\\left(x \\right)}}{x^{2}}$"
      ],
      "text/plain": [
       "(x*cos(x) - sin(x))/x**2"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dsinc = simplify(sinc.diff(x))\n",
    "display(dsinc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{x^{5}}{840} + \\frac{x^{3}}{30} - \\frac{x}{3}$"
      ],
      "text/plain": [
       "-x**5/840 + x**3/30 - x/3"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dsinc_ts = taylor_series(dsinc)\n",
    "display(limit(dsinc_ts, x0, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\frac{d^2}{dx^2} \\text{sinc}(x)$\n",
    "\\begin{align}\n",
    "\\text{sinc}''(x) &= \\begin{cases} \n",
    "      \\frac{-1}{3} & x = 0 \\\\\n",
    "      \\frac{-x^2\\sin(x) - 2x\\cos(x) + 2\\sin(x)}{x^3} & \\text{otherwise}\n",
    "   \\end{cases}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{- x^{2} \\sin{\\left(x \\right)} - 2 x \\cos{\\left(x \\right)} + 2 \\sin{\\left(x \\right)}}{x^{3}}$"
      ],
      "text/plain": [
       "(-x**2*sin(x) - 2*x*cos(x) + 2*sin(x))/x**3"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ddsinc = simplify(dsinc.diff(x))\n",
    "display(ddsinc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{x^{4}}{168} + \\frac{x^{2}}{10} - \\frac{1}{3}$"
      ],
      "text/plain": [
       "-x**4/168 + x**2/10 - 1/3"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ddsinc_ts = taylor_series(ddsinc)\n",
    "display(limit(ddsinc_ts, x0, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Derivatives of $sinc(\\|x\\|)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def norm(x):\n",
    "    return sqrt(sum(x**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\sin{\\left(\\sqrt{x^{2} + y^{2} + z^{2}} \\right)}}{\\sqrt{x^{2} + y^{2} + z^{2}}}$"
      ],
      "text/plain": [
       "sin(sqrt(x**2 + y**2 + z**2))/sqrt(x**2 + y**2 + z**2)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "X = numpy.array(symbols(\"x y z\"))\n",
    "sinc_normx = sinc.subs({\"x\": norm(X)})\n",
    "display(sinc_normx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient of $\\text{sinc}(\\|x\\|)$\n",
    "\\begin{align}\n",
    "\\nabla \\text{sinc}(\\|x\\|) &= \\frac{\\text{sinc}'(\\|x\\|)}{\\|x\\|}x \\\\\n",
    "\\frac{\\text{sinc}'(x)}{x} &= \\begin{cases}\n",
    "      -\\frac{x^4}{840} + \\frac{x^2}{30} - \\frac{1}{3}& |x| \\leq \\epsilon\\\\\n",
    "      \\frac{x\\cos(x)-\\sin(x)}{x^3} & \\text{otherwise}\n",
    "   \\end{cases}\n",
    "\\end{align}\n",
    "where $\\epsilon = 10^{-4}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dsinc_normx = numpy.array(sinc_normx.diff(X))\n",
    "my_dsinc_normx = dsinc.subs({\"x\": norm(X)}) / norm(X) * X\n",
    "for i in range(3):\n",
    "    assert (simplify((dsinc_normx - my_dsinc_normx)[i]) == 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{x \\cos{\\left(x \\right)} - \\sin{\\left(x \\right)}}{x^{3}}$"
      ],
      "text/plain": [
       "(x*cos(x) - sin(x))/x**3"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dsinc_x = dsinc / x\n",
    "display(dsinc_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{x^{4}}{840} + \\frac{x^{2}}{30} - \\frac{1}{3}$"
      ],
      "text/plain": [
       "-x**4/840 + x**2/30 - 1/3"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dsinc_x_ts = taylor_series(dsinc_x)\n",
    "display(limit(dsinc_x_ts, x0, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hessian of $\\text{sinc}(\\|x\\|)$\n",
    "\\begin{align}\n",
    "\\nabla^2 \\text{sinc}(\\|x\\|) &= \\left(\\frac{\\text{sinc}''(\\|x\\|)}{\\|x\\|^2} - \\frac{\\text{sinc}'(\\|x\\|)}{\\|x\\|^3}\\right)xx^T + \\frac{\\text{sinc}'(\\|x\\|)}{\\|x\\|}I \\\\\n",
    "\\frac{\\text{sinc}''(\\|x\\|)}{\\|x\\|^2} - \\frac{\\text{sinc}'(\\|x\\|)}{\\|x\\|^3} &= \\begin{cases} \n",
    "      \\frac{x^{4}}{7560}-\\frac{x^{2}}{210}+\\frac{1}{15} & |x| \\leq \\epsilon \\\\\n",
    "      \\frac{-x^{2} \\sin (x)-3 x \\cos (x)+3 \\sin (x)}{x^{5}} & \\text{otherwise}\n",
    "   \\end{cases}\n",
    "\\end{align}\n",
    "where $\\epsilon = 0.1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddsinc_normx = numpy.array(Matrix(Matrix(dsinc_normx).diff(X)))\n",
    "my_ddsinc_normx = ((ddsinc.subs({\"x\": norm(X)}) / norm(X)**2 -\n",
    "                    dsinc.subs({\"x\": norm(X)}) / norm(X)**3) *\n",
    "                   X.reshape(3, 1) @ X.reshape(1, 3) +\n",
    "                   dsinc.subs({\"x\": norm(X)}) / norm(X) * numpy.eye(3))\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        assert (simplify(ddsinc_normx[i, j][0] - my_ddsinc_normx[i, j]) == 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{- x^{2} \\sin{\\left(x \\right)} - 3 x \\cos{\\left(x \\right)} + 3 \\sin{\\left(x \\right)}}{x^{5}}$"
      ],
      "text/plain": [
       "(-x**2*sin(x) - 3*x*cos(x) + 3*sin(x))/x**5"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ddsinc_x = simplify(ddsinc.subs({\"x\": x}) / x**2 - dsinc.subs({\"x\": x}) / x**3)\n",
    "display(ddsinc_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{x^{4}}{7560} - \\frac{x^{2}}{210} + \\frac{1}{15}$"
      ],
      "text/plain": [
       "x**4/7560 - x**2/210 + 1/15"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ddsinc_x_ts = taylor_series(ddsinc_x)\n",
    "display(limit(ddsinc_x_ts, x0, 0.0))"
   ]
  }
 ],
 "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.5"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}