{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Physics 503: Scientific Computing\n",
    "### Lecture 15: Ordinary Differential Equations 2 - Vector equations\n",
    "### J.R. Gladden, University of Mississippi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### ODEs of Vector Quantities\n",
    "ODEs we often deal with ODEs involving vector quantities - forces, electromagnetic fields, fluid flow, ... In these cases, we typically break the vector quanitity up into orthogonal components and solve each equation seperately.  For example, a 2D force problem can be broken up into $x$ and $y$ components using Newton's 2nd law:\n",
    "$$ \\frac{d^2 \\vec{r(t)}}{dt^2}  =  \\vec{a}(t) = \\frac{1}{m} \\vec{F(t)} $$\n",
    "becomes\n",
    "$$ \\frac{d^2 x(t)}{dt^2} = \\frac{1}{m} F_x(t) $$ \n",
    "and\n",
    "$$ \\frac{d^2 y(t)}{dt^2} = \\frac{1}{m} F_y(t) $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Example 1: An oscillating system - the simple pendulum\n",
    "Many physical systems are bounded and oscillate about some equilibirum point.  A classic example is a pendulm.  A closed form solution is simple to obtain for small angle amplitudes (~<10 degrees), but larger angles require a numerical solution.  See the slides to setup the physics of the problem. The ODE we are solving here is 1 dimensional since it is a constrained system.  The ODE is:\n",
    "\n",
    "$$ \\frac{d^2 \\phi}{dt^2} = -\\frac{g}{L} \\sin{\\phi} $$\n",
    "where is the angle is bounded to small values, so that $\\sin{\\phi} \\simeq \\phi$, an closed (not exact) solution can be found: $\\phi(t) = \\phi_0 \\sin(\\sqrt{ \\frac{g}{L} } t) $ which is known as the small angle approximation.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Imports "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib wx\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "source": [
    "#### Problem setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "\n",
    "#Choose a method\n",
    "#method = 'euler'\n",
    "method = 'verlet'\n",
    "\n",
    " #Initial conditions\n",
    "thetai=5.*np.pi/180.\n",
    "vi=0.0\n",
    "\n",
    " #Pendulum Physics\n",
    "g=9.81\n",
    "L=0.3\n",
    "w=np.sqrt(g/L)\n",
    "T = 2*np.pi/w\n",
    "\n",
    " #Simulation parameters\n",
    "tau=T/1000.  # Base the time step on a characteristic time scale of the system.\n",
    "maxstep=10000\n",
    "\n",
    " #Convert to angular acceleration\n",
    "g=g/L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "#### Define the Euler and Verlet solvers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def euler(r,v,tau):\n",
    "    for k in range(maxstep-1):\n",
    "        vc=v[k]\n",
    "        rc=r[k]\n",
    "        t=k*tau\n",
    "        times[k]=t\n",
    "         #Theory with no drag for comaprison\n",
    "        theta_saa[k]=thetai*np.cos(w*t)\n",
    "         \n",
    "         #Compute accel. vector at each time k\n",
    "        acc=-g*np.sin(rc)\n",
    "         \n",
    "         #Compute r and v for next time k\n",
    "        v[k+1] = vc + tau* acc\n",
    "        r[k+1] = rc+tau*vc\n",
    "     \n",
    "    times[-1]=maxstep*tau \n",
    "    return r,v\n",
    " \n",
    "def verlet(r,v,tau):\n",
    "    v[1] = v[0] - g*np.sin(r[0])*tau\n",
    "    r[1]=r[0]-tau*v[1]\n",
    "    for k in range(1,maxstep-1):\n",
    "        t=k*tau\n",
    "        times[k]=t\n",
    "         #Theory with no drag for comaprison\n",
    "        theta_saa[k]=thetai*np.cos(w*t)\n",
    "         \n",
    "        acc=-g*np.sin(r[k])      \n",
    "        r[k+1] = 2*r[k] - r[k-1] +tau**2*acc        \n",
    "        v[k+1] = (r[k+1] - r[k-1])/(2*tau)\n",
    "        \n",
    "    times[-1]=maxstep*tau\n",
    "    return r,v\n",
    "\n",
    "# This is a pretty crude method of computing the period. A better method would be to do a Fourier Transform.\n",
    "def getPeriod(amps):\n",
    "    periods=[]\n",
    "    for amp in amps:\n",
    "        theta=np.zeros(maxstep)\n",
    "        v=np.zeros(maxstep)\n",
    "        theta[0]=amp*np.pi/180\n",
    "        theta,v = verlet(theta,v,tau)\n",
    "        offset = 100\n",
    "        period = (theta[offset:int(T/tau*1.5)].argmax()+offset)*tau\n",
    "        plt.plot(times[offset:int(T/tau*1.5)],theta[offset:int(T/tau*1.5)])\n",
    "        periods.append(period)\n",
    "    periods = np.array(periods)\n",
    "    return amps, periods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "#### Run the solvers, compare with small angle approximation (SAA), and plot solutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Create container 2-D arrays for position, velocity and theory position\n",
    "plt.close('all')\n",
    "theta=np.zeros(maxstep)\n",
    "v=np.zeros(maxstep)\n",
    "theta_saa=np.zeros(maxstep)\n",
    " \n",
    "times=np.zeros(maxstep)\n",
    "theta[0]=thetai\n",
    "v[0]=vi\n",
    " \n",
    "if method == 'euler': theta,v=euler(theta,v,tau)\n",
    "    \n",
    "elif method == 'verlet': theta,v = verlet(theta,v,tau)\n",
    "else: \n",
    "    print \"Method needs to be 'euler' or 'verlet'.\"\n",
    "    import sys\n",
    "    sys.exit()\n",
    "\n",
    "plt.subplot(221)\n",
    "theta=theta*180/np.pi\n",
    "theta_saa=theta_saa*180/np.pi\n",
    "plt.plot(times,theta,'b-',label='Numeric')\n",
    "plt.plot(times,theta_saa,'g--',label='Small Angle Approx.')\n",
    "plt.xlabel(\"Time (Sec)\")\n",
    "plt.ylabel(\"Theta\")\n",
    "plt.legend()\n",
    "plt.title(\"Pendulum Using %s Method\\nTime step: %2.3e seconds\" % (method, tau))\n",
    "\n",
    "plt.subplot(222)\n",
    "     \n",
    "#Exercise\n",
    "\n",
    "amps = np.arange(5,90,10)\n",
    "\n",
    "amps,periods = getPeriod(amps)\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Amplitude\")\n",
    "plt.title(\"Period Shifting with Amplitude\")\n",
    "\n",
    "plt.subplot(223)\n",
    "plt.plot(amps,periods/T)\n",
    "plt.xlabel(\"Amplitude (degrees)\")\n",
    "plt.ylabel(\"Period Relative to SAA\")\n",
    "\n",
    "taus = np.logspace(-5,-3,20)\n",
    "maxtheta = np.zeros(len(taus))\n",
    "i=0\n",
    "for t in taus:\n",
    "    theta=np.zeros(maxstep)\n",
    "    v=np.zeros(maxstep)\n",
    "    theta[0]=thetai\n",
    "    v[0]=vi\n",
    "    theta,v=euler(theta,v,t)\n",
    "    maxtheta[i] = max(theta)\n",
    "    i+=1\n",
    "\n",
    "plt.subplot(224)\n",
    "plt.plot(taus,maxtheta,'bo')\n",
    "plt.xlabel(\"Time step (tau)\")\n",
    "plt.ylabel(\"Maximum Amplitude\")\n",
    "plt.title(\"Amplitude Growth with Euler\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Example 2: Projectile motion with air drag\n",
    "Here is an example computing the trajectory of a ball thrown in the air near the earth's surface which takes air drag into account using classical Newtonian mechanics. The drag force is proportional to the speed of the ball and is always in the opposite direction as the velocity.  Note how we accomplish that in the code... "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Usage and visualization\n",
    "Take a close look at this code. Note the way we are storing the position and velocity data in arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Breaking at step: 182\n"
     ]
    }
   ],
   "source": [
    "\n",
    "def euler(rc,vc,acc,tau):\n",
    "    v_new = vc + tau* acc\n",
    "   #r_new = rc+tau*vc\n",
    "    r_new = rc+tau*v_new\n",
    "    return r_new,v_new\n",
    "\n",
    "def leapfrog(rc,vc,v_old,acc,tau):\n",
    "    v_new = v_old + 2*tau*acc\n",
    "    r_new = rc + 2*tau*v_new\n",
    "    return r_new, v_new\n",
    "   \n",
    "plt.close('all')\n",
    "\n",
    "#Initial conditions\n",
    "# NOTE: All units in the MKS system  - meters, kilograms, seconds\n",
    "yi=0.\n",
    "xi=0.\n",
    "speed=10.\n",
    "theta=65.\n",
    "vx=speed*np.cos(theta*np.pi/180.)\n",
    "vy=speed*np.sin(theta*np.pi/180.)\n",
    "\n",
    "#Simulation parameters: maxtime = tau*maxstep\n",
    "tau=0.01\n",
    "maxstep=1000\n",
    "\n",
    "#Create container 2-D arrays for position, velocity and theory position\n",
    "#1 row for each time step, each row: x and y component of position\n",
    "#Format: r=[[x1,y1],[x2,y2],[x3,y3],...,[xn,yn]]\n",
    "r=np.zeros((maxstep,2))\n",
    "v=np.zeros((maxstep,2))\n",
    "r_nodrag=np.zeros((maxstep,2))\n",
    "\n",
    "times=np.zeros(maxstep)\n",
    "r[0,0]=xi\n",
    "r[0,1]=yi\n",
    "v[0,0]=vx\n",
    "v[0,1]=vy\n",
    "\n",
    "#Some air drag physics\n",
    "Cd=0.35\n",
    "area=9.3E-3\n",
    "g=9.81\n",
    "m=0.45\n",
    "\n",
    "#density of air in Kg/m^3\n",
    "rho= 1.225\n",
    "#rho =0.\n",
    "drag=-0.5*Cd*rho*area/m\n",
    "\n",
    "for step in range(maxstep-1):\n",
    "    vc=v[step]\n",
    "    rc=r[step]\n",
    "    t=step*tau\n",
    "    times[step]=t\n",
    "   #Theory with no drag for comaprison\n",
    "    r_nodrag[step,0]=xi + vx*t\n",
    "    r_nodrag[step,1]=yi + vy*t - 0.5*g*t**2\n",
    "   \n",
    "   #Compute accel. vector at each time step\n",
    "    acc=drag*np.linalg.norm(vc)*vc\n",
    "    acc[1]=acc[1]-g\n",
    "   \n",
    "   #Compute r and v for next time step\n",
    "    r[step+1],v[step+1] = euler(rc,vc,acc,tau)\n",
    "   #Check to see if current height (y) is below the ground\n",
    "    if rc[1] <0.:\n",
    "        print(\"Breaking at step: %i\" % step)\n",
    "        v=v[:step]\n",
    "        r=r[:step]\n",
    "        r_nodrag=r_nodrag[:step]\n",
    "        break\n",
    "\n",
    "\n",
    "ground=np.zeros(len(r),dtype=float)\n",
    "plt.plot(r[:,0],r[:,1],'b-')\n",
    "plt.plot(r_nodrag[:,0],r_nodrag[:,1],'g--')\n",
    "plt.plot(r[:,0],ground,'k-')\n",
    "plt.xlabel(\"X (m)\")\n",
    "plt.ylabel(\"Y (m)\")\n",
    "plt.legend((\"With drag\", \"No drag\"))\n",
    "plt.title(\"Projectile with Air Drag \\n Drag Coeff. is: %2.3f\" % abs(drag))\n",
    "plt.show()\n",
    "   \n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}