Lab.07-Central-Forces-and-Potential-Energy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('default')\n",
    "import scipy.integrate as integrate\n",
    "import sympy as sp\n",
    "sp.init_printing()\n",
    "from ipywidgets import interact\n",
    "from IPython.display import display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Central Forces and Potential Energy\n",
    "## Lab 7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Name:\n",
    "### Student ID:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Instructions:\n",
    "Skim through the entire lab before beginning.  Read the text carefully and complete the steps as indicated.  Submit your completed lab through the Brightspace dropbox.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our discussion the solar system, we explicitly used Newton's Law of Gravitation to calculate the force between two celestial bodies.  Since gravity is a **conservative** force, we also could have considered the gravitional potential energy $U_g$ and used\n",
    "\n",
    "$$ F_g = - \\frac{d}{dr}U_g$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this lab we'll start from a potential energy and use that to calculate the motion of two particles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The potential energy between two point particles with masses $m_1 = 1$ and $m_2 = 2$ is \n",
    "\n",
    "$$U(r)=-\\frac{1}{r},$$\n",
    "\n",
    "where $r$ is the distance between them.  \n",
    "\n",
    "This results in an attractive force acting along the line connecting the two particles (much like gravity). \n",
    "\n",
    "At time t=0, the positions of the particles are $r_1 = (-10,-1)$ and $r_2 = (10,1)$, and  their velocities are $v_1 = (2,0)$ and $v_2 = (-2,0)$.  I.e. the particles initially are almost on a collision-course. \n",
    "\n",
    "Writing out $F = m a$ for each particle in component form, use `odeint` to solve the resulting set of differential equations from t=0 to t=10.  Plot the trajectories of both particles on a single graph.  \n",
    "\n",
    "Note that in this problem, all motion is confined to the x-y plane.  Assume MKS units for all quantities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### a) Setup Problem with SymPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to determine the force on particle 1.  First, we write an expression for the potential energy as stated in the problem. We express the distance $r$ between the two particles in terms of their coordinates.\n",
    "\n",
    "Using the following SymPy symbols, define an expression of the potential energy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1, y1, x2, y2 = sp.symbols('x1, y1, x2, y2')\n",
    "\n",
    "\n",
    "U = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The force in the $x$ direction on a particle is given by the negative of the derivative of the energy with respect to the $x$ coordinate of the particle.  Let's find the $x$ component of the force on particle 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, symbolic evaluate the derivative of $U$ with respect to $x1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BLANK.diff( BLANK )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is almost all we need to define a function for the force in the $x$ direction. The force is the negative of this partial deriviative.\n",
    "\n",
    "We will want to evaluate this force numerically, so *lambdify* this force function so that it it works with NumPy arrays. This function, called `Fx` needs to take four arguments: x1, y1, x2, y2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Fx = sp.lambdify([ BLANK, BLANK, BLANK, BLANK], BLANK, 'numpy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do the same for the force in the $y$ direction: create a *lambdified* function called `Fy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Fy = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### b) Numerical solution\n",
    "\n",
    "Let us start solving this problem by entering the parameters of the system (fill in the missing values)\n",
    "\n",
    "    m1 = 1  \n",
    "    m2 = \n",
    "\n",
    "    x1i = -10  \n",
    "    y1i =  \n",
    "    x2i =  \n",
    "    y2i = 1\n",
    "\n",
    "    v1x = 2  \n",
    "    v1y = 0  \n",
    "    v2x =  \n",
    "    v2y =  \n",
    "\n",
    "    tf = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m1 = 1  \n",
    "m2 = \n",
    "\n",
    "x1i = -10  \n",
    "y1i =  \n",
    "x2i =  \n",
    "y2i = 1\n",
    "\n",
    "v1x = 2  \n",
    "v1y = 0  \n",
    "v2x =  \n",
    "v2y =  \n",
    "\n",
    "tf = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enter the initial conditions for both particles (fill in the blanks)  \n",
    "    \n",
    "    Q0 = [x1i, ?, ?, v1y, ?, ?, ?, v2y]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Q0 = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we define our differential equation function that we pass to odeint. We need 4 equations (x, y of particles 1 and 2). The force on particle 2 should be equal in magnitude, but opposite in direction to the force in particle 1.\n",
    "\n",
    "eq1x = Fx(x1, y1, x2, y2) / m1  \n",
    "eq1y = Fy(x1, y1, x2, y2) / m1  \n",
    "eq2x = ?  \n",
    "eq2y = ?  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To be able to use `integrate.odeint` we need to define the right-hand-side of our system of equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def RHS(Q, t):\n",
    "    x1, vx1, y1, vy1, x2, vx2, y2, vy2 = Q\n",
    "    eq1x = BLANK\n",
    "    eq1y = BLANK\n",
    "    eq2x = BLANK\n",
    "    eq2y = BLANK\n",
    "\n",
    "    return [vx1, eq1x, vy1, eq1y, vx2, eq2x, vy2, eq2y]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having neatly obtained all our equations, we now employ `integrate.odeint` to solve the system of equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 2000\n",
    "t = np.linspace(0, tf, N)\n",
    "P = integrate.odeint(BLANK, BLANK, BLANK)\n",
    "x1, vx1, y1, vy1, x2, vx2, y2, vy2 = P.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### c) Plots "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the trajectories of the particles, the path of particle 1 in black, that of particle 2 in red."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel('x (m)')\n",
    "plt.ylabel('y (m)')\n",
    "plt.xlim(-10, 10)\n",
    "plt.ylim(-1.05, 1.05)\n",
    "plt.plot(BLANK, BLANK, BLANK)\n",
    "plt.plot(BLANK, BLANK, BLANK)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using @interact, make an *animation* showing the interaction between these two particles over time. At any one time, show the position of the particle with a large circle and only the *previous* trajectory up to that point.  \n",
    "\n",
    "You can use a counter variable like `i` but the title should be the current time in the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact(i=(0, N, 20))\n",
    "def explore(BLANK=BLANK):\n",
    "    plt.xlabel('x (m)')\n",
    "    plt.ylabel('y (m)')\n",
    "    plt.xlim(-10, 10)\n",
    "    plt.ylim(-1.05, 1.05)\n",
    "    plt.plot(BLANK, BLANK, BLANK)\n",
    "    plt.plot(BLANK, BLANK, BLANK)\n",
    "    plt.plot(BLANK, BLANK, BLANK)\n",
    "    plt.plot(BLANK, BLANK, BLANK)\n",
    "    \n",
    "    plt.title(BLANK)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### d) Discussion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which particle is deflected more?  Why does this make sense?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer**:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the speed of each particle as a function of time.  Which particle is moving faster at t=10?  (Pick the appropriate commands.)\n",
    "\n",
    "plt.plot(t,x1p + y1p,'k-');  \n",
    "plt.plot(t,x2p + y2p,'r-');\n",
    "\n",
    "or\n",
    "\n",
    "plt.plot(t,np.sqrt(vx1\\*\\*2 + vy1\\*\\*2),'k-');  \n",
    "plt.plot(t,np.sqrt(vx2\\*\\*2 + vy2\\*\\*2),'r-');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel('t (s)')\n",
    "plt.ylabel('Speed (m/s)')\n",
    "plt.plot(t, BLANK, 'k-');\n",
    "plt.plot(t, BLANK, 'r-');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer**: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know energy should be conserved.  But the plot above shows that the kinetic energy increases during the *collision*.  How can this be?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer**:"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}