{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Download\n",
"\n",
"## Notebook link: http://tinyurl.com/eng122-oct19\n",
"\n",
"# Submission\n",
"\n",
"## Step 1: Rename your notebook\n",
"\n",
"LAST1_NAME1_LAST2_NAME2\n",
"\n",
"for example:\n",
"\n",
"MOORE_JASON_HOWLETT_JAMES\n",
"\n",
"## Step 2: Turn in your notebook and a PDF version to **Canvas** by midnight Saturday."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Introduction\n",
"\n",
"The center of mass of a book (modeled as a homogeneous cuboid of mass $m$, length $l$, and height $h$) lies directly above the top of a cylinder of radius $r$. There is sufficient friction to prevent slipping of the two surfaces when the book is perturbed from this equilibrium position. Let $\\theta$ be the angle between the vertical and the radial line that passes through the contact point when the book oscillates in a vertical plane perpendicular to the axis of the cylinder."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.display import SVG, Latex\n",
"SVG('book-balance.svg')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Non-linear Equation of Motion\n",
"\n",
"The book oscillates at an angular rate, $\\dot{\\theta}$, and the magnitude of the velocity of the center of mass of the book can be shown to be $v = \\frac{1}{2} \\sqrt{\\left(h^{2} + 4 r^{2} \\theta^{2}\\right) \\dot{\\theta}^{2}}$. The moment of inertia of the book is approximately $\\frac{m}{12} \\left(h^{2} + l^{2}\\right)$.\n",
"\n",
"Thus, the total kinetic energy is:\n",
"\n",
"$$T = \\frac{m \\dot{\\theta}^{2}}{24} \\left(h^{2} + l^{2}\\right) + \\frac{m \\dot{\\theta}^{2}}{8} \\left(h^{2} + 4 r^{2} \\theta^{2}\\right)$$\n",
"\n",
"The potential energy is:\n",
"\n",
"$$U = - g m \\left(\\frac{h}{2} - r \\theta \\operatorname{sin}\\left(\\theta\\right) + r - \\left(\\frac{h}{2} + r\\right) \\operatorname{cos}\\left(\\theta\\right)\\right)$$\n",
"\n",
"The Langragian can be formed and then finally the equation of motion:\n",
"\n",
"$$ - \\frac{g h}{2} m \\operatorname{sin}\\left(\\theta\\right) + g m r \\theta \\operatorname{cos}\\left(\\theta\\right) + \\frac{h^{2} m}{3} \\ddot{\\theta} + \\frac{l^{2} m}{12} \\ddot{\\theta} + m r^{2} \\theta^{2} \\ddot{\\theta} + m r^{2} \\theta \\dot{\\theta}^{2} = 0$$\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from IPython.display import Latex"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib notebook"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Question\n",
"\n",
"Linearize the non-linear equation of motion about $\\theta=0$ and type the result below:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Type answer here*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Question\n",
"\n",
"Derive an approximate expression for the natural frequency of small oscillations by linearizing the above equation of motion about $\\theta=0$, the equilibrium point and compute the natural frequency in rad/s and Hz along with the period of osciallation. Use the parameters $m= 1.058\\textrm{ kg}$, $l = 0.238 \\textrm{ m}$, $g = 9.81 \\textrm{ ms}^{-1}$, $h = 0.029 \\textrm{ m}$, and $r=0.042 \\textrm{ m}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Type answer here*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m = 1.058 # kg\n",
"l = 0.238 # m\n",
"g = 9.81 # m/s^2\n",
"h = 0.029 # m\n",
"r = 0.042 # m"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"omega_n = \n",
"Latex('$\\omega_n = {:1.3f} \\\\textrm{{ rad/s}}$'.format(omega_n))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fn = \n",
"Latex('$f_n = {:1.3f} \\\\textrm{{ Hz}}$'.format(fn))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"Tn = \n",
"Latex('$T_n = {:1.3f} \\\\textrm{{ s}}$'.format(Tn))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Question\n",
"\n",
"Describe how the mass affects the natural frequency of the system."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Type answer here*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Question\n",
"\n",
"Are there any limits to the size of the book (the ratios $h/r$ or $l/r$) by a requirement of stability of the oscillations? Said another way, how high a pile of books can you balance? Hint; it is probably difficult to balance a very high pile of book. Recall that if the effective stiffness is negative you will get unstable behavior in an $m$-$k$ model."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Type answer here*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Question\n",
"\n",
"Simulate the linear and non-linear equations of motion that predict the free response of the book released from rest at an initial angle $\\theta_0=1\\textrm{ deg}$ for 10 seconds. Use `scipy.integrate.odeint` for the numerical integration. Plot the results of each simulation on the same graph."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from scipy.integrate import odeint"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def eval_nonlin_rhs(x, t):\n",
" \"\"\"Evaluates the right hand side of the non-linear differential equations.\n",
" \n",
" Parameters\n",
" ==========\n",
" x : array_like, shape(2, )\n",
" An array of the states: [theta, omega].\n",
" t : float\n",
" The value of time.\n",
" \n",
" Returns\n",
" =======\n",
" xdot : array_like, shape(2, )\n",
" An array of the derivatives of the states: [thetadot, omegadot].\n",
" \n",
" \"\"\"\n",
" \n",
" theta, omega = x\n",
" \n",
" # type your first order equations of motion here\n",
"\n",
" return thetadot, omegadot"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def eval_lin_rhs(x, t):\n",
" \"\"\"Evaluates the right hand side of the non-linear differential equations.\n",
" \n",
" Parameters\n",
" ==========\n",
" x : array_like, shape(2, )\n",
" An array of the states: [theta, omega].\n",
" t : float\n",
" The value of time.\n",
" \n",
" Returns\n",
" =======\n",
" xdot : array_like, shape(2, )\n",
" An array of the derivatives of the states: [thetadot, omegadot].\n",
" \n",
" \"\"\"\n",
" \n",
" theta, omega = x\n",
" \n",
" # type your first order linear equations of motion here\n",
" \n",
" return thetadot, omegadot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Simulate the system."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"t = np.linspace(0, 10, num=1000)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x_nonlin = odeint(eval_nonlin_rhs, [np.deg2rad(1), 0], t)\n",
"x_lin = odeint(eval_lin_rhs, [np.deg2rad(1), 0], t)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot the trajectory of $\\theta$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, ax = plt.subplots(1, 1)\n",
"\n",
"# type both plot commands here as `ax.plot(...)`\n",
"\n",
"ax.set_xlabel('Time [s]')\n",
"ax.set_ylabel('$\\\\theta$ [deg]')\n",
"ax.legend(['Non-linear', 'Linear']);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Question\n",
"\n",
"Using your program, create a function that calculates the period of the non-linear model to three significant figures of the 11 oscillations when $\\theta_0= 1, 10 \\textrm{ and } 20 \\textrm{ deg}$. Compare these results to the period predicted by the linear model. By how much and why do they differ?\n",
"\n",
"*Hint: Look for sign changes with `np.sign()`, use boolean indexing to extract important times, and finally `np.diff()` and `np.mean()` can be useful for finding the delta times and averaging. Note that `np.diff()` returns one fewer item in the array it operates on.*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def find_period(t, theta):\n",
" \"\"\"Computes the period of oscillation based on the trajectory of theta.\n",
" \n",
" Parameters\n",
" ==========\n",
" t : array_like, shape(n,)\n",
" An array of monotonically increasing time values.\n",
" theta : array_like, shape(n,)\n",
" An array of values for theta at each time in ``t``.\n",
" \n",
" Returns\n",
" =======\n",
" T : float\n",
" An estimate of the period of oscillation.\n",
" \n",
" \"\"\"\n",
" \n",
" # type your code here\n",
" \n",
" return T"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x_nonlin = odeint(eval_nonlin_rhs, [np.deg2rad(1), 0], t)\n",
"T_d = find_period(t, x_nonlin[:, 0])\n",
"T_d"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# repeat for 10 and 20 degrees"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Question\n",
"\n",
"Using your plot above plot the linear and non-linear time repsonses for $\\theta_0=20 \\textrm{ deg}$. What do you observe?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Type your code here"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"# Question\n",
"\n",
"Make a plot of the period vs $\\theta_0=1,2,..,25$ for the non-linear system. Also plot a horizontal line for the linear period for comparison using `ax.axhline()`.\n",
"\n",
"*Hint: Use a `for` loop to iterate through `np.arange(1, 25)` and collect your results in an initially empty list with `periods.append()`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"initial_thetas = np.arange(1, 25) # deg\n",
"\n",
"periods = []\n",
"\n",
"# type your for loop here"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, ax = plt.subplots(1, 1)\n",
"\n",
"# type your plot commands here using `ax.plot()` and `ax.axhline()`\n",
"\n",
"ax.set_xlabel(r'$\\theta_0$ [rad]')\n",
"ax.set_ylabel(r'$T$ Period [s]')\n",
"ax.legend(['Non-linear', 'Linear']);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Question\n",
"\n",
"Compare the period predicted by your model to the period measured in class. If it does not match, what are the possible explanations? Is the linear model a good model to use for predicting motion of the system?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Type your answer here.*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Question\n",
"\n",
"Derive the expression for the magnitude of the velocity of the center of mass of the book and the height value used in the potential energy expression."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Type answer here*"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [conda root]",
"language": "python",
"name": "conda-root-py"
},
"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.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 1
}