CE4 Path Integrals and Numerical Integration_Complete_student.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computational Exercise 4: Path Integrals and Numerical Integration\n",
    "\n",
    "You saw in Tutorial 4 that electric potential is path-independent: no matter how you integrate between two points, the potential difference will be the same as long as the points and electric field stay the same. For this exercise, we'll be using the same setup as in Part 3 of that Tutorial: an E-field of the form $\\mathbf{E} = y \\hat{x} + x \\hat{y}$. As you saw in that tutorial, this vector field looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "x = np.linspace(-2, 2, 10)\n",
    "y = np.linspace(-2, 2, 10)\n",
    "\n",
    "X, Y = np.meshgrid(x, y)\n",
    "Z = X+Y\n",
    "\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.quiver(X, Y, Y, X, Z)\n",
    "\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You also found that the potential difference between two points using two different paths was the same, as seen below. Since we're using python, we'll define these two points as $(1, 2)$ and $(4, 6)$, rather than working symbolically. If you're unconvinced, feel free to change these points and see what happens!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "p1 = (1, 2)\n",
    "p2 = (4, 6)\n",
    "\n",
    "###Path 1 (red) is a segment from x1 to x2 at y1 followed by a \n",
    "###   segment from y1 to y2 at x2\n",
    "\n",
    "path1 = -(p2[0]*p1[1] - p1[0]*p1[1] + p2[0]*p2[1] - p2[0]*p1[1])\n",
    "\n",
    "###Path 2 (blue) is a segment from y1 to y2 at x1 followed by a \n",
    "###   segment from x1 to x2 at y2\n",
    "\n",
    "path2 = -(p1[0]*p2[1] - p1[0]*p1[1] + p2[0]*p2[1] - p1[0]*p2[1])\n",
    "\n",
    "print(\"Path 1 potential difference is {}, path 2 potential \\\n",
    "difference is {}.\".format(path1, path2))\n",
    "\n",
    "#Make quiver plot for plotting\n",
    "x = np.linspace(p1[0]-1, p2[0]+1, 15)\n",
    "y = np.linspace(p1[1]-1, p2[1]+1, 15)\n",
    "\n",
    "X, Y = np.meshgrid(x, y)\n",
    "Z = X+Y\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.quiver(X, Y, Y, X, Z)\n",
    "\n",
    "#Make points for p1 and p2\n",
    "plt.scatter(*p1, c='black', zorder=10)\n",
    "plt.scatter(*p2, c='black', zorder=10)\n",
    "\n",
    "#Add path 1 arrows\n",
    "ax.arrow(*p1, p2[0]-p1[0], 0, color='r', shape='full', \n",
    "         length_includes_head=True, head_width=.1)\n",
    "ax.arrow(p2[0], p1[1], 0, p2[1]-p1[1], color='r', shape='full', \n",
    "         length_includes_head=True, head_width=.1)\n",
    "\n",
    "#Add path 2 arrows\n",
    "ax.arrow(*p1, 0, p2[1]-p1[1], color='b', shape='full', \n",
    "         length_includes_head=True, head_width=.1)\n",
    "ax.arrow(p1[0], p2[1], p2[0]-p1[0], 0, color='b', shape='full', \n",
    "         length_includes_head=True, head_width=.1)\n",
    "\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$')\n",
    "ax.set_xlim(p1[0]-1, p2[0]+1)\n",
    "ax.set_ylim(p1[1]-1, p2[1]+1)\n",
    "\n",
    "fig.tight_layout(pad=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This works well, but it doesn't tell us how path integrals really work. For this exercise, we'll be building a path integral function that lets you integrate over any path between two points. We're going to rely on a very simple method for integration, much like a Riemann sum, but with vectors. The basic parts of this are:\n",
    "* A function that returns the electric field vector at each point (this is given)\n",
    "* A function which returns the potential difference between two given points, using vector calculations\n",
    "* Code that iterates through a path and sums the differences\n",
    "\n",
    "First, we've included a function that gives you the electric field vector at any given point. You don't need to change this at all, but you're welcome to if you want to experiment with different electric fields.\n",
    "\n",
    "This function returns what is known as a <code>tuple</code>, which you'll be working with a lot here. A tuple is much like a list, as both are ordered lists of objects. The difference is tuples cannot be changed, which makes them great for coordinate points or vectors. You can see this difference as tuples are generated with parentheses (as in the cell below), while lists use square brackets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def E(pt):\n",
    "    \"\"\"This function returns the Electric field vector at a given\n",
    "        point, a tuple of the form (x, y)\"\"\"\n",
    "    Ex = pt[1]  # Ex = y\n",
    "    Ey = pt[0]  # Ey = x\n",
    "    return (Ex, Ey)\n",
    "\n",
    "print('Electric field at (0,1): ',E([0,1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# a) Calculate the potential difference\n",
    "\n",
    "Next, we need the function for finding potential differences. The basic scaffolding is laid out for you, but you will need to edit the lines to include correct calculations. The function should take as input start and end points and return the potential difference between those two points. \n",
    "\n",
    "As a reminder, Griffiths Equation 1.22 is $V(\\mathbf{b}) - V(\\mathbf{a}) = -\\int_{a}^{b} \\mathbf{E} \\cdot d\\mathbf{l}$\n",
    "\n",
    "Some functions and tips that may come in handy:\n",
    "* `np.subtract(a, b)` lets you subtract $b$ from $a$ with two arrays or vectors, element-wise. For example if $a = (4, 3)$ and $b = (2, 0)$, `np.subtract(a, b)` returns $(4-2,3-0)=(2, 3)$.\n",
    "* `np.dot(a, b)` returns the dot product of $a$ and $b$\n",
    "\n",
    "*Something to think about:* When calculating the E-field, should you calculate it at point $a$, point $b$, or find both field strengths and average them? How does that difference change things, and does it matter?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def V(a, b):\n",
    "    \"\"\"This function returns the potential difference between two\n",
    "    points a and b, which are both tuples\"\"\"\n",
    "    sep = 1 #The separation vector between pt a and b\n",
    "    e = 1 #The E-field\n",
    "    v = 1 #The potential difference between a and b\n",
    "    return v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# b) Numerically integrate the potential differences\n",
    "\n",
    "Now, to approximate the integral, we need to make the difference between points $a$ and $b$ very small and sum over each section with a for loop. If we were to make the distance between the points infinitely small (the resolution infinite), it would exactly equal the integral.\n",
    "\n",
    "_Things to try:_\n",
    "\n",
    "- You may want to try changing the resolution of this numerical integration by changing the number of steps in the integration `res`.\n",
    "    \n",
    "- You may also try integrating over a different path by altering the calculation for `y_arr`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "res = 5 #Try changing the resolution and see how and if it changes the results!\n",
    "\n",
    "#Make a nice curve with res points between p1 and p2 defined above\n",
    "x_arr  = np.linspace(p1[0], p2[0], res)\n",
    "a = (p2[1] - p1[1])/(np.cosh(p2[0]) - np.cosh(p1[0]))\n",
    "b = p1[1] - a*np.cosh(p1[0])\n",
    "y_arr = a*np.cosh(x_arr) + b\n",
    "\n",
    "v = 0 #Initialize the potential along the path\n",
    "\n",
    "###YOUR CODE TO SUM OVER THE LINE HERE\n",
    "### This code should calculate the full potential, v, along the path\n",
    "### by integrating over the potential differences between each points in your arrays\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.quiver(X, Y, Y, X, Z)\n",
    "\n",
    "plt.plot(x_arr, y_arr)\n",
    "\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$')\n",
    "ax.set_xlim(p1[0]-1, p2[0]+1)\n",
    "ax.set_ylim(p1[1]-1, p2[1]+1)\n",
    "    \n",
    "    \n",
    "print(\"The potential difference over the path is {:.2f}\".format(v))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# c) Analysis of your results\n",
    "In the cell below, compare your calculation for the potential difference between points $a$ and $b$ to your answer from Tutorial 4, part 3 (feel free to check the solutions if you didn't finish the tutorial). Do you find reasonable convergence? What resolution did you feel was sufficient?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}