Skip to content
Permalink
716aeff5c4
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

compmech-project03/03_IVPs_project.ipynb

Go to file
 
 
Cannot retrieve contributors at this time
975 lines (975 sloc) 92.4 KB
Raw Blame
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Initial Value Problems - Project\n",
"\n",
"![Initial condition of firework with FBD and sum of momentum](../images/firework.png)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are going to end this module with a __bang__ by looking at the flight path of a firework. Shown above is the initial condition of a firework, the _Freedom Flyer_ in (a), its final height where it detonates in (b), the applied forces in the __Free Body Diagram (FBD)__ in (c), and the __momentum__ of the firework $m\\mathbf{v}$ and the propellent $dm \\mathbf{u}$ in (d). \n",
"\n",
"The resulting equation of motion is that the acceleration is proportional to the speed of the propellent and the mass rate change $\\frac{dm}{dt}$ as such\n",
"\n",
"$$\\begin{equation}\n",
"m\\frac{dv}{dt} = u\\frac{dm}{dt} -mg - cv^2.~~~~~~~~(1)\n",
"\\end{equation}$$\n",
"\n",
"If we assume that the acceleration and the propellent momentum are much greater than the forces of gravity and drag, then the equation is simplified to the conservation of momentum. A further simplification is that the speed of the propellant is constant, $u=constant$, then the equation can be integrated to obtain an analytical rocket equation solution of [Tsiolkovsky](https://www.math24.net/rocket-motion/) [1,2], \n",
"\n",
"$$\\begin{equation}\n",
"m\\frac{dv}{dt} = u\\frac{dm}{dt}~~~~~(2.a)\n",
"\\end{equation}$$\n",
"\n",
"$$\\begin{equation}\n",
"\\frac{m_{f}}{m_{0}}=e^{-\\Delta v / u},~~~~~(2.b) \n",
"\\end{equation}$$\n",
"\n",
"where $m_f$ and $m_0$ are the mass at beginning and end of flight, $u$ is the speed of the propellent, and $\\Delta v=v_{final}-v_{initial}$ is the change in speed of the rocket from beginning to end of flight. Equation 2.b only relates the final velocity to the change in mass and propellent speed. When you integrate Eqn 2.a, you will have to compare the velocity as a function of mass loss. \n",
"\n",
"Your first objective is to integrate a numerical model that converges to equation (2.b), the Tsiolkovsky equation. Next, you will add drag and gravity and compare the results _between equations (1) and (2)_. Finally, you will vary the mass change rate to achieve the desired detonation height. \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"1. Create a `simplerocket` function that returns the velocity, $v$, the acceleration, $a$, and the mass rate change $\\frac{dm}{dt}$, as a function of the $state = [position,~velocity,~mass] = [y,~v,~m]$ using eqn (2.a). Where the mass rate change $\\frac{dm}{dt}$ and the propellent speed $u$ are constants. The average velocity of gun powder propellent used in firework rockets is $u=250$ m/s [3,4]. \n",
"\n",
"$\\frac{d~state}{dt} = f(state)$\n",
"\n",
"$\\left[\\begin{array}{c} v\\\\a\\\\ \\frac{dm}{dt} \\end{array}\\right] = \\left[\\begin{array}{c} v\\\\ \\frac{u}{m}\\frac{dm}{dt} \\\\ \\frac{dm}{dt} \\end{array}\\right]$\n",
"\n",
"Use [two integration methods](../notebooks/03_Get_Oscillations.ipynb) to integrate the `simplerocket` function, one explicit method and one implicit method. Demonstrate that the solutions converge to equation (2.b) the Tsiolkovsky equation. Use an initial state of y=0 m, v=0 m/s, and m=0.25 kg. \n",
"\n",
"Integrate the function until mass, $m_{f}=0.05~kg$, using a mass rate change of $\\frac{dm}{dt}=0.05$ kg/s. \n",
"\n",
"_Hint: your integrated solution will have a current mass that you can use to create $\\frac{m_{f}}{m_{0}}$ by dividing state[2]/(initial mass), then your plot of velocity(t) vs mass(t)/mass(0) should match Tsiolkovsky's_"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"plt.rcParams.update({'font.size': 22})\n",
"plt.rcParams['lines.linewidth'] = 3"
]
},
{
"cell_type": "code",
"execution_count": 93,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'velocity (m/s)')"
]
},
"execution_count": 93,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXyU9bX48c9JQgg7lh0ksoOyGgIIKoui1CpCXapURWwR0WqvWrBYrYa2Fi54r9aK4vJT8aLUUisgoigqgrJjABHZ90VZhJQlISQ5vz/mSZiQzGQmmckzy3m/XvOamWeb8xCdM99dVBVjjDHGW4LbARhjjIk8lhyMMcaUYMnBGGNMCZYcjDHGlGDJwRhjTAlJbgcQCvXr19cWLVq4HYYxxkSV1atXH1bVBqXti4nk0KJFC1atWuV2GMYYE1VEZJevfVatZIwxpgRLDsYYY0qw5GCMMaYESw7GGGNKiIkGaROYWZn7mDx/E/uPZdO0bjXGDmrP0IubuR2WiQJnzpxh79695OTkuB2KCVJiYiJ169alfv36JCQEXh6w5BAncia0YujpIwwFMpaNI6P/RJgNOR/VI+XR7W6HZyLc3r17qVWrFi1atEBE3A7HBEhVOXPmDD/88AN79+4lNTU14HMrpVpJRJ4UERWRTs77diKyVEQ2O89tvY71uc+UX8rpI0Wvx3/xaKnbjfElJyeHevXqWWKIMiJCcnIyzZo14+TJk0GdG/bkICJpwCXAbq/NU4EpqtoOmAK8FOA+Y4xLLDFEr2Cqk4rOCUMcRUSkKp4v+PsAdbY1BNKAGc5hM4A0EWngb18444wHGQvHIeOzkPFZAEWvMxaOczkyY0wkCnfJ4U/AdFXd4bWtObBPVfMBnOf9znZ/+4oRkVEiskpEVh06dCjMtxH9MvpPRJ+sgz5ZB6DodUb/iS5HZmLVrMx9XDrxM1qO+4BLJ37GrMx9Ibt2t27dyM7ODtn1du7cSf369YO6/siRI1m8eDEAs2bNYsWKFT6PHTFiBM8//3yxbWPGjCEjI6P8QYdZ2BqkRaQ30AMIy09TVX0ZeBkgPT3dlrMzJoLMytzHo//+huwz+QDsO5bNo//+BiAkPeTWrFlT4WtU9Pqvvvpq0etZs2aRnp5Oz549wxlWpQpnyaEf0AHYISI7gfOB+UBroJmIJAI4z02BPc7D1z5TETUaFr18st+EotcZC8dBRh3PY7K1/ZvQmDx/U1FiKJR9Jp/J8zeF5PoiwokTJwDP3GqPP/44vXv3JjU1lbfffptnn32Wnj170qZNm6Jf94WlgzFjxtCzZ086d+5ctM/f9b/77juuvvpqunTpQufOnZk2bRoA/fv3Z+7cucyfP585c+YwceJEunXrxptvvlmue5o0aRI9e/YkLS2NwYMH8/333wMlSx3e70eMGMHo0aO54ooraNu2LcOHDydUSz+HreSgqhOBojoLJ0Fcp6rrReQ+YBgw3XnOVNVDznFrfO0zFTB2C5n//VO04w1kMLZo8/gvHj1btXTyoEvBmWjTYtwH5Tpv37HsMs/dOfHaoK97+vRpli5dysqVK+nfvz+TJk1ixYoV/POf/+TRRx/lyy+/BODIkSN06dKFp59+mi+++IJhw4axbds2n9fNy8tjyJAhPPXUU9x8881F1/A2aNAgrr/+etLT07n//vt9XmvixInFShv79+/nvvvuA2D69Ols3bqVZcuWkZCQwIsvvsjvfvc73nrrrTLvff369SxYsICEhAQuvvhiFixYwFVXXVXmeWVxa5zDaGCaiDwBHAWGB7jPlFNO9knanFqD9p4Gq8aWfYIxfpT1BX7pxM/Yd6xknX2zutX4atwVIY/nlltuASAtLY1Tp04Vve/evTtbt24tOi45OZnbb78dgH79+lGtWjU2bdpE7dq1S73upk2byMvLK0oMAPXq1StXjOPGjSuWPMaMGVP0es6cOaxatYq0tDTAk5Tq1KkT0HWHDh1KSkoK4Ln/bdu2RVdyUNUWXq83Ar18HOdznym/Tcs/ompySzrUa0TGwnHFxjoU9mB6st8EMlyKz8SWsYPaF2tzAKhWJZGxg9qH5fMKvxwTExNLvM/Ly/N5nqr67aIbqiqasqgqjz/+OL/61a9K7EtKSqKgoKDo/bmj1AvvFcq+32DY3Epx4tS3H5J1/gDAei6Z8Bt6cTMm3NCZZnWrIXhKDBNu6Oz6dC25ubm8/fbbACxevJicnBzat/edsDp06EBSUhIzZ84s2nZutRJA7dq1ycrKKndc119/PS+88AJHjx4FPNVka9euBaB169asXLkSgAMHDvD555+X+3OCYdNnxAEtKKD54cXk9StfQ5kx5TH04mauJ4Nz1atXjy1bttCrVy9OnTrFjBkzSE5O9nl8UlISs2fP5v777+dPf/oTCQkJjBkzhjvuuKPYcXfccQcjRoxg5syZPPzwwwwfHlxt+B133MHhw4fp168fAAUFBdx333107dqVUaNGcdNNN9G1a1fatWtHr16VU7EilVVsCqf09HS1leB827lxDdX+cQMNn9iKJCR4eiU5jc8ZC8cVlRi8X1OjIYzd4lbIJsJ89913XHjhhW6HUSE7d+4kPT2dw4cPux2KK0r7G4rIalVNL+14KznEgf0rZ1Ol3qU0KhxC7/Wln5FxttHLei4ZYwpZm0McqL3nM5IvvMbtMIxxVYsWLeK21FAeVnKIcf/J+pGWpzeScEnpXQ+t55IxpjSWHGLcliXvk1ytI51rlt5nOqP/xKKqJBmfVdSDyePRUs8xxsQ+Sw4xLm/TR+S3uNLtMIwxUcbaHGJYQX4+rY8t4fyeP/d9kM25ZIwphZUcYtjWdUtISahBaquLfB9kPZeMMaWwkkMMO/z1HA40uNztMEw8mtz2bMnT+xGhpVDvmU6nTp3KM8884/f4VatWcdtttwFw7NgxJk2a5PPYnTt3kpSURLdu3ejSpQvp6eksWbKkaH9GRkaxeZZefvllWrduzbZt25g9ezbdu3enU6dOdOzYkf/5n/+pyG0GxUoOMazegS/I6/94wMdbzyUTMr5Km1FQCh09enSZx6SnpxfNmFqYHB555BGfx9etW7dojYgXX3yRX//613z33Xcljps8eTKvv/46ixYtolmzZhw+fJj333+fpk2bkpWVRffu3enZsyeXXx7+H31WcohRh37YS5O8PbTreXXA59icSyZaLF++nAEDBtC9e3e6d+/OBx94pgFfuHAhbdu2LZrnaMSIEYwbN67o9d13383AgQNp3749d999N7m5uSWufe4v+QkTJtC5c2e6du1Knz59KCgoYOHChaSnewYW/+Y3v+HYsWN069aNPn36lBl7//792b17d4ntjz/+OO+8805RYgDo1asXTZs2BaBOnTpceOGF7Nq1K5h/qnKzkkOM2r50Fsk1unNxckrZB5chY+G4s+0RNq2GAU8VUbjOzfA/gd2xY8cYPXo08+bNo0mTJhw4cIAePXqwfv16+vfvz/Dhw/n1r3/N9ddfz+bNm4utobB8+XKWLFlCSkoKP/vZz3j55Zf9rsEwbdo05syZw1dffUXt2rU5cuQICQnFf1NPmTKF9PT0gFene++997j11luLbXvjjTdo3LgxS5Ys8Tl9+MaNG1m2bBkvvfRSQJ9TUZYcYlTS1k/Iax3knO41GhYV+717LlnjtCmhjC9wvwmgrHPLsGTJEnbs2ME115wd9S8ibN26lfT0dB577DEGDhzI7373O1avXk1S0tmvuVtuuYWaNWsCcOedd/Luu+/6TQ5z587l3nvvLfrCLu9aDoUli0OHDpGXl8eyZcuK7e/ZsyeZmZl89NFH/OIXvyhx/oEDBxgyZAhTpkwpKkmEmyWHGHQm9zRtT6zkTO+/B3eij55LxkQSVaVLly4sWrSo1P1ZWVns3r2bqlWrcuTIEVJTU31ex99aDoXHhEJhm0N+fj4PPfQQt956K8uXLy/af9FFFzFx4kQGDRoEUCxBHDx4kIEDBzJ27NhSE0e4hLXNQURmichaEckUkcUi0s3ZvlNENorIGucxyOucdiKyVEQ2O8+R2b0hgm1a9SkHExtTr0np/1MEImPhOGR8VlGjdOHrjIXjQhWmiWVe42cC2h6EPn36sGXLlmLrGqxcubLoi/yuu+5i5MiRTJs2jWHDhnH8+PGi42bOnMnJkyfJy8tj+vTpDBgwwO9nDR48mBdffLHoGr7Wcjh16lRAi+wkJiYyadIkDhw4wOzZs4vt69KlC/Pnz+e3v/0t//znP4s+76qrruL+++9n5MiRZV4/lMJdcrhTVbMARGQI8BqQ5uy7SVXXl3LOVGCKqk4XkduBl4DQrysYw/6z7gNONO1Pmwpcw6bVMBUSxnap8847jzlz5jB27FgefPBBcnNzadWqFe+//z5/+9vfyM7O5ve//z0iws0338yoUaOYMWMGAH379mXo0KHs3r2bvn37MmrUKL+fNXz4cPbt28cll1xCUlIStWrVKlFi+clPfsJtt91G586dOe+884p1Uy1NSkoKf/nLXxg/fjxDhgwptq8wQQwaNAhVZdWqVWzevJmXXnqpqK3hv/7rv7jrrruC/WcLWqWt5yAiw4Hfqmq6iOwErjs3OYhIQ2AzUE9V80UkETgCtFXVQ76ubes5FLfjT53Ju+452qb5/1Xkl1e10rnJwdZ9iD+xsJ7DiBEjSE9P99vGEMuCXc8h7F1ZReRVEdkNPAXc6bXrLRFZJyIviEhdZ1tzYJ+q5gM4z/ud7SYA+3Zsok7BMVp37VuxC/mYVgMoNhbCGqiNiU1hb5BW1ZEAInIHMBn4GXC5qu4RkarAs8DzwO3BXFdERgGjAJ8NTvFoz/JZJNXpTbqz0Hq5eTdOA2TYWAcT3d544w23Q4gqlTYITlX/DxggIvVUdY+z7TTwAnCpc9geoJlTnYTz3NTZfu71XlbVdFVNb9CgQaXcQzRI2fkpie0HlX1gkKyB2pj4EraSg4jUBM4rTAQiMhj4EcgRkTqqmiWefmS3AmsAVPWgiKwBhgHTnedMf+0N5qxTJ4/TJnsdBX3eDvm1rYHamPgSzmqlGsBMEakB5ONJDIOBRsC7TqkgEdgA3Od13mhgmog8ARwFhocxxpiyadk8Uqq24cK69SvtM230tDGxKWzJQVV/AC7xsftiP+dtBHqFJagYd3rDh+Q2r0APJX9s9LQxccUm3osRWlDABUe+pEmPIWUfXB5jt3imPcjIson4TFAyMkJ7vV69etGtWzcuuuiioqmwu3Xr5rPv/5QpU3juuefKvO7tt9/O1KlTAc8keIUT9oXL1q1bady4cVg/oyJs+owYsWPjaqoJNG+XVvbBFWRTe5tgjB8f2gRROO3Ezp07A5rw7je/+U3oPjyOWMkhRny/6n321LsMSQj/n9Sm9jaR6LvvvuOSSy6ha9eudOrUiWeffRYoXgrIy8vjoYceomPHjnTq1IlHHnmE/Px8v9dds2YNnTt35ssvvwTg9ddfp3PnznTp0oUbb7yRQ4cOoaq0bNmSb7/9tui8Z555hrvvvpv8/HzuueceOnToQNeuXenbt+QYpJycHG666SYeeeQRvvrqK7p161Zsf8eOHVmxYkWF/n2CZckhRtTZ8xkpF11T9oFhZOtOm0IZGSDiecDZ16GuYvL2/PPPM2TIENauXcv69eu58847Sxzz4osvsmHDBtasWcPq1atZsWIFr732ms9rfvzxxwwfPpyZM2dy2WWXsXbtWv74xz/yySefsG7dOtq1a8eDDz6IiDB8+HCmTZtWdO60adO46667+Prrr1m0aBEbNmxg7dq1JeZUOnLkCFdffTUDBgxg0qRJXHrppVSpUoWvvvoKgM8//5zq1avTs2fPEP1LBcaSQwzI+vEwF+Rupd0lP6ucD/QxetpGTptCGRmg6nnA2dfhTA59+/bllVde4YknnuDzzz+nbt26JY5ZsGABd911F1WqVKFq1aqMGDGCBQsWlHq9Dz/8kDFjxvDxxx/ToUMHAD777DOuu+66oraCe+65p+j8ESNG8NZbb5Gfn09mZibZ2dn06dOHNm3acPr0ae6++26mT59ebCbYU6dOcemll/LQQw8Vq/564IEHeOGFFwBPm4kbVWOWHGLAlqWz2V69MynVa1XOB1rjtIlAt9xyC4sWLaJly5Y89dRTpTZQlzZNt69pu9u3b09ubi6rV68O6PyWLVvSpk0bPv74Y9544w1GjBgBeCYK3LBhAzfddBOZmZl07NiRgwc9P55SUlLo1asXs2fPLla9deutt7J48WIyMzP58ssvSywOVBksOcQA3TyfnJYDXflsGzltyvLkk5XzOVu2bKFp06bcdddd/PGPfyy1jv6qq67i9ddfJy8vj9zcXN58800GDiz9/51WrVrx8ccfM3bsWP71r38BcOWVV/L+++8Xfbm/8sorxc4fMWIEr7zyCu+88w7Dh3uGaB08eJCcnByuueYaJk2aRPXq1dm5cycACQkJvPHGG6SkpPDLX/6SM2fOAJCcnMydd97J9ddfz/Dhw0lJqfiKjsGy3kpRLj8/n9ZZS8m94c+ufL6NnDZlCWdVkrd//OMf/OMf/yA5ORkRKWqQ9nbvvfeyffv2ogbfa665xu/016mpqSxYsIBBgwaRnZ3NHXfcwZ///GeuvPJKRIQ2bdoUW7bz5ptv5oEHHuCyyy4rWgd6165d3HPPPeTl5ZGfn8/QoUPp0aMH27ZtAzwlj6lTp/Lggw9y4403MnPmTKpWrcrIkSP561//yujRo0P5zxSwSpuyO5ziecru71Z9TvV5D3DBE6UtjVEJfEztbdN6x5ZYmLI72rzxxhu89957JRqwyyvYKbut5BDljq2ZS1ajflzgVgA2ctqYkBs4cCC7du1i7ty5rsVgySHK1T+wkLyB7lQpAbbutDFh4KsHVWWy5BDFDu3fRcP8A9TofqXboQA2cjrWldZTx0SH8jQfWHKIYtuXzqJKrZ6kJVd1OxTAGqdjWWJiImfOnCE5OdntUEw5ZGdnU6VKlaDOsa6sUSxp+ydoG3e6sAbDRk5Hv7p16/LDDz9QUFDgdigmCKrKqVOn2LdvHw0bNiz7BC9WcohSp09n0/bEavIumep2KGdZ43TMql+/Pnv37mXTpk1uh2KCVKVKFRo1akTt2rWDOs+SQ5TavOITqlY5n3aNznc7lLOscTpmJSQk2FrtccaqlaLUiW/m8WPTMC3sEwI2ctqY6BbW5CAis0RkrYhkishiEenmbG8nIktFZLPz3NbrHJ/7zFnNDi2iftpgt8Pwyde03p6d1v5gTKQLd8nhTlXtqqoXA08DhXPjTgWmqGo7YArwktc5/vYZYM+2b6mhJ2ndpY/boQTNZm41JjqENTmoapbX2zpAgYg0BNKAGc72GUCaiDTwty+ccUabvctns71uHyQh0e1QfPMxrbcxJjqEvUFaRF4FrgYE+CnQHNinqvkAqpovIvud7eJn36FzrjsKGAXEXUNZjd2fkt9tuNth+Oc9l1L/CUVtD2CD44yJBmFvkFbVkaqaCvwBmBzC676squmqmt6gQfwULE4ez6JV9re07RO57Q3nsvYHY6JPpfVWUtX/AwYAe4FmIpII4Dw3BfY4D1/7DLBp2Vx2pbSnZu2fuB1KhVn7gzGRK2zJQURqikhzr/eDgR+Bg8AaYJizaxiQqaqHVNXnvnDFGW3OfPcRJ1IjYy6lgFn7gzFRJ5xtDjWAmSJSA8jHkxgGq6qKyGhgmog8ARwFvCvQ/e2La1pQQIsfvyL3p79zO5TgWPuDMVEnbMlBVX8ALvGxbyPQK9h98W7btyuoJkk0b9vV7VDKzdfkfEXzL4EtDmRMBLAR0lHk0Oo57K1/OcTgtMnW/mBMZLHkEEXq7vuc6p1+5nYYFWPtD8ZEBZt4L0ocPfw9zXN3UKXXT90OpWKs/cGYqGAlhyixdekcttboRtWUGm6HEjI2/sGYyGXJIVpsnk9uq8hf2CcUrP3BGPdZcogCeWfO0Ob4clpc8nO3Qwkta38wJmL5bXMQkZrALcCVwPlANrAWeFdVl4c/PAOwJXMhVRPr0er81m6HElrW/mBMxPJZchCRcXhGK6cBHwJ/Bp4HDgN/FZHPRaRDpUQZ546umcvBxv3dDiOsrP3BmMjir1opC+igqr9R1f9T1U9U9X1VnaSqV+IZyRxBa1TGroY/LKJu1yjvwlpO1v5gjDt8Viup6ov+TlTVTYCtNh5m3+/dTv38g9RKu8LtUMKrRsOiL39rfzDGfWWOcxCRh4H/p6pZIvJ/QA/gt6r6cdijM+xc9h5Vaveie1IVt0MJL2t/MCaiBNJbaYSTGAYADYFfAX8Nb1imUNXtC6Dt1W6HUams/cEY9wWSHPKd5wHAW6q6JMDzTAXlZJ+izclM2vQZ6nYoEcHaH4ypPIF8yWeLyGPAbcDHIiJAcnjDMgCbV3zEvuSW1KnX2O1QKpeNfzDGdYHMrTQCuA8Yq6rfi0hr4K2wRmUAOLX+Q0416+92GJXP2h+McZ2/cQ7/LSJ9VHWzqj6oqv8GUNVtqmo/58JMVTn/8GIadr/e7VBcZe0PxrjDX7XSOuAhEdkqIq+KyGARSamswOLd7i3rSNbTtOxo6x6VxtofjAkvn8lBVd9S1ZuBi4B/A9cB34rIeyIyQkTq+7uwiNQTkXkisklE1onIv0WkgbNvp4hsFJE1zmOQ13ntRGSpiGx2nuPyZ+H+FbPZ+ZNLkYQ4b/u39gdjXFFmm4Oq5gLznAcicgkwFBgDdPJ3KjBJVRc6500GJgK/dvbfpKrrSzlvKjBFVaeLyO3AS0CMjwArqdaeT8lPH+V2GO4LoP2h3wWLWYgtMWpMKAX8s1REkkWkOp7qpj+pqr/EgKr+WJgYHMuAC8r4jIZ45nKa4WyaAaQVljjixfGsH2mZs5G2va91O5SI4qv94Ytdl589yKqYjAmJMpODiNwoInvwzMh6HDjhPAdMRBKAe4E5XpvfcqqbXhCRus625sA+Vc0HcJ73O9vPveYoEVklIqsOHToUTDgRb8vSuWyvdhHVa9Yt+2BjjAmDQEoOk4EbgCqqmqiqCaqaGOTn/B1PUnneeX+5qnbFMxWHeG0PmKq+rKrpqpreoEFsFSzObPyIUxdc6XYYkcer/aHfBYuR8VlFVUuFrzMWjnMrOmNiSiDjHA6o6sryfoCIPA20BQaragGAqu5xnk+LyAucLVHsAZqJSKKq5otIItDU2R4XCvILaHVsCWcGP+Z2KJHHqy2hqI0BT2IorGrKWDjO070VrP3BmAoIpOTwdxH5s4ikichFhY9ALi4iTwHdgaGqetrZVkNE6jivBbgVz7oRqOpB5/Uw5xLDgExVja16Iz+2fbOEHKlG09Yd3Q4lKlkXV2NCI5CSQzPgYeBOzs6zpEArfyeJSEfgD8BmYIknD7AD+B3wrlMqSAQ24BmBXWg0ME1EngCOAsMDvZlYcPjr95GGfUs2spjibIpvY8IqkOTwW6CNqh4I5sKq+i2e9oTSXOznvI1A3I78Om//QvL7/cHtMCKfdXE1JqwCSQ67gk0MpnyOHNxH07zdpPQcVPbBpkhG/4lk9J8IFG9/8E4YVsVkTHACSQ4rRGQGMBPIKdyoqvPCFlWc2rZkNlVqpHFxVZulxBjjrkCSQ3fn+QGvbYozYtqETsLWj8lrfZXbYUQfr/aHwi6uhayKyZjyCWT6jAGVEUi8O3Mml7YnVnCm97NuhxJ9AujialVMxgTH35TdTcs6WUTibBWa8Nm86lMOJzamfpMWbodijDF+Sw7viMh64G1ghdc4hVTgp3gWAfof4N1wBxkP/rPuA4436UdrtwOJdlbFZExI+EsOfYGbgSeBPiJyGkgBfgDeA36pqjvDHmGcaPzDIs5ca1VKFWZVTMaEhM/koKoK/BP4p4gkAfWBbFXN8nWOKZ/9uzZTp+Aodbv2dTsUY4wBAuuthKrmAd+HOZa4tWfZLBLqXEKPpID+HCZQAVQx2VrUxpTOvo0iQNWdC8jv/Au3w4g9NlGfMeUW52tQui/75AnanFpHm95D3A4lLtlEfcaUzpKDyzYtn8eeqq2pc15srUkRcWwtamOCUma1kohsxLMYzzRVDWoFOFO2nA0fknN+3C2RXflsoj5jghJIyWEY0A3YKiIvikjnMMcUN7SggNQjX9I4/Xq3Q4krtha1MWUrMzmoaqaqjgTaA1uAeSKySERuCHt0MW7XpkwSVLmgQ/eyDzbGmEoUTG+lXkB/4BTwETBaRG5R1VvCEVg8OLBqDon1LqVxgjX9VCobRW1MmQJpc/gdntXZtgF/B+Y5A+T+KiJbwxxfTKu15zMKev3G7TDij42iNqZMgfxkbQUMVtWfquoHTmIo5LPUICL1RGSeiGwSkXUi8m8RaeDsayciS0Vks/Pc1us8n/tiSdbRI7Q4vYV2va91OxRThqKxEBl1YHJM/udoTAmBJIddztKdRUTkEQBVXe3nPAUmqWp7Ve2Cp+Qx0dk3FZiiqu2AKcBLXuf52xczti6dzfZqnUipXsvtUOKbVxfXwiqmwlJD4WsbC2HiUSBtDrcCkwLYVoyq/ggs9Nq0DLhXRBoCaUDhqjYzgOedUoX42qeqhwKINWrkb5pPXsuBbodhgq1iMiZO+EwOInIVcDXQVES8E0EdH6f4JCIJwL3AHKA5sE9V8wFUNV9E9jvbxc++Q+dccxQwCiA1NTXYkFxVkJ9P66yl5P58vNuhGD+sodrEM3/VSrnACTzVQye9HhuBYLux/t251vPliLFUqvqyqqaranqDBtE1unjLmsWcSKhNkxYd3A7FeDtnFLWNhTDxzN+U3V8AX4jIu6q6vrwfICJPA23xNGoXiMgeoJmIJDolg0SgKbAHT8nB176Y8WPm+0ijvlzgdiCmOK9SQEZG0AVkY2KKv2qlm1V1JtBXREosNKCqL5R1cRF5CugOXFu4kpyqHhSRNXhGXk93njML2xT87YsV9Q98wZkrrUopotlYCBPn/DVIdwJmAj1K2aelbCtGRDoCfwA2A0tEBGCHqv4cz7iJaSLyBHAUGO51qr99Ue/wgd00yt9PtXRrjI5oNhbCxDl/1UpPOs93lefCqvotnmqi0vZtxDPiOqh9sWD70lkk1fi890cAABQ7SURBVEwnLbmq26EYY4xPgYyQ/j3witM1FRGpB/xKVSeHO7hYlLjtEwraXu12GCYYAa4oR/8JZPSfePYcq2YyUSyQcQ7DVPW/C9+o6hER+SVgySFIuadP0/bkas70LrO5xkSSAKqYCt8XJQerZjJRLpAR0qVVDdnyouWweeXHfJ/UjHqNmrsdigkzm3LDRLtAvuS3iMjDwDN4EsVDgE24Vw7Hv/kAadrf7TBMRXhVMT3ZbwIZC8cVm17Du8rJShEmmgWSHH6Lp1vpX/H0UloC3B7OoGJV00OLOTP4RbfDMBXhPRYCIKNOURKwKTdMLCkzOajqfuAKEanhvD8Z9qhi0L7tG6hZcILzulzqdigmjGw8hIkVAbUdiMggYCCgIvKJqn4S3rBiz57ls0ms25seiYluh2JC6dxqprJKEVbFZKJEIF1ZH8EzEG2Gs+l/RWSaqj4d1shiTLVdn1LQ9Ta3wzChZlNumBgVSMnhdqC3qh4HEJHngK8ASw4BOnUii9bZ69He17sdigknm3LDxJBAkoMUJgYAVT0uzlwYJjCbls0juWo7Otat53YoJpxsyg0TQwJJDitF5HXgFTy9lUYCq8IaVYzJ3fAhualXuB2GiTAZC8edrYqyUoSJMIEkhweAJ4Dn8Ixz+AT4cziDiiVaUMAFP35F7qCH3A7FVKYAqpjAxkKYyBVIV9aTwO8rIZaYtGPDSlIkkeZtu7odiqlMtvyoiXL+1nO4z9+JgaznYOCH1XNIqH85TRMCmanExANrqDbRwF/JobR1HAqVuZ6D8ai793POXPqw22EYN9lYCBOF/K3nUK51HMxZxw7/wPm526nS6xq3QzFuCnIshDVUm0gQyCC46sCjQCtVvU1EOgAdVHVW2KOLcluXzSGpRje6VavhdigmUlhDtYkSgfRWehE4AHRz3u/FM1rab3IQkaeBG4EWQGdVXe9s3wnkOA+A36vqfGdfO2AaUA84AgxX1aj92aSb55Pb8kq3wzCRJMiGaitFGLcEkhw6q+qdzvxKqOoJEQmkdXUW8DdgcSn7bipMFueYCkxR1ekicjvwEhCVAwTy8/Jo859l5Paa4HYoJkpYKcJEkkC+5HO934hISiDnqeqXqron0EBEpCGQxtk5nGYAaSLSINBrRJItmQs5llCPRqm20IvxoUbDopdP9puAPlmnqPTg/doYNwRSclgkIn8AqopIf+BhYHYFP/ctZwqOL4E/qOoxoDmwT1XzAVQ1X0T2O9sPnXsBERkFjAJITU2tYDihd3TNXKRxX1q6HYiJXH4aqq27q3FbICWHx/CMjD4OTAJW4KxzUk6Xq2pXPF1lBXi+PBdR1ZdVNV1V0xs0iLzCRcPvv6Bu1+vcDsNEiwBKEV/suvzs8VbFZMIskJJDkqo+BTwVig8srGpS1dMi8gIwx9m1B2gmIolOqSERaOpsjyoH9+2gXv5BaqYNcDsUEy2su6uJMIEkhz0iMgt4TVWXVOTDnNXkklQ1y6lWuhVYA6CqB0VkDTAMz7Kkw4BMVS1RpRTpdiybRVKtnnSvkux2KCYaWXdXEwECSQ7tgF8Cz4lIbeAN4E1V3evvJGfdhxuAxsACETkCDAbedUoFicAGwHuajtHANBF5AjiKZ5GhqJO87RPyOwx2OwwTray7q4kAohr4TBgi0hEYA9yuqlXCFlWQ0tPTddWqyJhF/HTOKXIntCL//kzqNmjidjgm2mUUTw6+FOvZlGET+pnAiMhqVU0vbV9As8GJSIKIXAeMB67FU3owpdi0fD77kltYYjChYd1djUsCmT7jf/G0DXyLZ/TyHaqaHe7AotWpbz/kVNN+bodhYoV1dzUuCaTN4UegVzAD2uJZs0OLyf35q26HYWKRze5qKlEgi/38pTICiQW7t6wjRXM4v1Nvt0Mxsci6u5pKFEjJwQRo34rZJP6kDw1tYR8Tbtbd1YSZJYcQqrn7U/K6j3Q7DBMPrLurCTNLDiFy4vgxWuV8B71tygzjLitFmFCw+o8Q2bxkLjtSLqRGrbpuh2LiTZDdXTMWjvOMn8ioA5Nt1mBTOis5hEjexg/JvcAW9jEuCLK7K1gpwpTNSg4hoAUFtDi6hPN7DnU7FBPvrBRhQsRKDiGw7ZulVJUUmrfp7HYoJt5ZKcKEiJUcQuDQ13PY1+Dysg80pjJZKcJUgJUcQuC8fV+Q13ec22EYU5yVIkwFWMmhgn48uJ9meTtp13OQ26EY45uVIkyQrORQQduWziapRhoXp1RzOxRjfLNShAmSlRwqSLZ+TF6rq9wOw5jAWSnCBMBKDhWQdyaXNsdXcKb3M26HYkzgrBRhAhC2koOIPC0iO0RERaST1/Z2IrJURDY7z20D2ReJNq/+nMOJjWjQtIXboRhTPlaKMD6Es+QwC/gbsPic7VOBKao6XURuB14CrghgX8TJWjsXadKXNm4HYkx5WSnC+BC2koOqfnnuAkEi0hBIA2Y4m2YAaSLSwN++cMVYUY1+WMR53WyiPRMjrBRhvFR2m0NzYJ+q5gOoar6I7He2i599h869kIiMAkYBpKamVlL4Z32/ewvnFfxI7W79K/2zjQkLK0UYL1HbW0lVX1bVdFVNb9Cg8gsXu5bNYlvtXiQmWZu+iUFWioh7lf3NtgdoJiKJTskgEWjqbBc/+yJO1R0LyOt0k9thGBMeVoqIe5VaclDVg8AaYJizaRiQqaqH/O2rzBgDkZN9kjan1tLukiFuh2JM+FkpIi6JqobnwiLPATcAjYHDwBFV7SgiHYBpwHnAUWC4qm5yzvG5z5/09HRdtWpVWO6jNGs//xfJS/6XCx9bUmmfaUxE8CpFeJcazuWdMDIWjjtbqrAlSiOKiKxW1fTS9oWzt9JvVfV8VU1S1caq2tHZvlFVe6lqO+d5k9c5PvdFkuxv55F1fsT2sDUmfIIsRQCM/+LRs2+syilqWGtqkLSggNQjX3J6wHS3QzGm8gXZFvFkvwnFjslYOO7seVaKiGhR21vJLbs3ryFB82lxYaklMWPih59SRGFSKCw1yPgsZHyWlSKiiJUcgrR/1RwS6l1G4wTLqybO+SlFZPSfWNTOIOOzipJGsd5NVoqIaJYcglRr92fk97zX7TCMiSw1GhaVBM6tSgLr/hqNLDkE4T/HjtDi9GYSL7nW7VCMiSzepQiAya8VSxZWiog+lhyCsHXpHJKqdaJLzdpuh2JMZLNBdFHPKs6DkL9pPtktrnQ7DGOiiw2ii0pWcghQQX4+LY8tJXfIk26HYkx0qUApIuODX5Fx0qqb3GAlhwBtXfsVJxNq0rTlhW6HYkz0CrIUYV1f3WMlhwAdyXwfadiXC9wOxJhoVo5SRNHx1mhdqazkEKB6BxZSq7P1UjImZPyUIrzZADp3WMkhAIe/30PjvH1U63GV26EYEzv8lCK8u7ta11d3WHIIwPals0iq2Z205Kpuh2JMbPIziM4ard1h1UoBSNz2CXmtr3Y7DGNi19gtkJEFGVlkXPta0WZrtHaPlRzKcCb3NG1PrCK39/Nuh2JMfLBG64hgJYcybF65gB+SmlK/carboRgTf6zR2jVWcijD8W8+QJv0x8ZpGuMCa7R2jWvJQUR2AjnOA+D3qjpfRNrhWSq0HnAEz1Khrv1VmxxcxOnrprj18caYQtZoXancrla6SVW7OY/5zrapwBRVbQdMAV5yK7j9OzZSs+A/tOl6mVshGGMKWaN1pYqoaiURaQikAYUDCmYAz4tIA1U9VNnx7F4+i8S6vemRmFjZH22M8ccarcPO7ZLDWyKyTkReEJG6QHNgn6rmAzjP+53txYjIKBFZJSKrDh0KT96otnMBie1/GpZrG2NCpAKN1hkf/MpmgPXBzeRwuap2BXoAAgTVV1RVX1bVdFVNb9CgQciDyz55nNbZ62nT+/qQX9sYE0Le1U2Fa0I4rLqp/FyrVlLVPc7zaRF5AZgDPAw0E5FEVc0XkUSgKbCnsuPbtOwDqlRtS8e69Sr7o40x5VWORutCVt1UnCvJQURqAEmqmiUiAtwKrFHVgyKyBhgGTHeeM91obzi94UNyUq+o7I81xlSEd1vE5LZw0vP63KVKvVnvptK5VXJoBLzrlAwSgQ3Afc6+0cA0EXkCOAoMr+zgtKCAC458yemr/quyP9oYEyoVGCMx/otH437ZUleSg6puBy72sW8j0KtyIypu58bVVJUEUtt1czMMY0yoWHVT0CKqK2uk+H7VHBLqXUbTBLc7cxljQsKqm4Jm336lqLPnM1I6/sztMIwx4WC9mwJiJYdzZP14iNTcbST1usbtUIwx4WbVTT5ZcjjHlqWzSarelW7Va7odijEm3Ky6ySdLDufQzfM53fJKt8MwxlQ2691UjCUHL/l5ebTOWkbujX91OxRjjJususmSg7cta76gasJ5tEy1OVaMiWtW3WTJwdvRNXOhUT9auh2IMSZyhKi6KdoShXVl9dLgwBfU6Xqt22EYYyLVOTPAeiuc9fXc14WirRuslRwch/bvpEH+99TobvMpGWN8qGB1k7eM/mfPicSShCUHx46ls0is1ZPuVZLdDsUYEw0CrG7yVvj+yX4TIr7KyaqVHFW2f0JBm6vdDsMYE438VDeVNur63JHZkVjlZCUH4HTOKVqf/Jq83i+7HYoxJhr5qW7ylrFwXLFEEMldYS05AJtXfEJyUirtGzZzOxRjTLQro10iWrrCWnIATn47jxPN+rsdhjEm1ngnCoCMkhP9QWR2hbU2B6DZwUU0SBvsdhjGmFgXoq6wGR/8CjLqeB6TwzNoN+6Tw96t60nRbFp37uN2KMaYWOc9Xfi1rxVtfrLfhFKnCAf3xkxEZLWSiLQDpgH1gCPAcFUNXRnKq9vZ+YUv/nSes6/0PsnGGBNSFegKW3zbMZTQVzdFZHIApgJTVHW6iNwOvATY6DRjTGzyM9FfMIkilKWIiEsOItIQSAOucjbNAJ4XkQaqesi9yIwxJkwC7ApbVgN2KImqhuXC5SUi3YE3VbWj17YNwO2q+rXXtlHAKOdte2BToJ/RvUlCd1/7Vh8oWB100NGnPnDY7SAqmd1zfIipe+7aKKFrUoLnR/z+401oWusAAKsPXOznrNXBfIddoKoNStsRcSWHQKnqy0CFR62JyCpVTQ9BSFHD7jk+2D3Hh3DdcyT2VtoDNBORRADnuamz3RhjTCWIuOSgqgeBNcAwZ9MwINPaG4wxpvJEarXSaGCaiDwBHAWGh/Gz4nFCJbvn+GD3HB/Ccs8R1yBtjDHGfRFXrWSMMcZ9lhyMMcaUEBfJQUTaichSEdnsPJeYqUpEEkVkiohsE5GtIjLSjVhDJcB7/qOIfCsia0VktYgMciPWUAnknr2ObS8ip0Tk6cqMMdQCvWcR+YWIfCMi653nRpUda6gE+N92QxH5QETWichGEXlBRCK1jdUvEXlaRHaIiIpIJx/HhP77S1Vj/gF8hmcQHcDtwGelHDMcmI8nYTYA9gIt3I49zPc8CKjuvO4KHAOquR17OO/Z2ZcILATeBp52O+5K+DunAxuAxs77OkCK27GH+Z6fLfzbAlWA5cAv3I69nPd7GdAc2Al08nFMyL+/Yr5B2pmOYzNQT1XznXETR4C26tU9VkQ+AF5X1X85758HdqnqZDfirohA7/mccwRPcuioqnsrL9rQCOaeReQx4DRQE6ipqmMqPeAQCOK/7beAT1X1NR+XihpB3PMzQHXgXud5MXC/qn7lQtghISI7getUdX0p+0L+/RUP1UrNgX2qmg/gPO93tntLBXZ5vd9dyjHRItB79jYc2BaNicER0D2LSBc8JaZnKj3C0Av073wR0EpEFonI1yLyuPNjIBoFes9/BtoBB4DvgfnRnBgCEPLvr3hIDqYMItIPz/9Mw8o6NpqJSBXgFWB04ZdLnEgCuuCZzLIfcA1wh6sRhd/NwDqgCdAM6CsiN7kbUnSJh+QQ6HQcu4ELvN6nlnJMtAh4ChIR6Q1MB4aqasCTF0agQO65CdAamOcU0R8E7haRaB04FejfeRfwL1U9rarHgdlAz0qNNHQCvecHgLdUtUBVs/Dc84BKjbRyhfz7K+aTgwY+HcdMPF8UCSLSABgKvFt5kYZOoPcsIj2Ad4Cb1GvG22gUyD2r6m5Vra+qLVS1BZ5Gy1dUdVSJC0aBIP7bfhu4WjyqAFcCaysv0tAJ4p53AD8FEJFkYCBQoq4+hoT++8vtlvjKeAAd8PRW2Ow8t3e2zwPSndeJwIvANucxyu24K+GeVwKH8PzPVvjo7Hbs4bznc47PIPp7KwXyd04A/hf4DvjWeZ3gduxhvufWwCfAN3h6ak0BktyOvZz3+xye3kd5eNpPvi3lfkP+/RXzvZWMMcYEL+arlYwxxgTPkoMxxpgSLDkYY4wpwZKDMcaYEiw5GGOMKcGSgzEhJiKvisi1QZ5TXUTecWbU3Cgi14UrPmMCYcnBmBASkQQ8I3E/DfLUMcBxVW0DDAZeFZGaoY7PmEBZcjCmDM48+o+JyEoR2S4iV4rIBBHJdNZHuNDr8N7AGlXNEZEMEfmHiMxzSgTviMjFIvKZM+++94yZtwBTAVR1C7AKzxxIxrjCkoMxgTmmqj2A3+OZp+dLVb0YeBN4zOu4oc7+Qt3xTPHQHs/I3ol4vvS7AHd6LVQTS7MCmxhgycGYwLzjPH8NqKp+4LxfDbTxOu5a4AOv9/NVNUs9s8CuAz5RzwR4J4FNeKZ5MCbiWHIwJjA5znM+noWC8HqfBCAiFwEHVfVIKecVHnvu+8KlK2NpVmATAyw5GBM6QyhepRSMmcA9AE5VUw/goxDFZUzQLDkYEzoVSQ6TgboishWYi2dWzeMhi8yYINmsrMaEgIg0wdO+0MXtWIwJBUsOxhhjSrBqJWOMMSVYcjDGGFOCJQdjjDElWHIwxhhTgiUHY4wxJVhyMMYYU8L/BwiKsuYW/awqAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"plt.rcParams.update({'font.size': 11})\n",
"plt.rcParams['lines.linewidth'] = 1\n",
"\n",
"def simplerocket(state,dmdt=0.05, u=250):\n",
" '''Computes the right-hand side of the differential equation\n",
" for the acceleration of a rocket, without drag or gravity, in SI units.\n",
" \n",
" Arguments\n",
" ---------- \n",
" state : array of three dependent variables [y v m]^T\n",
" dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n",
" u : speed of propellent expelled (default is 250 m/s)\n",
" \n",
" Returns\n",
" -------\n",
" derivs: array of three derivatives [v (u/m*dmdt-g-c/mv^2) -dmdt]^T\n",
" '''\n",
" dstate = np.zeros(np.shape(state))\n",
" \n",
" dstate[0] = state[1]\n",
" dstate[1] = u*dmdt/state[2]\n",
" dstate[2] = - dmdt\n",
" return dstate\n",
"\n",
"def rk2_step(state, rhs, dt):\n",
" '''Update a state to the next time increment using modified Euler's method.\n",
" \n",
" Arguments\n",
" ---------\n",
" state : array of dependent variables\n",
" rhs : function that computes the RHS of the DiffEq\n",
" dt : float, time increment\n",
" \n",
" Returns\n",
" -------\n",
" next_state : array, updated after one time increment'''\n",
" \n",
" mid_state = state + rhs(state) * dt*0.5 \n",
" next_state = state + rhs(mid_state)*dt\n",
" \n",
" return next_state\n",
"\n",
"def heun_step(state,rhs,dt,etol=0.000001,maxiters = 100):\n",
" '''Update a state to the next time increment using the implicit Heun's method.\n",
" \n",
" Arguments\n",
" ---------\n",
" state : array of dependent variables\n",
" rhs : function that computes the RHS of the DiffEq\n",
" dt : float, time increment\n",
" etol : tolerance in error for each time step corrector\n",
" maxiters: maximum number of iterations each time step can take\n",
" \n",
" Returns\n",
" -------\n",
" next_state : array, updated after one time increment'''\n",
" e=1\n",
" eps=np.finfo('float64').eps\n",
" next_state = state + rhs(state)*dt\n",
" ################### New iterative correction #########################\n",
" for n in range(0,maxiters):\n",
" next_state_old = next_state\n",
" next_state = state + (rhs(state)+rhs(next_state))/2*dt\n",
" e=np.sum(np.abs(next_state-next_state_old)/np.abs(next_state+eps))\n",
" if e<etol:\n",
" break\n",
" ############### end of iterative correction #########################\n",
" return next_state\n",
"\n",
"#initialize solution array\n",
"N = 300 \n",
"\n",
"num_rk2 = np.zeros([N,3])\n",
"num_heun = np.zeros([N,3])\n",
" \n",
"#Set intial conditions\n",
"y0 = 0\n",
"v0 = 0\n",
"m0 = 0.25\n",
"dt = 0.05\n",
"\n",
"num_rk2[0,0] = y0\n",
"num_rk2[0,1] = v0\n",
"num_rk2[0,2] = m0\n",
"num_heun[0,0] = y0\n",
"num_heun[0,1] = v0\n",
"num_heun[0,2] = m0\n",
"\n",
" \n",
"for i in range(N-1):\n",
" num_rk2[i+1] = rk2_step(num_rk2[i], simplerocket, dt)\n",
" if num_rk2[i+1,2] <= 0.05 :\n",
" break\n",
"for i in range(N-1):\n",
" num_heun[i+1] = heun_step(num_heun[i], simplerocket, dt)\n",
" if num_heun[i+1,2] <= 0.05 :\n",
" break\n",
"\n",
"#for i in range(N-1):\n",
"# if num_rk2[i,0] == 0.0 and num_rk2[i,1] == 0.0 and num_rk2[i,2] == 0.0:\n",
"# plt.plot(num_rk2[i,0],num_rk2[i,2]/m0,'s-')#,label='explicit RK2')\n",
"# num_rk2 = np.delete(num_rk2, i, axis=0)\n",
"# print(len(num_rk2), i)\n",
" \n",
"#print(num_rk2)\n",
"#print(num_heun)\n",
"\n",
"#print(np.nonzero(num_rk2))\n",
"\n",
"plt.plot(num_heun[:,2]/m0, num_heun[:,1], 'o-',label='implicit Heun')\n",
"plt.plot(num_rk2[:,2]/m0, num_rk2[:,1], 's-',label='explicit RK2')\n",
"plt.plot(np.exp(-(num_rk2[:,1]-v0)/250),num_rk2[:,1],'b+',label='Tsiolkovsky')\n",
"plt.ylim(0,420)\n",
"plt.legend();\n",
"plt.xlabel('m/m0')\n",
"plt.ylabel('velocity (m/s)')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"2. You should have a converged solution for integrating `simplerocket`. Now, create a more relastic function, `rocket` that incorporates gravity and drag and returns the velocity, $v$, the acceleration, $a$, and the mass rate change $\\frac{dm}{dt}$, as a function of the $state = [position,~velocity,~mass] = [y,~v,~m]$ using eqn (1). Where the mass rate change $\\frac{dm}{dt}$ and the propellent speed $u$ are constants. The average velocity of gun powder propellent used in firework rockets is $u=250$ m/s [3,4]. \n",
"\n",
"$\\frac{d~state}{dt} = f(state)$\n",
"\n",
"$\\left[\\begin{array}{c} v\\\\a\\\\ \\frac{dm}{dt} \\end{array}\\right] = \n",
"\\left[\\begin{array}{c} v\\\\ \\frac{u}{m}\\frac{dm}{dt}-g-\\frac{c}{m}v^2 \\\\ \\frac{dm}{dt} \\end{array}\\right]$\n",
"\n",
"Use [two integration methods](../notebooks/03_Get_Oscillations.ipynb) to integrate the `rocket` function, one explicit method and one implicit method. Demonstrate that the solutions converge to equation (2.b) the Tsiolkovsky equation. Use an initial state of y=0 m, v=0 m/s, and m=0.25 kg. \n",
"\n",
"Integrate the function until mass, $m_{f}=0.05~kg$, using a mass rate change of $\\frac{dm}{dt}=0.05$ kg/s, . \n",
"\n",
"Compare solutions between the `simplerocket` and `rocket` integration, what is the height reached when the mass reaches $m_{f} = 0.05~kg?$\n"
]
},
{
"cell_type": "code",
"execution_count": 94,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"When the mass reaches m = 0.05kg, the height is 425.5m\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3iUVfbA8e9J6EhRigqSVaQpPQlFRIqi2CiuWFAIoIDoimtDsRLWAou6NhDEXQEVUVkUbLsqq4gKqGBCEZEiIE1BFH5IC0nO7495Z5y0ySRT3pnkfJ4nz2TeNudFnMM99773iqpijDHG+EtwOwBjjDGxx5KDMcaYAiw5GGOMKcCSgzHGmAIsORhjjCmggtsBhEPdunX11FNPdTsMY4yJKytWrPhFVesVtq9MJIdTTz2V5cuXux2GMcbEFRHZWtQ+KysZY4wpwJKDMcaYAiw5GGOMKcCSgzHGmALKRIe0MfHu2LFjbN++nSNHjrgdiiljEhMTqV27NnXr1iUhIfj2gCUHY2LA9u3bqVGjBqeeeioi4nY4poxQVY4dO8bPP//M9u3bSUpKCvrcqJSVRGSciKiItHLeNxORpSKy3nlt6ndskfuMKauOHDlCnTp1LDGYsBIRKlWqRMOGDTl48GCJzo14chCRZKAz8KPf5mnAFFVtBkwBng9ynzFlliUGEyklKSf5zolAHD4iUhnPF/xNgDrb6gPJwBznsDlAsojUC7QvknEaY4zJK9Ith78Br6jqZr9tjYAdqpoD4LzudLYH2peHiIwUkeUisnzPnj0Rvg1jYs/8jB2cPfFjThv7HmdP/Jj5GTvCdu127dpx+PDhsF1vy5Yt1K1bt0TXHz58OJ999hkA8+fP56uvviry2KFDhzJ58uQ82+68807S09NLH3QYFRZfrItYh7SInAV0AMZG4vqqOh2YDpCammrL2ZlyZX7GDu55czWHj+UAsGPfYe55czUA/ds3DPn6mZmZIV8j1Ov/85//9P0+f/58UlNT6dixYyTDirrs7GwqVIjNcUGRbDl0B1oAm0VkC3AK8AFwOtBQRBIBnNcGwDbnp6h9xhjHYx9870sMXoeP5fDYB9+H5foiwu+//w545i67//77Oeuss0hKSuLVV1/lqaeeomPHjjRp0sT3r3tv6+DOO++kY8eOtG7d2rcv0PW/++47LrjgAtq0aUPr1q2ZNWsWAD169ODdd9/lgw8+4O2332bixIm0a9eOl156qVT3NGnSJDp27EhycjJ9+vThp59+Agr+q97//dChQxk1ahTnnnsuTZs2JS0tjWCWVt6xYwfnnXcebdu2pX///vzyyy95rj969GguvPBCOnXqBMC1115LamoqrVu35rLLLuO3337zHX/ffffRpEkTOnXqxN13301qamqp7r+kIpayVHUiMNH73kkQl6rqGhG5CRgIvOK8ZqjqHue4zKL2GVNenDr2vVKdt2Pf4WLP3TLxkhJf9+jRoyxdupSvv/6aHj16MGnSJL766iveeOMN7rnnHj7//HMA9u7dS5s2bXj88cf59NNPGThwIJs2bSryutnZ2fTr149HHnmEK664wncNf71796Zv376kpqZy8803F3mtiRMn5mlt7Ny5k5tuugmAV155hY0bN7Js2TISEhKYOnUqd9xxB7Nnzy723tesWcPChQtJSEigffv2LFy4kPPPPz/gObfccgvdunVj3Lhx/PDDD7Rt25YLL7zQt3/p0qV8+umnVK9eHYCnn37aV3a7//77+fvf/87EiRN55513ePfdd1m5ciVVq1ZlwIABxcYbLm61Z0YBs0TkQeA3IC3IfcaUC8V9gZ898WN27CtYs29YuypfjD037PFcddVVACQnJ3Po0CHf+5SUFDZu3Og7rlKlSgwaNAiA7t27U7VqVb7//ntq1qxZ6HW///57srOzfYkBoE6dOqWKcezYsXmSx5133un7/e2332b58uUkJycDnqRUq1atoK7bv39/qlSpAnjuf9OmTcUmh08++YRnnnkGgMaNG3Peeefl2T9gwABfYgB46aWXmD17NllZWRw8eJBmzZr5rnPllVf6jh0yZAgPPfRQUHGHKmrJQVVP9ft9HdCpiOOK3GeM8RjTu3mePgeAqhUTGdO7eUQ+z/vlmJiYWOB9dnZ2keepasAhusGUaMJBVbn//vu57rrrCuyrUKECubm5vvf5n1L33isUf7/BOu6443y/f/bZZ0ydOpUlS5ZQr149Xn31VaZPn+6L260hzja3kjFxqH/7hkz4c2sa1q6K4GkxTPhz67B0RociKyuLV199FfB86R05coTmzYtOWC1atKBChQrMnTvXty1/WQmgZs2a7N+/v9Rx9e3bl+eee85Xyz969CgrV64E4PTTT+frr78GYNeuXXzyySdBXfOtt94iLa3wwsa5557LjBkzANi8eTP/+9//irzOvn37qFWrFnXq1OHo0aO8+OKLvn09e/Zk7ty5HDp0iNzcXF5++eWgYguH2OwmN8YUq3/7hq4ng/zq1KnDhg0b6NSpE4cOHWLOnDlUqlSpyOMrVKjAggULuPnmm/nb3/5GQkICd955J4MHD85z3ODBgxk6dChz587l9ttvL/JLuSiDBw/ml19+oXv37gDk5uZy00030bZtW0aOHMmAAQNo27YtzZo183USF2fjxo1Flsuefvpp0tLSmDt3Ls2bNw9Yhrrooot45ZVXaNGiBaeccgqpqam+Ybt9+/ZlyZIltG3bloYNG9K5c+c8ndWRJNFq1kVSamqq2kpwJp599913nHHGGW6HEZItW7aQmpqaZ2ROWXb55ZfzxBNPEOklig8cOECNGjXIzc1l+PDhNGjQgIcffrjE1yns75iIrFDVQoc/WcvBGGNKYd68eVH5nLS0NLZs2cLhw4dJSUnhrrvuisrnWnIwxoTFqaeeWm5aDdH01ltvufK51iFtjDGmAEsOxhhjCrDkYIwxpgBLDsYYYwqw5GCMMaYASw7GxKPHmkJ6rYI/j8Xmqrr+M51OmzaNJ598MuDxy5cv59prrwU8TxBPmjSpyGO3bNlChQoVaNeuHW3atCE1NZUlS5b49qenp+eZZ2n69OmcfvrpbNq0iQULFpCSkkKrVq1o2bIlTzzxRKnub+bMmVGdFC8abCirMfHo4O6SbY8ho0aNKvaY1NRU34yp3uQQaHx/7dq1fWtETJ06leuvv57vvvuuwHGPPfYYM2bMYPHixTRs2JBffvmFd955hwYNGrB//35SUlLo2LEj55xzTinvrqCcnBzfnFTxxFoOxpgCvvzyS3r27ElKSgopKSm8955nGvBFixbRtGlT3zxHQ4cOZezYsb7fR4wYQa9evWjevDkjRowgKyurwLXz/0t+woQJtG7dmrZt29KlSxdyc3NZtGiRb92Cv/zlL+zbt4927drRpUuXYmPv0aMHP/74Y4Ht999/P6+//rovMQB06tSJBg0aAFCrVi3OOOMMtm7dWuxnZGVlccMNN9C8eXPOPffcPKvUzZw5kwsvvJDBgweTkpLC6tWreeKJJ+jQoQPt27fnrLPOyrPY0bx582jRogXt27fn0UcfzbPWhZus5WBMLEoPbjrpUp2bHngCu3379jFq1Cjef/99Tj75ZHbt2kWHDh1Ys2YNPXr0IC0tjeuvv56+ffuyfv36PGsofPnllyxZsoQqVapw8cUXM3369IBrMMyaNYu3336bL774gpo1a7J3714SEvL+m3XKlCmkpqYGvTrdW2+9xdVXX51n28yZMznppJNYsmRJkfMhrVu3jmXLlvH8888X+xnPP/88mzdvZs2aNRw7doxu3brlmUbj888/Z+XKlZx++ukANGzYkDvuuAOAhQsXMmrUKJYtW8bu3bsZOXIky5Yto2nTpsWW26LJkoMxsaiYL/CACaC4c4uxZMkSNm/ezEUXXeTbJiJs3LiR1NRU7rvvPnr16sUdd9zBihUr8ixzedVVV/mmox4yZAjz5s0LmBzeffddbrzxRt8XdmnXcvC2LPbs2UN2djbLli3Ls79jx45kZGTw3//+lyuvvLLA+bt27aJfv35MmTLF15II5JNPPmHIkCFUrFiRihUrMmjQIN+CRwBdu3b1JQaAFStW8Oijj/Lrr7+SkJDA+vXrAVi2bBnJyck0berpK7ruuuu4/fbbS/VnEG6WHIwxeagqbdq0YfHixYXu379/Pz/++COVK1dm7969JCUlFXmd4tYiCNfEn94+h5ycHG677TauvvpqvvzyS9/+M888k4kTJ9K7d2+APAli9+7d9OrVizFjxhSaOEoTt/96DVlZWQwYMIDFixeTnJzMzp07fWUtN9drKE5E+xxEZL6IrBSRDBH5TETaOdu3iMg6Ecl0fnr7ndNMRJaKyHrnNTaHXxjjpur1S7a9BLp06cKGDRvyrGvw9ddf+74Qhw0bxvDhw5k1axYDBw7kwIEDvuPmzp3LwYMHyc7O5pVXXqFnz54BP6tPnz5MnTrVd42i1nI4dOhQUIvsJCYmMmnSJHbt2sWCBQvy7GvTpg0ffPABt9xyC2+88Ybv884//3xuvvlmhg8fnuf4HTt20KJFi0I/57zzzuPll18mOzubw4cP+9awKMyRI0fIzs6mUaNGADz33HO+fZ07d2bFihW+1fRmzpxZ7D1GS6RbDkNUdT+AiPQDXgSSnX0DVHVNIedMA6ao6isiMgh4Hgj/uofGxLMxGyJ26eOPP563336bMWPGcOutt5KVlUXjxo155513ePrppzl8+DB33303IsIVV1zByJEjmTNnDgDdunWjf//+/Pjjj3Tr1o2RI0cG/Ky0tDR27NhB586dqVChAjVq1CjQYjnhhBO49tprad26Nccff3yeYaqFqVKlCg8//DDjx4+nX79+efZ5E0Tv3r1RVZYvX8769et5/vnnfX0Nf/3rXxk2bBg7d+7MUzLzN3LkSFatWkXLli055ZRT6N69O5s3by702Jo1a/K3v/2NDh06kJSUlKdcd+KJJzJt2jQuueQS6tatS58+fahYsSLVqlULeI/RELX1HEQkDbhFVVNFZAtwaf7kICL1gfVAHVXNEZFEYC/QVFX3FHVtW8/BxLuysJ7D0KFDSU1NDdjHEE/+8Y9/UL9+fd+a2JHiXa8BYMaMGfzrX//K038RLjG3noOI/BO4ABDgQr9ds8VTbPscuFdV9wGNgB2qmgPgJIidzvYik4MxxoRbtDqGn3nmGebOnUt2djYnnHACL7zwQlQ+tzjRbDkMBgaq6sUi0khVt4lIZeApoIaqDhKRFOAlVW3pd95aYJCqfpPveiOBkQBJSUkpwYxNNiZWlYWWg4ltJW05RO0hOFV9GegpInVUdZuz7SjwHHC2c9g2oKFTTsJ5beBsz3+96aqaqqqp9erVi8o9GGNMeRGx5CAix4lII7/3fYBfgSMiUsvZJsDVQCaAqu52fh/onDYQyAjU32CMMSb8ItnnUB2YKyLVgRw8iaEPcCIwz2kVJAJrgZv8zhsFzBKRB4HfgLQIxmiMMaYQEUsOqvoz0LmI3e0DnLcO6BSRoIwxxgTFJt4zJs6lp4f3ep06daJdu3aceeaZvqmw27Vrx7Bhw4o8Z8qUKTzzzDPFXnvQoEFMmzYN8EyE5520L1I2btzISSedFNHPCFV2djYiwpEjR9wOJQ+bPsOYODd+fHgThHfaiS1btgQ94d1f/vKX8AVgfLKzs4t8EC/SrOVgjAnad999R+fOnWnbti2tWrXiqaeeAvK2ArKzs7ntttto2bIlrVq14q677iInJyfgdTMzM2ndurXv4a8ZM2bQunVr2rRpw+WXX86ePXtQVU477TS+/fZb33lPPvkkI0aMICcnhxtuuIEWLVrQtm1bunXrVuAzjhw5woABA7jrrrv44osvaNeuXZ79LVu2zDP1dmEWLlxISkoKI0aMoE2bNrRv3943iR7Ao48+SqtWrWjVqhXXX389hw4dKvQ6c+fOpXnz5nTp0oUJEyb4tntbEY8//jjdu3fnkUceITMzk65du5KcnEzLli159tlnfcdv27aNnj170rJlS/r168cll1zia5mFypKDMXEoPR1EPD/wx+/hLjHlN3nyZPr168fKlStZs2YNQ4YMKXDM1KlTWbt2LZmZmaxYsYKvvvqKF198schrfvjhh6SlpTF37ly6du3KypUreeCBB/joo49YtWoVzZo149Zbb0VESEtLY9asWb5zZ82axbBhw/jmm29YvHgxa9euZeXKlQXmVdq7dy8XXHABPXv2ZNKkSZx99tlUrFiRL774AvDMslqtWjU6duxY7J/B6tWrGT16NKtWraJ///488sgjALzzzju89tprLF26lNWrV3PkyBHfPn+7du1i1KhRvPvuuyxZsqTQhYBEhE8//ZRx48bRuHFjPv74Y7755huWLVvG5MmTfQnp5ptvpnfv3nz77bc89dRTRU6WWBqWHIyJQ+npoOr5gT9+j3Ry6NatGy+88AIPPvggn3zyCbVr1y5wzMKFCxk2bBgVK1akcuXKDB06lIULFxZ6vf/85z/ceeedfPjhh75J7j7++GMuvfRSX1/BDTfc4Dt/6NChzJ49m5ycHDIyMjh8+DBdunShSZMmHD16lBEjRvDKK6/kmen00KFDnH322dx22215yl+jR4/2TYI3ZcqUoEtjZ555Jm3atAE8E+dt2rTJd9/XXHMNNWrUQEQYMWJEofe9dOlSOnbs6Jumu7D5p/yT7sGDBxk6dCitW7ema9eu/PTTT6xatQrwJDVvX9Bpp51Gjx49grqHYFhyMMYE7aqrrmLx4sWcdtppPPLII4V2Uhc2DXVR01I3b96crKwsVqxYEdT5p512Gk2aNOHDDz9k5syZDB06FPBMFrh27VoGDBhARkYGLVu2ZPduz5KpVapUoVOnTixYsCBPeevqq6/ms88+IyMjg88//7zAAkFFqVKliu/3xMRE32yxwd53MLNS+E/5PXbsWJKSksjMzGTlypWkpKTk6byO1JTflhyMiXPjxkXvszZs2ECDBg0YNmwYDzzwQKE1+vPPP58ZM2aQnZ1NVlYWL730Er169Sr0eo0bN+bDDz9kzJgx/Pvf/wY802G/8847vi/3F154Ic/5Q4cO5YUXXuD1118nLc3zGNTu3bs5cuQIF110EZMmTaJatWps2bIFgISEBGbOnEmVKlW45pprOHbsGACVKlViyJAh9O3bl7S0NN+Xfk5ODi1atODnn38u0Z/N+eefz5w5c/j9999RVf71r38Vet9dunTh66+/9rU4/FfSK8y+ffto1KgRiYmJrFy50lcKA8+SqN5pvrdu3cqiRYtKFHMgNlrJmDgX6VKSv9dee43XXnuNSpUqISK+Dml/N954Iz/88IOvw/eiiy4KOAw2KSmJhQsX0rt3bw4fPszgwYN56KGHOO+88xARmjRpkmfpziuuuILRo0fTtWtX36I5W7du5YYbbiA7O5ucnBz69+9Phw4dfF/AIsK0adO49dZbufzyy5k7dy6VK1dm+PDhPProo4waNcp3/d27d/Pbb78VWjILpE+fPqxevZrOnT2Pd3Xq1Il77723wHEnn3wyU6dO5eKLL6Zu3boMGDAg4HUffPBBX19LkyZNOOecc3z7Jk+eTFpaGq+++iotWrSgS5cu1KoVwhKzfqI28V4k2ZTdJt7ZxHvumDlzJm+99VaeDuw33niDH374IeLPYITD4cOHqVSpEomJiezYsYMOHTqwePFimjRpUuDYmJuy2xhjYlGvXr3YunUr7777bp7twS4VGgvWrVvHsGHDUFWys7N5+OGHC00MpWHJwRhTLhU1giqetG/fPqiHFEvDOqSNiRFlocRrYlNp/m5ZcjAmBiQmJvpG0RgTbocPH6ZixYolOseSgzExoHbt2vz888/k5ua6HYopQ1SVQ4cOsWPHDurXr1+ic63PwZgYULduXbZv387333/vdiimjKlYsSInnngiNWvWLNF5lhyMiQEJCQkkJSW5HYYxPlZWMsYYU0BEk4OIzBeRlSKSISKfiUg7Z3szEVkqIuud16Z+5xS5zxhjTHREuuUwRFXbqmp74HHAO2/vNGCKqjYDpgDP+50TaJ8xxpgoiGhyUNX9fm9rAbkiUh9IBuY42+cAySJSL9C+SMZpjDEmr4h3SIvIP4ELAAEuBBoBO1Q1B0BVc0Rkp7NdAuzbk++6I4GRgHXkGWNMmEW8Q1pVh6tqEnAv8FgYrztdVVNVNbVePWtYGGNMOEVttJKqvgz0BLYDDUUkEcB5bQBsc36K2meMMSZKIpYcROQ4EWnk974P8CuwG8gEBjq7BgIZqrpHVYvcF6k4jTHGFBTJPofqwFwRqQ7k4EkMfVRVRWQUMEtEHgR+A9L8zgu0zxhjTBRELDmo6s9A5yL2rQM6lXSfMcaY6LAnpI0xxhRgycEYY0wBlhyMMcYUYMnBGGNMAZYcjDHGFGDJwRhjTAEBh7KKyHHAVcB5wCnAYWAlME9Vv4x8eMYYY9xQZHIQkbHAcOAD4D/AT0AV4AzgURFJAG50nkswxhhThgRqOewHWqhqdr7t7wCTRKQ5ntlSLTkYY0wZU2RyUNWpgU5U1e8BWw3dGGPKoGI7pEXkdhGp5fz+soisE5ELIh+aMcYYtwQzWmmoqu4XkZ5AfeA64NHIhmWMMcZNwSSHHOe1JzBbVZcEeZ4xxpg4FcysrIdF5D7gWuBsERGgUmTDMsYY46agykpAPWCMqv4ENAZmRzIoY4wx7gr0nMPfgQVOGelW73ZV3QRMiEJsxhhjXBKo5bAKuE1ENorIP0Wkj4hUiVZgxhhj3FNkclDV2ap6BXAm8CZwKfCtiLwlIkNFpG6gC4tIHRF5X0S+F5FVIvKmiNRz9m1xhsRmOj+9/c5rJiJLRWS989o0PLdqjDEmWMX2Oahqlqq+r6o3qOrpwN+BFsCi4k4FJqlqc1VtA2wCJvrtH6Cq7ZyfD/y2TwOmqGozYArwfAnuxxhjTBgEPSRVRCqJSDU85aa/qWqrQMer6q+qushv0zLgT8V8Rn0gGZjjbJoDJHtbHMYYY6IjmCekLxeRbXhmZD0A/O68Bs07SR/wtt/m2U656TkRqe1sawTsUNUcAOd1p7M9/zVHishyEVm+Z8+ekoRjjDGmGMG0HB4D/gxUVNVEVU1Q1cQSfs6zeJLKZOf9OaraFugAiN/2oKnqdFVNVdXUevWsYWGMMeEUTHLYpapfq2puaT5ARB4HmgJXea+hqtuc16PAc8DZzuHbgIYikuicmwg0cLYbY4yJkmCSw7Mi8pCIJIvImd6fYC4uIo8AKUB/JxEgItX9JvIT4GogE0BVdzu/D3QuMRDIUFWrGxljTBQFM31GQ+B2YAh/zLOkeJ6ULpKItATuBdYDSzx5gM3AHcA8p1WQCKwFbvI7dRQwS0QeBH4D0oK9GWOMMeERTHK4BWiiqrtKcmFV/RZPf0Jh2gc4bx3QqSSfZYwxJryCKSttLWliMMYYE9+CaTl8JSJzgLnAEe9GVX0/YlEZY4xxVTDJIcV5He23TQFLDsYYU0YVmxxUtWc0AjER9lhTOLi74Pbq9WHMhujH44L0dM+PMaZ4gabsbqCqOwOdLCInOWs8mFhXWGLwbk+v5fm9jCeK8eM9r5YgjCleoA7p10VkqoicIyKVvRtFJMmZumIJfzy8ZsqCohJIGeJNEMaYwAIlh27AJ8A44DcR+U1EDgOLgTOAa1R1XhRiNNGUXsvz81jZmCk9PR1EPD9eItCjh1sRGRMfiiwrqaoCbwBviEgFoC5wWFX3Rys446IyUm7ylpDytxg+/dSTIBYtinJAxsSJoKbsVtVsVf3JEkM5FeflpvR0UP3jvff3Tz91JRxj4kLQ6zmY+JZTtU5oFygD5abu3T2vVmIypniWHMqJ5VXPZmmDNEjf7ykTlZa33BSHyWLRoj8ShD9vickY84dgHoIzcW7HD9/S7NdPSBi9wrPBv//A269QWnFWcvL2MXhbD6qe363EZExewawEt05EbhaRGtEIyITfzgXjWZd0DbXqnFhwZyitCK84bEVYicmYwIIpKw0E2gEbneceWkc4JhNGW9d9w+n7l9Lq8rGFHzBmg6fUFGq5CfKWnGI8UViJyZjAik0OqpqhqsOB5sAG4H0RWSwif454dCZke99NZ33jodSodULxB/snilDFQblp0SIbxWRMUUrSId0J6AEcAv4LjBKR1yMRlAmPTauW0Oj3VbS9/K6Sn1yOyk1WYjKmoGD6HO4QkQ3AbcDzQAtVfVRVL+CPGVtNDDrwn/Fsaj6SqtVL0V3k34oIxwinGBaoxGTzMJnyKpjRSo2BPs4KbfldVdRJIlIHeBk4HTgKbARuUNU9ItIMmAXUAfYCaaq6wTmvyH0meOuW/4+TDm/i+MveCs8FQx3hFONPWxc1ismY8irYleDyJAYRuQtAVVcEOE+BSaraXFXbAJuAic6+acAUVW0GTMHTIiGIfSZIxz56iB9b/4XKVaqF/+Lhek4iBstN+UtM48dbicmUT8Ekh6uD3JaHqv6qqov8Ni0D/iQi9YFkYI6zfQ6QLCL1Au0LIk7j+PaL9zg+axft+9wUmQ8I1winGCw3LVoE48YV7Ki2TmpT3gRaz+F84AKggYhM8ttV4pqCiCQANwJvA42AHaqaA6CqOSKy09kuAfbtyXfNkcBIgKSkpJKGVGZpbi4Jix7mp/a3ckqlysWfEKoyWm7K3zkNf/Q/WD+EKQ8C9TlkAb/jKQ8d9Nu+C5hQws951rnWZKB9Cc8tlKpOB6YDpKamajGHlxurP32TmjkHaHbxiOh/ePX6pW8NxNAssN4V43r0yNtiyD+zqyUJU5aJauDvVRFppaprSv0BIo8DbfB0ah91SkfrgTpOyyART8dzUzwth0L3qeqeIj6C1NRUXb58eWlDLDM0N5eNj3TgQIfRJF841N1gQp2WIxzPWoSByB+d04W9GhPPRGSFqqYWti9QWekKVZ0LdBORbvn3q+pzQXzwI3iGu16iqked83aLSCaeJ69fcV4zvF/+gfaZwDI+mk1NlHbnD3Y7lNBaERAzrYhx4/5oIXjLS1ZmMuVBoLJSK2Au0KGQfcX+m0lEWgL34mkJLBHP/1GbVfUyYBQwS0QeBPglk8oAABeYSURBVH4D0vxODbTPFCEnO5sTvpzE/q73kZCY6HY4eb/QH2saWrnJRfm/+P1LS/6/W4IwZU2xZaV4YGUlWP7udI7LfJHm9y5BEmJ4JvZQyk0x0GldWHnJv3VhTDwJVFYK5gnpu0XkBL/3dURkTDgDNKHJPpbFiSueJLvHvbGdGCDun5EorMxkz0KYsiiYJ6QHqurfvW9Uda+IXAM8FrmwTElkvDOVahXr0aprX7dDKV6cl5v8Wwjp6fZEtSm7gkkOhf2Vt0WCYsTRI4c4ZdWz7L94qtuhlFycPyNhz0KYsiyYGsQGEbldPBJE5A488ySZGJC54Fl2Vz2NFh3PdzuU0MTZk9bp6Z7WQv4J+8aP9/xYcjDxLpjnHBrgGVbaBc8opSXAIFXdFfnwglNeO6QPHzzAgcfasL//SzRtd47b4YRPnHVaF/UMhPdhOmNiVamec/BS1Z3AuSJS3Xl/sJhTTJSsfOsJKlVvSXJZSgwQ+pPWURboWQiwBGHiU1BDWUWkN9ALT8vhI1X9KNKBlUR5bDn8/n+/cfQfbTlw1Zucekahib9sCKXTOsqtCG9LwVoQJl6E1HJwpudO44+ZUv8hIrNU9fEwxmhKaPW8iVSs2YHUspwYILRO6yjP1+SfBAprQXiPMSYeBDPqaBBwlqoeABCRZ4AvAEsOLtn/6x5abJ3N74P+43Yo0RUH5Sb/0Ur+LQjve0sOJl4ENZTVmxgAVPWAiI3odtPaeY+QeHw3OjZp7XYo0RUnQ18DtSBsqKuJF8Ekh69FZAbwAp4+h+FA+Srwx5Bfd+/gjB1zOXLdIrdDcVeMtyL8v/xtPiYTj4IZylodeBA4D88DcR8BD8XSqKXy1CG9bOooJCeLTje/6HYosSPGh77aUFcTq0IdynoQuDvsUZkS271jM2f8/DbHRi51O5TYEuOLDBU31BUsSZjYU2TLQUQCLkAczHoO0VJeWg5fTh6GVqhK51Ex80cfe0IZ+goRXWSoqKGutnCQcUtpWw6FrePgZX+Vo2znlu9p/suH6F/KfhIMSQx3WltHtYknRSYHVR0WzUBMYNvnp5N7ylV0rney26HEjxjstLaOahMvgumQrgbcAzRW1WtFpAXQQlXnRyPAYJT1stK2DSs5bvalJPw1g1rH13U7nPgUg53WgcpL1lltoiGkxX6AqUBFoJ3zfjswLogPfVxENouIikgrv+1bRGSdiGQ6P7399jUTkaUist55dWdFlxjz89vjWXfaYEsMoYjBRYaK6qgW+aMlYQnCuCWY5xxaq+oQ75e4qv4uIsEklfnA08BnhewboKprCtk+DZiiqq+IyCDgeeDcID6rzNr87ZecemAFVUfa0NWQxOAiQ4UtHORtOfivMmcJwrghmOSQ5f9GRKoQRItDVT93jg8qEBGpDyQD3oUJ5gCTRaSequ4J6iJl0L73xvNz0+vpXKO226GUHTHWae3/5W8LCJlYEUwLYLGI3AtUFpEewBvAghA/d7aIrBKR50TE+63XCNihqjkAzutOZ3sBIjJSRJaLyPI9e8pm7tiQsZiGh76j3WW3ux1K2RVDiwyNG+dpOYzLV7S1BYSMG4LpkK4I3AX0xfOE9NvARFXNDuoDRLYAl3rLSCLSSFW3iUhl4CmghqoOEpEU4CVVbel37lo8Cwt9E+gzymqH9KqJ53G48YV0unKM26GUDzHUaW1PVZtoCOkJaaCCqj4CPBKOYFR1m/N6VESew5NsALYBDUUkUVVzRCQRaOBsL3fWLvsvdY/+SN1+o90OpfyIoaGvtoCQcVswyWGbiMwHXlTVJaF8mDNPUwVV3e/M7Ho1kAmgqrtFJBMYiGdZ0oFARnnsb9DcXPj4Yba3uYUGlau4HU75EWqntX/LI8SWRFGd1f7TflsrwkRSMGWlE4BrgKFATWAmnvLP9mLOewb4M3AS8AuwF+gDzAMSnZ+1wC3e9aidZyhmAccDvwFpqvp9cTdR1spKqxcvoNaie2hwTyYVKlZyOxwTSrkJwpIo/B+S82elJhOKQGWloJYJ9btQS+BOPP0AFcMUX8jKUnLQ3Fw2PNqJA+1vIOWS4W6HYyD0+Zog5Dmb8s/LBDY3kwldqH0OOM81XIyn9dANT+vBRMDKj1+nVm4W7S+02UtiRqhDX/3PK2Urwoa7mmgrdiiriPwDz1PRf8XzYNufVHVEpAMrj3Jzcqix9O/s6zyGhMREt8MxhQll6CuE3AKx4a4mWoLpc7gfmOUdZRSLykpZacX7/6LGN9Noeu+XSEIwj6AY17k4/NWGu5pQhbrYz8PhD8nkl30si3rLn2B/94ctMcQTFxcaCma4qyUKU1ol6pCOVWWh5fD1/ClU//ZVzrjnM0sO8cqlhYaKWkQIrLPaBBbqrKwmwo5lHaXhyqfh3PstMcSzMRs8X/ClHZlUyplfi1pEyB6aM6EIarSSiaxvFjxL1cqn0Oasi9wOxYRLlJ+29h+tlH+uS0sSpjSsrOSyI4cPsv/vrdnf90WaJfdwOxwTCVHutC7uWQjrhzBeVlaKYZlvPcnOas0tMZRlUV5oyDvM1RYSMqGwloOLDv2+n0OPt+H/BrxB41ad3A7HREOUO63tyWoTSNimz4hV8Zocls66j0q/fEvKHTGzHLeJpiiWm4pacyt/K8OUL1ZWikH/t28vzTe/RN1L090OxbgliuWmYJ+stiRhvKzl4JKl/7qTCge20+HW19wOxcSCUGd+LUG5qagnq/33mfLBWg4xZt8vP9Fi22s07JfudigmVoQ6Z1MJWxH2XIQpjj3n4ILv5j1EQp3z6HRaC7dDMbEi1IWGvII4r7CFhPxZkjBgZaWo++WnH6kw7SyODl/Miaec7nY4JtZFodPanosov6ysFEM2zRvPuhMvtcRgghOFTmt7LsIUJmLJQUQeF5HNIqIi0spvezMRWSoi653XpsHsKwt++nEDzff8l6Z/fsDtUEy88J+vKdREUQT/kUreUU3gefUvLBS1VKkpmyLZcpiPZ9W4rfm2TwOmqGozYArwfJD74t6P88fzXYPLqXPiKW6HYuJRFCb2y7/iXP6pwG3Ia/kRseSgqp/nXyBIROoDycAcZ9McIFlE6gXaF6kYo2n7xjU0/XURZw6wVoMJgwiXmwI9F2GlpvIh2n0OjYAdqpoD4LzudLYH2leAiIwUkeUisnzPnj1RCT4UP72dzro/XUutE8pErjNuC7UV4VVEuSl/C8FKTeVP3HZIq+p0VU1V1dR69WL7C3frdyto/H9f0urPd7sdiimLIvyMhH/rwUpN5Ue0k8M2oKGIJAI4rw2c7YH2xbW976Wz/vRh1Kh1gtuhmLIonJ3WhSQK75e+lZrKl6gmB1XdDWQCA51NA4EMVd0TaF80Ywy3jSu/oNHvq2n75zFuh2LKgwiWm6zUVL5EcijrMyKyHTgFWCgi3zq7RgGjRWQ9MNp5TxD74tLB/47nhxY3ULV6DbdDMeVNBMtNVmoq++wJ6Qha9/VCar83iuPHrqJylWpuh2PKs1DXkfDK99S19+np9PSiWwyq9pR1rLL1HFyyZkJ3DjW7jI6X3+p2KMb8IUIzwAaa5dVme41NNn2GC9Z8/ja1s36mfZ8b3Q7FmLwiVG6yUlPZYskhAjQ3lwqfTuCn9rdSsVJlt8MxJq8IjW6yUU1li5WVImDlJ3Op9dl4Gt2bSWIFmxXdxIkIlJus1BTbrKwURZqbS/XPJ/BbpzGWGEx8iUC5yUpN8cuSQ5hlfvQyAO3OH+xyJMaUUATKTenVPYnCSk3xx8pKYZSTnc22R9uzv+sDtD33SrfDMSY8wlxuslJT7LCyUpRkvP9PjiRUp02PAW6HYkz4hKvc5JScrNQUHyw5hMmxrKOclPEU2T3vQxLsj9WUIeEqNwEc3E06nkQxrtczQZWawBKFG6ysFCZfzXuKauvfotU9n7oahzFRE2q5CXwlp0Clpvy/m/CxslKEHT1yiKTVk6nQyxbyMeVIqK0I8JWbxvV6xrcpf6nJ+ztYCyKaLDmEQeb8p/m5amNadOjldijGRE8Yy03pZz/gSRLdJ6CTmhbaSvCWmqxPIjqsrBSiwwcP8Ptjrdl/2Ss0advVlRiMiTlhGOFU3GgmKzWFzspKEbTyrcfZVr21JQZj/IVhhNO47hPyTMthpaboskd4Q3Bg/6803TiDA1e96XYoxsQWv2m9SztdeHqPiXAQz+im7mNJv+RF5K4NeY7JnyRsavDwsbJSCJbOuJuKv20i9fZ/R/2zjYlLIZabZPx+dFwtqF4fuWtDkaUmSxLBsbJSBOzf+zMtts7mpL7pbodiTPwIsdw0rvsEANLfuw4outRkS5WGzrXkICJbRGSdiGQ6P72d7c1EZKmIrHdeC65RGAPWvvkoG47vzilNWrkdijHxI8QRTuk9Jvpex3Wf4GlF5GNPWoeH2y2HAarazvn5wNk2DZiiqs2AKcDz7oVXuL0/b+eMHf8m6bJ0t0MxJn75J4pS8CYKAB1Xy9eq8LInrUPjdnLIQ0TqA8nAHGfTHCBZROq5F1VBG958iO/rXchJSTHZqDEm/oRQbvKVmpxk4W1NFPa8hM3+GjzXOqRFZAuwHxDgc+Be4HTgJVVt6XfcWmCQqn6T7/yRwEiApKSklK1bt0Yl7t07NlP5hbM5NnIpdRv8KSqfaUy5UsrRTQDpi8aS3mMiMr7o1ojN/vqHWO2QPkdV2wId8CSIySU5WVWnq2qqqqbWqxe9hsXmN9P57qT+lhiMiZQQ+iW8rQdvf4T1SZSea8lBVbc5r0eB54CzgW1AQxFJBHBeGzjbXbdz8zqa711IiwEPuh2KMeVDKfsl/PsjILg+CUsSebmSHESkuojUcn4X4GogU1V3A5nAQOfQgUCGqu5xI878dixI57tGV1O77kluh2JM+VPKfoki+yTG1UIn/dFvaB3XebnVcjgRWCQiq4A1QDPgJmffKGC0iKwHRjvvXffj+kxO3/cFLS+/1+1QjCmfSllu8m9F+LceZPx+3xPX+UtNlijsCemgrXiiP1l1W3LWkEci+jnGmBIqRQe2t+M6fdFYxn96T6HHlIcJ/mK1Qzpu/LDmS/50IIO2l9/ldijGmPxK0S/h/zAdUGTHdXme4M+SQxD2v5/OxmbDqXZcGFa+MsZETin6JfxLTYFGOJW3tSQsORRj/TeLaHDoe9pddrvboRhjilOKfgn/4a/+/BOF9zW9uqcDuzyMcLLkUIwjH/6NLS1vokrV6m6HYowpiRImisI6rtMXjQXwPVQnd234o9T0WNMy3XFtHdIBrF36H2p/+Ffqjl1FpcpVwn59Y4wLStiBXdxT1+PGeRJDPHZcW4d0KWhuLnz8MDva3mKJwZiypIQd2MFO8Oc7Pj0cQbrPkkMR1nw2n+o5+0i+NCYeszDGREIJOrCLepjOX1nquLbkUAjNzaXyZxP4JfV2EivYSqrGlFkl6Jco6mG6QjuuqQXpteK649qSQyFW/u81KuRm0b73ULdDMcZESykSRZEd1+P3/9FH4U0Sj3lGOsVLorDkkE9uTg41l/6d/zvrbhISE90OxxjjhiATRVAr0zlJIv296/IkilhPEpYc8sn47wyyEyrR9ryr3Q7FGBMLguzALrbj+tN7/mhNHNztSRLpnvJT+vnPhDvqkFly8JN9LIt6K/7B0W73Ign2R2OMySeIDuygOq69rQmnHDV+4S2eJNFjgq/85Db7BvST8d50fq9wAq269nM7FGNMLAqi3BSo4zpga8J5X9pV8MLNkoMj6+gRGq58Gs6931oNxpjilSBRlLg1EQOtCPsWdGQseJZfKidxZucL3Q7FGBNvikkUpWlNeFsR6T0mePomopwoLDkARw79zmlrn6PKBbb8pzEmREEmioAr1Pm1LHzrTUQ5UVhyADLf+gc7qrWgWXJ3t0MxxpQl/okiX7IIuEKd3/MS4Ndx7SQK77DYSCaKmJx4T0SaAbOAOsBeIE1VNxR1fIkn3ksPsC5DCRcyN8aYUilkAsCgVqgbVwsZv9/XupDx+9BxtT07q9f3JKQgxePEe9OAKaraDJgCPO9yPMYYE16FPD8R1Ap1+VoT4DfrXxhHOsVcchCR+kAyMMfZNAdIFpF67kVljDERVEjfREmHwYZbzJWVRCQFeElVW/ptWwsMUtVv/LaNBEY6b5sD3wf7GSknJ6QUtW/FrtwVJQ46/tQFfnE7iCizey4f4v6e256Y0LZCAr4ZP3ceOJkGNXYBsGJXe1JOzmDFrvYBrrCiJN9hf1LVQv/hHbdTjqrqdGB6qNcRkeVF1dzKKrvn8sHuuXyI1D3HXFkJ2AY0FJFEAOe1gbPdGGNMFMRcclDV3UAmMNDZNBDIUNU97kVljDHlS6yWlUYBs0TkQeA3IC2CnxVyaSoO2T2XD3bP5UNE7jnmOqSNMca4L+bKSsYYY9xnycEYY0wB5SI5iEgzEVkqIuud1wKTkYhIoohMEZFNIrJRRIa7EWu4BHnPD4jItyKyUkRWiEhvN2INl2Du2e/Y5iJySEQej2aM4RbsPYvIlSKyWkTWOK8nRjvWcAny73Z9EXlPRFaJyDoReU5EYrWPNSAReVxENouIikirIo4J//eXqpb5H+BjPA/RAQwCPi7kmDTgAzwJsx6wHTjV7dgjfM+9gWrO722BfUBVt2OP5D07+xKBRcCrwONuxx2F/86pwFrgJOd9LaCK27FH+J6f8v63BSoCXwJXuh17Ke+3K9AI2AK0KuKYsH9/lfkOaWc6jvVAHVXNcZ6b2As0Vb/hsSLyHjBDVf/tvJ8MbFXVx9yIOxTB3nO+cwRPcmipqtujF214lOSeReQ+4ChwHHCcqt4Z9YDDoAR/t2cD/1PVF10KNWxKcM9PAtWAG53Xz4CbVfULF8IOCxHZAlyqqmsK2Rf276/yUFZqBOxQ1RwA53Wns91fErDV7/2PhRwTL4K9Z39pwKZ4TAyOoO5ZRNrgaTE9GfUIwy/Y/85nAo1FZLGIfCMi9zv/GIhHwd7zQ0AzYBfwE/BBPCeGIIT9+6s8JAdTDBHpjud/poHFHRvPRKQi8AIwyvvlUk5UANoA5wPdgYuAwa5GFHlXAKuAk4GGQDcRGeBuSPGlPCSHYKfj+BH4k9/7pEKOiRdBT0EiImcBrwD9VTXoyQtjUDD3fDJwOvC+00S/FRghIvH64FSw/523Av9W1aOqegBYAHSMaqThE+w9jwZmq2ququ7Hc889oxppdIX9+6vMJwcNfjqOuXi+KBKc6cH7A/OiF2n4BHvPItIBeB0YoH4z3sajYO5ZVX9U1bqqeqqqnoqn0/IFVR1Z4IJxoAR/t18FLhCPisB5wMroRRo+JbjnzcCFACJSCegFFKjVlyHh//5yuyc+Gj9ACzyjFdY7r82d7e8Dqc7vicBUYJPzM9LtuKNwz18De/D8z+b9ae127JG853zHpxP/o5WC+e+cAPwD+A741vk9we3YI3zPpwMfAavxjNSaAlRwO/ZS3u8zeEYfZePpP/m2kPsN+/dXmR+tZIwxpuTKfFnJGGNMyVlyMMYYU4AlB2OMMQVYcjDGGFOAJQdjjDEFWHIwJsxE5J8ickkJz6kmIq87M2quE5FLIxWfMcGw5GBMGIlIAp4ncf9XwlPvBA6oahOgD/BPETku3PEZEyxLDsYUw5lH/z4R+VpEfhCR80RkgohkOOsjnOF3+FlApqoeEZF0EXlNRN53WgSvi0h7EfnYmXfff8bMq4BpAKq6AViOZw4kY1xhycGY4OxT1Q7A3Xjm6flcVdsDLwH3+R3X39nvlYJniofmeJ7snYjnS78NMMRvoZqyNCuwKQMsORgTnNed128AVdX3nPcrgCZ+x10CvOf3/gNV3a+eWWBXAR+pZwK8g8D3eKZ5MCbmWHIwJjhHnNccPAsF4fe+AoCInAnsVtW9hZznPTb/e+/SlWVpVmBTBlhyMCZ8+pG3pFQSc4EbAJxSUwfgv2GKy5gSs+RgTPiEkhweA2qLyEbgXTyzah4IW2TGlJDNympMGIjIyXj6F9q4HYsx4WDJwRhjTAFWVjLGGFOAJQdjjDEFWHIwxhhTgCUHY4wxBVhyMMYYU4AlB2OMMQX8P4jIyhW1qky6AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"plt.rcParams.update({'font.size': 11})\n",
"plt.rcParams['lines.linewidth'] = 1\n",
"\n",
"def rocket(state,dmdt=0.05, u=250,c=0.18e-3):\n",
" '''Computes the right-hand side of the differential equation\n",
" for the acceleration of a rocket, with drag, in SI units.\n",
" \n",
" Arguments\n",
" ---------- \n",
" state : array of three dependent variables [y v m]^T\n",
" dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n",
" u : speed of propellent expelled (default is 250 m/s)\n",
" c : drag constant for a rocket set to 0.18e-3 kg/m\n",
" Returns\n",
" -------\n",
" derivs: array of three derivatives [v (u/m*dmdt-g-c/mv^2) -dmdt]^T\n",
" '''\n",
"\n",
" dstate = np.zeros(np.shape(state))\n",
" \n",
" dstate[0] = state[1]\n",
" dstate[1] = u*dmdt/state[2] - 9.8 - c*state[1]**2/state[2]\n",
" dstate[2] = - dmdt\n",
" return dstate\n",
"\n",
"\n",
"def simplerocket(state,dmdt=0.05, u=250):\n",
" '''Computes the right-hand side of the differential equation\n",
" for the acceleration of a rocket, without drag or gravity, in SI units.\n",
" \n",
" Arguments\n",
" ---------- \n",
" state : array of three dependent variables [y v m]^T\n",
" dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n",
" u : speed of propellent expelled (default is 250 m/s)\n",
" \n",
" Returns\n",
" -------\n",
" derivs: array of three derivatives [v (u/m*dmdt-g-c/mv^2) -dmdt]^T\n",
" '''\n",
" dstate = np.zeros(np.shape(state))\n",
" \n",
" dstate[0] = state[1]\n",
" dstate[1] = u*dmdt/state[2]\n",
" dstate[2] = - dmdt\n",
" return dstate\n",
"\n",
"def rk2_step(state, rhs, dt):\n",
" '''Update a state to the next time increment using modified Euler's method.\n",
" \n",
" Arguments\n",
" ---------\n",
" state : array of dependent variables\n",
" rhs : function that computes the RHS of the DiffEq\n",
" dt : float, time increment\n",
" \n",
" Returns\n",
" -------\n",
" next_state : array, updated after one time increment'''\n",
" \n",
" mid_state = state + rhs(state) * dt*0.5 \n",
" next_state = state + rhs(mid_state)*dt\n",
" \n",
" return next_state\n",
"\n",
"def heun_step(state,rhs,dt,etol=0.000001,maxiters = 100):\n",
" '''Update a state to the next time increment using the implicit Heun's method.\n",
" \n",
" Arguments\n",
" ---------\n",
" state : array of dependent variables\n",
" rhs : function that computes the RHS of the DiffEq\n",
" dt : float, time increment\n",
" etol : tolerance in error for each time step corrector\n",
" maxiters: maximum number of iterations each time step can take\n",
" \n",
" Returns\n",
" -------\n",
" next_state : array, updated after one time increment'''\n",
" e=1\n",
" eps=np.finfo('float64').eps\n",
" next_state = state + rhs(state)*dt\n",
" ################### New iterative correction #########################\n",
" for n in range(0,maxiters):\n",
" next_state_old = next_state\n",
" next_state = state + (rhs(state)+rhs(next_state))/2*dt\n",
" e=np.sum(np.abs(next_state-next_state_old)/np.abs(next_state+eps))\n",
" if e<etol:\n",
" break\n",
" ############### end of iterative correction #########################\n",
" return next_state\n",
"\n",
"#initialize solution array\n",
"N = 300 \n",
"\n",
"num_rk2 = np.zeros([N,3])\n",
"num_heun = np.zeros([N,3])\n",
" \n",
"#Set intial conditions\n",
"y0 = 0\n",
"v0 = 0\n",
"m0 = 0.25\n",
"dt = 0.05\n",
"\n",
"num_rk2[0,0] = y0\n",
"num_rk2[0,1] = v0\n",
"num_rk2[0,2] = m0\n",
"num_heun[0,0] = y0\n",
"num_heun[0,1] = v0\n",
"num_heun[0,2] = m0\n",
"\n",
" \n",
"for i in range(N-1):\n",
" num_rk2[i+1] = rk2_step(num_rk2[i], rocket, dt)\n",
" if num_rk2[i+1,2] <= 0.05 :\n",
" break\n",
"for i in range(N-1):\n",
" num_heun[i+1] = heun_step(num_heun[i], rocket, dt)\n",
" if num_heun[i+1,2] <= 0.05 :\n",
" break\n",
"\n",
"#for i in range(N-1):\n",
"# if num_rk2[i,0] == 0.0 and num_rk2[i,1] == 0.0 and num_rk2[i,2] == 0.0:\n",
"# plt.plot(num_rk2[i,0],num_rk2[i,2]/m0,'s-')#,label='explicit RK2')\n",
"# num_rk2 = np.delete(num_rk2, i, axis=0)\n",
"# print(len(num_rk2), i)\n",
" \n",
"#print(num_rk2)\n",
"#print(num_heun)\n",
"\n",
"\n",
"#print(np.nonzero(num_rk2))\n",
"\n",
"plt.plot(num_heun[:,2]/m0, num_heun[:,1], 'o-',label='implicit Heun, drag')\n",
"plt.plot(num_rk2[:,2]/m0, num_rk2[:,1], 's-',label='explicit RK2, drag')\n",
"plt.plot(np.exp(-(num_rk2[:,1]-v0)/250),num_rk2[:,1],'b+',label='Tsiolkovsky, no drag')\n",
"plt.ylim(0,420)\n",
"plt.legend();\n",
"plt.xlabel('m/m0')\n",
"plt.ylabel('velocity (m/s)')\n",
"\n",
"print('When the mass reaches m = 0.05kg, the height is 425.5m')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"3. Solve for the mass change rate that results in detonation at a height of 300 meters. Create a function `f_dm` that returns the final height of the firework when it reaches $m_{f}=0.05~kg$. The inputs should be \n",
"\n",
"$f_{m}= f_{m}(\\frac{dm}{dt},~parameters)$\n",
"\n",
"where $\\frac{dm}{dt}$ is the variable we are using to find a root and $parameters$ are the known values, `m0=0.25, c=0.18e-3, u=250`. When $f_{m}(\\frac{dm}{dt}) = 0$, we have found the correct root. \n",
"\n",
"Plot the height as a function of time and use a star to denote detonation at the correct height with a `'*'`-marker\n",
"\n",
"Approach the solution in two steps, use the incremental search [`incsearch`](../notebooks/04_Getting_to_the_root.ipynb) with 5-10 sub-intervals _we want to limit the number of times we call the function_. Then, use the modified secant method to find the true root of the function.\n",
"\n",
"a. Use the incremental search to find the two closest mass change rates within the interval $\\frac{dm}{dt}=0.05-0.4~kg/s.$\n",
"\n",
"b. Use the modified secant method to find the root of the function $f_{m}$.\n",
"\n",
"c. Plot your solution for the height as a function of time and indicate the detonation with a `*`-marker."
]
},
{
"cell_type": "code",
"execution_count": 235,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 0.3424583 0.34616337 1.07120461 ... 222.60168834 222.69262015\n",
" 222.78345394]\n",
"The root of the function f_m is 240.00979999999996\n",
"[[0.00000000e+00 0.00000000e+00 2.50000000e-01]\n",
" [5.05637971e-02 2.02255189e+00 2.47500000e-01]\n",
" [2.02885336e-01 4.07030968e+00 2.45000000e-01]\n",
" [4.58229953e-01 6.14347500e+00 2.42500000e-01]\n",
" [8.17872928e-01 8.24224398e+00 2.40000000e-01]\n",
" [1.28309920e+00 1.03668067e+01 2.37500000e-01]\n",
" [1.85520303e+00 1.25173468e+01 2.35000000e-01]\n",
" [2.53548772e+00 1.46940407e+01 2.32500000e-01]\n",
" [3.32526518e+00 1.68970573e+01 2.30000000e-01]\n",
" [4.22585554e+00 1.91265572e+01 2.27500000e-01]\n",
" [5.23858677e+00 2.13826920e+01 2.25000000e-01]\n",
" [6.36479416e+00 2.36656038e+01 2.22500000e-01]\n",
" [7.60581987e+00 2.59754245e+01 2.20000000e-01]\n",
" [8.96301235e+00 2.83122750e+01 2.17500000e-01]\n",
" [1.04377258e+01 3.06762644e+01 2.15000000e-01]\n",
" [1.20313197e+01 3.30674894e+01 2.12500000e-01]\n",
" [1.37451577e+01 3.54860331e+01 2.10000000e-01]\n",
" [1.55806077e+01 3.79319647e+01 2.07500000e-01]\n",
" [1.75390403e+01 4.04053381e+01 2.05000000e-01]\n",
" [1.96218285e+01 4.29061910e+01 2.02500000e-01]\n",
" [2.18303469e+01 4.54345444e+01 2.00000000e-01]\n",
" [2.41659705e+01 4.79904009e+01 1.97500000e-01]\n",
" [2.66300741e+01 5.05737444e+01 1.95000000e-01]\n",
" [2.92240312e+01 5.31845382e+01 1.92500000e-01]\n",
" [3.19492128e+01 5.58227245e+01 1.90000000e-01]\n",
" [3.48069864e+01 5.84882229e+01 1.87500000e-01]\n",
" [3.77987152e+01 6.11809293e+01 1.85000000e-01]\n",
" [4.09257563e+01 6.39007145e+01 1.82500000e-01]\n",
" [4.41894598e+01 6.66474233e+01 1.80000000e-01]\n",
" [4.75911672e+01 6.94208725e+01 1.77500000e-01]\n",
" [5.11322102e+01 7.22208499e+01 1.75000000e-01]\n",
" [5.48139093e+01 7.50471129e+01 1.72500000e-01]\n",
" [5.86375718e+01 7.78993870e+01 1.70000000e-01]\n",
" [6.26044905e+01 8.07773641e+01 1.67500000e-01]\n",
" [6.67159422e+01 8.36807009e+01 1.65000000e-01]\n",
" [7.09731851e+01 8.66090177e+01 1.62500000e-01]\n",
" [7.53774580e+01 8.95618962e+01 1.60000000e-01]\n",
" [7.99299773e+01 9.25388784e+01 1.57500000e-01]\n",
" [8.46319359e+01 9.55394642e+01 1.55000000e-01]\n",
" [8.94845002e+01 9.85631104e+01 1.52500000e-01]\n",
" [9.44888087e+01 1.01609228e+02 1.50000000e-01]\n",
" [9.96459689e+01 1.04677182e+02 1.47500000e-01]\n",
" [1.04957056e+02 1.07766286e+02 1.45000000e-01]\n",
" [1.10423108e+02 1.10875806e+02 1.42500000e-01]\n",
" [1.16045127e+02 1.14004952e+02 1.40000000e-01]\n",
" [1.21824073e+02 1.17152882e+02 1.37500000e-01]\n",
" [1.27760862e+02 1.20318696e+02 1.35000000e-01]\n",
" [1.33856365e+02 1.23501435e+02 1.32500000e-01]\n",
" [1.40111403e+02 1.26700082e+02 1.30000000e-01]\n",
" [1.46526744e+02 1.29913553e+02 1.27500000e-01]\n",
" [1.53103100e+02 1.33140704e+02 1.25000000e-01]\n",
" [1.59841126e+02 1.36380322e+02 1.22500000e-01]\n",
" [1.66741412e+02 1.39631126e+02 1.20000000e-01]\n",
" [1.73804485e+02 1.42891767e+02 1.17500000e-01]\n",
" [1.81030800e+02 1.46160822e+02 1.15000000e-01]\n",
" [1.88420742e+02 1.49436793e+02 1.12500000e-01]\n",
" [1.95974615e+02 1.52718110e+02 1.10000000e-01]\n",
" [2.03692646e+02 1.56003123e+02 1.07500000e-01]\n",
" [2.11574977e+02 1.59290104e+02 1.05000000e-01]\n",
" [2.19621661e+02 1.62577246e+02 1.02500000e-01]\n",
" [2.27832658e+02 1.65862656e+02 1.00000000e-01]\n",
" [2.36207832e+02 1.69144363e+02 9.75000000e-02]\n",
" [2.44746948e+02 1.72420308e+02 9.50000000e-02]\n",
" [2.53449662e+02 1.75688349e+02 9.25000000e-02]\n",
" [2.62315527e+02 1.78946253e+02 9.00000000e-02]\n",
" [2.71343976e+02 1.82191705e+02 8.75000000e-02]\n",
" [2.80534327e+02 1.85422300e+02 8.50000000e-02]\n",
" [2.89885773e+02 1.88635545e+02 8.25000000e-02]\n",
" [2.99397383e+02 1.91828858e+02 8.00000000e-02]\n",
" [3.09068094e+02 1.94999567e+02 7.75000000e-02]\n",
" [3.18896706e+02 1.98144914e+02 7.50000000e-02]\n",
" [3.28881880e+02 2.01262050e+02 7.25000000e-02]\n",
" [3.39022133e+02 2.04348038e+02 7.00000000e-02]\n",
" [3.49315830e+02 2.07399856e+02 6.75000000e-02]\n",
" [3.59761187e+02 2.10414392e+02 6.50000000e-02]\n",
" [3.70356259e+02 2.13388453e+02 6.25000000e-02]\n",
" [3.81098940e+02 2.16318756e+02 6.00000000e-02]\n",
" [3.91986958e+02 2.19201941e+02 5.75000000e-02]\n",
" [4.03017871e+02 2.22034563e+02 5.50000000e-02]\n",
" [4.14189064e+02 2.24813097e+02 5.25000000e-02]\n",
" [4.25497741e+02 2.27533944e+02 5.00000000e-02]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00]]\n",
"80\n"
]
},
{
"data": {
"text/plain": [
"Text(0, 0.5, 'height (m)')"
]
},
"execution_count": 235,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de3xU9ZnH8c9DANmgoNwUEI1QLoqQQBAVahHRWrTuanW1gsVda0FqZbu6u121S72ttWrdlcqi1CreoGhF23opLYqCIGgQUFRE5SIXCaBBQYFAePaPMwlDkplMkrmefN+v17wmc67PyYF58nt+55yfuTsiIiLRmmU6ABERyT5KDiIiUoOSg4iI1KDkICIiNSg5iIhIDc0zHUAydOjQwQsKCjIdhohITlmyZMk2d+9Y27xQJIeCggJKSkoyHYaISE4xs3Wx5qmsJCIiNSg5iIhIDUoOIiJSg5KDiIjUEIoOaZFk2L9/P9u2bWP79u1UVFRkOhyRpMjLy+Pwww+nQ4cONGuWeHtAyUEa7dNP4fvfh5kz4aijMh1Nw23YsAEzo6CggBYtWmBmmQ5JpFHcnb1791JaWsqGDRs45phjEl5XZSVptFtvhddeg1tuyXQkjfPVV1/RtWtXWrZsqcQgoWBmtGzZkq5du/LVV1/Va10lB2mwv/s7MIMpU2D//uDdLJieq+rT7BbJFQ35d63/CdJgq1fDqFGQ32w3APn5MHo0rFmT4cBEpNGUHKTBOneGNm1g9/4WtGqxj927g8+53O+QbcyMnTt31mudTZs2MXz48EZv/6abbqK8vDzmulOmTKFPnz4MGDCAHTt21CvGumzfvp0777zzoGlXXnkl8+fPT+p+6uucc87h448/BmDatGmsWrWqat60adO46KKLGrTdgoICVqxYkZQYk0XJQRqldN0uruJ+Fv3iL1x1FWzenOmIpEuXLsydO7fR27n55pvjJodJkybx2GOPsXTpUg477LBG7y9abcnhwQcf5LTTTkvqfurrhRdeoEePHkDN5JAK+/btS+n241FykEaZdfXLTOYnFHb7nMmTYdasTEcUPpMmTeKkk06ie/fuPP3001XTFy9ezPDhwykuLqa4uJjnn38egLVr19KhQ4eq5Z5++umqv/Bvv/32Gq2F2rZ/9dVXAzBkyBCKiorYvn37QTFdcsklfPzxx/zgBz9g9OjRvPLKKwwaNKhqfvTnV155haKiIsaNG0f//v0pLCzk/fffr1r2oYceorCwkMLCQk466SRKS0u5+uqr2b59O0VFRQwZMgSA008/neeeew6A0tJSLrjgAvr370+/fv149NFHq7ZXUFDAxIkTOfXUUykoKOC+++6r83c8e/Zszj33XAC2bNmCmfHUU08BcOedd3LDDTdUbXvFihU8/PDDlJSUMGHCBIqKipgzZw4AX375JZdccgl9+/Zl6NChbI7x19L8+fPp168fgwcPZsKECUQP11xQUMCtt97K8OHDGTduHJs3b646z3379uU//uM/qpb94osvuPDCC+nTpw8jRoxgzJgx/Nu//Vudx5sQd8/5V3FxsUuG3HKLO7hPm5bpSBrtvffeqzkRUvdKAOC/+c1v3N39tdde8y5duri7e1lZmRcVFfmmTZvc3X3Tpk3etWtXLysr8zVr1nj79u3d3b20tNTbtWvnq1atcnf3e+65xwHfsWNH3O1XzqtcrjbHHnusv/POO+7uPnfuXI/+fxj9ee7cud68eXN/66233N39tttu81GjRlXN69Gjh3/66afu7r5jxw7ftWvXQcdQadiwYf7nP//Z3d0vvvhi//nPf1517EcddVRVLMcee6xfd9117u6+Zs0ab926ddzjcHf/6quvvF27dl5eXu7Tp0/3U0891ceNG+fu7t/+9rd9zpw5NY45Oh5394cfftgPP/xw/+STT9zd/corr/Qbbrihxr52797tXbp08blz57q7+8yZMx04KP7x48dXLb9r166q+MvLy3348OH+4osvurv7tdde6z/84Q/d3f2zzz7zgoKCqmOvrrZ/30CJx/heVctBGuett6B16+BypTBKZXpI0Pe//30ATjnlFDZt2sTu3btZuHAha9asYeTIkRQVFTFy5EjMjI8++uigdRctWsTAgQPp2bMnAFdccUVC20+23r17M2DAgKr9VNbtn3/+ecaMGcNRkY6qQw89lFatWtW5vTlz5jBu3DgAOnfuzLnnnntQKa3ymAoKCjjiiCPYsGFD3O3l5+fTt29fFi9ezJw5c5g4cSILFiygvLyckpIShg4dmtBxDh06lG7dutU4zmgffPAB+fn5nH766QBcfPHFtG3b9qBlxowZU/VzRUUF//7v/05hYSHFxcWsWLGCZcuWATB37lz++Z//GYB27dpx/vnnJxRnInQTnDTO0qUwcGC9vuykfiq/LPPy8oCgDu3u9O/fn3nz5tVYfu3atVU/u3ud92zUtv36at68Ofuj/kConmCiv/Dz8vKq9uGN+HdT/biiP8faXzwjRozgpZdeYtGiRUyZMoUjjzyS6dOnU1hYmFDCSnS/iRzzoYceWvXzPffcQ1lZGYsXL6ZVq1aMHTu26vebyPltKLUcpOE++wzKyuAb31BySLMhQ4bw4YcfHvTX8ptvvlnji+eUU05hyZIlVS2KadOmJbyPww47jC+++CKhZY877jhWr15NWVkZ7s6MGTMSWu+8887j0UcfpbS0FICdO3eyZ88e2rRpw9dffx3zS/3MM89k6tSpAGzevJkXXnghoSu0rr/++ph9ECNGjODhhx+mW7dutGzZkhEjRnDTTTcxYsSIWpdv06ZNwr+faH369GHXrl1Vif0Pf/hD3O1s376dzp0706pVKzZu3Mgf//jHqnnDhw/nkUceAaCsrOygeY2l5CANt3QpDBgAzZqFt6yUpY444gj+9Kc/cfPNN1NYWMjxxx/PTTfdVCM5HHnkkdx///2ce+65DB06lF27dtGiRQvy8/Pr3Md1113HGWecUWuHdHVdu3bluuuuo7i4mDPPPJPOnTsndBzDhg3j+uuv58wzz6SwsJAzzjiD7du3065dO0aPHk2/fv2qOqSjTZo0ieXLl9O/f3/OOuss7rjjDvr27Vvn/t5+++2qElZ1J598Mtu2batKBiNGjGDdunWcccYZtS4/duxYbr31VgYMGFDVIZ2IQw45hBkzZnD11VczePBgSkpK4j7WYsKECSxYsIABAwYwfvz4g5LVxIkT2bJlC3379uWyyy5j6NChNUpUDRarMyKZL+AXgAMnRj73Al4HVkXee0YtG3NerJc6pDPkV79y/+lP3X/0I/cHHsh0NI1Wa4d0CHz55ZdVPz/00EM+dOjQDEaTORUVFT548GCvqKjIdChJU15e7rt27XJ39y+++ML79evnf/vb32pdtr4d0invczCzgcApwCdRk+8HJrv742Z2GfAAcEYC8ySbvPUWnHMOLFigslIWmzRpEk899RT79u2jXbt2/Pa3v810SBnRrFkzFi9enOkwkqqsrIyRI0dSUVHB7t27GTVqFGeeeWZStm2ewv/UZnYI8AowCpgLfBfYQtAqaO/uFWaWB3wG9AQs1jx33xprP4MGDXKNIZ0B3/gG/PnPcO+9UFgI48dnOqJGef/99zn++OMzHYZIStT279vMlrj7oNqWT3Wfwy3A4+4e/bSdbsBGd68AiLxvikyPN+8gZjbWzErMrGTr1ph5Q1KlrAxKS6FXr6DPQS0HkVBJWXIws1OBk4D/S8X23X2quw9y90EdO3ZMxS4kniVLgs7ovLzgUaxKDiKhksqWwzCgD7DGzNYCRwOzgR5A10jJiMh7F2B95BVrnmSTJUug8nEJZrpaSSRkUpYc3P0Od+/i7gXuXgBsAM529yeBZcClkUUvBZa6+1Z33xJrXqrilAYqKTmQHFRWEgmdTN3ncBVwjZmtAq6JfE5knmSL6OSgspJI6KQtOURaECsiP69095PdvVfk/YOo5WLOkyyxbRt8/nlwtRKorEQwjvawYal/ZHld4zvU9qjrZFu7dm3V3cmVosc5kHDQHdJSf0uWQHFxUE4ClZXInnG0M5Ucosc5kHBQcpD6iy4pQZMuK6V6HO1Zs2bRp08fhgwZwm233VY1PdZYDrWNg/DRRx8xYsQI+vfvz8CBA/nLX/5StR0z4/bbb691vIjRo0czaNAg+vXrxwUXXEBZWVnVPt577z2KioqqRj6LHsmsofuTLBPr1ulceunxGWl2/vnuM2ce+Hzdde533pm5eJKkIY/P2LTJfdQo9/z84Dnc+fnuo0e7R4YnaJTKsRhWrlzp7u6/+tWvHPD169cnNJZDpcGDB/uDDz7o7u7vvvuut2/f3rds2eLuHnc8h61bt1b9fOONN/rPfvYzd685doP7weMcNHR/kloaz0FSr7KsVKkJl5WqxtHeDa1akdRxtCvHYujduzcQPOgN4K233kpoLAeAHTt2sGzZsqpn/p9wwgkUFRWxaNGiqmVijefw6KOPUlxcTL9+/Zg+fXrVGALxNGZ/kl00noPUT2kp7NgB3bsfmNaEy0oQ/EquugrGjoWpU4PO6WTwGL9TT3Ash3jbqG3sg+jxHObPn8+UKVNYuHAhHTt2ZPr06TX6GeoTc137k+yjloPUT+XNb9EDjDTxq5VmzYLJk4PHSyVzHO1TTz2VpUuX8uGHHwLw4IMPAjBw4MCYYzlUHwehTZs2FBUVVT3zf+XKlSxfvpyTTz457r63b99O27Ztad++PXv27OGhhx6qmhdvHIOG7k+yj5KD1E/1zmho8i2HVOnUqRNTp07lvPPOY8iQITRvHjT0443lUNs4CE888QSPP/44/fv3Z9SoUTz22GPU9ciZkSNH0qNHD/r06cPIkSMZOHBg1bz+/fvTu3dvTjzxxKoO6WgN2Z9kn5Q+lTVd9FTWNPr7v4fLL4cLLzww7cYbIT8/eM9heiqrhFm2PZVVwkYtB5EmQclBErdpE+zdC9WHNGzifQ4iYaTkIImrbDVEd0ZDqC5l3a8kJyHUkH/XSg6SuOr3N1QKSVmpdevWbNy4kfLy8piXZIrkEnenvLycjRs30rp163qtq/scJHElJfCjH9WcbgYVFemPJ8mOPvpotm3bxrp163TtvYRG8+bNadu2LR06dKjfeimKR8LGHRYvhtoGp2/WDELwZdqsWTM6depEp06dMh2KSMaprCSJWb06eJpcly4154WkrCQiByg5SGIWL4ZYd7nqaiWR0FFykMTESw4hulpJRAJKDpKYxYth8ODa56msJBI6Sg5St/JyeOed2i9jBZWVREJIyUHqtnx5MF70oYfWPl9lJZHQUXKQusXrbwCVlURCSMlB6pZIclBZSSRUlBykbnUlB5WVREJHyUHi+/xz2LwZ4o1zoLKSSOgoOUh8b7wRPIk1Mt5vrVRWEgkdJQeJr66SEqisJBJCSg4SXyLJQWUlkdBRcpDY3IOyUiLJQWUlkVBRcpDYVq+G/Hzo3Dn+cmo5iISOkoPElkhJCdTnIBJCSg4SW6LJQS0HkdBRcpDY4j2JNZr6HERCR8lBardrF6xYEdzjUBeVlURCR8lBardkCZxwQtAhXReVlURCR8lBardgAQwdmtiyKiuJhI6Sg9Ru4UIYMiSxZVVWEgkdJQepyT1IDvVpOSg5iIRKSpODmT1rZsvNbKmZzTezosj0Xmb2upmtirz3jFon5jxJkw8/hNatoUuXxJZXWUkkdFLdcrjc3QvdfQBwN/BQZPr9wGR37wVMBh6IWifePEmHBQsSLymBykoiIZTS5ODuX0R9bAvsN7NOwEBgRmT6DGCgmXWMNy+VcUo19SkpgcpKIiHUPNU7MLMHgW8DBnwH6AZsdPcKAHevMLNNkekWZ97WVMcqEQsWwI9/nPjyKiuJhE7KO6Td/Up3Pwa4AbgrWds1s7FmVmJmJVu3Km8kzeefw4YN0K9f4uuorCQSOmm7WsndHwOGAxuArmaWBxB57wKsj7xizau+vanuPsjdB3XsqKpT0rz+evDIjOb1aFSqrCQSOilLDmZ2qJl1i/p8HvA5sAVYBlwamXUpsNTdt7p7zHmpilOqqW9/A6isJBJCqexzaA08ZWatgQqCxHCeu7uZXQU8YmYTgTJgTNR68eZJqi1YANdfX791VFYSCZ2UJQd3LwVOiTFvJVDrs6DjzZMU27s3eKbSKbWetthUVhIJHd0hLQcsWwbHHQdt29ZvPZWVREJHyUEOaEh/A6isJBJCSg5yQH3vjK6kspJI6Cg5SMC9cclBZSWRUFFykMDq1cF79+71X1dlJZHQUXKQwKuvwrBhQSugvlRWEgkdJQcJVCaHhlByEAkdJQcJzJsH3/pWw9ZVn4NI6CR0E5yZ9QKOBnYBK9x9R0qjkvT65BP4+mvo06dh66vPQSR0YiYHMzsMuA74IbAHKAVaAd3NbBFwl7u/nJYoJbVefTVoNTSkvwFUVhIJoXgth5eBx4DiyAPxADCzZsA3gXFm9g13n5riGCXVGlNSApWVREIoXnIY6u7l1Se6+35gHjDPzFqmLDJJn1dfhWuuafj6KiuJhE7M5FA9MUQSQfOo+V/Xljwkx3z6KXz2GZx4YsO3obKSSOjUebWSmV1oZusJOqN3ADsj7xIG8+bBN78Z/PXfUCoriYROIt8IdwHfA1q4e567N3P3vBTHJenSmPsbKqmsJBI6iSSHT939zUhfg4RNMpKDykoioZPIfQ6/MbNbgWeA3ZUT3f29lEUl6bF1K2zcCEVFjduOykoioZNIcugKXAtcTjDcJ4ADDXhCm2SV+fODp7DmNbJKqLKSSOgkkhwmAN9w909THYykWTJKSqCykkgIJdLnsE6JIaSSmRxUVhIJlURaDm+Y2QzgKQ7uc3ghZVFJ6n3+eTCGQ3Fx47elspJI6CSSHCq/PaJvoXVAySGXzZ0b3N/QokXjt6Wykkjo1Jkc3H14OgKRNHvpJRgxIjnbUllJJHRi9jmYWZe6Vjazo5IbjqRNMpODykoioROvQ3qmmU0xs9PM7JDKiWZ2jJmNNbOFwNDUhyhJt2FD0OfQv39ytqeykkjoxCsrfQv4R+AXwBAz20MwnkMpwQ1xo9x9bcojlOR76SUYPrxxz1OKprKSSOjEeyqrA08CT5pZc6ADsMvdv0hXcJIiySwpgcpKIiGU0J+O7r7P3TcrMYSAO8yZk9zkoLKSSOgkqa4gOWPlyuDy1R49krdNJQeR0FFyaGoqS0oNHS+6NupzEAmdRAb7aZPINMkRye5vAPU5iIRQIi2HVxKcJtmuoiJ4ntLwJN/XqLKSSOjEvFopcoVSS6CZmf0dUFmHaAvkpyE2SbY334QuXYJXMqmsJBI68VoONxKMF90P+Cry807gfeCJ1IcmSTd7NnznO8nfrspKIqETMzm4+83u3gyYEhk3uvJ1uLvfmsYYJVlmz4azz07+dlVWEgmdRB689xMAM2sZvby7f53CuCTZysrgnXfgtNOSv22VlURCJ5GrlS4wsw0EYznsICgt7Uh1YJJkL70UPKK7Vavkb1tlJZHQSeRqpbuAi4Hm7p4XKS01ctBhSbtU9TeAykoiIZRIcvjc3Re6e73qBmbW3sxeMLMPzOxtM5tlZh0j83qZ2etmtiry3jNqvZjzpIHc4S9/SU1/A6isJBJC8cZzyDezfOAZMxtvZu0qp0Wm18WBO929t7v3Bz4G7ojMux+Y7O69gMnAA1HrxZsnDfH++0Hpp3fv1GxfZSWR0InXIb2T4Au+8v6GyVGfHYhbWnL3zzn4ZrlFwHgz6wQMBM6KTJ8B3BdpVVisee6+NcFjkuoqr1JK5iMzoqmsJBI68S5lbRbVx9Cs2ud69TmYWTNgPPAnoBuw0d0rIvupADZFpsebV32bY82sxMxKtm5V3ogrVZewVlJZSSR0ErlaKb/6qwH7+Q1BS+S+BqxbK3ef6u6D3H1Qx44dk7XZ8Nm1CxYsSP7zlKKprCQSOol0SFdeulr1MrPdZjbPzOosYpvZ3UBP4JJIp/Z6oKuZ5UXm5wFdItPjzZOGmDsXBgyAww9P3T5UVhIJnTpvgiN4jMYu4CGCPoF/4sBwoQ8Ap8da0cz+GygGznX3PQDuvsXMlgGXAo9H3pdW9inEmycN8Nxz8N3vpnYfKiuJhE4iyeEidy+O+nyvmc1399PM7LpYK5lZX+AGYBWw0ILO0DXufgFwFfCImU0EyoAxUavGmyf14Q7PPw8vvpja/aisJBI6iSSHfDPr7u6rAcysO8F40gD7Yq3k7u9y4Eqn6vNWAifXd57U04oVkJcHxx+f2v2orCQSOokkh58Db5jZEoJLWIuBq8zsUOCpVAYnjVRZUkrVJayVVFYSCZ1EHrz3tJnNJ/hr3oBF7r4lMvv2VAYnjfTcc/CLX6R+PyoriYROIi0HIsngzymORZJp27agrDRsWOr3pbKSSOjEGwnuJXcfYWZbCcpJVbMAd/dOKY9OGu7FF+GMM+CQQ1K/L5WVREInXsvhssj7oHQEIkmWjktYK6nlIBI68R6f8WnkfR3BJaXt3X1d5StdAUoD7N0Lf/0rnHNOevanPgeR0Enk8RnnAO8CsyKfB5mZ+h+y2bx50LMndO6cnv2p5SASOok8PuNm4CSC1gPuXgL0SGVQ0kizZsEFF6Rvf+pzEAmdRK9W2mwHXyu/JzXhSKPt3w/PPhsMC5ouKiuJhE4iLYcdZnYkkSuWzOx0YHsqg5JGePNNaNsW+vRJ3z5VVhIJnURaDv8JvAgcZ2avEDxh9e9TGZQ0QrpLSqCykkgIJXKH9BtmNhwYQnCPw0J3V8shG7nDM8/AjBnp3a/KSiKhk2ifwxdm9lLl8maW7+5fpzQyqb/33oM9e2DgwPTuV2UlkdBJ5FLW75nZBoIxHXZwYPAfyTbPPAPnn5/6B+1Vp7KSSOgk0iF9J3Ax0KKhY0hLmsyaBd/7Xvr3q7KSSOgkUlb63N0XpjwSaZy1a2H9ehg6NP37rmypuKe/1SIiKRGz5WBm+WaWDzxjZuPNrF3ltMh0ySZPPhm0Gpon1I2UGmo9iIRGvG+SnQT3NlT+KTg56rMDKi1lk5kz4a67Mrd/lZZEQiVmcnD3RPojJBt8+CFs3JiesRti0RVLIqGiBBAGM2fCP/5jMF50puiKJZFQUXIIg5kz4ZJLMhuDykoioaLkkOvefRfKymDIkMzGobKSSKgoOeS6ylZDswyfSpWVREJFySGXucPvf5/5khKorCQSMkoOuWzpUti3D046KdORqKwkEjJKDrnsscdg9OjsuCtZZSWRUMng7bTSKHv3wvTp8NprmY4koJaDSKio5ZCrZs+GHj2gZ89MRxJQn4NIqCg55KpHH4XLL890FAeo5SASKkoOuaisDP76V7j44kxHcoD6HERCRckhFz35JHz723DEEZmO5ACVlURCRckhFz3yCIwZk+koDqaykkioKDnkmg8/hI8/hrPPznQkB1NZSSRUlBxyze9+B5ddBi1aZDqSg6msJBIqus8hl5SXw7Rp8OqrmY6kJpWVREJFLYdc8uyzcPzx0Lt3piOpSWUlkVBRcsglU6fCuHGZjqJ2KiuJhErKkoOZ3W1ma8zMzezEqOm9zOx1M1sVee+ZyLwm76OP4O234YILMh1J7VRWEgmVVLYcngW+BayrNv1+YLK79wImAw8kOK9pmzo1uCP6kEMyHUntVFYSCZWUdUi7+2sAFvXEUDPrBAwEzopMmgHcZ2YdAYs1z923pirOnLBnT3BvQ7Y8ZK82KiuJhEq6+xy6ARvdvQIg8r4pMj3evBrMbKyZlZhZydatIc8dzzwDJ56YPQ/Zq43KSiKhkrMd0u4+1d0Hufugjh07Zjqc1Lr3XrjmmkxHEZ/KSiKhku77HNYDXc0sz90rzCwP6BKZbnHmNV2LFkFpKZx3XqYjiU9lJZFQSWvLwd23AMuASyOTLgWWuvvWePPSGWPW+d//hQkTIC8v05HEp7KSSKik8lLWSWa2ATgamGNm70ZmXQVcY2argGsin0lgXtOzfn3waO4rrsh0JHVTWUkkVFJ5tdIEYEIt01cCJ8dYJ+a8Jmny5ODpq23aZDqSuqmsJBIqerZStvrqK3jwQVi8ONORJEZlJZFQydmrlULvkUdg6NBgnOhcoLKSSKio5ZCN9u6FO++EGTMyHUniVFYSCRW1HLLRE09A9+5w6qmZjiRxKiuJhIpaDtmmogJ++UuYMiXTkdSPkoNIqKjlkG2efhratYPhwzMdSf2oz0EkVJQcsok73H473Hhj8GWbS9TnIBIqSg7Z5Pnng/dzz81sHA2hspJIqCg5ZIv9++HnP4eJE3Ov1QAqK4mEjJJDtvj976FVq+wd6a0uKiuJhIquVsoG5eXwX/8Fv/tdbrYaQGUlkZBRyyEb/Pa3wUA+p5+e6UgaTmUlkVBRyyHTdu6E226DF17IdCSNo7KSSKio5ZBpv/41DBsGAwZkOpLGUVlJJFTUcsiktWth0iRYsiTTkTSeykoioaKWQyZdey389KdQUJDpSBpPZSWRUFHLIVNmz4bly2H69ExHkhwqK4mEiloOmVBeHowLfe+9wb0NYaCykkioKDlkwq9/Db16wXe/m+lIkkdlJZFQUVkp3d57D+65B958M9ORJJfKSiKhopZDOu3bB5dfHtzXEIZO6GgqK4mEipJDOt11Fxx+OIwdm+lIkk9lJZFQUVkpXd55JygnLVmSu89PikdlJZFQUcshHXbtgssugzvugGOOyXQ0qaGykkioKDmkw09/CiecAFdckelIUkdlJZFQUVkp1aZPh7lzoaQknOWkSioriYSKkkMqvfUW/Mu/wJw50KZNpqNJLSUHkVBRWSlVNm+G88+H+++HwsJMR5N66nMQCRUlh1TYuRP+4R/ghz+ECy/MdDTpoT4HkVBRcki28nK46CI48USYODHT0aSPykoioaLkkEz79sGYMXDIIfDAA+HugK5OZSWRUFGHdLLs3QujRsGOHfDMM9C8if1qVVYSCZUm9g2WIl9/DZdeGvzl/Oyz4XkMd32orCQSKiorNVZpKQwfDm3bwtNPN83EACoriYSMkkNjLF4MJ58M3/kOPPIItGyZ6YgyR2UlkVBRWakh9u8PRnH75S9h6tTgfq1UCcMAAAb0SURBVIamTmUlkVBRcqivlSvhyiuDnxctgu7dMxtPtlBZSSRUVFZKVGkp/OQn8M1vBlclzZunxBBNZSWRUMnK5GBmvczsdTNbFXnvmbFgVq6EH/8Yjj8+uDy18nOzrPzVZY7KSiKhkq1lpfuBye7+uJldBjwAnJGWPZeXBw/MmzsX/vAH2LgxGLnt3Xehc+e0hJCTVFYSCZWsSw5m1gkYCJwVmTQDuM/MOrr71qTu7Msv4Wc/C56FtH07rF4Na9ZAr17wrW/B3XcH73l5Sd1tKDVrBtOmBf0wIpI+V18djBeTZFmXHIBuwEZ3rwBw9woz2xSZXpUczGwsMBbgmIaOrtayJfTvD4ceGjxS+7jjoEcPaN26scfQ9Pzrv8LChZmOQqTpSdH3VTYmh4S4+1RgKsCgQYMaVuxu1QrGj09mWE3XgAHBS0RCIRt7VdcDXc0sDyDy3iUyXURE0iDrkoO7bwGWAZdGJl0KLE16f4OIiMSUrWWlq4BHzGwiUAaMyXA8IiJNSlYmB3dfCZyc6ThERJqqrCsriYhI5ik5iIhIDUoOIiJSg5KDiIjUYB6Ch6WZ2VZgXQNX7wBsS2I4maRjyU5hOZawHAfoWCod6+4da5sRiuTQGGZW4u6DMh1HMuhYslNYjiUsxwE6lkSorCQiIjUoOYiISA1KDpGH94WEjiU7heVYwnIcoGOpU5PvcxARkZrUchARkRqUHEREpIYmnRzMrJeZvW5mqyLvPTMdU6LMbK2ZrTSzZZHX2ZHpWX9MZna3ma0xMzezE6Omx4w9G48rznHUem4i87LuOADMrL2ZvWBmH5jZ22Y2y8w6RublzHmp4zhy8bw8a2bLzWypmc03s6LI9NSfE3dvsi/gZeCyyM+XAS9nOqZ6xL4WODEXjwn4JsGwrwcdQ7zYs/G44hxHrecmW48jEks74PSoz3cBv8u181LHceTieWkb9fM/AG+l65xk/OAz+EvvBGwH8iKf8yKfO2Y6tgTjr/EPPdeOKfoY4sWe7ceVaHLI9uOoFuuFwJxcPi/RxxGG80Iwrk1Jus5JUy4rdQM2unsFQOR9U2R6rngi0nT+PzM7nNw+pnix5+JxVT83kCPHYWbNgPHAn8jh81LtOCrl3HkxswfN7BPgv4HLSdM5acrJIded5u6FwEmAAfdlOB45INfPzW+AneRe3NVVP46cPC/ufqW7HwPcQFAmS4umnBzWA13NLA8g8t4lMj3rufv6yPse4P+AoeT2McWLPaeOK8a5gRw4DjO7G+gJXOLu+8nR81LLceT0eQFw98eA4cAG0nBOmmxycPctwDLg0sikS4Gl7r41c1Elxsxam1nbyM8GfB9YlsvHFC/2XDquWOcGsv/fnJn9N1AMnB/5As3J81LbceTieTGzQ82sW9Tn84DPgfSck0x3smTyBfQBFgOrIu+9Mx1TgnF3B5YCbwPvAk8BnXPlmIBJBH/97AM2A+/WFXs2HldtxxHv3GTrcUTi6gs48EHky2UZ8EyunZdYx5GL5wU4ElgEvBM5jpeBgek6J3p8hoiI1NBky0oiIhKbkoOIiNSg5CAiIjUoOYiISA1KDiIiUoOSgzRpZnaTmbWM+nyLmV2Swv0VmdnFCSzzxzqWOcfMHkhudCIH6FJWadLMzIHD3H1nmvb3T8B33f2iOMu8CNzi7q/Xsa0lwMXu/nFyoxRRy0GaMDObHPlxYeT5/oeb2TQz+0lk/k1m9vvI+AAfmdlMMxtgZi+b2cdmdlfUtjqb2R/M7A0ze8fMbqhlf+2BW4AzI/ubVMsyxxDctPR65HMnM5sT2eY7ZvY/UYs/CVyRvN+IyAFKDtJkufvVkR+HuHuRu2+vZbFigkcQ9Ca48/QOYCTQH7g8aiCVR4FJ7j44ss5IMzur2v4+AyYSPEK6yN0n1LK/YcAbUZ9HA+vcvZ+79yNILpVeB0YkfsQiiWue6QBEstxsd/8CwMzeBpZ78LyePWb2AdDDzDYBpwMdg8f2AHAYcDzwt3ru72igNOrzIuDaSCvlVWB21LzNkeVFkk7JQSS+3VE/V9TyuTlBC9yBk9x9byP3twtoVfnB3V+PDA15FvAD4D8JRqAjstyuRu5PpFYqK0lTtwNo25gNuPsOYD7BFzcAZtbNzI6qZfEv69jfOwQlrMrtHAd86e6/B64FiiOD2EDQMlnemNhFYlFykKbu18DLlR3SjdjOaOCEyo5jYCZQ2/ZeAlpHBo2v0SENvAYcV/l4aYJy1VIzWwa8CFzlkfEJgLOBpxsRs0hMupRVJMuY2fXAbnf/nzjLtCd4hPNJ7l6etuCkyVDLQST73EPdfQndgfFKDJIqajmIiEgNajmIiEgNSg4iIlKDkoOIiNSg5CAiIjUoOYiISA3/D29Mr0WvojA2AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def rocket1(state,dmdt=0.05, u=250,c=0.18e-3):\n",
" '''Computes the right-hand side of the differential equation\n",
" for the acceleration of a rocket, with drag, in SI units.\n",
" \n",
" Arguments\n",
" ---------- \n",
" state : array of three dependent variables [y v m]^T\n",
" dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n",
" u : speed of propellent expelled (default is 250 m/s)\n",
" c : drag constant for a rocket set to 0.18e-3 kg/m\n",
" Returns\n",
" -------\n",
" derivs: array of three derivatives [v (u/m*dmdt-g-c/mv^2) -dmdt]^T\n",
" '''\n",
"# print(dmdt)\n",
" dstate = np.zeros(np.shape(state))\n",
" \n",
" dstate[0] = state[1]\n",
" dstate[1] = u*dmdt/state[2] - 9.8 - c*state[1]**2/state[2]\n",
" dstate[2] = - dmdt\n",
" return dstate\n",
"\n",
"def f_m(dmdt,m0=0.25, c=0.18e-3, u=250):\n",
" ''' define a function f_m(dmdt) that returns \n",
" height_desired-height_predicted[-1]\n",
" here, the time span is based upon the value of dmdt\n",
" \n",
" arguments:\n",
" ---------\n",
" dmdt: the unknown mass change rate\n",
" m0: the known initial mass\n",
" c: the known drag in kg/m\n",
" u: the known speed of the propellent\n",
" \n",
" returns:\n",
" --------\n",
" error: the difference between height_desired and height_predicted[-1]\n",
" when f_m(dmdt)= 0, the correct mass change rate was chosen\n",
" '''\n",
" \n",
" mid_state = np.zeros(3)\n",
" for i in range(300):\n",
" mid_state = num_rk2[i] + rocket1(num_rk2[i],dmdt) * dt*0.5 \n",
" num_rk2[i+1] = num_rk2[i] + rocket1(mid_state,dmdt)*dt\n",
" if num_rk2[i+1,2] <= 0.05 :\n",
" break\n",
" max_height = num_rk2[i+1,0]\n",
"# print(dmdt)\n",
" error = 300 - max_height\n",
"# print(max_height)\n",
" return error\n",
"\n",
"def incsearch(func,xmin,xmax,ns=50):\n",
" '''incsearch: incremental search root locator\n",
" xb = incsearch(func,xmin,xmax,ns):\n",
" finds brackets of x that contain sign changes\n",
" of a function on an interval\n",
" arguments:\n",
" ---------\n",
" func = name of function\n",
" xmin, xmax = endpoints of interval\n",
" ns = number of subintervals (default = 50)\n",
" returns:\n",
" ---------\n",
" xb(k,1) is the lower bound of the kth sign change\n",
" xb(k,2) is the upper bound of the kth sign change\n",
" If no brackets found, xb = [].'''\n",
" x = np.linspace(xmin,xmax,ns)\n",
" f = func(x)\n",
" sign_f = np.sign(f)\n",
" delta_sign_f = sign_f[1:]-sign_f[0:-1]\n",
" i_zeros = np.nonzero(delta_sign_f!=0)\n",
" nb = len(i_zeros[0])\n",
" xb = np.block([[ x[i_zeros[0]+1]],[x[i_zeros[0]] ]] )\n",
"\n",
" if nb==0:\n",
" print('no brackets found\\n')\n",
" print('check interval or increase ns\\n')\n",
" else:\n",
" print('number of brackets: {}\\n'.format(nb))\n",
" return xb\n",
" \n",
"N = 1300\n",
"step = (0.4 - 0.05)/N\n",
"temp_array = np.zeros(N)\n",
"\n",
"for i in range(N):\n",
" temp_array[i] = f_m(0.05+i*step, 0.25, 0.18e-3, 250)\n",
"# print(f_m(0.05+i*step, 0.25, 0.18e-3, 250))\n",
"# if abs(f_m(0.05+i*step, 0.25, 0.18e-3, 250)) < 1 :\n",
"# print(f_m(0.05+i*step, 0.25, 0.18e-3, 250))\n",
" \n",
"temp_array = np.sort(abs(temp_array))\n",
"print(temp_array)\n",
"\n",
"def mod_secant(func,dx,x0,es=0.0001,maxit=50):\n",
" '''mod_secant: Modified secant root location zeroes\n",
" root,[fx,ea,iter]=mod_secant(func,dfunc,xr,es,maxit,p1,p2,...):\n",
" uses modified secant method to find the root of func\n",
" arguments:\n",
" ----------\n",
" func = name of function\n",
" dx = perturbation fraction\n",
" xr = initial guess\n",
" es = desired relative error (default = 0.0001 )\n",
" maxit = maximum allowable iterations (default = 50)\n",
" p1,p2,... = additional parameters used by function\n",
" returns:\n",
" --------\n",
" root = real root\n",
" fx = func evaluated at root\n",
" ea = approximate relative error ( )\n",
" iter = number of iterations'''\n",
"\n",
" iter = 0;\n",
" xr=x0\n",
" for iter in range(0,maxit):\n",
" xrold = xr;\n",
" dfunc=(func(xr+dx)-func(xr))/dx;\n",
" xr = xr - func(xr)/dfunc;\n",
" if xr != 0:\n",
" ea = abs((xr - xrold)/xr) * 100;\n",
" else:\n",
" ea = abs((xr - xrold)/1) * 100;\n",
" if ea <= es:\n",
" break\n",
" return xr,[func(xr),ea,iter]\n",
"\n",
"\n",
"root,out = mod_secant(f_m,0.001,1,es=0,maxit=100)\n",
"\n",
"print('The root of the function f_m is', root)\n",
"print(num_heun)\n",
"\n",
"t = np.linspace(0,len(num_heun),len(num_heun))\n",
"for i in range(len(num_heun)):\n",
" if num_heun[i,2] < 0.05 :\n",
" xd = i\n",
" yd = num_heun[i,0]\n",
" print(i)\n",
" break\n",
"\n",
"plt.rcParams.update({'font.size': 11})\n",
"plt.rcParams['lines.linewidth'] = 1\n",
"\n",
"plt.plot(t, num_heun[:,0], 'r-',label='height function, with drag')\n",
"plt.plot(xd, yd, 'b*',label='detonation') \n",
"\n",
"plt.legend();\n",
"plt.xlabel('time t (s)')\n",
"plt.ylabel('height (m)')\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n",
"\n",
"1. Math 24 _Rocket Motion_. <https://www.math24.net/rocket-motion/\\>\n",
"\n",
"2. Kasdin and Paley. _Engineering Dynamics_. [ch 6-Linear Momentum of a Multiparticle System pp234-235](https://www.jstor.org/stable/j.ctvcm4ggj.9) Princeton University Press \n",
"\n",
"3. <https://en.wikipedia.org/wiki/Specific_impulse>\n",
"\n",
"4. <https://www.apogeerockets.com/Rocket_Motors/Estes_Motors/13mm_Motors/Estes_13mm_1_4A3-3T>"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}