Lab.02-Skydiver-Problem.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Skydiver Problem\n",
    "## Lab 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Name:\n",
    "### Student ID:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Instructions:\n",
    "*Complete the notebook by reading all of the markdown cells and filling in the *blanks*. In the code cells, you will need to replace with the correct code everywhere you see a `___`. You will also need to add new code as appropirate. Please complete this lab during the lab slot and ask for help from the instructor or the TA as often as needed.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This lab involves modelling a skydiver falling through the air.  Initially, air resistance is very small, but upon deployment of a parachute, air drag suddenly increases to a large value.  The mass of the jumper is 60 kg, acceleration due to gravity is 9.81 m/s^2 , the jumper starts from a height of 1000m and opens the parachute at about 200 m above the ground.  The jumper feels a drag force $F=-k v |v|$, where $k$ changes to a high value after the parachute is opened."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Model the action of the parachute with a hyperbolic tangent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " \n",
    " To model the action of the parachute we will use the hyperbolic tangent function.  Plot `np.tanh(x-5)` from x=0 to 10.  Don't worry about making the graph look perfect, but make sure the range of the graph goes from -1.2 to 1.2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(___, ___, 0.01)\n",
    "y = ___\n",
    "\n",
    "fig, axes = plt.subplots()\n",
    "plt.ylim(___, ___)\n",
    "plt.plot(x, y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see from the graph, to a good approximation the function goes from one constant value to another.  Now plot the function 6 - 2 tanh( (x-5)/0.1 ) from $x$ = 0 to 10. Choose an appropriate range so that you can see both the top and bottom of your plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code to make plot goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how we have transformed the function so that it changes more abruptly, has a wider range and different end values.  Let's now choose appropriate parameters to model our system, and store them in variables\n",
    "\n",
    "* m = 60  - mass of jumper in kg\n",
    "* g = 9.81  - acceleration due to gravity\n",
    "* h = 1000   - initial height in m\n",
    "* y0 = 200  -  height at which the jumper opens the parachute\n",
    "* kd = 0.48  -  drag coefficient without parachute in kg /m\n",
    "* kp = 30.24 -  drag coefficient with parachute fully deployed\n",
    "* d = 20  -  roughly the distance (m) the jumper falls during the parachute opening"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize all the variables here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a new Python function that implements the the formula\n",
    "$$k(y) = kp - (kp-kd)/2 (1+tanh( 2 (y-y0)/d  )).$$  \n",
    "\n",
    "k(y) is our mathematical model for the drag coefficient as a function of the jumper's height above the ground."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def k(y):\n",
    "    return ___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the function k(y), with y ranging from 0 to 1000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.arange(0, 1000, 0.1)\n",
    "\n",
    "fig, axes = plt.subplots()\n",
    "plt.xlabel(\"height (m)\")\n",
    "plt.ylabel(\"drag coefficient (kg/m)\")\n",
    "plt.plot(y, k(y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check your values for k(0) and k(1000).  They should be 30.24 and 0.48 respectively.  \n",
    "\n",
    "k(0)  \n",
    "k(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(k(0), k(1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Solving the problem numerically"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Now that we've modelled the skydiver problem with a height-varying drag coefficient, let's find the height and velocity as functions of time using Euler's method, assuming initial velocity of zero. \n",
    "\n",
    "The equation that we want to solve is\n",
    "$$ \\frac{d^2y}{dt^2} = -g - k(y)/m \\frac{dy}{dt} \\left| \\frac{dy}{dt} \\right| \\;, $$\n",
    "\n",
    "along with appropriate initial conditions.\n",
    "\n",
    "The first task is to rewrite this second-order differential equation as pair of two first-order differential equations by introducing $v = \\frac{dy}{dt}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "\\frac{dy}{dt} &= \\mbox{___} \\\\\n",
    "\\frac{dv}{dt} &= \\mbox{___}\n",
    "\\end{align}\n",
    "\n",
    "*Determine what should be on the right hand side and edit this Markdown cell to replace the `\\mbox{___}` accordingly.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implement your numerical solution below using Euler's method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computation routine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute():\n",
    "    tmax = 400\n",
    "    dt = ___   # you should check that your choice of dt gives a good solution\n",
    "    \n",
    "    N = round(tmax/dt)\n",
    "    \n",
    "    # allocate memory for y, v, and t arrays\n",
    "    y = np.zeros(N)\n",
    "    v = np.zeros(N)\n",
    "    t = np.zeros(N)\n",
    "    \n",
    "    # initialize variables\n",
    "    y[0] = ___\n",
    "    v[0] = ___\n",
    "    t[0] = ___\n",
    "    \n",
    "    # solve using Euler's method\n",
    "    for i in range(N-1):\n",
    "        y[i+1] = ___\n",
    "        v[i+1] = ___\n",
    "        t[i+1] = ___\n",
    "        \n",
    "        if y[i+1] <= 0:\n",
    "            break\n",
    "            \n",
    "    # truncate arrays to when skydiver has landed (y<=0)\n",
    "    N = i+2\n",
    "    y = y[:N]\n",
    "    v = v[:N]\n",
    "    t = t[:N]\n",
    "    \n",
    "    return (y, v, t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run your numerical model and plot the height $y$ vs time $t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y, v, t = compute()\n",
    "\n",
    "# plot the solution\n",
    "fig, axes = plt.subplots()\n",
    "\n",
    "plt.plot(___, ___)\n",
    "plt.xlabel(___)\n",
    "plt.ylabel(___)\n",
    "plt.title('Skydiver height vs time')\n",
    "plt.ylim(0, h)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Analysis of solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimate the total time the skydiver spends in the air.\n",
    "\n",
    "* Easier (approximate) method: Use the plot to estimate time when when y = 0.\n",
    "* Harder (exact) method: Use linear interpolation to estimate the time when time y = 0 (we did this is a previous lecture).\n",
    "* *Try either approach (or both!) as you prefer.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_landing = ___\n",
    "print(\"The skydiver lands after\", t_landing, \"s.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the velocity $v$ versus time $t$.  Remember to use good axes labels with units and add a title!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the velocity vs time\n",
    "fig, axes = plt.subplots()\n",
    "\n",
    "plt.plot(___, ___)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the approximate terminal speed of the skydiver before the parachute is deployed?  What about after?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can look at the actually values of v (using `print(v[i])` at some appropirately chosen values of `i`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Answer: Terminal speed without the parachute is `___` m/s, while with the parachute the it is `___` m/s."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate and plot the acceleration of the skydiver (include axes labels with units).\n",
    "\n",
    "* *Hint: you will need to calculate acceleration as the change in velocity divided by the change in time. You have already computed both velocity and time.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.zeros_like(v)\n",
    "N = len(a)\n",
    "for i in range(N-1):\n",
    "    a[i+1] = ___\n",
    "    \n",
    "fig, axes = plt.subplots()\n",
    "plt.plot(t,a)\n",
    "plt.xlabel(___)\n",
    "plt.ylabel(___)\n",
    "plt.ylim(___, ___)\n",
    "plt.title(___)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the largest acceleration felt by the skydiver during his descent (express your answer as a multiple of g)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the max() function returns the maximum of an array\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Answer:  The largest acceleration felt is approximately `___` g."
   ]
  }
 ],
 "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": 1
}