{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fdd8f5a3",
   "metadata": {},
   "source": [
    "This notebook uses Julia to do symbolic computation: *it must run with a Julia kernel!*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9c9b91b",
   "metadata": {},
   "source": [
    "Consider the problem of computing the slope of the tangent line to a point on a circle."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd956f59",
   "metadata": {},
   "source": [
    "# 1. Implicit Differentiation with SymPy.jl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "873a6cc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "using SymPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8452b87f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$x$"
      ],
      "text/plain": [
       "x"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = Sym(\"x\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8f6ff717",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$y$"
      ],
      "text/plain": [
       "y"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = SymFunction(\"y\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fe3e51ea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\frac{d}{d x} y{\\left(x \\right)}$"
      ],
      "text/plain": [
       "d       \n",
       "──(y(x))\n",
       "dx      "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dy = diff(y(x), x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7472ffdb-4968-4394-9aa0-a0f177fe52ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Sym{PyCall.PyObject}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(dy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4377e244",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$x^{2} + y^{2}{\\left(x \\right)} - 1$"
      ],
      "text/plain": [
       " 2    2       \n",
       "x  + y (x) - 1"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = x^2 + y(x)^2 - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "aa774319-2625-4bee-b023-149820e2f659",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Sym{PyCall.PyObject}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "122a1fe2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$2 x + 2 y{\\left(x \\right)} \\frac{d}{d x} y{\\left(x \\right)}$"
      ],
      "text/plain": [
       "             d       \n",
       "2⋅x + 2⋅y(x)⋅──(y(x))\n",
       "             dx      "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dc = diff(c, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a925a3e0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\left[\\begin{smallmatrix}- \\frac{x}{y{\\left(x \\right)}}\\end{smallmatrix}\\right]$"
      ],
      "text/plain": [
       "1-element Vector{Sym{PyCall.PyObject}}:\n",
       " -x/y(x)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dydx = solve(dc, dy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a8828f0e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Vector{Sym{PyObject}}\u001b[90m (alias for \u001b[39m\u001b[90mArray{Sym{PyCall.PyObject}, 1}\u001b[39m\u001b[90m)\u001b[39m"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(dydx)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2f7d3e5",
   "metadata": {},
   "source": [
    "From the output above, we see that ``dydx`` is an array.  Arrays in Julia start at one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "26a0ff2a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$- \\frac{x}{y{\\left(x \\right)}}$"
      ],
      "text/plain": [
       "-x  \n",
       "────\n",
       "y(x)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "symbolicslope = dydx[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31976b58",
   "metadata": {},
   "source": [
    "Let us now take a random point on the circle, which is generated by a random angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "65781c43",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.0739308808121364"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "angle = 2*pi*rand()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe1104d8",
   "metadata": {},
   "source": [
    "Because Julia is oriented towards numerical computing, the ``pi`` is a 64-bit floating point constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "1f3356c6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.9977118154165967, 0.06761015735907423)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(a, b) = (cos(angle), sin(angle))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c4a3b0b",
   "metadata": {},
   "source": [
    "We store the coordinates of the point in a dictionary, which is a convenient data structure for substitution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "261efe32",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}\\begin{cases}y{\\left(x \\right)} & \\text{=>} &0.0676102\\\\x & \\text{=>} &-0.997712\\\\\\end{cases}\\end{equation*}"
      ],
      "text/plain": [
       "Dict{Sym{PyCall.PyObject}, Float64} with 2 entries:\n",
       "  y(x) => 0.0676102\n",
       "  x    => -0.997712"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sd = Dict(x => a, y(x) => b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "427a49d9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$14.7568332095102$"
      ],
      "text/plain": [
       "14.7568332095102"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "slope = symbolicslope.subs(sd)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cc0b89c",
   "metadata": {},
   "source": [
    "We see that the ``subs`` executed a simultaneous substitution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "023f1ed0",
   "metadata": {},
   "source": [
    "Let us now do the slope computation via a Taylor series.  The ``series`` of ``SymPy`` works only for one variable, but it works well on purely symbolic expressions.  We declare new symbols ``X`` and ``Y`` for the variables in the circle equation and ``A`` and ``B`` for the coordinates of the point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "35b99ea1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(X, Y, A, B)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X, Y, A, B = Sym(\"X, Y, A, B\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "afc3f89e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$X^{2} + Y^{2} - 1$"
      ],
      "text/plain": [
       " 2    2    \n",
       "X  + Y  - 1"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nc = X^2 + Y^2  - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "82fca3c6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$-1 + Y^{2} + 2 A \\left(- A + X\\right) + A^{2} + O\\left(\\left(- A + X\\right)^{2}; X\\rightarrow A\\right)$"
      ],
      "text/plain": [
       "      2                   2    ⎛        2       ⎞\n",
       "-1 + Y  + 2⋅A⋅(-A + X) + A  + O⎝(-A + X) ; X → A⎠"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cx = SymPy.series(nc, X, A, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03da9838",
   "metadata": {},
   "source": [
    "To convert the ``cx`` into a symbolic expression on which we can apply the ``series`` again, we select the coefficients and make the expression ``px``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "30b2b406",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$2 A$"
      ],
      "text/plain": [
       "2⋅A"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cx.coeff(X-A,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "cd7a4f89",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$A^{2} + Y^{2} - 1$"
      ],
      "text/plain": [
       " 2    2    \n",
       "A  + Y  - 1"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cx.coeff(X, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "93eb2ac2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$A^{2} + 2 A \\left(- A + X\\right) + Y^{2} - 1$"
      ],
      "text/plain": [
       " 2                   2    \n",
       "A  + 2⋅A⋅(-A + X) + Y  - 1"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "px = cx.coeff(X, 0) + cx.coeff(X-A, 1)*(X-A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "e7d2f46c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$-1 + 2 B \\left(- B + Y\\right) + B^{2} + 2 A \\left(- A + X\\right) + A^{2} + O\\left(\\left(- B + Y\\right)^{2}; Y\\rightarrow B\\right)$"
      ],
      "text/plain": [
       "                     2                   2    ⎛        2       ⎞\n",
       "-1 + 2⋅B⋅(-B + Y) + B  + 2⋅A⋅(-A + X) + A  + O⎝(-B + Y) ; Y → B⎠"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cy = series(px, Y, B, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "4e57b255-b31e-42b9-b8f0-9cda52787dc5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$A^{2} + 2 A \\left(- A + X\\right) + B^{2} + 2 B \\left(- B + Y\\right) - 1$"
      ],
      "text/plain": [
       " 2                   2                   \n",
       "A  + 2⋅A⋅(-A + X) + B  + 2⋅B⋅(-B + Y) - 1"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "py = cy.coeff(Y, 0) + cy.coeff(Y-B, 1)*(Y-B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "2537262d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$- 1.99542363083319 X + 0.135220314718148 Y - 2.0$"
      ],
      "text/plain": [
       "-1.99542363083319⋅X + 0.135220314718148⋅Y - 2.0"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cY = py.subs(Dict(A=>a, B=>b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "81b73227",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$14.7568332095102$"
      ],
      "text/plain": [
       "14.7568332095102"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "- cY.coeff(X, 1)/cY.coeff(Y, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "67541133",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$14.7568332095102$"
      ],
      "text/plain": [
       "14.7568332095102"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "slope"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95e64f9c",
   "metadata": {},
   "source": [
    "# 2. Implicit Differentiation via Symbolics.jl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "de8e9df6",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Symbolics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "358d5a93",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$ \\begin{equation}\n",
       "\\left[\n",
       "\\begin{array}{c}\n",
       "x \\\\\n",
       "y\\left( x \\right) \\\\\n",
       "\\end{array}\n",
       "\\right]\n",
       "\\end{equation}\n",
       " $$"
      ],
      "text/plain": [
       "2-element Vector{Num}:\n",
       "    x\n",
       " y(x)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@variables x, y(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5aa6d61a",
   "metadata": {},
   "source": [
    "Also ``SymPy`` has a ``Differential`` and in order not to confuse, because we have done ``using SymPy`` earlier in this notebook, we have have to declare that we are using the ``Differential`` from ``Symbolics``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "838389bc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Differential(x)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D = Symbolics.Differential(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "7285b976",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$ \\begin{equation}\n",
       "\\frac{\\mathrm{d} y\\left( x \\right)}{\\mathrm{d}x}\n",
       "\\end{equation}\n",
       " $$"
      ],
      "text/plain": [
       "Differential(x)(y(x))"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dy = D(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "0fe0777f-3f3f-45ef-a67b-dccd241a4ffd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Num"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(dy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "9963d0c2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$ \\begin{equation}\n",
       "-1 + x^{2} + \\left( y\\left( x \\right) \\right)^{2}\n",
       "\\end{equation}\n",
       " $$"
      ],
      "text/plain": [
       "-1 + x^2 + y(x)^2"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = x^2 + y^2 - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "433470bd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$ \\begin{equation}\n",
       "\\frac{\\mathrm{d}}{\\mathrm{d}x} \\left( -1 + x^{2} + \\left( y\\left( x \\right) \\right)^{2} \\right)\n",
       "\\end{equation}\n",
       " $$"
      ],
      "text/plain": [
       "Differential(x)(-1 + x^2 + y(x)^2)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dc = D(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "bf158b5e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$ \\begin{equation}\n",
       "2 x + 2 \\frac{\\mathrm{d} y\\left( x \\right)}{\\mathrm{d}x} y\\left( x \\right)\n",
       "\\end{equation}\n",
       " $$"
      ],
      "text/plain": [
       "2x + 2Differential(x)(y(x))*y(x)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq = expand_derivatives(dc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "71a16762",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$ \\begin{equation}\n",
       "\\frac{ - x}{y\\left( x \\right)}\n",
       "\\end{equation}\n",
       " $$"
      ],
      "text/plain": [
       "(-x) / y(x)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SBslope = Symbolics.solve_for(eq, dy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a754c3da",
   "metadata": {},
   "source": [
    "We have to make a new dictionary with the coordinates for the points, mapping them to the Symbolics variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "c6d23b26",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dict{Num, Float64} with 2 entries:\n",
       "  y(x) => 0.0676102\n",
       "  x    => -0.997712"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nsd = Dict(x => a, y => b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "e6ebfe66",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$ \\begin{equation}\n",
       "14.757\n",
       "\\end{equation}\n",
       " $$"
      ],
      "text/plain": [
       "14.756833209510194"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "substitute(SBslope, nsd)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8d4f51d",
   "metadata": {},
   "source": [
    "Let us compare with the slope we computed earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "3b87a118",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$14.7568332095102$"
      ],
      "text/plain": [
       "14.7568332095102"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "slope"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.10.4",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}