phys512/derivatives_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "77869faf",
   "metadata": {},
   "source": [
    "## Derivatives"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d468e8a",
   "metadata": {},
   "source": [
    "#### Second order finite differences\n",
    "\n",
    "To keep track of the error terms, write out the Taylor expansions up to the 4th order terms:\n",
    "\n",
    "$$f(x+\\Delta x) \\approx f(x) + \\Delta x {df\\over dx} + {(\\Delta x)^2\\over 2}{d^2f\\over dx^2} + {(\\Delta x)^3\\over 6}{d^3f\\over dx^3} + {(\\Delta x)^4\\over 24}{d^4f\\over dx^4}$$\n",
    "\n",
    "$$f(x-\\Delta x) \\approx f(x) - \\Delta x {df\\over dx} + {(\\Delta x)^2\\over 2}{d^2f\\over dx^2}-{(\\Delta x)^3\\over 6}{d^3f\\over dx^3} + {(\\Delta x)^4\\over 24}{d^4f\\over dx^4}$$\n",
    "\n",
    "Then add and subtract:\n",
    "\n",
    "$$f(x+\\Delta x) + f(x-\\Delta x) \\approx 2f(x) + (\\Delta x)^2{d^2f\\over dx^2}  + {(\\Delta x)^4\\over 12}{d^4f\\over dx^4}$$\n",
    "\n",
    "$$f(x+\\Delta x) - f(x-\\Delta x) \\approx 2\\Delta x {df\\over dx} + {(\\Delta x)^3\\over 3}{d^3f\\over dx^3}$$\n",
    "\n",
    "Therefore\n",
    "\n",
    "$${d^2f\\over dx^2}\\approx {f(x+\\Delta x) + f(x-\\Delta x) -2f(x)\\over (\\Delta x)^2} - {(\\Delta x)^2\\over 12}{d^4f\\over dx^4}$$\n",
    "\n",
    "$${df\\over dx} \\approx {f(x+\\Delta x) - f(x-\\Delta x)\\over 2\\Delta x} - {(\\Delta x)^2\\over 6}{d^3f\\over dx^3}$$\n",
    "\n",
    "For comparison, the first order derivative with the error term kept in is\n",
    "\n",
    "$${df\\over dx} \\approx {f(x+\\Delta x) - f(x)\\over \\Delta x} - {\\Delta x\\over 2}{d^2f\\over dx^2}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2da5f3eb",
   "metadata": {},
   "source": [
    "#### Optimal step size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "05710fd2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAHNCAYAAADli4RZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAADQ7klEQVR4nOzdd3hTZRvA4V9W9160BTooe7SlZciUvbeyNyI4UBCRoX6iIqLIUhRBLaCoCIogWwSRvaGsAqVQoNC9d9Mk5/sjEKgUaLAhLX3v68pFcs57znlOaZMn75RJkiQhCIIgCIIgGMjNHYAgCIIgCEJZIxIkQRAEQRCEfxEJkiAIgiAIwr+IBEkQBEEQBOFfRIIkCIIgCILwLyJBEgRBEARB+BeRIAmCIAiCIPyLSJAEQRAEQRD+RSRIgiAIgiAI/yISJEEQBEEQhH8RCZIgCIIgCMK/iARJECoonU7HRx99REBAACqVioCAAObOnUvt2rXR6XRGn2/p0qX4+PhQUFBggmiNs3LlSmQyGdeuXTN3KIIglFMiQRKECmrJkiW899579OvXj+XLl7Nw4UI+/fRTpk2bhlxu/FvDqFGjUKvVLFu2zATRCuaSnZ3NzJkz6dKlCy4uLshkMlauXGnusATB5GSSJEnmDkIQhCcvNDQUNzc3/vzzTwAWLVrEzJkzSUhIwMrK6rHOOW3aNNasWUN0dDQymaw0wzWKVqulsLAQS0tLs8bxNLh27Rr+/v74+PhQrVo1/vnnH1asWMGoUaPMHZogmJSoQRKECig/P5/Tp0/TunVrw7YVK1bQq1evx06OAAYMGMD169fZvXt3aYT52BQKBVZWViI5KgVeXl7ExcVx/fp1PvvsM3OHIwhPjEiQBKGCeeGFF7C2tkar1fLuu+8ik8nw8vLizJkzdOjQ4b7yt27dwsrKijFjxhTZvnPnTlQqFW+88YZhW2hoKC4uLvzxxx8mvYesrCwmTZqEn58flpaWeHh40LFjR06ePAkU3wfp/fffRyaTERUVxahRo3BycsLR0ZHRo0eTm5tb7HWMuXdTM1cslpaWeHp6muTcglCWiQRJECqYoUOHMn78eAA+//xzVq1axUsvvQRASEjIfeUrV67M2LFj+fHHH7l+/ToAFy9epH///nTt2pX58+cXKR8SEsKBAweKvXZhYSHJycklejyso/hLL73E119/zXPPPceSJUuYMmUK1tbWXLhw4ZH3P2DAALKyspgzZw4DBgxg5cqVfPDBB8WWNfbeTelxYimtn7cgVEiSIAgVzttvvy3Z2tpKWq1WkiRJevfddyVAysrKKrb8zZs3JUtLS+nll1+WkpOTpYCAACk4OFjKzs6+r+y4ceMka2vrYs+ze/duCSjRIzo6+oHxOzo6Sq+++uoD969YseK+c8ycOVMCpDFjxhQp27dvX8nV1fWB5zLm3k3N2FhK6+d9x7FjxyRAWrFiRenemCCUQconl4oJglBWnDlzhnr16hlGq6WkpKBUKrGzsyu2fOXKlXnxxRf59ttvOXnyJHl5eezZswdbW9v7yjo7O5OXl0dubi42NjZF9gUFBfHXX3+VKMaHNes4OTlx5MgRYmNj8fb2LtH57rhTW3ZHq1atWL9+PZmZmTg4ONxX3ph7NzVjYymtn7cgVEQiQRKECuj06dN07tzZqGOmTJnCl19+yZkzZ9i3bx+VK1cutpx0e2BscR2knZ2di+3nZKy5c+cycuRIqlatSmhoKN26dWPEiBFUq1btkcf6+PjcFxNAWlpasQkSlPzenwRjYimtn7cgVEQiQRKECiY9PZ2YmBgaNGhg2Obq6opGoyErKwt7e/tij5s9ezYAGo0GFxeXB54/LS0NGxsbrK2t79unVqtJTU0tUZzu7u4oFIpi9w0YMMBQ87Njxw4+++wzPv30U37//Xe6du360PM+6JzSQ2Y8Kem9PwnGxFJaP29BqIhEJ21BqGDOnDkDQGBgoGFb7dq1AYiOji72mM8++4zvvvuOL7/8EqVSafiQLk50dDR16tQpdt/Bgwfx8vIq0SMmJuah9+Hl5cUrr7zChg0biI6OxtXV9aFxPS5j7t3UjI2lNH/eglDRiBokQahgTp8+DRRNkJo1awbA8ePHi2wH2LBhA9OnT2fWrFm8+uqrXL58mSVLlvDOO+/g7+9/3/lPnjzJ0KFDi712afSJ0Wq1ZGdn4+joaNjm4eGBt7d3qS9zYuy9m9LjxCL6IAnCf2DuXuKCIDxZY8eOlSpXrnzf9vr160uDBw8usu348eOSjY2NNHz4cMO2W7duSZaWltILL7xw3zmOHz8uAdLOnTtLP/Db0tLSJFtbW2nkyJHSggULpG+++UYaMGCABEjz58+XJOnho9iSkpKKnK+4snfuxZh7B6Rnn3221O7zv8RS2hYvXizNmjVLevnllyVA6tevnzRr1ixp1qxZUnp6usmvLwjmIBIkQahgmjRpInXt2vW+7QsWLJDs7Oyk3NxcSZIkKSYmRvLy8pJatGgh5efnFyn78ssvSyqVSrp69WqR7dOmTZN8fHwknU5nsvgLCgqkt956SwoKCpLs7e0lW1tbKSgoSFqyZImhzH9NkIy996ysLAmQBg0aVLo3+xixmIKvr+9/mh5AEMojsRabIAgAZGRkUK1aNebOncsLL7xg9PEFBQX4+fkxffp0Jk6caIIIy66tW7fSo0cPTp8+XaTzuyAI5ZfopC0IAgCOjo5MnTqVzz777LFmVV6xYgUqleq+eYYqgt27dzNo0CCRHAnCU6RC1yBt3ryZN998E51Ox7Rp0xg7dqy5QxIEQRAEoQyosAmSRqOhbt267N69G0dHR0JDQzl48CCurq7mDk0QBEEQBDOrsE1sR48epV69elSuXBk7Ozu6du3Kjh07zB2WIAiCIAhlQLlNkPbu3UvPnj3x9vZGJpOxYcOG+8p89dVX+Pn5YWVlRdOmTTl69KhhX2xsbJEp+itXrsytW7eeROiCIAiCIJRx5TZBysnJISgoiK+++qrY/WvWrGHy5MnMnDmTkydPEhQUROfOnUlMTHzCkQqCIAiCUN6U25m0u3bt+tA1lxYsWMCLL77I6NGjAVi6dClbtmxh+fLlTJ8+HW9v7yI1Rrdu3aJJkyYPPF9BQUGRWXp1Oh2pqam4uroWuyinIAiCIAhljyRJZGVl4e3tjVz+kHoiM87BVGoAaf369YbXBQUFkkKhKLJNkiRpxIgRUq9evSRJkqTCwkKpevXq0s2bN6WsrCypZs2aUnJy8gOvcWeSOfEQD/EQD/EQD/Eo/4+YmJiH5hbltgbpYZKTk9FqtVSqVKnI9kqVKnHx4kUAlEol8+fPp23btuh0OqZOnfrQEWwzZsxg8uTJhtcZGRn4+PgQExODg4ODaW5EEARBEIRSlZmZSdWqVbG3t39ouacyQSqpXr160atXrxKVtbS0xNLS8r7tDg4OIkESBEEQhHLmUd1jym0n7Ydxc3NDoVCQkJBQZHtCQoJYsVoQBEEQhEd6KhMkCwsLQkND2bVrl2GbTqdj165dNGvWzIyRCYIgCIJQHpTbJrbs7GyioqIMr6OjowkPD8fFxQUfHx8mT57MyJEjadSoEU2aNGHRokXk5OQYRrUJgqCn0+lQq9XmDkMoB1QqFQqFwtxhCMITUW4TpOPHj9O2bVvD6zsdqEeOHMnKlSsZOHAgSUlJvPfee8THxxMcHMz27dvv67gtCBWZWq0mOjr6sRanFSomJycnPD09xfQmwlOvwq7F9l9lZmbi6OhIRkaG6KQtlEuSJHHjxg0KCwsfPR+IUOFJkkRubi6JiYk4OTnh5eVl7pAE4bGU9PO73NYgCYLw32g0GnJzc/H29sbGxsbc4QjlgLW1NQCJiYl4eHiI5jbhqSa+MgpCBaXVagH9oAZBKKk7yXRhYaGZIxEE0xIJkiBUcKIviWAM8fsiVBQiQRIEQRAEQfgXkSAJgiAIgiD8i0iQBEF4arRp04ZJkyaZOwyjXLt2DZlMRnh4uLlDEQThHiJBEgShXBk1ahQymey+R1RUFL///juzZs36T+eXyWRs2LChdIIVBKHcEsP8BUEod7p
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def f(x):\n",
    "    return np.sin(x)\n",
    "\n",
    "dx_vals = 10**np.linspace(-13, -1, 100)\n",
    "\n",
    "x = 1.0 # Evaluate the derivatives at this value of x\n",
    "\n",
    "err = np.zeros_like(dx_vals)   # to store the fractional error from the 1st order difference\n",
    "err2 = np.zeros_like(dx_vals)  # same for the second order\n",
    "\n",
    "for i, dx in enumerate(dx_vals):\n",
    "    dfdx1 = (f(x + dx) - f(x)) / dx\n",
    "    err[i] = np.abs((dfdx1 - np.cos(x))/np.cos(x))\n",
    "\n",
    "    dfdx2 = (f(x + dx) - f(x - dx)) / (2*dx)\n",
    "    err2[i] = np.abs((dfdx2 - np.cos(x))/np.cos(x))\n",
    "\n",
    "# Plot the results\n",
    "plt.plot(dx_vals, err, label=\"First order\")\n",
    "plt.plot(dx_vals, err2, label=\"Second order\")\n",
    "\n",
    "plt.ylim((1e-2*min(err2), 1.0))\n",
    "\n",
    "# and the analytic expressions\n",
    "plt.plot(dx_vals, 1e-16 * f(x) / dx_vals, \":\", label=r'$\\epsilon f(x)/\\Delta x$')\n",
    "plt.plot(dx_vals, dx_vals, \":\", label=r'$\\Delta x$')\n",
    "plt.plot(dx_vals, dx_vals**2, \":\", label=r'$(\\Delta x)^2$')\n",
    "\n",
    "plt.title(r'$f(x) = \\sin x, \\ \\ \\ x=%lg$' % (x,))\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.xlabel(r'$\\Delta x$')\n",
    "plt.ylabel('Fractional error in derivative')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "127d61b8",
   "metadata": {},
   "source": [
    "The optimal step size for the 1st order expression is $\\sim \\epsilon^{1/2}\\sim 10^{-8}$ and the associated error is also $\\sim \\epsilon^{1/2}$. \n",
    "\n",
    "For the second order expression, the finite difference error is $\\propto (\\Delta x)^2$ rather than $\\Delta x$. Following the same argument as for the first order derivative leads to an optimal step size of order $\\epsilon^{1/3}\\sim 10^{-16/3}\\sim 10^{-5.3}$, and an error $\\epsilon^{2/3}\\sim 10^{-10.6}$.\n",
    "\n",
    "Notice that a higher order derivative allows us to get a better error with a larger step size."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5969ebb",
   "metadata": {},
   "source": [
    "#### Automatic derivatives\n",
    "\n",
    "Inspecting the code, we can see that the operators that are implemented inside the class return a new instance of Var that is initialized with (1) the value obtained from the operation, and (2) a set of children that give the partial derivatives with respect to the variables involved in the operation. \n",
    "\n",
    "The multiplication operator is instructive to look at. For the operation `x*y`, there are two variables `x` and `y`, so we need to return (1) the operation value `x*y`, (2) the partial derivative with respect to `x`, which is $\\partial(xy)/\\partial x = y$, and (3) the partial derivative with respect to `y`, which is $\\partial(xy)/\\partial y = x$. The corresponding function definition is\n",
    "\n",
    "```\n",
    "def __mul__(self, other):\n",
    "        return Var(self.value * other.value, [(other.value, self), (self.value, other)])\n",
    "```\n",
    "\n",
    "For $exp(x)$, there is only one argument to the function so we just have to return one partial derivative $(\\partial/\\partial x)\\exp(x) = \\exp(x)$. So we can add \n",
    "```\n",
    "def exp(self):\n",
    "        return Var(math.exp(self.value), [(math.exp(self.value), self)])\n",
    "```\n",
    "\n",
    "which we use for example as \n",
    "\n",
    "`f= x * y + x.exp()`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cbe50bd",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}