{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "92b5fa74",
   "metadata": {},
   "source": [
    "# Derivatives\n",
    "\n",
    " Börge Göbel, Jan 2022"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0711bcd",
   "metadata": {},
   "source": [
    "## 1. Mathematical definition of a derivative"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27ddc9e7",
   "metadata": {},
   "source": [
    "![Derivative](figure_04_derivatives.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ad10ced",
   "metadata": {},
   "source": [
    "\\\\( f'(x)=\\lim_{h\\rightarrow 0}\\frac{f(x+h)-f(x)}{h}\\\\)\n",
    "\n",
    "\\\\( f'(x)=\\lim_{h\\rightarrow 0}\\frac{f(x)-f(x-h)}{h}\\\\)\n",
    "\n",
    "\\\\( f'(x)=\\lim_{h\\rightarrow 0}\\frac{f(x+h)-f(x-h)}{2h}\\\\)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f91acf34",
   "metadata": {},
   "source": [
    "All three definitions are equivalent in case of a continuous function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f64d0ed",
   "metadata": {},
   "source": [
    "## 2. Numerical implementation of first-order derivatives"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d1a296b",
   "metadata": {},
   "source": [
    "Forward differences\n",
    "\\\\( f'(x_n)\\approx\\frac{f(x_{n+1})-f(x_n)}{x_{n+1}-x_n}\\\\)\n",
    "\n",
    "Backward differences\n",
    "\\\\( f'(x)\\approx\\frac{f(x_n)-f(x_{n-1})}{x_n-x_{n-1}}\\\\)\n",
    "\n",
    "Central differences\n",
    "\\\\( f'(x)\\approx\\frac{f(x_{n+1})-f(x_{n-1})}{x_{n+1}-x_{n-1}}\\\\)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6c934ba",
   "metadata": {},
   "source": [
    "### Example function: \\\\( f(x)=\\sin(x)x-\\frac{1}{100}x^3 \\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d336e105",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db71be64",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "51dcb403",
   "metadata": {},
   "source": [
    "### Why is central differences (typically) better than forward and backward differences?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "186b1b11",
   "metadata": {},
   "source": [
    "Forward differences\n",
    "\\\\( f'(x)=\\frac{f(x+h)-f(x)}{h}+\\mathcal{O}(h)\\\\)\n",
    "\n",
    "Backward differences\n",
    "\\\\( f'(x)=\\frac{f(x)-f(x-h)}{h}+\\mathcal{O}(h)\\\\)\n",
    "\n",
    "Central differences\n",
    "\\\\( f'(x)=\\frac{f(x+h)-f(x-h)}{2h}+\\mathcal{O}(h^2)\\\\)\n",
    "\n",
    "- \\\\(\\mathcal{O}(h^n)\\\\) means that the error is proportional to h^n.\n",
    "- Since \\\\(h\\\\) is small, the central differences method is more accurate. Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0b27259",
   "metadata": {},
   "source": [
    "Taylor expansion: \n",
    "\n",
    "\\\\(f(x+h)=f(x)+f'(x)h+\\frac{1}{2}f''(x)h^2+\\frac{1}{6}f'''(x)h^3+\\dots\\\\)\n",
    "\n",
    "\\\\(f(x-h)=f(x)-f'(x)h+\\frac{1}{2}f''(x)h^2-\\frac{1}{6}f'''(x)h^3\\pm\\dots\\\\)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9eadb72",
   "metadata": {},
   "source": [
    "- From the first and second line we can imediately see the  \\\\(\\mathcal{O}(h)\\\\) dependence of the forward and backward differences methods"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f04d750",
   "metadata": {},
   "source": [
    "\\\\(f'(x)=\\frac{1}{h}\\left[f(x+h)-f(x)-\\frac{1}{2}f''(x)h^2-\\frac{1}{6}f'''(x)h^3+\\dots\\right]=\\frac{f(x+h)-f(x)}{h}-\\frac{1}{2}f''(x)h-\\frac{1}{6}f'''(x)h^2+\\dots\\\\)\n",
    "\n",
    "\\\\(f'(x)=\\frac{1}{h}\\left[f(x)-f(x-h)+\\frac{1}{2}f''(x)h^2-\\frac{1}{6}f'''(x)h^3\\pm\\dots\\right]=\\frac{f(x)-f(x-h)}{h}+\\frac{1}{2}f''(x)h-\\frac{1}{6}f'''(x)h^2\\pm\\dots\\\\)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2cfa98f",
   "metadata": {},
   "source": [
    "- To find the \\\\(\\mathcal{O}(h^2)\\\\) dependence of the central differences method, we have to subtract the two terms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "547ef3fd",
   "metadata": {},
   "source": [
    "\\\\(f(x+h)-f(x-h)=2f'(x)h+\\frac{1}{3}f'''(x)h^3+\\dots\\\\)\n",
    "\n",
    "\\\\(\\frac{1}{2h}\\left[f(x+h)-f(x-h)\\right]=f'(x)+\\frac{1}{6}f'''(x)h^2+\\dots\\\\)\n",
    "\n",
    "\\\\( f'(x)=\\frac{f(x+h)-f(x-h)}{2h}+\\mathcal{O}(h^2)\\\\)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c7bd5d2",
   "metadata": {},
   "source": [
    "### Higher accuracy:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8548fd7",
   "metadata": {},
   "source": [
    "Richardson: \\\\(f'(x)=\\frac{1}{12h}\\left[f(x-2h)-8f(x-h)+8f(x+h)-f(x+2h)\\right]+\\mathcal{O}(h^4)\\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4c9a19b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "6b57d06c",
   "metadata": {},
   "source": [
    "### Even higher accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01d07a1a",
   "metadata": {},
   "source": [
    "Iteration formula:\n",
    "\n",
    "\\\\(D_{n+1}=\\frac{2^{2n}D_n(h)-D_n(2h)}{2^{2n}-1}\\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "310265c4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "480468c9",
   "metadata": {},
   "source": [
    "Calculate f'(x) for \\\\( f(x)=\\sin(x)x-\\frac{1}{100}x^3 \\\\) at \\\\(x = 3\\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e50a3b7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c21fa72e",
   "metadata": {},
   "source": [
    "Comparison to D1Richardson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81bd5c3d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "8eee3925",
   "metadata": {},
   "source": [
    "Analytical result: \\\\( f'(3)=3\\cos(3)+\\sin(3)-\\frac{3}{100}\\cdot 3^2 \\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "957819f6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "89a93179",
   "metadata": {},
   "source": [
    "## 3. Second derivatives"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaa2c1de",
   "metadata": {},
   "source": [
    "We derive \\\\(f'(x)\\\\) another time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "048cdc0b",
   "metadata": {},
   "source": [
    "\\\\( f''(x)=\\lim_{h\\rightarrow 0}\\frac{f'(x+h)-f'(x)}{h}\\\\)\n",
    "\n",
    "\\\\( f''(x)=\\lim_{h\\rightarrow 0}\\frac{f'(x)-f'(x-h)}{h}\\\\)\n",
    "\n",
    "\\\\( f''(x)=\\lim_{h\\rightarrow 0}\\frac{f'(x+h)-f'(x-h)}{2h}\\\\)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "818a0d68",
   "metadata": {},
   "source": [
    "This gives us many possibilities for the definition of \\\\(f''(x)\\\\) based on \\\\(f(x)\\\\), e. g."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30c84d6c",
   "metadata": {},
   "source": [
    "- Double forward & double backward\n",
    "\n",
    "\\\\(f''(x)=\\lim_{h\\rightarrow 0}\\frac{\\left[f(x+2h)-f(x+h)\\right]-\\left[f(x+h)-f(x)\\right]}{h^2}=\\lim_{h\\rightarrow 0}\\frac{f(x+2h)-2f(x+h)+f(x)}{h^2}\\\\)\n",
    "\n",
    "\\\\(f''(x)=\\lim_{h\\rightarrow 0}\\frac{\\left[f(x)-f(x-h)\\right]-\\left[f(x-h)-f(x-2h)\\right]}{h^2}=\\lim_{h\\rightarrow 0}\\frac{f(x)-2f(x-h)+f(x-2h)}{h^2}\\\\)\n",
    "\n",
    "- Forward and backward\n",
    "\n",
    "\\\\(f''(x)=\\lim_{h\\rightarrow 0}\\frac{\\left[f(x+h)-f(x)\\right]-\\left[f(x)-f(x-h)\\right]}{h^2}=\\lim_{h\\rightarrow 0}\\frac{f(x+h)-2f(x)+f(x-h)}{h^2}\\\\)\n",
    "\n",
    "- Double central (same result as forward and backward for \\\\(2h=g\\\\))\n",
    "\n",
    "\\\\(f''(x)=\\lim_{h\\rightarrow 0}\\frac{\\left[f(x+2h)-f(x)\\right]-\\left[f(x)-f(x-2h)\\right]}{(2h)^2}=\\lim_{h\\rightarrow 0}\\frac{f(x+2h)-2f(x)+f(x-2h)}{4h^2}=\\lim_{g\\rightarrow 0}\\frac{f(x+g)-2f(x)+f(x-g)}{g^2}\\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4aac9a16",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "389a0ca6",
   "metadata": {},
   "source": [
    "### Higher accuracy:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fc006fd",
   "metadata": {},
   "source": [
    "Richardson: \\\\(f''(x)=\\frac{1}{12h^2}\\left[-f(x-2h)+16f(x-h)-30f(x)+16f(x+h)-f(x+2h)\\right]+\\mathcal{O}(h^4)\\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00cb7a9a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "58915d60",
   "metadata": {},
   "source": [
    "## 4. Gradient, Divergence & Curl"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ab706cf",
   "metadata": {},
   "source": [
    "Now we consider a multidimensional function which means, the function depends on multiple variables \n",
    "\n",
    "\\\\( f(x,y,z)\\\\)\n",
    "\n",
    "or it is a function that has multiple dimensions itself\n",
    "\n",
    "\\\\( \\vec{g}(x,y,z)=\\begin{pmatrix}g_x(x,y,z)\\\\g_y(x,y,z)\\\\g_z(x,y,z)\\end{pmatrix}\\\\)\n",
    "\n",
    "With the nabla operator \\\\( \\nabla = \\begin{pmatrix}\\frac{\\partial}{\\partial x}\\\\\\frac{\\partial}{\\partial y}\\\\\\frac{\\partial}{\\partial z}\\end{pmatrix}\\\\) we can calculate:\n",
    "\n",
    "- gradient \\\\( \\nabla f(x,y,z) = \\begin{pmatrix}\\frac{\\partial}{\\partial x}f(x,y,z)\\\\\\frac{\\partial}{\\partial y}f(x,y,z)\\\\\\frac{\\partial}{\\partial z}f(x,y,z)\\end{pmatrix}\\\\)\n",
    "\n",
    "- curl \\\\(\\nabla\\times \\vec{g}(x,y,z) = \\begin{pmatrix}\n",
    "\\frac{\\partial}{\\partial y}g_z(x,y,z) - \\frac{\\partial}{\\partial z}g_y(x,y,z)\\\\\n",
    "\\frac{\\partial}{\\partial z}g_x(x,y,z) - \\frac{\\partial}{\\partial x}g_z(x,y,z)\\\\\n",
    "\\frac{\\partial}{\\partial x}g_y(x,y,z) - \\frac{\\partial}{\\partial y}g_x(x,y,z)\\\\\n",
    "\\end{pmatrix}\\\\)\n",
    "\n",
    "- divergence \\\\(\\nabla\\cdot \\vec{g}(x,y,z) = \\frac{\\partial}{\\partial x}g_x(x,y,z)+\\frac{\\partial}{\\partial y}g_y(x,y,z)+\\frac{\\partial}{\\partial z}g_z(x,y,z)\\\\)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15291943",
   "metadata": {},
   "source": [
    "### Example \n",
    "\n",
    "\\\\(\n",
    "f(\\vec{r})=f(x,y,z) = \\exp(-x^2-y^4)\n",
    "\\\\)\n",
    "\n",
    "\\\\(\n",
    "\\vec{g}(\\vec{r})=\\frac{\\vec{r}}{r}=\\begin{pmatrix}x/\\sqrt{x^2+y^2+z^2}\\\\y/\\sqrt{x^2+y^2+z^2}\\\\z/\\sqrt{x^2+y^2+z^2}\\end{pmatrix}\n",
    "\\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cb65689",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(r):\n",
    "    return np.exp(-r[0]**2-r[1]**4)\n",
    "\n",
    "def g(r):\n",
    "    return r / np.linalg.norm(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba1b7773",
   "metadata": {},
   "outputs": [],
   "source": [
    "x3, y3 = np.meshgrid(np.linspace(-2,2,201),np.linspace(-2,2,201))\n",
    "z3 = f( np.array([ x3, y3 ]) )\n",
    "\n",
    "plotproj = plt.axes(projection='3d')\n",
    "plotproj.contour3D(x3,y3,z3,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d7a7995",
   "metadata": {},
   "outputs": [],
   "source": [
    "x3, y3, z3 = np.meshgrid(np.linspace(-2,2,11),np.linspace(-2,2,11),np.linspace(-2,2,11))\n",
    "values = g( np.array([ x3, y3, z3 ]) )\n",
    "\n",
    "arrowplot = plt.axes(projection='3d')\n",
    "arrowplot.axis(False)\n",
    "\n",
    "scale=7\n",
    "arrowplot.quiver(\n",
    "    x3, y3, z3,\n",
    "    values[0]*scale,values[1]*scale,values[2]*scale\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30f53c1c",
   "metadata": {},
   "source": [
    "### Gradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38debaa9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "68e78a58",
   "metadata": {},
   "source": [
    "- analytical solution \n",
    "\n",
    "\\\\( \\nabla f(x,y,z) = \\begin{pmatrix}-2x\\exp(-x^2-y^4)\\\\-4y^3\\exp(-x^2-y^4)\\\\0\\end{pmatrix} \\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "913eda64",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "982455eb",
   "metadata": {},
   "source": [
    "### Divergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed356fe0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3b2b9195",
   "metadata": {},
   "source": [
    "- analytical solution \n",
    "\n",
    "\\\\( \\nabla \\cdot \\vec{g}(\\vec{r}) = \\frac{2}{r} \\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8513209a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "cd90f9c9",
   "metadata": {},
   "source": [
    "### Curl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5a17aa9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "941807c3",
   "metadata": {},
   "source": [
    "- analytical solution \n",
    "\n",
    "\\\\( \\nabla \\times \\vec{g}(\\vec{r}) = 0 \\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cc55f33",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}