Lecture3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 3, Solving Systems of ODE's and Newton's 2nd law\n",
    "\n",
    "\n",
    "**Overview:**\n",
    "* Reducing nth order ODE's to 1st order.\n",
    "* Solving Newton's 2nd law.\n",
    "    * A simple fix for the Euler method: Reading Euler-Cromer method Euler_Cromer.pdf\n",
    "* Intro to object oriented programing in Python, beginning to develop our \"particle\" class.\n",
    "\n",
    "**Next Lecture:** Developing a physical model and more object oriented practice, continuing development of our particle class.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this lecture workbook, we will begin to explore object oriented programming in python. We will do this by developing a class called `particle.py` that will be used throughout most of PHYS 1600/2600.  \n",
    "\n",
    "In this workbook, we will begin developing our class within the Jupyter Notebook. This is a good workflow for rapid testing and modifying our class, but is not efficient in the long-run. In the next lecture, we will move our class outside of the Jupyter notebook and into a separate .py file. We will then use an `import` statement inside of a separate notebook any time we would like to use our class. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tasks:\n",
    "\n",
    "**Understanding the class structure**\n",
    "\n",
    "---\n",
    "\n",
    "1. Create an instance of the Particle class called \"free_ball\". Explore the attributes of free_ball, do you understand what `__init__` is doing? Do you understand what the purpose of the argument `self` is ?\n",
    "2. Try calling the Euler_step method of free_ball a number of times. After each call to Euler_step use the results method to display the current position and velocity. Do you understand how this is working?\n",
    "2.  Create an instance of the FallingParticle subclass called \"fall_ball\". Explore the attributes of fall_ball, how do they compare with free_ball? Are any attributes the same? Are any different?\n",
    "3. Calculate the trajectory of a particle falling under gravity using the Euler method.\n",
    "4. Plot the trajectory.\n",
    "\n",
    "**Extending the class and exploring limitations of our numerical method**\n",
    "\n",
    "---\n",
    "\n",
    "5. Create a sub-class of Particle for a 1-Dimensional particle in a simple harmonic potential.\n",
    "6. Create an instance of your new simple harmonic potential subclass called SHO.\n",
    "6. Generate phase space and trajectory plots for SHO, comparing the exact and  Euler solutions. You might want to extend your previous plot method to do this. What do you notice?\n",
    "6. In the SHO sub-class, include a method to calculate the total energy of the SHO, plot the relative error in energy for the Euler solution.\n",
    "7. Reading: Newman 8.1 - 8.4 and Euler_Cromer.pdf.\n",
    " \n",
    " \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The Particle class is defined below**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt # for plotting          \n",
    "import numpy as np\n",
    "%matplotlib notebook  \n",
    "\n",
    "class Particle (object):\n",
    "\n",
    "    \"\"\"Class that describes particle\"\"\"\n",
    "    def __init__(self, m = 1.0, y0=1.0, v0=0.0,  tf = 10.0, dt = 0.01):\n",
    "        '''\n",
    "        This is the initialization method. It is run automatically everytime a new instance of our particle class\n",
    "        is created. \n",
    "        '''\n",
    "        \n",
    "        print(\"Particle init'ed\")\n",
    "        self.m = m\n",
    "        self.y = y0\n",
    "        self.v = v0\n",
    "        self.t = 0.0\n",
    "        self.tf = tf\n",
    "        self.dt = dt\n",
    "        npoints = int(tf/dt) # always starting at t = 0.0\n",
    "        self.npoints = npoints\n",
    "        self.tarray = np.linspace(0.0, tf,npoints, endpoint = True) # include final timepoint\n",
    "        \n",
    "    def F(self, y, v, t):\n",
    "        '''\n",
    "        force on a free particle\n",
    "        F = 0.0\n",
    "        '''\n",
    "        return 0.0\n",
    "\n",
    "    def Euler_step(self): \n",
    "        '''\n",
    "        Method to increment a single time step.\n",
    "        '''\n",
    "        a = self.F(self.y, self.v, self.t) / self.m\n",
    "        self.y += self.v * self.dt\n",
    "        self.v += a * self.dt\n",
    "        self.t += self.dt\n",
    "\n",
    "    def trajectory(self):\n",
    "        '''\n",
    "        Loop over time steps to generate a trajectory\n",
    "        \n",
    "        '''\n",
    "        y_array = []\n",
    "        v_array = []\n",
    "        \n",
    "        for ii in range(self.npoints):\n",
    "            v_array.append(self.v)\n",
    "            y_array.append(self.y)\n",
    "            \n",
    "            #propagate in time using an Euler Method\n",
    "            self.Euler_step()\n",
    "        \n",
    "        self.y_array = np.array(y_array)\n",
    "        self.v_array = np.array(v_array)\n",
    "\n",
    "    def plot(self):\n",
    "        fig1 = plt.figure(figsize = [8,4])\n",
    "        ax1 = fig1.add_subplot(121)\n",
    "        ax2 = fig1.add_subplot(122)\n",
    "        \n",
    "        ax1.plot(self.tarray, self.y_array, \"k\", label = 'euler')\n",
    "        ax2.plot(self.y_array, self.v_array, \"k\", label = 'euler')\n",
    "    \n",
    "    \n",
    "        ax1.set_title('Trajectory')\n",
    "        ax1.set_xlabel(\"t\")\n",
    "        ax1.set_ylabel(\"x\")\n",
    "        \n",
    "        ax2.set_title('Phase space')\n",
    "        ax2.set_xlabel(\"v\")\n",
    "        ax2.set_ylabel(\"x\")\n",
    "\n",
    "        ax1.legend()\n",
    "        ax2.legend()\n",
    "\n",
    "    def results(self):\n",
    "        \"\"\"\" \n",
    "        Method to display the results at a given final time\n",
    "        \"\"\"\n",
    "        \n",
    "        print('\\nPosition and Velocity at Final Time:')\n",
    "        print('Euler:')\n",
    "        print('t = {0:0.2f} s | y = {1:0.3f} m | v = {2:0.3f} m/s'.format(self.t, self.y , self.v))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we create a new class spawned from our Particle class. This FallingParticle class inherits the methods of the parent class: Particle. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class FallingParticle (Particle):\n",
    "    \"\"\"\n",
    "    Subclass of Particle Class that describes a falling particle\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self,m = 1.0, x0 = 1.0 , v0 = 0.0, tf = 10.0,  dt = 0.1):\n",
    "        \n",
    "        print(\"Falling particle init'd\")\n",
    "        super().__init__(m,x0,v0,tf,dt)   # call initialization method of the super (parent) class\n",
    "        \n",
    "        # use super() instead of \n",
    "        # Particle.__init__(self)\n",
    "        \n",
    "        #what does super do?\n",
    "        #From Stack Overflow:\n",
    "        #def __init__(self):\n",
    "        #    mro = type(self).mro()             # Get the Method Resolution Order.\n",
    "        #    check_next = mro.index(ChildB) + 1 # Start looking after *this* class.\n",
    "        #    while check_next < len(mro):\n",
    "        #        next_class = mro[check_next]\n",
    "        #    if '__init__' in next_class.__dict__:\n",
    "        #        next_class.__init__(self)\n",
    "        #        break\n",
    "        #    check_next += 1\n",
    "    \n",
    "    def F(self, x, v, t):\n",
    "        '''\n",
    "        Force on a particle moving under a gravitational field near the Earth's surface.\n",
    "        F = mg\n",
    "        '''\n",
    "        \n",
    "        g = 9.8\n",
    "        return  -self.g*self.m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "freeball = Particle(m = 10)\n",
    "freeball.results()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}