{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Physics 503: Scientific Computing\n",
    "### Lecture 14: Ordinary Differential Equations 1\n",
    "### J.R. Gladden, University of Mississippi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### ODE equations\n",
    "See accompanying slides which present the theory behind the methods used here to compute numerical approximations to ordinary differential equations.  The following is one of the simplest schemes using the Euler method for solving a ODE of the form: $$ \\frac{\\mathrm{d} u(t)}{\\mathrm{d}t} = f(u(t)$$.  Given an initital condition $u(0) = u_0$, the solution at time step $k+1$ is approximated by $$u_{k+1} = u_k + \\Delta t f(u_k)$$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib wx\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "def forwardeuler(f,u0,T,n):\n",
    "     dt=T/float(n)\n",
    "     t=np.linspace(0,T,n+1)\n",
    "     #Use arrays rather than lists - faster!\n",
    "     # ForwardEuler needs extra point at the end to compute values at t=T\n",
    "     u=np.zeros(n+1)\n",
    "     u[0]=u0\n",
    "     for k in range(n):\n",
    "         u[k+1]=u[k]+dt*f(u[k])\n",
    "     return u, t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Usage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def f(u):\n",
    "     return -b*u         #exponential decay\n",
    "     #return 0.5*u*(1.0-u/800.)  #logistic population growth\n",
    "\n",
    " #set up parameters\n",
    "nsteps=3000\n",
    "u0=5.0\n",
    "Tfinal=10.0\n",
    "b=2.0\n",
    "\n",
    " #compute solution and exact for comparison\n",
    "u,t = forwardeuler(f,u0,Tfinal,nsteps)\n",
    "uexact=u0*np.exp(-b*t)\n",
    "udiff = abs(uexact - u)/max(uexact)*100\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "plt.plot(t,u,'o',label='Euler')\n",
    "plt.plot(t,uexact,'-',lw=2,label='Exact')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('u(t)')\n",
    "plt.legend()\n",
    "plt.title('Euler solution with %i time steps'%nsteps)\n",
    "plt.twinx()\n",
    "plt.plot(t,udiff,'r-',label='Error')\n",
    "plt.ylabel('Error (%)')\n",
    "plt.legend(loc=4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Projectile Motion with air drag\n",
    "Classic mechanics problem based on Newton's Laws.  Here we use both Euler method and the Leapfrog method to solve for both the velocity and position as a function of time given a known force. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "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",
    "   "
   ]
  },
  {
   "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": 15,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Breaking at step: 577\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",
    "yi=0.\n",
    "xi=0.\n",
    "speed=35.\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "np.linalg.norm?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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
}