From 8474b2321021d77b718dae98d6d7e681ae75ad55 Mon Sep 17 00:00:00 2001 From: Andrew Date: Thu, 14 Sep 2023 15:37:30 +0100 Subject: [PATCH] WIP --- .../notebooks/01_vib/vib_app.ipynb | 5486 +++++++++++++++++ .../notebooks/01_vib/vib_gen.ipynb | 2250 +++++++ 2 files changed, 7736 insertions(+) create mode 100644 fdm-jupyter-book/notebooks/01_vib/vib_app.ipynb create mode 100644 fdm-jupyter-book/notebooks/01_vib/vib_gen.ipynb diff --git a/fdm-jupyter-book/notebooks/01_vib/vib_app.ipynb b/fdm-jupyter-book/notebooks/01_vib/vib_app.ipynb new file mode 100644 index 00000000..485e0ea0 --- /dev/null +++ b/fdm-jupyter-book/notebooks/01_vib/vib_app.ipynb @@ -0,0 +1,5486 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "ename": "ModuleNotFoundError", + "evalue": "No module named 'mayavi'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[1], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mmayavi\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m mlab\n\u001b[1;32m 2\u001b[0m mlab\u001b[38;5;241m.\u001b[39minit_notebook()\n\u001b[1;32m 3\u001b[0m mlab\u001b[38;5;241m.\u001b[39mtest_plot3d()\n", + "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'mayavi'" + ] + } + ], + "source": [ + "from mayavi import mlab\n", + "mlab.init_notebook()\n", + "mlab.test_plot3d()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "# Applications of vibration models\n", + "
\n", + "\n", + "The following text derives some of the most well-known physical\n", + "problems that lead to second-order ODE models of the type addressed in\n", + "this book. We consider a simple spring-mass system; thereafter\n", + "extended with nonlinear spring, damping, and external excitation; a\n", + "spring-mass system with sliding friction; a simple and a physical\n", + "(classical) pendulum; and an elastic pendulum.\n", + "\n", + "## Oscillating mass attached to a spring\n", + "
\n", + "\n", + "\n", + "\n", + "
\n", + "\n", + "

Simple oscillating mass.

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "The most fundamental mechanical vibration system is depicted in [Figure](#vib:app:mass_spring:fig). A body with mass $m$ is attached to a\n", + "spring and can move horizontally without friction (in the wheels). The\n", + "position of the body is given by the vector $\\textbf{r}(t) = u(t)\\textbf{i}$, where\n", + "$\\textbf{i}$ is a unit vector in $x$ direction.\n", + "There is\n", + "only one force acting on the body: a spring force $\\textbf{F}_s =-ku\\textbf{i}$, where\n", + "$k$ is a constant. The point $x=0$, where $u=0$, must therefore\n", + "correspond to the body's position\n", + "where the spring is neither extended nor compressed, so the force\n", + "vanishes.\n", + "\n", + "The basic physical principle that governs the motion of the body is\n", + "Newton's second law of motion: $\\textbf{F}=m\\textbf{a}$, where\n", + "$\\textbf{F}$ is the sum of forces on the body, $m$ is its mass, and $\\textbf{a}=\\ddot{r}$\n", + "is the acceleration. We use the dot for differentiation with respect\n", + "to time, which is\n", + "usual in mechanics. Newton's second law simplifies here\n", + "to $-\\textbf{F}_s=m\\ddot u\\textbf{i}$, which translates to" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "-ku = m\\ddot u\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Two initial conditions are needed: $u(0)=I$, $\\dot u(0)=V$.\n", + "The ODE problem is normally written as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "m\\ddot u + ku = 0,\\quad u(0)=I,\\ \\dot u(0)=V\\thinspace .\n", + "\\label{vib:app:mass_spring:eqx} \\tag{1}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It is\n", + "not uncommon to divide by $m$\n", + "and introduce the frequency $\\omega = \\sqrt{k/m}$:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\ddot u + \\omega^2 u = 0,\\quad u(0)=I,\\ \\dot u(0)=V\\thinspace .\n", + "\\label{vib:app:mass_spring:equ} \\tag{2}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is the model problem in the first part of this chapter, with the\n", + "small difference that we write the time derivative of $u$ with a dot\n", + "above, while we used $u^{\\prime}$ and $u^{\\prime\\prime}$ in previous\n", + "parts of the book.\n", + "\n", + "\n", + "Since only one scalar mathematical quantity, $u(t)$, describes the\n", + "complete motion, we say that the mechanical system has one degree of freedom\n", + "(DOF).\n", + "\n", + "### Scaling\n", + "\n", + "For numerical simulations it is very convenient to scale\n", + "([2](#vib:app:mass_spring:equ)) and thereby get rid of the problem of\n", + "finding relevant values for all the parameters $m$, $k$, $I$, and $V$.\n", + "Since the amplitude of the oscillations are dictated by $I$ and $V$\n", + "(or more precisely, $V/\\omega$), we scale $u$ by $I$ (or $V/\\omega$ if\n", + "$I=0$):" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\bar u = \\frac{u}{I},\\quad \\bar t = \\frac{t}{t_c}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The time scale $t_c$ is normally chosen as the inverse period $2\\pi/\\omega$ or\n", + "angular frequency $1/\\omega$, most often as $t_c=1/\\omega$.\n", + "Inserting the dimensionless quantities $\\bar u$ and $\\bar t$ in\n", + "([2](#vib:app:mass_spring:equ)) results in the scaled problem" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\frac{d^2\\bar u}{d\\bar t^2} + \\bar u = 0,\\quad \\bar u(0)=1,\\ \\frac{\\bar u}{\\bar t}(0)=\\beta = \\frac{V}{I\\omega},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $\\beta$ is a dimensionless number. Any motion that starts from rest\n", + "($V=0$) is free of parameters in the scaled model!\n", + "\n", + "### The physics\n", + "\n", + "The typical physics of the system in [Figure](#vib:app:mass_spring:fig) can be described as follows. Initially,\n", + "we displace the body to some position $I$, say at rest ($V=0$). After\n", + "releasing the body, the spring, which is extended, will act with a\n", + "force $-kI\\ii$ and pull the body to the left. This force causes an\n", + "acceleration and therefore increases velocity. The body passes the\n", + "point $x=0$, where $u=0$, and the spring will then be compressed and\n", + "act with a force $kx\\ii$ against the motion and cause retardation. At\n", + "some point, the motion stops and the velocity is zero, before the\n", + "spring force $kx\\ii$ has worked long enough to push the body in\n", + "positive direction. The result is that the body accelerates back and\n", + "forth. As long as there is no friction forces to damp the motion, the\n", + "oscillations will continue forever.\n", + "\n", + "## General mechanical vibrating system\n", + "
\n", + "\n", + "\n", + "\n", + "
\n", + "\n", + "

General oscillating system.

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "The mechanical system in [Figure](#vib:app:mass_spring:fig) can easily be\n", + "extended to the more general system in [Figure](#vib:app:mass_gen:fig),\n", + "where the body is attached to a spring and a dashpot, and also subject\n", + "to an environmental force $F(t)\\textbf{i}$. The system has still only one\n", + "degree of freedom since the body can only move back and forth parallel to\n", + "the $x$ axis. The spring force was linear, $\\textbf{F}_s=-ku\\textbf{i}$,\n", + "in the section [Oscillating mass attached to a spring](#vib:app:mass_spring), but in more general cases it can\n", + "depend nonlinearly on the position. We therefore set $\\textbf{F}_s=s(u)\\textbf{i}$.\n", + "The dashpot, which acts\n", + "as a damper, results in a force $\\textbf{F}_d$ that depends on the body's\n", + "velocity $\\dot{u}$ and that always acts against the motion.\n", + "The mathematical model of the force is written $\\textbf{F}_d =f(\\dot u)\\textbf{i}$.\n", + "A positive $\\dot{u}$ must result in a force acting in the positive $x$\n", + "direction.\n", + "Finally, we have the external environmental force $\\textbf{F}_e = F(t)\\textbf{i}$.\n", + "\n", + "Newton's second law of motion now involves three forces:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "F(t)\\textbf{i} - f(\\dot u)\\textbf{i} - s(u)\\textbf{i} = m\\ddot{u} \\textbf{i} \\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The common mathematical form of the ODE problem is" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "m\\ddot u + f(\\dot u) + s(u) = F(t),\\quad u(0)=I,\\ \\dot u(0)=V\\thinspace .\n", + "\\label{vib:app:mass_gen:equ} \\tag{3}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is the generalized problem treated in the last part of the\n", + "present chapter, but with prime denoting the derivative instead of the dot.\n", + "\n", + "The most common models for the spring and dashpot are linear: $f(\\dot u)\n", + "=b\\dot u$ with a constant $b\\geq 0$, and $s(u)=ku$ for a constant $k$.\n", + "\n", + "### Scaling\n", + "\n", + "A specific scaling requires specific choices of $f$, $s$, and $F$.\n", + "Suppose we have" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "f(\\dot u) = b|\\dot u|\\dot u,\\quad s(u)=ku,\\quad F(t)=A\\sin(\\phi t)\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We introduce dimensionless variables as usual, $\\bar u = u/u_c$ and\n", + "$\\bar t = t/t_c$. The scale $u_c$ depends both on the initial conditions\n", + "and $F$, but as time grows, the effect of the initial conditions die out\n", + "and $F$ will drive the motion. Inserting $\\bar u$ and $\\bar t$ in the\n", + "ODE gives" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "m\\frac{u_c}{t_c^2}\\frac{d^2\\bar u}{d\\bar t^2}\n", + "+ b\\frac{u_c^2}{t_c^2}\\left\\vert\\frac{d\\bar u}{d\\bar t}\\right\\vert\n", + "\\frac{d\\bar u}{d\\bar t} + ku_c\\bar u = A\\sin(\\phi t_c\\bar t)\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We divide by $u_c/t_c^2$ and demand the coefficients of the\n", + "$\\bar u$ and the forcing term from $F(t)$ to have unit coefficients.\n", + "This leads to the scales" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "t_c = \\sqrt{\\frac{m}{k}},\\quad u_c = \\frac{A}{k}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The scaled ODE becomes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\frac{d^2\\bar u}{d\\bar t^2}\n", + "+ 2\\beta\\left\\vert\\frac{d\\bar u}{d\\bar t}\\right\\vert\n", + "\\frac{d\\bar u}{d\\bar t} + \\bar u = \\sin(\\gamma\\bar t),\n", + "\\label{vib:app:mass_gen:scaled} \\tag{4}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where there are two dimensionless numbers:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\beta = \\frac{Ab}{2mk},\\quad\\gamma =\\phi\\sqrt{\\frac{m}{k}}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The $\\beta$ number measures the size of the damping term (relative to unity)\n", + "and is assumed to be small, basically because $b$ is small. The $\\phi$\n", + "number is the ratio of the time scale of free vibrations and the time scale\n", + "of the forcing.\n", + "The scaled initial conditions have two other dimensionless numbers\n", + "as values:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\bar u(0) = \\frac{Ik}{A},\\quad \\frac{d\\bar u}{d\\bar t}=\\frac{t_c}{u_c}V = \\frac{V}{A}\\sqrt{mk}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## A sliding mass attached to a spring\n", + "
\n", + "\n", + "Consider a variant of the oscillating body in the section [Oscillating mass attached to a spring](#vib:app:mass_spring)\n", + "and [Figure](#vib:app:mass_spring:fig): the body rests on a flat\n", + "surface, and there is sliding friction between the body and the surface.\n", + "[Figure](#vib:app:mass_sliding:fig) depicts the problem.\n", + "\n", + "\n", + "\n", + "
\n", + "\n", + "

Sketch of a body sliding on a surface.

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "The body is attached to a spring with spring force $-s(u)\\ii$.\n", + "The friction force is proportional to the normal force on the surface,\n", + "$-mg\\jj$, and given by $-f(\\dot u)\\ii$, where" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "f(\\dot u) = \\left\\lbrace\\begin{array}{ll}\n", + "-\\mu mg,& \\dot u < 0,\\\\ \n", + "\\mu mg, & \\dot u > 0,\\\\ \n", + "0, & \\dot u=0\n", + "\\end{array}\\right.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here, $\\mu$ is a friction coefficient. With the signum function" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\mbox{sign(x)} = \\left\\lbrace\\begin{array}{ll}\n", + "-1,& x < 0,\\\\ \n", + "1, & x > 0,\\\\ \n", + "0, & x=0\n", + "\\end{array}\\right.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "we can simply write $f(\\dot u) = \\mu mg\\,\\hbox{sign}(\\dot u)$\n", + "(the sign function is implemented by `numpy.sign`).\n", + "\n", + "The equation of motion becomes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "m\\ddot u + \\mu mg\\hbox{sign}(\\dot u) + s(u) = 0,\\quad u(0)=I,\\ \\dot u(0)=V\\thinspace .\n", + "\\label{vib:app:mass_sliding:equ} \\tag{5}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## A jumping washing machine\n", + "
\n", + "\n", + "A washing machine is placed on four springs with efficient dampers.\n", + "If the machine contains just a few clothes, the circular motion of\n", + "the machine induces a sinusoidal external force from the floor and the machine will\n", + "jump up and down if the frequency of the external force is close to\n", + "the natural frequency of the machine and its spring-damper system.\n", + "\n", + "[hpl 1: Not finished. This is a good example on resonance.]\n", + "\n", + "\n", + "\n", + "## Motion of a pendulum\n", + "
\n", + "\n", + "### Simple pendulum\n", + "\n", + "A classical problem in mechanics is the motion of a pendulum. We first\n", + "consider a [simplified pendulum](https://en.wikipedia.org/wiki/Pendulum) (sometimes also called a\n", + "mathematical pendulum): a small body of mass $m$ is\n", + "attached to a massless wire and can oscillate back and forth in the\n", + "gravity field. [Figure](#vib:app:pendulum:fig_problem) shows a sketch\n", + "of the problem.\n", + "\n", + "\n", + "\n", + "
\n", + "\n", + "

Sketch of a simple pendulum.

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "The motion is governed by Newton's 2nd law, so we need to find\n", + "expressions for the forces and the acceleration. Three forces on the\n", + "body are considered: an unknown force $S$ from the wire, the gravity\n", + "force $mg$, and an air resistance force, $\\frac{1}{2}C_D\\varrho A|v|v$,\n", + "hereafter called the drag force, directed against the velocity\n", + "of the body. Here, $C_D$ is a drag coefficient, $\\varrho$ is the\n", + "density of air, $A$ is the cross section area of the body, and $v$ is\n", + "the magnitude of the velocity.\n", + "\n", + "We introduce a coordinate system with polar coordinates and unit\n", + "vectors $\\ir$ and $\\ith$ as shown in [Figure](#vib:app:pendulum:fig_forces). The position of the center of mass\n", + "of the body is" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\rpos(t) = x_0\\ii + y_0\\jj + L\\ir,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $\\ii$ and $\\jj$ are unit vectors in the corresponding Cartesian\n", + "coordinate system in the $x$ and $y$ directions, respectively. We have\n", + "that $\\ir = \\cos\\theta\\ii +\\sin\\theta\\jj$.\n", + "\n", + "\n", + "\n", + "
\n", + "\n", + "

Forces acting on a simple pendulum.

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "The forces are now expressed as follows.\n", + "\n", + " * Wire force: $-S\\ir$\n", + "\n", + " * Gravity force: $-mg\\jj = mg(-\\sin\\theta\\,\\ith + \\cos\\theta\\,\\ir)$\n", + "\n", + " * Drag force: $-\\frac{1}{2}C_D\\varrho A |v|v\\,\\ith$\n", + "\n", + "Since a positive velocity means movement in the direction of $\\ith$,\n", + "the drag force must be directed along $-\\ith$ so it works against the\n", + "motion. We assume motion in air so that the added mass effect can\n", + "be neglected (for a spherical body, the added mass is $\\frac{1}{2}\\varrho V$,\n", + "where $V$ is the volume of the body). Also the buoyancy effect\n", + "can be neglected for motion in the air when the density difference\n", + "between the fluid and the body is so significant.\n", + "\n", + "The velocity of the body is found from $\\rpos$:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\v(t) = \\dot\\rpos (t) = \\frac{d}{d\\theta}(x_0\\ii + y_0\\jj + L\\ir)\\frac{d\\theta}{dt} = L\\dot\\theta\\ith,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "since $\\frac{d}{d\\theta}\\ir = \\ith$. mathcal{I}_t follows that $v=|\\v|=L\\dot\\theta$.\n", + "The acceleration is" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\acc(t) = \\dot\\v(r) = \\frac{d}{dt}(L\\dot\\theta\\ith)\n", + "= L\\ddot\\theta\\ith + L\\dot\\theta\\frac{d\\ith}{d\\theta}\\dot\\theta =\n", + "= L\\ddot\\theta\\ith - L\\dot\\theta^2\\ir,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "since $\\frac{d}{d\\theta}\\ith = -\\ir$.\n", + "\n", + "Newton's 2nd law of motion becomes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "-S\\ir + mg(-\\sin\\theta\\,\\ith + \\cos\\theta\\,\\ir) -\n", + "\\frac{1}{2}C_D\\varrho AL^2|\\dot\\theta|\\dot\\theta\\,\\ith\n", + "= mL\\ddot\\theta\\dot\\theta\\,\\ith - L\\dot\\theta^2\\ir,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "leading to two component equations" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "-S + mg\\cos\\theta = -L\\dot\\theta^2,\n", + "\\label{vib:app:pendulum:ir} \\tag{6}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "-mg\\sin\\theta - \\frac{1}{2}C_D\\varrho AL^2|\\dot\\theta|\\dot\\theta\n", + "= mL\\ddot\\theta\\thinspace .\n", + "\\label{vib:app:pendulum:ith} \\tag{7}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "From ([6](#vib:app:pendulum:ir)) we get an expression for\n", + "$S=mg\\cos\\theta + L\\dot\\theta^2$, and from ([7](#vib:app:pendulum:ith))\n", + "we get a differential equation for the angle $\\theta(t)$. This latter\n", + "equation is ordered as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "m\\ddot\\theta + \\frac{1}{2}C_D\\varrho AL|\\dot\\theta|\\dot\\theta\n", + "+ \\frac{mg}{L}\\sin\\theta = 0\\thinspace .\n", + "\\label{vib:app:pendulum:thetaeq} \\tag{8}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Two initial conditions are needed: $\\theta=\\Theta$ and $\\dot\\theta = \\Omega$.\n", + "Normally, the pendulum motion is started from rest, which means $\\Omega =0$.\n", + "\n", + "Equation ([8](#vib:app:pendulum:thetaeq)) fits the general model\n", + "used in ([vib:ode2](#vib:ode2)) in the section [vib:model2](#vib:model2) if we define\n", + "$u=\\theta$, $f(u^{\\prime}) = \\frac{1}{2}C_D\\varrho AL|\\dot u|\\dot u$,\n", + "$s(u) = L^{-1}mg\\sin u$, and $F=0$.\n", + "If the body is a sphere with radius $R$, we can take $C_D=0.4$ and $A=\\pi R^2$.\n", + "[Exercise 4: Simulate a simple pendulum](#vib:exer:pendulum_simple) asks you to scale the equations\n", + "and carry out specific simulations with this model.\n", + "\n", + "### Physical pendulum\n", + "\n", + "The motion of a compound or physical pendulum where the wire is a rod with\n", + "mass, can be modeled very similarly. The governing equation is\n", + "$I\\acc = \\boldsymbol{T}$ where $I$ is the moment of inertia of the entire body about\n", + "the point $(x_0,y_0)$, and $\\boldsymbol{T}$ is the sum of moments of the forces\n", + "with respect to $(x_0,y_0)$. The vector equation reads" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\rpos\\times(-S\\ir + mg(-\\sin\\theta\\ith + \\cos\\theta\\ir) -\n", + "\\frac{1}{2}C_D\\varrho AL^2|\\dot\\theta|\\dot\\theta\\ith)\n", + "= I(L\\ddot\\theta\\dot\\theta\\ith - L\\dot\\theta^2\\ir)\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The component equation in $\\ith$ direction gives the equation of motion\n", + "for $\\theta(t)$:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "I\\ddot\\theta + \\frac{1}{2}C_D\\varrho AL^3|\\dot\\theta|\\dot\\theta\n", + "+ mgL\\sin\\theta = 0\\thinspace .\n", + "\\label{vib:app:pendulum:thetaeq_physical} \\tag{9}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Dynamic free body diagram during pendulum motion\n", + "
\n", + "\n", + "\n", + "Usually one plots the mathematical quantities as functions of time to\n", + "visualize the solution of ODE models. [Exercise 4: Simulate a simple pendulum](#vib:exer:pendulum_simple) asks you to do this for the motion of a\n", + "pendulum in the previous section. However, sometimes it is more\n", + "instructive to look at other types of visualizations. For example, we\n", + "have the pendulum and the free body diagram in Figures\n", + "[vib:app:pendulum:fig_problem](#vib:app:pendulum:fig_problem) and\n", + "[vib:app:pendulum:fig_forces](#vib:app:pendulum:fig_forces). We may think of these figures as\n", + "animations in time instead. Especially the free body diagram will show both the\n", + "motion of the pendulum *and* the size of the forces during the motion.\n", + "The present section exemplifies how to make such a dynamic body\n", + "diagram.\n", + "Two typical snapshots of free body diagrams are displayed below\n", + "(the drag force is magnified 5 times to become more visual!).\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "
\n", + "\n", + "
\n", + "

The drag force is magnified 5 times!

\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from IPython.display import HTML\n", + "_s = \"\"\"\n", + "
\n", + "\n", + "
\n", + "

The drag force is magnified 5 times!

\n", + "\n", + "\n", + "\n", + "\n", + "\"\"\"\n", + "HTML(_s)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "\n", + "Dynamic physical sketches, coupled to the numerical solution of\n", + "differential equations, requires a program to produce a sketch for\n", + "the situation at each time level.\n", + "[Pysketcher](https://github.com/hplgit/pysketcher) is such a tool.\n", + "In fact (and not surprising!) Figures [vib:app:pendulum:fig_problem](#vib:app:pendulum:fig_problem) and\n", + "[vib:app:pendulum:fig_forces](#vib:app:pendulum:fig_forces) were drawn using Pysketcher.\n", + "The details of the drawings are explained in the\n", + "[Pysketcher tutorial](http://hplgit.github.io/pysketcher/doc/web/index.html).\n", + "Here, we outline how this type of sketch can be used to create an animated\n", + "free body diagram during the motion of a pendulum.\n", + "\n", + "Pysketcher is actually a layer of useful abstractions on top of\n", + "standard plotting packages. This means that we in fact apply Matplotlib\n", + "to make the animated free body diagram, but instead of dealing with a wealth\n", + "of detailed Matplotlib commands, we can express the drawing in terms of\n", + "more high-level objects, e.g., objects for the wire, angle $\\theta$,\n", + "body with mass $m$, arrows for forces, etc. When the position of these\n", + "objects are given through variables, we can just couple those variables\n", + "to the dynamic solution of our ODE and thereby make a unique drawing\n", + "for each $\\theta$ value in a simulation.\n", + "\n", + "### Writing the solver\n", + "\n", + "Let us start with the most familiar part of the current problem:\n", + "writing the solver function. We use Odespy for this purpose.\n", + "We also work with dimensionless equations. Since $\\theta$ can be\n", + "viewed as dimensionless, we only need to introduce a dimensionless time,\n", + "here taken as $\\bar t = t/\\sqrt{L/g}$.\n", + "The resulting dimensionless mathematical model for $\\theta$,\n", + "the dimensionless angular velocity $\\omega$, the\n", + "dimensionless wire force $\\bar S$, and the dimensionless\n", + "drag force $\\bar D$ is then" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\frac{d\\omega}{d\\bar t} = - \\alpha|\\omega|\\omega - \\sin\\theta,\n", + "\\label{vib:app:pendulum_bodydia:eqth} \\tag{10}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\frac{d\\theta}{d\\bar t} = \\omega,\n", + "\\label{vib:app:pendulum_bodydia:eqomega} \\tag{11}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\bar S = \\omega^2 + \\cos\\theta,\n", + "\\label{vib:app:pendulum_bodydia:eqS} \\tag{12}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\bar D = -\\alpha |\\omega|\\omega,\n", + "\\label{vib:app:pendulum_bodydia:eqD} \\tag{13}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "with" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\alpha = \\frac{C_D\\varrho\\pi R^2L}{2m}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "as a dimensionless parameter expressing the ratio of the drag force and\n", + "the gravity force. The dimensionless $\\omega$ is made non-dimensional\n", + "by the time, so $\\omega\\sqrt{L/g}$ is the corresponding angular\n", + "frequency with dimensions.\n", + "\n", + "\n", + "\n", + "A suitable function for computing\n", + "([10](#vib:app:pendulum_bodydia:eqth))-([13](#vib:app:pendulum_bodydia:eqD))\n", + "is listed below." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "! Should I convert this to Devito?\n", + "\n", + "def simulate(alpha, Theta, dt, T):\n", + " import odespy\n", + "\n", + " def f(u, t, alpha):\n", + " omega, theta = u\n", + " return [-alpha*omega*abs(omega) - sin(theta),\n", + " omega]\n", + "\n", + " import numpy as np\n", + " Nt = int(round(T/float(dt)))\n", + " t = np.linspace(0, Nt*dt, Nt+1)\n", + " solver = odespy.RK4(f, f_args=[alpha])\n", + " solver.set_initial_condition([0, Theta])\n", + " u, t = solver.solve(\n", + " t, terminate=lambda u, t, n: abs(u[n,1]) < 1E-3)\n", + " omega = u[:,0]\n", + " theta = u[:,1]\n", + " S = omega**2 + np.cos(theta)\n", + " drag = -alpha*np.abs(omega)*omega\n", + " return t, theta, omega, S, drag" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Drawing the free body diagram\n", + "\n", + "The `sketch` function below applies Pysketcher objects to build\n", + "a diagram like that in [Figure](#vib:app:pendulum:fig_forces),\n", + "except that we have removed the rotation point $(x_0,y_0)$ and\n", + "the unit vectors in polar coordinates as these objects are not\n", + "important for an animated free body diagram." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "try:\n", + " from pysketcher import *\n", + "except ImportError:\n", + " print 'Pysketcher must be installed from'\n", + " print 'https://github.com/hplgit/pysketcher'\n", + " sys.exit(1)\n", + "\n", + "# Overall dimensions of sketch\n", + "H = 15.\n", + "W = 17.\n", + "\n", + "drawing_tool.set_coordinate_system(\n", + " xmin=0, xmax=W, ymin=0, ymax=H,\n", + " axis=False)\n", + "\n", + "def sketch(theta, S, mg, drag, t, time_level):\n", + " \"\"\"\n", + " Draw pendulum sketch with body forces at a time level\n", + " corresponding to time t. The drag force is in\n", + " drag[time_level], the force in the wire is S[time_level],\n", + " the angle is theta[time_level].\n", + " \"\"\"\n", + " import math\n", + " a = math.degrees(theta[time_level]) # angle in degrees\n", + " L = 0.4*H # Length of pendulum\n", + " P = (W/2, 0.8*H) # Fixed rotation point\n", + "\n", + " mass_pt = path.geometric_features()['end']\n", + " rod = Line(P, mass_pt)\n", + "\n", + " mass = Circle(center=mass_pt, radius=L/20.)\n", + " mass.set_filled_curves(color='blue')\n", + " rod_vec = rod.geometric_features()['end'] - \\\n", + " rod.geometric_features()['start']\n", + " unit_rod_vec = unit_vec(rod_vec)\n", + " mass_symbol = Text('$m$', mass_pt + L/10*unit_rod_vec)\n", + "\n", + " rod_start = rod.geometric_features()['start'] # Point P\n", + " vertical = Line(rod_start, rod_start + point(0,-L/3))\n", + "\n", + " def set_dashed_thin_blackline(*objects):\n", + " \"\"\"Set linestyle of objects to dashed, black, width=1.\"\"\"\n", + " for obj in objects:\n", + " obj.set_linestyle('dashed')\n", + " obj.set_linecolor('black')\n", + " obj.set_linewidth(1)\n", + "\n", + " set_dashed_thin_blackline(vertical)\n", + " set_dashed_thin_blackline(rod)\n", + " angle = Arc_wText(r'$\\theta$', rod_start, L/6, -90, a,\n", + " text_spacing=1/30.)\n", + "\n", + " magnitude = 1.2*L/2 # length of a unit force in figure\n", + " force = mg[time_level] # constant (scaled eq: about 1)\n", + " force *= magnitude\n", + " mg_force = Force(mass_pt, mass_pt + force*point(0,-1),\n", + " '', text_pos='end')\n", + " force = S[time_level]\n", + " force *= magnitude\n", + " rod_force = Force(mass_pt, mass_pt - force*unit_vec(rod_vec),\n", + " '', text_pos='end',\n", + " text_spacing=(0.03, 0.01))\n", + " force = drag[time_level]\n", + " force *= magnitude\n", + " air_force = Force(mass_pt, mass_pt -\n", + " force*unit_vec((rod_vec[1], -rod_vec[0])),\n", + " '', text_pos='end',\n", + " text_spacing=(0.04,0.005))\n", + "\n", + " body_diagram = Composition(\n", + " {'mg': mg_force, 'S': rod_force, 'air': air_force,\n", + " 'rod': rod, 'body': mass\n", + " 'vertical': vertical, 'theta': angle,})\n", + "\n", + " body_diagram.draw(verbose=0)\n", + " drawing_tool.savefig('tmp_%04d.png' % time_level, crop=False)\n", + " # (No cropping: otherwise movies will be very strange!)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Making the animated free body diagram\n", + "\n", + "mathcal{I}_t now remains to couple the `simulate` and `sketch` functions.\n", + "We first run `simulate`:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "from math import pi, radians, degrees\n", + "import numpy as np\n", + "alpha = 0.4\n", + "period = 2*pi # Use small theta approximation\n", + "T = 12*period # Simulate for 12 periods\n", + "dt = period/40 # 40 time steps per period\n", + "a = 70 # Initial amplitude in degrees\n", + "Theta = radians(a)\n", + "\n", + "t, theta, omega, S, drag = simulate(alpha, Theta, dt, T)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The next step is to run through the time levels in the simulation and\n", + "make a sketch at each level:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "for time_level, t_ in enumerate(t):\n", + " sketch(theta, S, mg, drag, t_, time_level)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The individual sketches are (by the `sketch` function) saved in files\n", + "with names `tmp_%04d.png`. These can be combined to videos using\n", + "(e.g.) `ffmpeg`. A complete function `animate` for running the\n", + "simulation and creating video files is\n", + "listed below." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "def animate():\n", + " # Clean up old plot files\n", + " import os, glob\n", + " for filename in glob.glob('tmp_*.png') + glob.glob('movie.*'):\n", + " os.remove(filename)\n", + " # Solve problem\n", + " from math import pi, radians, degrees\n", + " import numpy as np\n", + " alpha = 0.4\n", + " period = 2*pi # Use small theta approximation\n", + " T = 12*period # Simulate for 12 periods\n", + " dt = period/40 # 40 time steps per period\n", + " a = 70 # Initial amplitude in degrees\n", + " Theta = radians(a)\n", + "\n", + " t, theta, omega, S, drag = simulate(alpha, Theta, dt, T)\n", + "\n", + " # Visualize drag force 5 times as large\n", + " drag *= 5\n", + " mg = np.ones(S.size) # Gravity force (needed in sketch)\n", + "\n", + " # Draw animation\n", + " import time\n", + " for time_level, t_ in enumerate(t):\n", + " sketch(theta, S, mg, drag, t_, time_level)\n", + " time.sleep(0.2) # Pause between each frame on the screen\n", + "\n", + " # Make videos\n", + " prog = 'ffmpeg'\n", + " filename = 'tmp_%04d.png'\n", + " fps = 6\n", + " codecs = {'flv': 'flv', 'mp4': 'libx264',\n", + " 'webm': 'libvpx', 'ogg': 'libtheora'}\n", + " for ext in codecs:\n", + " lib = codecs[ext]\n", + " cmd = '%(prog)s -i %(filename)s -r %(fps)s ' % vars()\n", + " cmd += '-vcodec %(lib)s movie.%(ext)s' % vars()\n", + " print(cmd)\n", + " os.system(cmd)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Motion of an elastic pendulum\n", + "
\n", + "\n", + "\n", + "Consider a pendulum as in [Figure](#vib:app:pendulum:fig_problem), but\n", + "this time the wire is elastic. The length of the wire when it is not\n", + "stretched is $L_0$, while $L(t)$ is the stretched\n", + "length at time $t$ during the motion.\n", + "\n", + "Stretching the elastic wire a distance $\\Delta L$ gives rise to a\n", + "spring force $k\\Delta L$ in the opposite direction of the\n", + "stretching. Let $\\boldsymbol{n}$ be a unit normal vector along the wire\n", + "from the point $\\rpos_0=(x_0,y_0)$ and in the direction of $\\ith$, see\n", + "[Figure](#vib:app:pendulum:fig_forces) for definition of $(x_0,y_0)$\n", + "and $\\ith$. Obviously, we have $\\boldsymbol{n}=\\ith$, but in this modeling\n", + "of an elastic pendulum we do not need polar coordinates. Instead, it\n", + "is more straightforward to develop the equation in Cartesian\n", + "coordinates.\n", + "\n", + "A mathematical expression for $\\boldsymbol{n}$ is" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{n} = \\frac{\\rpos-\\rpos_0}{L(t)},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $L(t)=||\\rpos-\\rpos_0||$ is the current length of the elastic wire.\n", + "The position vector $\\rpos$ in Cartesian coordinates reads\n", + "$\\rpos(t) = x(t)\\ii + y(t)\\jj$, where $\\ii$ and $\\jj$ are unit vectors\n", + "in the $x$ and $y$ directions, respectively.\n", + "It is convenient to introduce the Cartesian components $n_x$ and $n_y$\n", + "of the normal vector:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{n} = \\frac{\\rpos-\\rpos_0}{L(t)} = \\frac{x(t)-x_0}{L(t)}\\ii + \\frac{y(t)-y_0}{L(t)}\\jj = n_x\\ii + n_y\\jj\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The stretch $\\Delta L$ in the wire is" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\Delta t = L(t) - L_0\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The force in the wire is then $-S\\boldsymbol{n}=-k\\Delta L\\boldsymbol{n}$.\n", + "\n", + "The other forces are the gravity and the air resistance, just as in\n", + "[Figure](#vib:app:pendulum:fig_forces). For motion in air we can\n", + "neglect the added mass and buoyancy effects. The main difference is\n", + "that we have a *model* for $S$ in terms of the motion (as soon as we\n", + "have expressed $\\Delta L$ by $\\rpos$). For simplicity, we drop the air\n", + "resistance term (but [Exercise 6: Simulate an elastic pendulum with air resistance](#vib:exer:pendulum_elastic_drag) asks\n", + "you to include it).\n", + "\n", + "Newton's second law of motion applied to the body now results in" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "m\\ddot\\rpos = -k(L-L_0)\\boldsymbol{n} - mg\\jj\n", + "\\label{vib:app:pendulum_elastic:eq1} \\tag{14}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The two components of\n", + "([14](#vib:app:pendulum_elastic:eq1)) are" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\ddot x = -\\frac{k}{m}(L-L_0)n_x,\n", + "\\label{_auto1} \\tag{15}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\label{vib:app:pendulum_elastic:eq2a} \\tag{16} \n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\ddot y = - \\frac{k}{m}(L-L_0)n_y - g\n", + "\\label{vib:app:pendulum_elastic:eq2b} \\tag{17}\\thinspace .\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Remarks about an elastic vs a non-elastic pendulum\n", + "\n", + "Note that the derivation of the ODEs for an elastic pendulum is more\n", + "straightforward than for a classical, non-elastic pendulum,\n", + "since we avoid the details\n", + "with polar coordinates, but instead work with Newton's second law\n", + "directly in Cartesian coordinates. The reason why we can do this is that\n", + "the elastic pendulum undergoes a general two-dimensional motion where\n", + "all the forces are known or expressed as functions of $x(t)$ and $y(t)$,\n", + "such that we get two ordinary differential equations.\n", + "The motion of the non-elastic pendulum, on the other hand, is constrained:\n", + "the body has to move along a circular path, and the force $S$ in the\n", + "wire is unknown.\n", + "\n", + "The non-elastic pendulum therefore leads to\n", + "a *differential-algebraic* equation, i.e., ODEs for $x(t)$ and $y(t)$\n", + "combined with an extra constraint $(x-x_0)^2 + (y-y_0)^2 = L^2$\n", + "ensuring that the motion takes place along a circular path.\n", + "The extra constraint (equation) is compensated by an extra unknown force\n", + "$-S\\boldsymbol{n}$. Differential-algebraic equations are normally hard\n", + "to solve, especially with pen and paper.\n", + "Fortunately, for the non-elastic pendulum we can do a\n", + "trick: in polar coordinates the unknown force $S$ appears only in the\n", + "radial component of Newton's second law, while the unknown\n", + "degree of freedom for describing the motion, the angle $\\theta(t)$,\n", + "is completely governed by the asimuthal component. This allows us to\n", + "decouple the unknowns $S$ and $\\theta$. But this is a kind of trick and\n", + "not a widely applicable method. With an elastic pendulum we use straightforward\n", + "reasoning with Newton's 2nd law and arrive at a standard ODE problem that\n", + "(after scaling) is easy to solve on a computer.\n", + "\n", + "### Initial conditions\n", + "\n", + "What is the initial position of the body? We imagine that first the\n", + "pendulum hangs in equilibrium in its vertical position, and then it is\n", + "displaced an angle $\\Theta$. The equilibrium position is governed\n", + "by the ODEs with the accelerations set to zero.\n", + "The $x$ component leads to $x(t)=x_0$, while the $y$ component gives" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "0 = - \\frac{k}{m}(L-L_0)n_y - g = \\frac{k}{m}(L(0)-L_0) - g\\quad\\Rightarrow\\quad\n", + "L(0) = L_0 + mg/k,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "since $n_y=-11$ in this position. The corresponding $y$ value is then\n", + "from $n_y=-1$:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "y(t) = y_0 - L(0) = y_0 - (L_0 + mg/k)\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let us now choose $(x_0,y_0)$ such that the body is at the origin\n", + "in the equilibrium position:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "x_0 =0,\\quad y_0 = L_0 + mg/k\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Displacing the body an angle $\\Theta$ to the right leads to the\n", + "initial position" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "x(0)=(L_0+mg/k)\\sin\\Theta,\\quad y(0)=(L_0+mg/k)(1-\\cos\\Theta)\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The initial velocities can be set to zero: $x'(0)=y'(0)=0$.\n", + "\n", + "### The complete ODE problem\n", + "\n", + "We can summarize all the equations as follows:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{align*}\n", + "\\ddot x &= -\\frac{k}{m}(L-L_0)n_x,\n", + "\\\\ \n", + "\\ddot y &= -\\frac{k}{m}(L-L_0)n_y - g,\n", + "\\\\ \n", + "L &= \\sqrt{(x-x_0)^2 + (y-y_0)^2},\n", + "\\\\ \n", + "n_x &= \\frac{x-x_0}{L},\n", + "\\\\ \n", + "n_y &= \\frac{y-y_0}{L},\n", + "\\\\ \n", + "x(0) &= (L_0+mg/k)\\sin\\Theta,\n", + "\\\\ \n", + "x'(0) &= 0,\n", + "\\\\ \n", + "y(0) & =(L_0+mg/k)(1-\\cos\\Theta),\n", + "\\\\ \n", + "y'(0) &= 0\\thinspace .\n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We insert $n_x$ and $n_y$ in the ODEs:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\ddot x = -\\frac{k}{m}\\left(1 -\\frac{L_0}{L}\\right)(x-x_0),\n", + "\\label{vib:app:pendulum_elastic:x} \\tag{18}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\ddot y = -\\frac{k}{m}\\left(1 -\\frac{L_0}{L}\\right)(y-y_0) - g,\n", + "\\label{vib:app:pendulum_elastic:y} \\tag{19}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "L = \\sqrt{(x-x_0)^2 + (y-y_0)^2},\n", + "\\label{vib:app:pendulum_elastic:L} \\tag{20}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "x(0) = (L_0+mg/k)\\sin\\Theta,\n", + "\\label{vib:app:pendulum_elastic:x0} \\tag{21}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "x'(0) = 0,\n", + "\\label{vib:app:pendulum_elastic:vx0} \\tag{22}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "y(0) =(L_0+mg/k)(1-\\cos\\Theta),\n", + "\\label{vib:app:pendulum_elastic:y0} \\tag{23}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "y'(0) = 0\\thinspace .\n", + "\\label{vib:app:pendulum_elastic:vy0} \\tag{24}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Scaling\n", + "\n", + "The elastic pendulum model can be used to study both an elastic pendulum\n", + "and a classic, non-elastic pendulum. The latter problem is obtained\n", + "by letting $k\\rightarrow\\infty$. Unfortunately,\n", + "a serious problem with the ODEs\n", + "([18](#vib:app:pendulum_elastic:x))-([19](#vib:app:pendulum_elastic:y)) is that for large $k$, we have a very large factor $k/m$ multiplied by a\n", + "very small number $1-L_0/L$, since for large $k$, $L\\approx L_0$ (very\n", + "small deformations of the wire). The product is subject to\n", + "significant round-off errors for many relevant physical values of\n", + "the parameters. To circumvent the problem, we introduce a scaling. This\n", + "will also remove physical parameters from the problem such that we end\n", + "up with only one dimensionless parameter,\n", + "closely related to the elasticity of the wire. Simulations can then be\n", + "done by setting just this dimensionless parameter.\n", + "\n", + "The characteristic length can be taken such that in equilibrium, the\n", + "scaled length is unity, i.e., the characteristic length is $L_0+mg/k$:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\bar x = \\frac{x}{L_0+mg/k},\\quad \\bar y = \\frac{y}{L_0+mg/k}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We must then also work with the scaled length $\\bar L = L/(L_0+mg/k)$.\n", + "\n", + "Introducing $\\bar t=t/t_c$, where $t_c$ is a characteristic time we\n", + "have to decide upon later, one gets" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{align*}\n", + "\\frac{d^2\\bar x}{d\\bar t^2} &=\n", + "-t_c^2\\frac{k}{m}\\left(1 -\\frac{L_0}{L_0+mg/k}\\frac{1}{\\bar L}\\right)\\bar x,\\\\ \n", + "\\frac{d^2\\bar y}{d\\bar t^2} &=\n", + "-t_c^2\\frac{k}{m}\\left(1 -\\frac{L_0}{L_0+mg/k}\\frac{1}{\\bar L}\\right)(\\bar y-1)\n", + "-t_c^2\\frac{g}{L_0 + mg/k},\\\\ \n", + "\\bar L &= \\sqrt{\\bar x^2 + (\\bar y-1)^2},\\\\ \n", + "\\bar x(0) &= \\sin\\Theta,\\\\ \n", + "\\bar x'(0) &= 0,\\\\ \n", + "\\bar y(0) & = 1 - \\cos\\Theta,\\\\ \n", + "\\bar y'(0) &= 0\\thinspace .\n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For a non-elastic pendulum with small angles, we know that the\n", + "frequency of the oscillations are $\\omega = \\sqrt{L/g}$. mathcal{I}_t is therefore\n", + "natural to choose a similar expression here, either the length in\n", + "the equilibrium position," + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "t_c^2 = \\frac{L_0+mg/k}{g}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "or simply the unstretched length," + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "t_c^2 = \\frac{L_0}{g}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These quantities are not very different (since the elastic model\n", + "is valid only for quite small elongations), so we take the latter as it is\n", + "the simplest one.\n", + "\n", + "The ODEs become" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{align*}\n", + "\\frac{d^2\\bar x}{d\\bar t^2} &=\n", + "-\\frac{L_0k}{mg}\\left(1 -\\frac{L_0}{L_0+mg/k}\\frac{1}{\\bar L}\\right)\\bar x,\\\\ \n", + "\\frac{d^2\\bar y}{d\\bar t^2} &=\n", + "-\\frac{L_0k}{mg}\\left(1 -\\frac{L_0}{L_0+mg/k}\\frac{1}{\\bar L}\\right)(\\bar y-1)\n", + "-\\frac{L_0}{L_0 + mg/k},\\\\ \n", + "\\bar L &= \\sqrt{\\bar x^2 + (\\bar y-1)^2}\\thinspace .\n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can now identify a dimensionless number" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\beta = \\frac{L_0}{L_0 + mg/k} = \\frac{1}{1+\\frac{mg}{L_0k}},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "which is the ratio of the unstretched length and the\n", + "stretched length in equilibrium. The non-elastic pendulum will have\n", + "$\\beta =1$ ($k\\rightarrow\\infty$).\n", + "With $\\beta$ the ODEs read" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\frac{d^2\\bar x}{d\\bar t^2} =\n", + "-\\frac{\\beta}{1-\\beta}\\left(1- \\frac{\\beta}{\\bar L}\\right)\\bar x,\n", + "\\label{vib:app:pendulum_elastic:x:s} \\tag{25}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\frac{d^2\\bar y}{d\\bar t^2} =\n", + "-\\frac{\\beta}{1-\\beta}\\left(1- \\frac{\\beta}{\\bar L}\\right)(\\bar y-1)\n", + "-\\beta,\n", + "\\label{vib:app:pendulum_elastic:y:s} \\tag{26}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\bar L = \\sqrt{\\bar x^2 + (\\bar y-1)^2},\n", + "\\label{vib:app:pendulum_elastic:L:s} \\tag{27}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\bar x(0) = (1+\\epsilon)\\sin\\Theta,\n", + "\\label{vib:app:pendulum_elastic:x0:s} \\tag{28}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\frac{d\\bar x}{d\\bar t}(0) = 0,\n", + "\\label{vib:app:pendulum_elastic:vx0:s} \\tag{29}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\bar y(0) = 1 - (1+\\epsilon)\\cos\\Theta,\n", + "\\label{vib:app:pendulum_elastic:y0:s} \\tag{30}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\frac{d\\bar y}{d\\bar t}(0) = 0,\n", + "\\label{vib:app:pendulum_elastic:vy0:s} \\tag{31}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have here added a parameter $\\epsilon$, which is an additional\n", + "downward stretch of the wire at $t=0$. This parameter makes it possible\n", + "to do a desired test: vertical oscillations of the pendulum. Without\n", + "$\\epsilon$, starting the motion from $(0,0)$ with zero velocity will\n", + "result in $x=y=0$ for all times (also a good test!), but with\n", + "an initial stretch so the body's position is $(0,\\epsilon)$, we\n", + "will have oscillatory vertical motion with amplitude $\\epsilon$ (see\n", + "[Exercise 5: Simulate an elastic pendulum](#vib:exer:pendulum_elastic)).\n", + "\n", + "### Remark on the non-elastic limit\n", + "\n", + "We immediately see that as $k\\rightarrow\\infty$ (i.e., we obtain a non-elastic\n", + "pendulum), $\\beta\\rightarrow 1$, $\\bar L\\rightarrow 1$, and we have\n", + "very small values $1-\\beta\\bar L^{-1}$ divided by very small values\n", + "$1-\\beta$ in the ODEs. However, it turns out that we can set $\\beta$\n", + "very close to one and obtain a path of the body that within the visual\n", + "accuracy of a plot does not show any elastic oscillations.\n", + "(Should the division of very small values become a problem, one can\n", + "study the limit by L'Hospital's rule:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\lim_{\\beta\\rightarrow 1}\\frac{1 - \\beta \\bar L^{-1}}{1-\\beta}\n", + "= \\frac{1}{\\bar L},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "and use the limit $\\bar L^{-1}$ in the ODEs for $\\beta$ values very\n", + "close to 1.)\n", + "\n", + "## Vehicle on a bumpy road\n", + "
\n", + "\n", + "\n", + "\n", + "
\n", + "\n", + "

Sketch of one-wheel vehicle on a bumpy road.

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "We consider a very simplistic vehicle, on one wheel, rolling along a\n", + "bumpy road. The oscillatory nature of the road will induce an external\n", + "forcing on the spring system in the vehicle and cause vibrations.\n", + "[Figure](#vib:app:bumpy:fig:sketch) outlines the situation.\n", + "\n", + "To derive the equation that governs the motion, we must first establish\n", + "the position vector of the black mass at the top of the spring.\n", + "Suppose the spring has length $L$ without any elongation or compression,\n", + "suppose the radius of the wheel is $R$, and suppose the height of the\n", + "black mass at the top is $H$. With the aid of the $\\rpos_0$ vector\n", + "in [Figure](#vib:app:bumpy:fig:sketch), the position $\\rpos$ of\n", + "the center point of the mass is" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\rpos = \\rpos_0 + 2R\\jj + L\\jj + u\\jj + \\frac{1}{2} H\\jj,\\ \n", + "\\label{_auto2} \\tag{32}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $u$ is the elongation or compression in the spring according to\n", + "the (unknown and to be computed) vertical displacement $u$ relative to the\n", + "road. If the vehicle travels\n", + "with constant horizontal velocity $v$ and $h(x)$ is the shape of the\n", + "road, then the vector $\\rpos_0$ is" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\rpos_0 = vt\\ii + h(vt)\\jj,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "if the motion starts from $x=0$ at time $t=0$.\n", + "\n", + "The forces on the mass is the gravity, the spring force, and an optional\n", + "damping force that is proportional to the vertical velocity $\\dot u$. Newton's\n", + "second law of motion then tells that" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "m\\ddot\\rpos = -mg\\jj - s(u) - b\\dot u\\jj\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This leads to" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "m\\ddot u = - s(u) - b\\dot u - mg -mh''(vt)v^2\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To simplify a little bit, we omit the gravity force $mg$ in comparison with\n", + "the other terms. Introducing $u'$ for $\\dot u$ then gives a standard\n", + "damped, vibration equation with external forcing:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "mu'' + bu' + s(u) = -mh''(vt)v^2\\thinspace .\n", + "\\label{_auto3} \\tag{33}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since the road is normally known just as a set of array values, $h''$ must\n", + "be computed by finite differences. Let $\\Delta x$ be the spacing between\n", + "measured values $h_i= h(i\\Delta x)$ on the road. The discrete second-order\n", + "derivative $h''$ reads" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "q_i = \\frac{h_{i-1} - 2h_i + h_{i+1}}{\\Delta x^2}, \\quad i=1,\\ldots,N_x-1\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We may for maximum simplicity set\n", + "the end points as $q_0=q_1$ and $q_{N_x}=q_{N_x-1}$.\n", + "The term $-mh''(vt)v^2$ corresponds to a force with discrete time values" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "F^n = -mq_n v^2,\\quad \\Delta t = v^{-1}\\Delta x\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This force can be directly used in a numerical model" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "[mD_tD_t u + bD_{2t} u + s(u) = F]^n\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Software for computing $u$ and also making an animated sketch of\n", + "the motion like we did in the section [Dynamic free body diagram during pendulum motion](#vib:app:pendulum_bodydia)\n", + "is found in a separate project on the web:\n", + ". You may start looking at the\n", + "\"tutorial\":\n", + "\"http://hplgit.github.io/bumpy/doc/pub/bumpy.pdf\".\n", + "\n", + "## Bouncing ball\n", + "
\n", + "\n", + "A bouncing ball is a ball in free vertically fall until it impacts the\n", + "ground, but during the impact, some kinetic energy is lost, and a new\n", + "motion upwards with reduced velocity starts. After the motion is retarded,\n", + "a new free fall starts, and the process is repeated. At some point the\n", + "velocity close to the ground is so small that the ball is considered\n", + "to be finally at rest.\n", + "\n", + "The motion of the ball falling in air is governed by Newton's second\n", + "law $F=ma$, where $a$ is the acceleration of the body, $m$ is the mass,\n", + "and $F$ is the sum of all forces. Here, we neglect the air resistance\n", + "so that gravity $-mg$ is the only force. The height of the ball is\n", + "denoted by $h$ and $v$ is the velocity. The relations between $h$, $v$, and\n", + "$a$," + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "h'(t)= v(t),\\quad v'(t) = a(t),\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "combined with Newton's second law gives the ODE model" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "h^{\\prime\\prime}(t) = -g,\n", + "\\label{vib:app:bouncing:ball:h2eq} \\tag{34}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "or expressed alternatively as a system of first-order equations:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "v'(t) = -g,\n", + "\\label{vib:app:bouncing:ball:veq} \\tag{35} \n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "h'(t) = v(t)\\thinspace .\n", + "\\label{vib:app:bouncing:ball:heq} \\tag{36}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These equations govern the motion as long as the ball is away from\n", + "the ground by a small distance $\\epsilon_h > 0$. When $h<\\epsilon_h$,\n", + "we have two cases.\n", + "\n", + "1. The ball impacts the ground, recognized by a sufficiently large negative\n", + " velocity ($v<-\\epsilon_v$). The velocity then changes sign and is\n", + " reduced by a factor $C_R$, known as the [coefficient of restitution](http://en.wikipedia.org/wiki/Coefficient_of_restitution).\n", + " For plotting purposes, one may set $h=0$.\n", + "\n", + "2. The motion stops, recognized by a sufficiently small velocity\n", + " ($|v|<\\epsilon_v$) close to the ground.\n", + "\n", + "## Two-body gravitational problem\n", + "
\n", + "\n", + "Consider two astronomical objects $A$ and $B$ that attract each other\n", + "by gravitational forces. $A$ and $B$ could be two stars in a binary\n", + "system, a planet orbiting a star, or a moon orbiting a planet.\n", + "Each object is acted upon by the\n", + "gravitational force due to the other object. Consider motion in a plane\n", + "(for simplicity) and let $(x_A,y_A)$ and $(x_B,y_B)$ be the\n", + "positions of object $A$ and $B$, respectively.\n", + "\n", + "### The governing equations\n", + "\n", + "Newton's second law of motion applied to each object is all we need\n", + "to set up a mathematical model for this physical problem:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "m_A\\ddot\\x_A = \\F,\n", + "\\label{_auto4} \\tag{37}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "m_B\\ddot\\x_B = -\\F,\n", + "\\label{_auto5} \\tag{38}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $F$ is the gravitational force" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\F = \\frac{Gm_Am_B}{||\\rpos||^3}\\rpos,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\rpos(t) = \\x_B(t) - \\x_A(t),\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "and $G$ is the gravitational constant:\n", + "$G=6.674\\cdot 10^{-11}\\hbox{ Nm}^2/\\hbox{kg}^2$.\n", + "\n", + "### Scaling\n", + "\n", + "A problem with these equations is that the parameters are very large\n", + "($m_A$, $m_B$, $||\\rpos||$) or very small ($G$). The rotation time\n", + "for binary stars can be very small and large as well. mathcal{I}_t is therefore\n", + "advantageous to scale the equations.\n", + "A natural length scale could be the initial distance between the objects:\n", + "$L=\\rpos(0)$. We write the dimensionless quantities as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\bar\\x_A = \\frac{\\x_A}{L},\\quad\\bar\\x_B = \\frac{\\x_B}{L},\\quad\n", + "\\bar t = \\frac{t}{t_c}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The gravity force is transformed to" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\F = \\frac{Gm_Am_B}{L^2||\\bar\\rpos||^3}\\bar\\rpos,\\quad \\bar\\rpos = \\bar\\x_B - \\bar\\x_A,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "so the first ODE for $\\x_A$ becomes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\frac{d^2 \\bar\\x_A}{d\\bar t^2} =\n", + "\\frac{Gm_Bt_c^2}{L^3}\\frac{\\bar\\rpos}{||\\bar\\rpos||^3}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Assuming that quantities with a bar and their derivatives are around unity\n", + "in size, it is natural to choose $t_c$ such that the fraction $Gm_Bt_c/L^2=1$:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "t_c = \\sqrt{\\frac{L^3}{Gm_B}}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "From the other equation for $\\x_B$ we get another candidate for $t_c$ with\n", + "$m_A$ instead of $m_B$. Which mass we choose play a role if $m_A\\ll m_B$ or\n", + "$m_B\\ll m_A$. One solution is to use the sum of the masses:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "t_c = \\sqrt{\\frac{L^3}{G(m_A+m_B)}}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Taking a look at [Kepler's laws](https://en.wikipedia.org/wiki/Kepler%27s_laws_of_planetary_motion) of planetary motion, the orbital period for a planet around the star is given by the $t_c$ above, except for a missing factor of $2\\pi$,\n", + "but that means that $t_c^{-1}$ is just the angular frequency of the motion.\n", + "Our characteristic time $t_c$ is therefore highly relevant.\n", + "Introducing the dimensionless number" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\alpha = \\frac{m_A}{m_B},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "we can write the dimensionless ODE as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\frac{d^2 \\bar\\x_A}{d\\bar t^2} =\n", + "\\frac{1}{1+\\alpha}\\frac{\\bar\\rpos}{||\\bar\\rpos||^3},\n", + "\\label{_auto6} \\tag{39}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\frac{d^2 \\bar\\x_B}{d\\bar t^2} =\n", + "\\frac{1}{1+\\alpha^{-1}}\\frac{\\bar\\rpos}{||\\bar\\rpos||^3}\\thinspace .\n", + "\\label{_auto7} \\tag{40}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the limit $m_A\\ll m_B$, i.e., $\\alpha\\ll 1$,\n", + "object B stands still, say $\\bar\\x_B=0$, and object\n", + "A orbits according to" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\frac{d^2 \\bar\\x_A}{d\\bar t^2} = -\\frac{\\bar\\x_A}{||\\bar \\x_A||^3}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Solution in a special case: planet orbiting a star\n", + "\n", + "To better see the motion, and that our scaling is reasonable,\n", + "we introduce polar coordinates $r$ and $\\theta$:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\bar\\x_A = r\\cos\\theta\\ii + r\\sin\\theta\\jj,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "which means $\\bar\\x_A$ can be written as $\\bar\\x_A =r\\ir$. Since" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\frac{d}{dt}\\ir = \\dot\\theta\\ith,\\quad \\frac{d}{dt}\\ith = -\\dot\\theta\\ir,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "we have" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\frac{d^2 \\bar\\x_A}{d\\bar t^2} =\n", + "(\\ddot r - r\\dot\\theta^2)\\ir + (r\\ddot\\theta + 2\\dot r\\dot\\theta)\\ith\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The equation of motion for mass A is then" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{align*}\n", + "\\ddot r - r\\dot\\theta^2 &= -\\frac{1}{r^2},\\\\ \n", + "r\\ddot\\theta + 2\\dot r\\dot\\theta &= 0\\thinspace .\n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The special case of circular motion, $r=1$, fulfills the equations, since\n", + "the latter equation then gives $\\dot\\theta =\\hbox{const}$ and\n", + "the former then gives $\\dot\\theta = 1$, i.e., the motion is\n", + "$r(t)=1$, $\\theta(t)=t$, with unit angular frequency as expected and\n", + "period $2\\pi$ as expected.\n", + "\n", + "\n", + "## Electric circuits\n", + "\n", + "Although the term \"mechanical vibrations\" is used in the present\n", + "book, we must mention that the same type of equations arise\n", + "when modeling electric circuits.\n", + "The current $I(t)$ in a\n", + "circuit with an inductor with inductance $L$, a capacitor with\n", + "capacitance $C$, and overall resistance $R$, is governed by" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\ddot I + \\frac{R}{L}\\dot I + \\frac{1}{LC}I = \\dot V(t),\n", + "\\label{_auto8} \\tag{41}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $V(t)$ is the voltage source powering the circuit.\n", + "This equation has the same form as the general model considered in\n", + "the section [vib:model2](#vib:model2) if we set $u=I$, $f(u^{\\prime})=bu^{\\prime}$\n", + "and define $b=R/L$, $s(u) = L^{-1}C^{-1}u$, and $F(t)=\\dot V(t)$.\n", + "\n", + "\n", + "# Exercises\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "## Exercise 1: Simulate resonance\n", + "
\n", + "\n", + "\n", + "We consider the scaled ODE model\n", + "([4](#vib:app:mass_gen:scaled)) from the section [General mechanical vibrating system](#vib:app:mass_gen).\n", + "After scaling, the amplitude of $u$ will have a size about unity\n", + "as time grows and the effect of the initial conditions die out due\n", + "to damping. However, as $\\gamma\\rightarrow 1$, the amplitude of $u$\n", + "increases, especially if $\\beta$ is small. This effect is called\n", + "*resonance*. The purpose of this exercise is to explore resonance.\n", + "\n", + "\n", + "**a)**\n", + "Figure out how the `solver` function in `vib.py` can be called\n", + "for the scaled ODE ([4](#vib:app:mass_gen:scaled)).\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "Comparing the scaled ODE ([4](#vib:app:mass_gen:scaled))\n", + "with the ODE ([3](#vib:app:mass_gen:equ)) with dimensions, we\n", + "realize that the parameters in the latter must be set as\n", + "\n", + " * $m=1$\n", + "\n", + " * $f(\\dot u) = 2\\beta |\\dot u|\\dot u$\n", + "\n", + " * $s(u)=ku$\n", + "\n", + " * $F(t)=\\sin(\\gamma t)$\n", + "\n", + " * $I=Ik/A$\n", + "\n", + " * $V=\\sqrt{mk}V/A$\n", + "\n", + "The expected period is $2\\pi$, so simulating for $N$ periods means\n", + "$T=2\\pi N$. Having $m$ time steps per period means $\\Delta t = 2\\pi/m$.\n", + "\n", + "Suppose we just choose $I=1$ and $V=0$. Simulating for 20 periods with\n", + "60 time steps per period, implies the following\n", + "`solver` call to run the scaled model:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "u, t = solver(I=1, V=0, m=1, b=2*beta, s=lambda u: u,\n", + " F=lambda t: sin(gamma*t), dt=2*pi/60,\n", + " T=2*pi*20, damping='quadratic')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "**b)**\n", + "Run $\\gamma =5, 1.5, 1.1, 1$ for $\\beta=0.005, 0.05, 0.2$.\n", + "For each $\\beta$ value, present an image with plots of $u(t)$ for\n", + "the four $\\gamma$ values.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "An appropriate program is" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "from vib import solver, visualize, plt\n", + "from math import pi, sin\n", + "import numpy as np\n", + "\n", + "beta_values = [0.005, 0.05, 0.2]\n", + "beta_values = [0.00005]\n", + "gamma_values = [5, 1.5, 1.1, 1]\n", + "for i, beta in enumerate(beta_values):\n", + " for gamma in gamma_values:\n", + " u, t = solver(I=1, V=0, m=1, b=2*beta, s=lambda u: u,\n", + " F=lambda t: sin(gamma*t), dt=2*pi/60,\n", + " T=2*pi*20, damping='quadratic')\n", + " visualize(u, t, title='gamma=%g' %\n", + " gamma, filename='tmp_%s' % gamma)\n", + " print gamma, 'max u amplitude:', np.abs(u).max()\n", + " for ext in 'png', 'pdf':\n", + " cmd = 'doconce combine_images '\n", + " cmd += ' '.join(['tmp_%s.' % gamma + ext\n", + " for gamma in gamma_values])\n", + " cmd += ' resonance%d.' % (i+1) + ext\n", + " os.system(cmd)\n", + "raw_input()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For $\\beta = 0.2$ we see that the amplitude is not far from unity:\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "For $\\beta =0.05$ we see that as $\\gamma\\rightarrow 1$, the amplitude grows:\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Finally, a small damping ($\\beta = 0.005$) amplifies the amplitude significantly (by a factor of 10) for $\\gamma=1$:\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "For a very small $\\beta=0.00005$, the amplitude grows linearly up to\n", + "about 60 for $\\bar t\\in [0,120]$.\n", + "\n", + "\n", + "\n", + "Filename: `resonance`.\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "## Exercise 2: Simulate oscillations of a sliding box\n", + "
\n", + "\n", + "Consider a sliding box on a flat surface as modeled in the section [A sliding mass attached to a spring](#vib:app:mass_sliding). As spring force we choose the nonlinear\n", + "formula" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "s(u) = \\frac{k}{\\alpha}\\tanh(\\alpha u) = ku + \\frac{1}{3}\\alpha^2 ku^3 + \\frac{2}{15}\\alpha^4 k u^5 + \\mathcal{O}({u^6})\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**a)**\n", + "Plot $g(u)=\\alpha^{-1}\\tanh(\\alpha u)$ for various values of $\\alpha$.\n", + "Assume $u\\in [-1,1]$.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "Here is a function that does the plotting:" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "import scitools.std as plt\n", + "import numpy as np\n", + "\n", + "def plot_spring():\n", + " alpha_values = [1, 2, 3, 10]\n", + " s = lambda u: 1.0/alpha*np.tanh(alpha*u)\n", + " u = np.linspace(-1, 1, 1001)\n", + " for alpha in alpha_values:\n", + " print alpha, s(u)\n", + " plt.plot(u, s(u))\n", + " plt.hold('on')\n", + " plt.legend([r'$\\alpha=%g$' % alpha for alpha in alpha_values])\n", + " plt.xlabel('u'); plt.ylabel('Spring response $s(u)$')\n", + " plt.savefig('tmp_s.png'); plt.savefig('tmp_s.pdf')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "**b)**\n", + "Scale the equations using $I$ as scale for $u$ and $\\sqrt{m/k}$ as\n", + "time scale.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "Inserting the dimensionless dependent and independent variables," + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\bar u = \\frac{u}{I},\\quad \\bar t = \\frac{t}{\\sqrt{m/k}},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "in the problem" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "m\\ddot u + \\mu mg\\hbox{sign}(\\dot u) + s(u) = 0,\\quad u(0)=I,\\ \\dot u(0)=V,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "gives" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\frac{d^2\\bar u}{d\\bar t^2} + \\frac{\\mu mg}{kI}\\hbox{sign}\\left(\n", + "\\frac{d\\bar u}{d\\bar t}\\right) + \\frac{1}{\\alpha I}\\tanh(\\alpha I\\bar u)\n", + "= 0,\\quad \\bar u(0)=1,\\ \\frac{d\\bar u}{d\\bar t}(0)=\\frac{V\\sqrt{mk}}{kI}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can now identify three dimensionless parameters," + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\beta = \\frac{\\mu mg}{kI},\\quad\n", + "\\gamma = \\alpha I,\\quad \\delta = \\frac{V\\sqrt{mk}}{kI}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The scaled problem can then be written" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\frac{d^2\\bar u}{d\\bar t^2} + \\beta\\hbox{sign}\\left(\n", + "\\frac{d\\bar u}{d\\bar t}\\right) + \\gamma^{-1}\\tanh(\\gamma \\bar u)\n", + "= 0,\\quad \\bar u(0)=1,\\ \\frac{d\\bar u}{d\\bar t}(0)=\\delta\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The initial set of 7 parameters $(\\mu, m, g, k, \\alpha, I, V)$ are\n", + "reduced to 3 dimensionless combinations.\n", + "\n", + "\n", + "\n", + "**c)**\n", + "Implement the scaled model in b). Run it for some values of\n", + "the dimensionless parameters.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "We use Odespy to solve the ODE, which requires rewriting the ODE as a\n", + "system of two first-order ODEs:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{align*}\n", + "v' &= - \\beta\\hbox{sign}(v) - \\gamma^{-1}\\tanh(\\gamma \\bar u),\\\\ \n", + "u' &= v,\n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "with initial conditions $v(0)=\\delta$ and $u(0)=1$. Here, $u(t)$ corresponds\n", + "to the previous $\\bar u(\\bar t)$, while $v(t)$ corresponds to\n", + "$d\\bar u/d\\bar t (\\bar t)$. The code can be like this:" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "! Should I change this to Devito?\n", + "\n", + "def simulate(beta, gamma, delta=0,\n", + " num_periods=8, time_steps_per_period=60):\n", + " # Use oscillations without friction to set dt and T\n", + " P = 2*np.pi\n", + " dt = P/time_steps_per_period\n", + " T = num_periods*P\n", + " t = np.linspace(0, T, time_steps_per_period*num_periods+1)\n", + " import odespy\n", + " def f(u, t, beta, gamma):\n", + " # Note the sequence of unknowns: v, u (v=du/dt)\n", + " v, u = u\n", + " return [-beta*np.sign(v) - 1.0/gamma*np.tanh(gamma*u), v]\n", + " #return [-beta*np.sign(v) - u, v]\n", + "\n", + " solver = odespy.RK4(f, f_args=(beta, gamma))\n", + " solver.set_initial_condition([delta, 1]) # sequence must match f\n", + " uv, t = solver.solve(t)\n", + " u = uv[:,1] # recall sequence in f: v, u\n", + " v = uv[:,0]\n", + " return u, t" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We simulate for an almost linear spring in the regime of $\\bar u$ (recall\n", + "that $\\bar u\\in [0,1]$ since $u$ is scaled with $I$), which corresponds\n", + "to $\\alpha = 1$ in a) and therefore $\\gamma =1$. Then we can try a\n", + "spring whose force quickly flattens out like $\\alpha=5$ in a), which\n", + "corresponds to $\\gamma = 5$ in the scaled model. A third option is\n", + "to have a truly linear spring, e.g., $\\gamma =0.1$. After some\n", + "experimentation we realize that $\\beta=0,0.05, 0.1$ are relevant values.\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Filename: `sliding_box`.\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "## Exercise 3: Simulate a bouncing ball\n", + "
\n", + "\n", + "The section [Bouncing ball](#vib:app:bouncing_ball) presents a model for a bouncing\n", + "ball.\n", + "Choose one of the two ODE formulation, ([34](#vib:app:bouncing:ball:h2eq)) or\n", + "([35](#vib:app:bouncing:ball:veq))-([36](#vib:app:bouncing:ball:heq)),\n", + "and simulate the motion of a bouncing ball. Plot $h(t)$. Think about how to\n", + "plot $v(t)$.\n", + "\n", + "\n", + "\n", + "**Hint.**\n", + "A naive implementation may get stuck in repeated impacts for large time\n", + "step sizes. To avoid this situation, one can introduce a state\n", + "variable that holds the mode of the motion: free fall, impact, or rest.\n", + "Two consecutive impacts imply that the motion has stopped.\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "A tailored `solver` function and some plotting statements go like" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "ename": "SyntaxError", + "evalue": "Missing parentheses in call to 'print'. Did you mean print(...)? (4188502967.py, line 42)", + "output_type": "error", + "traceback": [ + "\u001b[0;36m Cell \u001b[0;32mIn[6], line 42\u001b[0;36m\u001b[0m\n\u001b[0;31m print '%4d v=%8.5f h=%8.5f %s' % (n, v[n+1], h[n+1], mode)\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m Missing parentheses in call to 'print'. Did you mean print(...)?\n" + ] + } + ], + "source": [ + "! Devito-ize this\n", + "\n", + "import numpy as np\n", + "\n", + "def solver(H, C_R, dt, T, eps_v=0.01, eps_h=0.01):\n", + " \"\"\"\n", + " Simulate bouncing ball until it comes to rest. Time step dt.\n", + " h(0)=H (initial height). T: maximum simulation time.\n", + " Method: Euler-Cromer.\n", + " \"\"\"\n", + " dt = float(dt)\n", + " Nt = int(round(T/dt))\n", + " h = np.zeros(Nt+1)\n", + " v = np.zeros(Nt+1)\n", + " t = np.linspace(0, Nt*dt, Nt+1)\n", + " g = 0.81\n", + "\n", + " v[0] = 0\n", + " h[0] = H\n", + " mode = 'free fall'\n", + " for n in range(Nt):\n", + " v[n+1] = v[n] - dt*g\n", + " h[n+1] = h[n] + dt*v[n+1]\n", + "\n", + " if h[n+1] < eps_h:\n", + " #if abs(v[n+1]) > eps_v: # handles large dt, but is wrong\n", + " if v[n+1] < -eps_v:\n", + " # Impact\n", + " v[n+1] = -C_R*v[n+1]\n", + " h[n+1] = 0\n", + " if mode == 'impact':\n", + " # impact twice\n", + " return h[:n+2], v[:n+2], t[:n+2]\n", + " mode = 'impact'\n", + " elif abs(v[n+1]) < eps_v:\n", + " mode = 'rest'\n", + " v[n+1] = 0\n", + " h[n+1] = 0\n", + " return h[:n+2], v[:n+2], t[:n+2]\n", + " else:\n", + " mode = 'free fall'\n", + " else:\n", + " mode = 'free fall'\n", + " print '%4d v=%8.5f h=%8.5f %s' % (n, v[n+1], h[n+1], mode)\n", + " raise ValueError('T=%g is too short simulation time' % T)\n", + "\n", + "import matplotlib.pyplot as plt\n", + "h, v, t = solver(\n", + " H=1, C_R=0.8, T=100, dt=0.0001, eps_v=0.01, eps_h=0.01)\n", + "plt.plot(t, h)\n", + "plt.legend('h')\n", + "plt.savefig('tmp_h.png'); plt.savefig('tmp_h.pdf')\n", + "plt.figure()\n", + "plt.plot(t, v)\n", + "plt.legend('v')\n", + "plt.savefig('tmp_v.png'); plt.savefig('tmp_v.pdf')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Filename: `bouncing_ball`.\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "## Exercise 4: Simulate a simple pendulum\n", + "
\n", + "\n", + "Simulation of simple pendulum can be carried out by using\n", + "the mathematical model derived in the section [Motion of a pendulum](#vib:app:pendulum)\n", + "and calling up functionality in the [`vib.py`](${src_vib}/vib.py)\n", + "file (i.e., solve the second-order ODE by centered finite differences).\n", + "\n", + "\n", + "**a)**\n", + "Scale the model. Set up the dimensionless governing equation for $\\theta$\n", + "and expressions for dimensionless drag and wire forces.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "The angle is measured in radians so we may think of this quantity as\n", + "dimensionless, or we may scale it by the initial condition to obtain\n", + "a primary unknown that lies in $[-1,1]$. We go for the former strategy here.\n", + "\n", + "Dimensionless time $\\bar t$ is introduced as $t/t_c$ for some suitable\n", + "time scale $t_c$.\n", + "\n", + "Inserted in the two governing equations\n", + "([8](#vib:app:pendulum:thetaeq)) and ([6](#vib:app:pendulum:ir)),\n", + "for the\n", + "two unknowns $\\theta$ and $S$, respectively, we achieve" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{align*}\n", + "-S + mg\\cos\\theta &= -\\frac{1}{t_c}L\\frac{d\\theta}{d\\bar t},\\\\ \n", + "\\frac{1}{t_c^2}m\\frac{d^2\\theta}{d\\bar t^2} +\n", + "\\frac{1}{2}C_D\\varrho AL \\frac{1}{t_c^2}\\left\\vert\n", + "\\frac{d\\theta}{d\\bar t}\\right\\vert\n", + "\\frac{d\\theta}{d\\bar t}\n", + "+ \\frac{mg}{L}\\sin\\theta &= 0\\thinspace .\n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We multiply the latter equation by $t_c^2/m$ to make each term\n", + "dimensionless:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\frac{d^2\\theta}{d\\bar t^2} +\n", + "\\frac{1}{2m}C_D\\varrho AL \\left\\vert\n", + "\\frac{d\\theta}{d\\bar t}\\right\\vert\n", + "\\frac{d\\theta}{d\\bar t}\n", + "+ \\frac{t_c^2g}{L}\\sin\\theta = 0\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Assuming that the acceleration term and the\n", + "gravity term to be the dominating terms, these should balance, so\n", + "$t_c^2g/L=1$, giving $t_c = \\sqrt{g/L}$. With $A=\\pi R^2$ we get the\n", + "dimensionless ODEs" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\frac{d^2\\theta}{d\\bar t^2} +\n", + "\\alpha\\left\\vert\\frac{d\\theta}{d\\bar t}\\right\\vert\\frac{d\\theta}{d\\bar t} +\n", + "\\sin\\theta = 0,\n", + "\\label{vib:exer:pendulum_simple:eq:ith:s} \\tag{42}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\frac{S}{mg} = \\left(\\frac{d\\theta}{d\\bar t}\\right)^2 + \\cos\\theta,\n", + "\\label{vib:exer:pendulum_simple:eq:ir:s} \\tag{43}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $\\alpha$ is a dimensionless drag coefficient" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\alpha = \\frac{C_D\\varrho\\pi R^2L}{2m}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that in ([43](#vib:exer:pendulum_simple:eq:ir:s)) we have divided by\n", + "$mg$, which is in fact a force scale, making the gravity force unity\n", + "and also $S/mg=1$ in the equilibrium position $\\theta=0$. We may\n", + "introduce" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\bar S = S/mg\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "as a dimensionless drag force.\n", + "\n", + "The parameter $\\alpha$ is about\n", + "the ratio of the drag force and the gravity force:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\frac{|\\frac{1}{2} C_D\\varrho \\pi R^2 |v|v|}{|mg|}\\sim\n", + "\\frac{C_D\\varrho \\pi R^2 L^2 t_c^{-2}}{mg}\n", + "\\left|\\frac{d\\bar\\theta}{d\\bar t}\\right|\\frac{d\\bar\\theta}{d\\bar t}\n", + "\\sim \\frac{C_D\\varrho \\pi R^2 L}{2m}\\Theta^2 = \\alpha \\Theta^2\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "(We have that $\\theta(t)/d\\Theta$ is in $[-1,1]$, so we expect\n", + "since $\\Theta^{-1}d\\bar\\theta/d\\bar t$ to be around unity. Here,\n", + "$\\Theta=\\theta(0)$.)\n", + "\n", + "Let us introduce $\\omega$ for the dimensionless angular velocity," + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\omega = \\frac{d\\theta}{d\\bar t}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "When $\\theta$ is computed, the dimensionless wire and drag forces\n", + "are computed by" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{align*}\n", + "\\bar S &= \\omega^2 + \\cos\\theta,\\\\ \n", + "\\bar D &= -\\alpha |\\omega|\\omega\\thinspace .\n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "**b)**\n", + "Write a function for computing\n", + "$\\theta$ and the dimensionless drag force and the force in the wire,\n", + "using the `solver` function in\n", + "the `vib.py` file. Plot these three quantities\n", + "below each other (in subplots) so the graphs can be compared.\n", + "Run two cases, first one in the limit of $\\Theta$ small and\n", + "no drag, and then a second one with $\\Theta=40$ degrees and $\\alpha=0.8$.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "The first step is to realize how to utilize the `solver` function for\n", + "our dimensionless model. Introducing `Theta` for $\\Theta$, the\n", + "arguments to `solver` must be set as" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "I = Theta\n", + "V = 0\n", + "m = 1\n", + "b = alpha\n", + "s = lambda u: sin(u)\n", + "F = lambda t: 0\n", + "damping = 'quadratic'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After computing $\\theta$, we need to find $\\omega$ by finite differences:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\omega^n = \\frac{\\theta^{n+1}-\\theta^{n-1}}{2\\Delta t},\n", + "\\ n=1,\\ldots,N_t-1,\\quad \\omega^0=\\frac{\\theta^1-\\theta^0}{\\Delta t},\n", + "\\ \\omega^{N_t}=\\frac{\\theta^{N_t}-\\theta^{N_t-1}}{\\Delta t}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The duration of the simulation and the time step can be computed on\n", + "basis of the analytical insight we have for small $\\theta$\n", + "($\\theta\\approx \\Theta\\cos(t)$). A complete function then reads" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "def simulate(Theta, alpha, num_periods=10):\n", + " # Dimensionless model requires the following parameters:\n", + " from math import sin, pi\n", + "\n", + " I = Theta\n", + " V = 0\n", + " m = 1\n", + " b = alpha\n", + " s = lambda u: sin(u)\n", + " F = lambda t: 0\n", + " damping = 'quadratic'\n", + "\n", + " # Estimate T and dt from the small angle solution\n", + " P = 2*pi # One period (theta small, no drag)\n", + " dt = P/40 # 40 intervals per period\n", + " T = num_periods*P\n", + "\n", + " theta, t = solver(I, V, m, b, s, F, dt, T, damping)\n", + " omega = np.zeros(theta.size)\n", + " omega[1:-1] = (theta[2:] - theta[:-2])/(2*dt)\n", + " omega[0] = (theta[1] - theta[0])/dt\n", + " omega[-1] = (theta[-1] - theta[-2])/dt\n", + "\n", + " S = omega**2 + np.cos(theta)\n", + " D = alpha*np.abs(omega)*omega\n", + " return t, theta, S, D" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Assuming imports like" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "the following function visualizes $\\theta$, $\\bar S$, and $\\bar D$\n", + "with three subplots:" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "def visualize(t, theta, S, D, filename='tmp'):\n", + " f, (ax1, ax2, ax3) = plt.subplots(3, sharex=True, sharey=False)\n", + " ax1.plot(t, theta)\n", + " ax1.set_title(r'$\\theta(t)$')\n", + " ax2.plot(t, S)\n", + " ax2.set_title(r'Dimensonless force in the wire')\n", + " ax3.plot(t, D)\n", + " ax3.set_title(r'Dimensionless drag force')\n", + " plt.savefig('%s.png' % filename)\n", + " plt.savefig('%s.pdf' % filename)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A suitable main program is" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "import math\n", + "# Rough verification that small theta and no drag gives cos(t)\n", + "Theta = 1.0\n", + "alpha = 0\n", + "t, theta, S, D = simulate(Theta, alpha, num_periods=4)\n", + "# Scale theta by Theta (easier to compare with cos(t))\n", + "theta /= Theta\n", + "visualize(t, theta, S, D, filename='pendulum_verify')\n", + "\n", + "Theta = math.radians(40)\n", + "alpha = 0.8\n", + "t, theta, S, D = simulate(Theta, alpha)\n", + "visualize(t, theta, S, D, filename='pendulum_alpha0.8_Theta40')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The \"verification\" case looks good (at least when the `solver` function\n", + "has been thoroughly verified in other circumstances):\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "The \"real case\" shows how quickly the drag force is reduced, even when\n", + "we set $\\alpha$ to a significant value (0.8):\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Filename: `simple_pendulum`.\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "## Exercise 5: Simulate an elastic pendulum\n", + "
\n", + "\n", + "The section [Motion of an elastic pendulum](#vib:app:pendulum_elastic) describes a model for an elastic\n", + "pendulum, resulting in a system of two ODEs. The purpose of this\n", + "exercise is to implement the scaled model, test the software, and\n", + "generalize the model.\n", + "\n", + "\n", + "**a)**\n", + "Write a function `simulate`\n", + "that can simulate an elastic pendulum using the scaled model.\n", + "The function should have the following arguments:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "def simulate(\n", + " beta=0.9, # dimensionless parameter\n", + " Theta=30, # initial angle in degrees\n", + " epsilon=0, # initial stretch of wire\n", + " num_periods=6, # simulate for num_periods\n", + " time_steps_per_period=60, # time step resolution\n", + " plot=True, # make plots or not\n", + " ):" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To set the total simulation time and the time step, we\n", + "use our knowledge of the scaled, classical, non-elastic pendulum:\n", + "$u^{\\prime\\prime} + u = 0$, with solution\n", + "$u = \\Theta\\cos \\bar t$.\n", + "The period of these oscillations is $P=2\\pi$\n", + "and the frequency is unity. The time\n", + "for simulation is taken as `num_periods` times $P$. The time step\n", + "is set as $P$ divided by `time_steps_per_period`.\n", + "\n", + "The `simulate` function should return the arrays of\n", + "$x$, $y$, $\\theta$, and $t$, where $\\theta = \\tan^{-1}(x/(1-y))$ is\n", + "the angular displacement of the elastic pendulum corresponding to the\n", + "position $(x,y)$.\n", + "\n", + "If `plot` is `True`, make a plot of $\\bar y(\\bar t)$\n", + "versus $\\bar x(\\bar t)$, i.e., the physical motion\n", + "of the mass at $(\\bar x,\\bar y)$. Use the equal aspect ratio on the axis\n", + "such that we get a physically correct picture of the motion. Also\n", + "make a plot of $\\theta(\\bar t)$, where $\\theta$ is measured in degrees.\n", + "If $\\Theta < 10$ degrees, add a plot that compares the solutions of\n", + "the scaled, classical, non-elastic pendulum and the elastic pendulum\n", + "($\\theta(t)$).\n", + "\n", + "Although the mathematics here employs a bar over scaled quantities, the\n", + "code should feature plain names `x` for $\\bar x$, `y` for $\\bar y$, and\n", + "`t` for $\\bar t$ (rather than `x_bar`, etc.). These variable names make\n", + "the code easier to read and compare with the mathematics.\n", + "\n", + "\n", + "\n", + "**Hint 1.**\n", + "Equal aspect ratio is set by `plt.gca().set_aspect('equal')` in\n", + "Matplotlib (`import matplotlib.pyplot as plt`)\n", + "and in SciTools by the command\n", + "`plt.plot(..., daspect=[1,1,1], daspectmode='equal')`\n", + "(provided you have done `import scitools.std as plt`).\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "**Hint 2.**\n", + "If you want to use Odespy to solve the equations, order the ODEs\n", + "like $\\dot \\bar x, \\bar x, \\dot\\bar y,\\bar y$ such that\n", + "`odespy.EulerCromer` can be applied.\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "Here is a suggested `simulate` function:" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [], + "source": [ + "! Devito-ize the odespy part?\n", + "\n", + "import odespy\n", + "import numpy as np\n", + "import scitools.std as plt\n", + "\n", + "def simulate(\n", + " beta=0.9, # dimensionless parameter\n", + " Theta=30, # initial angle in degrees\n", + " epsilon=0, # initial stretch of wire\n", + " num_periods=6, # simulate for num_periods\n", + " time_steps_per_period=60, # time step resolution\n", + " plot=True, # make plots or not\n", + " ):\n", + " from math import sin, cos, pi\n", + " Theta = Theta*np.pi/180 # convert to radians\n", + " # Initial position and velocity\n", + " # (we order the equations such that Euler-Cromer in odespy\n", + " # can be used, i.e., vx, x, vy, y)\n", + " ic = [0, # x'=vx\n", + " (1 + epsilon)*sin(Theta), # x\n", + " 0, # y'=vy\n", + " 1 - (1 + epsilon)*cos(Theta), # y\n", + " ]\n", + "\n", + " def f(u, t, beta):\n", + " vx, x, vy, y = u\n", + " L = np.sqrt(x**2 + (y-1)**2)\n", + " h = beta/(1-beta)*(1 - beta/L) # help factor\n", + " return [-h*x, vx, -h*(y-1) - beta, vy]\n", + "\n", + " # Non-elastic pendulum (scaled similarly in the limit beta=1)\n", + " # solution Theta*cos(t)\n", + " P = 2*pi\n", + " dt = P/time_steps_per_period\n", + " T = num_periods*P\n", + " omega = 2*pi/P\n", + "\n", + " time_points = np.linspace(\n", + " 0, T, num_periods*time_steps_per_period+1)\n", + "\n", + " solver = odespy.EulerCromer(f, f_args=(beta,))\n", + " solver.set_initial_condition(ic)\n", + " u, t = solver.solve(time_points)\n", + " x = u[:,1]\n", + " y = u[:,3]\n", + " theta = np.arctan(x/(1-y))\n", + "\n", + " if plot:\n", + " plt.figure()\n", + " plt.plot(x, y, 'b-', title='Pendulum motion',\n", + " daspect=[1,1,1], daspectmode='equal',\n", + " axis=[x.min(), x.max(), 1.3*y.min(), 1])\n", + " plt.savefig('tmp_xy.png')\n", + " plt.savefig('tmp_xy.pdf')\n", + " # Plot theta in degrees\n", + " plt.figure()\n", + " plt.plot(t, theta*180/np.pi, 'b-',\n", + " title='Angular displacement in degrees')\n", + " plt.savefig('tmp_theta.png')\n", + " plt.savefig('tmp_theta.pdf')\n", + " if abs(Theta) < 10*pi/180:\n", + " # Compare theta and theta_e for small angles (<10 degrees)\n", + " theta_e = Theta*np.cos(omega*t) # non-elastic scaled sol.\n", + " plt.figure()\n", + " plt.plot(t, theta, t, theta_e,\n", + " legend=['theta elastic', 'theta non-elastic'],\n", + " title='Elastic vs non-elastic pendulum, '\\\n", + " 'beta=%g' % beta)\n", + " plt.savefig('tmp_compare.png')\n", + " plt.savefig('tmp_compare.pdf')\n", + " # Plot y vs x (the real physical motion)\n", + " return x, y, theta, t" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "**b)**\n", + "Write a test function for testing that $\\Theta=0$ and $\\epsilon=0$\n", + "gives $x=y=0$ for all times.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "Here is the code:" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "def test_equilibrium():\n", + " \"\"\"Test that starting from rest makes x=y=theta=0.\"\"\"\n", + " x, y, theta, t = simulate(\n", + " beta=0.9, Theta=0, epsilon=0,\n", + " num_periods=6, time_steps_per_period=10, plot=False)\n", + " tol = 1E-14\n", + " assert np.abs(x.max()) < tol\n", + " assert np.abs(y.max()) < tol\n", + " assert np.abs(theta.max()) < tol" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "**c)**\n", + "Write another test function for checking that the pure vertical\n", + "motion of the elastic pendulum is correct.\n", + "Start with simplifying the ODEs for pure vertical motion and show that\n", + "$\\bar y(\\bar t)$ fulfills a vibration equation with\n", + "frequency $\\sqrt{\\beta/(1-\\beta)}$. Set up the exact solution.\n", + "\n", + "Write a test function that\n", + "uses this special case to verify the `simulate` function. There will\n", + "be numerical approximation errors present in the results from\n", + "`simulate` so you have to believe in correct results and set a\n", + "(low) tolerance that corresponds to the computed maximum error.\n", + "Use a small $\\Delta t$ to obtain a small numerical approximation error.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "For purely vertical motion, the ODEs reduce to $\\ddot x = 0$ and" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\frac{d^2\\bar y}{d\\bar t^2} = -\\frac{\\beta}{1-\\beta}(1-\\beta\\frac{1}{\\sqrt{(\\bar y - 1)^2}})(\\bar y-1) - \\beta = -\\frac{\\beta}{1-\\beta}(\\bar y-1 + \\beta) - \\beta\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have here used that $(\\bar y -1)/\\sqrt{(\\bar y -1)^2}=-1$ since\n", + "$\\bar y$ cannot exceed 1 (the pendulum's wire is fixed at the scaled\n", + "point $(0,1)$). In fact, $\\bar y$ will be around zero.\n", + "(As a consistency check, we realize that in equilibrium, $\\ddot\\bar y =0$,\n", + "and multiplying by $(1-\\beta)/\\beta$ leads to the expected $\\bar y=0$.)\n", + "Further calculations easily lead to" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\frac{d^2\\bar y}{d\\bar t^2} = -\\frac{\\beta}{1-\\beta}\\bar y = -\\omega^2\\bar y,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where we have introduced the frequency\n", + "$\\omega = \\sqrt{\\beta/(1-\\beta)}$.\n", + "Solving this standard ODE, with an initial stretching $\\bar y(0)=\\epsilon$\n", + "and no velocity, results in" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\bar y(\\bar t) = \\epsilon\\cos(\\omega\\bar t)\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that the oscillations we describe here are very different from\n", + "the oscillations used to set the period and time step in function\n", + "`simulate`. The latter type of oscillations are due to gravity when\n", + "a classical, non-elastic pendulum oscillates back and forth, while\n", + "$\\bar y(\\bar t)$ above refers to vertical *elastic* oscillations in the wire\n", + "around the equilibrium point in the gravity field. The angular frequency\n", + "of the vertical oscillations are given by $\\omega$ and the corresponding\n", + "period is $\\hat P = 2\\pi/\\omega$. Suppose we want to simulate for\n", + "$T=N\\hat P = N2\\pi/\\omega$ and use $n$ time steps per period,\n", + "$\\Delta\\bar t = \\hat P/n$. The `simulate` function operates with\n", + "a simulation time of `num_periods` times $2\\pi$. This means that we must set\n", + "`num_periods=N/omega` if we want to simulate to time $T=N\\hat P$.\n", + "The parameter `time_steps_per_period` must be set to $\\omega n$\n", + "since `simulate` has $\\Delta t$ as $2\\pi$ divided by `time_steps_per_period`\n", + "and we want $\\Delta t = 2\\pi\\omega^{-1}n^{-1}$.\n", + "\n", + "The corresponding test function can be written as follows." + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "def test_vertical_motion():\n", + " beta = 0.9\n", + " omega = np.sqrt(beta/(1-beta))\n", + " # Find num_periods. Recall that P=2*pi for scaled pendulum\n", + " # oscillations, while here we don't have gravity driven\n", + " # oscillations, but elastic oscillations with frequency omega.\n", + " period = 2*np.pi/omega\n", + " # We want T = N*period\n", + " N = 5\n", + " # simulate function has T = 2*pi*num_periods\n", + " num_periods = 5/omega\n", + " n = 600\n", + " time_steps_per_period = omega*n\n", + "\n", + " y_exact = lambda t: -0.1*np.cos(omega*t)\n", + " x, y, theta, t = simulate(\n", + " beta=beta, Theta=0, epsilon=0.1,\n", + " num_periods=num_periods,\n", + " time_steps_per_period=time_steps_per_period,\n", + " plot=False)\n", + "\n", + " tol = 0.00055 # ok tolerance for the above resolution\n", + " # No motion in x direction is epxected\n", + " assert np.abs(x.max()) < tol\n", + " # Check motion in y direction\n", + " y_e = y_exact(t)\n", + " diff = np.abs(y_e - y).max()\n", + " if diff > tol: # plot\n", + " plt.plot(t, y, t, y_e, legend=['y', 'exact'])\n", + " raw_input('Error in test_vertical_motion; type CR:')\n", + " assert diff < tol, 'diff=%g' % diff" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "**d)**\n", + "Make a function `demo(beta, Theta)` for simulating an elastic pendulum with a\n", + "given $\\beta$ parameter and initial angle $\\Theta$. Use 600 time steps\n", + "per period to get every accurate results, and simulate for 3 periods.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "The `demo` function is just" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "def demo(beta=0.999, Theta=40, num_periods=3):\n", + " x, y, theta, t = simulate(\n", + " beta=beta, Theta=Theta, epsilon=0,\n", + " num_periods=num_periods, time_steps_per_period=600,\n", + " plot=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Below are plots corresponding to $\\beta = 0.999$ (3 periods) and\n", + "$\\beta = 0.93$ (one period):\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Filename: `elastic_pendulum`.\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "## Exercise 6: Simulate an elastic pendulum with air resistance\n", + "
\n", + "\n", + "This is a continuation [Exercise 5: Simulate an elastic pendulum](#vib:exer:pendulum_elastic).\n", + "Air resistance on the body with mass $m$ can be modeled by the\n", + "force $-\\frac{1}{2}\\varrho C_D A|\\v|\\v$,\n", + "where $C_D$ is a drag coefficient (0.2 for a sphere), $\\varrho$\n", + "is the density of air (1.2 $\\hbox{kg }\\,{\\hbox{m}}^{-3}$), $A$ is the\n", + "cross section area ($A=\\pi R^2$ for a sphere, where $R$ is the radius),\n", + "and $\\v$ is the velocity of the body.\n", + "Include air resistance in the original model, scale the model,\n", + "write a function `simulate_drag` that is a copy of the `simulate`\n", + "function from [Exercise 5: Simulate an elastic pendulum](#vib:exer:pendulum_elastic), but with the\n", + "new ODEs included, and show plots of how air resistance\n", + "influences the motion.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "We start with the model\n", + "([18](#vib:app:pendulum_elastic:x))-([24](#vib:app:pendulum_elastic:vy0)).\n", + "Since $\\v = \\dot x\\ii + \\dot y\\jj$, the air resistance term\n", + "can be written" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "-q(\\dot x\\ii + \\dot y\\jj),\\quad q=\\frac{1}{2}\\varrho C_D A\\sqrt{\\dot x^2 + \\dot y^2}\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that for positive velocities, the pendulum is moving to the right\n", + "and the air resistance works against the motion, i.e., in direction of\n", + "$-\\v = -\\dot x\\ii - \\dot y\\jj$.\n", + "\n", + "We can easily include the terms in the ODEs:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\ddot x = -\\frac{q}{m}\\dot x -\\frac{k}{m}\\left(1 -\\frac{L_0}{L}\\right)(x-x_0),\n", + "\\label{vib:app:pendulum_elastic_drag:x} \\tag{44}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\ddot y = -\\frac{q}{m}\\dot y -\\frac{k}{m}\\left(1 -\\frac{L_0}{L}\\right)(y-y_0) - g,\n", + "\\label{vib:app:pendulum_elastic_drag:y} \\tag{45}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "L = \\sqrt{(x-x_0)^2 + (y-y_0)^2},\n", + "\\label{vib:app:pendulum_elastic_drag:L} \\tag{46}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\label{_auto9} \\tag{47}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The initial conditions are not affected.\n", + "\n", + "The next step is to scale the model. We use the same scales as in\n", + "[Exercise 5: Simulate an elastic pendulum](#vib:exer:pendulum_elastic), introduce $\\beta$, and $A=\\pi R^2$\n", + "to simplify the $-q\\dot x/m$ term to" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\frac{L_0}{2m}\\varrho C_D R^2\\beta^{-1}\n", + "\\sqrt{\\left(\\frac{d\\bar x}{d\\bar t}\\right)^2 +\n", + "\\left(\\frac{d\\bar y}{d\\bar t}\\right)^2}\n", + "= \\gamma \\beta^{-1}\n", + "\\sqrt{\\left(\\frac{d\\bar x}{d\\bar t}\\right)^2 +\n", + "\\left(\\frac{d\\bar y}{d\\bar t}\\right)^2},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $\\gamma$ is a second dimensionless parameter:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\gamma = \\frac{L_0}{2m}\\varrho C_D R^2\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The final set of scaled equations is then" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\frac{d^2\\bar x}{d\\bar t^2} = -\\gamma\\beta^{-1}\n", + "\\sqrt{\\left(\\frac{d\\bar x}{d\\bar t}\\right)^2 +\n", + "\\left(\\frac{d\\bar y}{d\\bar t}\\right)^2}\\frac{d\\bar x}{d\\bar t}\n", + "-\\frac{\\beta}{1-\\beta}\\left(1- \\frac{\\beta}{\\bar L}\\right)\\bar x,\n", + "\\label{vib:app:pendulum_elastic_drag:x:s} \\tag{48}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\frac{d^2\\bar y}{d\\bar t^2} =\n", + "-\\gamma\\beta^{-1}\n", + "\\sqrt{\\left(\\frac{d\\bar x}{d\\bar t}\\right)^2 +\n", + "\\left(\\frac{d\\bar y}{d\\bar t}\\right)^2}\\frac{d\\bar y}{d\\bar t}\n", + "-\\frac{\\beta}{1-\\beta}\\left(1- \\frac{\\beta}{\\bar L}\\right)(\\bar y-1)\n", + "-\\beta,\n", + "\\label{vib:app:pendulum_elastic_drag:y:s} \\tag{49}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\bar L = \\sqrt{\\bar x^2 + (\\bar y-1)^2},\n", + "\\label{vib:app:pendulum_elastic_drag:L:s} \\tag{50}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\bar x(0) = (1+\\epsilon)\\sin\\Theta,\n", + "\\label{vib:app:pendulum_elastic_drag:x0:s} \\tag{51}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\frac{d\\bar x}{d\\bar t}(0) = 0,\n", + "\\label{vib:app:pendulum_elastic_drag:vx0:s} \\tag{52}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\bar y(0) = 1 - (1+\\epsilon)\\cos\\Theta,\n", + "\\label{vib:app:pendulum_elastic_drag:y0:s} \\tag{53}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\frac{d\\bar y}{d\\bar t}(0) = 0,\n", + "\\label{vib:app:pendulum_elastic_drag:vy0:s} \\tag{54}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The new `simulate_drag` function is implemented below." + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "! Devito-ize odespy\n", + "\n", + "def simulate_drag(\n", + " beta=0.9, # dimensionless elasticity parameter\n", + " gamma=0, # dimensionless drag parameter\n", + " Theta=30, # initial angle in degrees\n", + " epsilon=0, # initial stretch of wire\n", + " num_periods=6, # simulate for num_periods\n", + " time_steps_per_period=60, # time step resolution\n", + " plot=True, # make plots or not\n", + " ):\n", + " from math import sin, cos, pi\n", + " Theta = Theta*np.pi/180 # convert to radians\n", + " # Initial position and velocity\n", + " # (we order the equations such that Euler-Cromer in odespy\n", + " # can be used, i.e., vx, x, vy, y)\n", + " ic = [0, # x'=vx\n", + " (1 + epsilon)*sin(Theta), # x\n", + " 0, # y'=vy\n", + " 1 - (1 + epsilon)*cos(Theta), # y\n", + " ]\n", + "\n", + " def f(u, t, beta, gamma):\n", + " vx, x, vy, y = u\n", + " L = np.sqrt(x**2 + (y-1)**2)\n", + " v = np.sqrt(vx**2 + vy**2)\n", + " h1 = beta/(1-beta)*(1 - beta/L) # help factor\n", + " h2 = gamma/beta*v\n", + " return [-h2*vx - h1*x, vx, -h2*vy - h1*(y-1) - beta, vy]\n", + "\n", + " # Non-elastic pendulum (scaled similarly in the limit beta=1)\n", + " # solution Theta*cos(t)\n", + " P = 2*pi\n", + " dt = P/time_steps_per_period\n", + " T = num_periods*P\n", + " omega = 2*pi/P\n", + "\n", + " time_points = np.linspace(\n", + " 0, T, num_periods*time_steps_per_period+1)\n", + "\n", + " solver = odespy.EulerCromer(f, f_args=(beta, gamma))\n", + " solver.set_initial_condition(ic)\n", + " u, t = solver.solve(time_points)\n", + " x = u[:,1]\n", + " y = u[:,3]\n", + " theta = np.arctan(x/(1-y))\n", + "\n", + " if plot:\n", + " plt.figure()\n", + " plt.plot(x, y, 'b-', title='Pendulum motion',\n", + " daspect=[1,1,1], daspectmode='equal',\n", + " axis=[x.min(), x.max(), 1.3*y.min(), 1])\n", + " plt.savefig('tmp_xy.png')\n", + " plt.savefig('tmp_xy.pdf')\n", + " # Plot theta in degrees\n", + " plt.figure()\n", + " plt.plot(t, theta*180/np.pi, 'b-',\n", + " title='Angular displacement in degrees')\n", + " plt.savefig('tmp_theta.png')\n", + " plt.savefig('tmp_theta.pdf')\n", + " if abs(Theta) < 10*pi/180:\n", + " # Compare theta and theta_e for small angles (<10 degrees)\n", + " theta_e = Theta*np.cos(omega*t) # non-elastic scaled sol.\n", + " plt.figure()\n", + " plt.plot(t, theta, t, theta_e,\n", + " legend=['theta elastic', 'theta non-elastic'],\n", + " title='Elastic vs non-elastic pendulum, '\\\n", + " 'beta=%g' % beta)\n", + " plt.savefig('tmp_compare.png')\n", + " plt.savefig('tmp_compare.pdf')\n", + " # Plot y vs x (the real physical motion)\n", + " return x, y, theta, t" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The plot of $\\theta$ shows the damping ($\\beta = 0.999$):\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Test functions for equilibrium and vertical motion are also included. These\n", + "are as in [Exercise 6: Simulate an elastic pendulum with air resistance](#vib:exer:pendulum_elastic_drag), except that\n", + "they call `simulate_drag` instead of `simulate`.\n", + "\n", + "\n", + "Filename: `elastic_pendulum_drag`.\n", + "\n", + "\n", + "\n", + "### Remarks\n", + "\n", + "Test functions are challenging to construct for the problem with\n", + "air resistance. You can reuse the tests from\n", + "[Exercise 6: Simulate an elastic pendulum with air resistance](#vib:exer:pendulum_elastic_drag) for `simulate_drag`,\n", + "but these tests does not verify the new terms arising from air\n", + "resistance.\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "## Exercise 7: Implement the PEFRL algorithm\n", + "
\n", + "\n", + "We consider the motion of a planet around a star (the section [Two-body gravitational problem](#vib:app:gravitation)).\n", + "The simplified case where one\n", + "mass is very much bigger than the other and one object is at rest,\n", + "results in the scaled ODE model" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{align*}\n", + "\\ddot x + (x^2 + y^2)^{-3/2}x & = 0,\\\\ \n", + "\\ddot y + (x^2 + y^2)^{-3/2}y & = 0\\thinspace .\n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**a)**\n", + "It is easy to show that $x(t)$ and $y(t)$ go like sine and cosine\n", + "functions. Use this idea to derive the exact solution.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "We may assume $x=C_x\\cos(\\omega t)$ and $y=C_y\\sin(\\omega t)$ for\n", + "constants $C_x,$, $C_y$, and $\\omega$. Inserted in the equations, we\n", + "see that $\\omega =1$. The initial conditions determine the other\n", + "constants, which we may choose as $C_x=C_y=1$ (the object starts\n", + "at $(1,0)$ with a velocity $(0,1)$). The motion is a perfect circle,\n", + "which should last forever.\n", + "\n", + "\n", + "\n", + "**b)**\n", + "One believes that a planet may orbit a star for billions of years.\n", + "We are now interested\n", + "in how accurate methods we actually need for such calculations.\n", + "A first task is to determine what the time interval of interest is in\n", + "scaled units. Take the earth and sun as typical objects and find\n", + "the characteristic time used in the scaling of the equations\n", + "($t_c = \\sqrt{L^3/(mG)}$), where $m$ is the mass of the sun, $L$ is the\n", + "distance between the sun and the earth, and $G$ is the gravitational\n", + "constant. Find the scaled time interval corresponding to one billion years.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "According to [Wikipedia](https://en.wikipedia.org/wiki/Solar_mass),\n", + "the mass of the sun is approximately $2\\cdot 10^{30}$ kg. This\n", + "is 332946 times the mass of the earth, implying that the\n", + "dimensionless constant $\\alpha \\approx 3\\cdot 10^{-6}$. With\n", + "$G=6.674\\cdot 10^{-11}\\hbox{ Nm}^2/\\hbox{kg}^2$, and the\n", + "[sun-earth distance](https://en.wikipedia.org/wiki/Astronomical_unit)\n", + "as (approximately) 150 million km, we have $t_c \\approx 5 028 388$ s.\n", + "This is about 58 days, which is the characteristic time, chosen as the\n", + "angular frequency of the oscillations. To get the period of one orbit we therefore must multiply by $2\\pi$. This gives about 1 year (and demonstrates the\n", + "fact mentioned about the scaling: the natural time scale is consistent with\n", + "Kepler's law about the period).\n", + "\n", + "Thus, one billion years correspond to 62,715,924,070 time units (dividing\n", + "one billion years by $t_c$), which corresponds to about 2000\n", + "\"time unit years\".\n", + "\n", + "\n", + "\n", + "**c)**\n", + "Solve the equations using 4th-order Runge-Kutta and the Euler-Cromer\n", + "methods. You may benefit from applying Odespy for this purpose. With\n", + "each solver, simulate 10000 orbits and print the maximum position\n", + "error and CPU time as a function of time step. Note that the maximum\n", + "position error does not necessarily occur at the end of the\n", + "simulation. The position error achieved with each solver will depend\n", + "heavily on the size of the time step. Let the time step correspond to\n", + "200, 400, 800 and 1600 steps per orbit, respectively. Are the results\n", + "as expected? Explain briefly. When you develop your program, have in\n", + "mind that it will be extended with an implementation of the other\n", + "algorithms (as requested in d) and e) later) and experiments with this\n", + "algorithm as well.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "The first task is to implement the right-hand side function for the\n", + "system of ODEs such that we can call up Odespy solvers (or make use of\n", + "other types of ODE software, e.g., from SciPy). The $2\\times 2$ system of\n", + "second-order ODEs must be expressed as a $4\\times 4$ system of first-order\n", + "ODEs. We have three different cases of right-hand sides:\n", + "\n", + "1. Common numbering of unknowns: $x$, $v_x$, $y$, $y_x$\n", + "\n", + "2. Numbering required by Euler-Cromer: $v_x$, $x$, $v_y$, $y$\n", + "\n", + "3. Numbering required by the PEFRL method: same as Euler-Cromer\n", + "\n", + "Most Odespy solvers can handle any convention for numbering of the unknowns.\n", + "The important point is that initial conditions and new values at the end of\n", + "the time step are filled in the right positions of a one-dimensional array\n", + "containing the unknowns.\n", + "Using Odespy to solve the system by the Euler-Cromer method, however, requires\n", + "the unknowns to appear as velocity 1st degree-of-freedom, displacement\n", + "1st degree-of-freedom, velocity 2nd degree-of-freedom, displacement\n", + "2nd degree-of-freedom, and so forth. Two alternative right-hand side\n", + "functions `f(u, t)` for Odespy solvers is then" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "def f_EC(u, t):\n", + " '''\n", + " Return derivatives for the 1st order system as\n", + " required by Euler-Cromer.\n", + " '''\n", + " vx, x, vy, y = u # u: array holding vx, x, vy, y\n", + " d = -(x**2 + y**2)**(-3.0/2)\n", + " return [d*x, vx, d*y, vy ]\n", + "\n", + "def f_RK4(u, t):\n", + " '''\n", + " Return derivatives for the 1st order system as\n", + " required by ordinary solvers in Odespy.\n", + " '''\n", + " x, vx, y, vy = u # u: array holding x, vx, y, vy\n", + " d = -(x**2 + y**2)**(-3.0/2)\n", + " return [vx, d*x, vy, d*y ]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In addition, we shall later in d) implement the PEFRL method and just\n", + "give the $g$ function as input to a system of the form $dv_x = g_x$,\n", + "$dv_y = g_y$, and $g$ becomes the vector $(g_x,g_y)$:\n", + "\n", + "Some prefer to number the unknowns differently, and with the RK4 method we\n", + "are free to use any numbering, including this one:" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "def g(u, v):\n", + " return np.array([-u])\n", + "def u_exact(t):\n", + " return np.array([3*np.cos(t)]).transpose()\n", + "I = u_exact(0)\n", + "V = np.array([0])\n", + "print 'V:', V, 'I:', I\n", + "\n", + "# Numerical parameters\n", + "w = 1\n", + "P = 2*np.pi/w\n", + "dt_values = [P/20, P/40, P/80, P/160, P/320]\n", + "T = 8*P\n", + "error_vs_dt = []\n", + "for n, dt in enumerate(dt_values):\n", + " u, v, t = solver_PEFRL(I, V, g, dt, T)\n", + " error = np.abs(u - u_exact(t)).max()\n", + " print 'error:', error\n", + " if n > 0:\n", + " error_vs_dt.append(error/dt**4)\n", + "for i in range(1, len(error_vs_dt)):\n", + " #print abs(error_vs_dt[i]- error_vs_dt[0])\n", + " assert abs(error_vs_dt[i]-\n", + " error_vs_dt[0]) < 0.1\n", + "\n", + "\n", + "s PEFRL(odespy.Solver):\n", + "\"\"\"Class wrapper for Odespy.\"\"\" # Not used!\n", + "quick_desctiption = \"Explicit 4th-order method for v'=-f, u=v.\"\"\"\n", + "\n", + "def advance(self):\n", + " u, f, n, t = self.u, self.f, self.n, self.t\n", + " dt = t[n+1] - t[n]\n", + " I = np.array([u[1], u[3]])\n", + " V = np.array([u[0], u[2]])\n", + " u, v, t = solver_PFFRL(I, V, f, dt, t+dt)\n", + " return np.array([v[-1], u[-1]])\n", + "\n", + "compute_orbit_and_error(\n", + "f,\n", + "solver_ID,\n", + "timesteps_per_period=20,\n", + "N_orbit_groups=1000,\n", + "orbit_group_size=10):\n", + "'''\n", + "For one particular solver:\n", + "Calculte the orbits for a multiple of grouped orbits, i.e.\n", + "number of orbits = orbit_group_size*N_orbit_groups.\n", + "Returns: time step dt, and, for each N_orbit_groups cycle,\n", + "the 2D position error and cpu time (as lists).\n", + "'''\n", + "def u_exact(t):\n", + " return np.array([np.cos(t), np.sin(t)])\n", + "\n", + "w = 1\n", + "P = 2*np.pi/w # scaled period (1 year becomes 2*pi)\n", + "dt = P/timesteps_per_period\n", + "Nt = orbit_group_size*N_orbit_groups*timesteps_per_period\n", + "T = Nt*dt\n", + "t_mesh = np.linspace(0, T, Nt+1)\n", + "E_orbit = []\n", + "\n", + "#print ' dt:', dt\n", + "T_interval = P*orbit_group_size\n", + "N = int(round(T_interval/dt))\n", + "\n", + "# set initial conditions\n", + "if solver_ID == 'EC':\n", + " A = [0,1,1,0]\n", + "elif solver_ID == 'PEFRL':\n", + " I = np.array([1, 0])\n", + " V = np.array([0, 1])\n", + "else:\n", + " A = [1,0,0,1]\n", + "\n", + "t1 = time.clock()\n", + "for i in range(N_orbit_groups):\n", + " time_points = np.linspace(i*T_interval, (i+1)*T_interval,N+1)\n", + " u_e = u_exact(time_points).transpose()\n", + " if solver_ID == 'EC':\n", + " solver = odespy.EulerCromer(f)\n", + " solver.set_initial_condition(A)\n", + " ui, ti = solver.solve(time_points)\n", + " # Find error (correct final pos: x=1, y=0)\n", + " orbit_error = np.sqrt(\n", + " (ui[:,1]-u_e[:,0])**2 + (ui[:,3]-u_e[:,1])**2).max()\n", + " elif solver_ID == 'PEFRL':\n", + " # Note: every T_inverval is here counted from time 0\n", + " ui, vi, ti = solver_PEFRL(I, V, f, dt, T_interval)\n", + " # Find error (correct final pos: x=1, y=0)\n", + " orbit_error = np.sqrt(\n", + " (ui[:,0]-u_e[:,0])**2 + (ui[:,1]-u_e[:,1])**2).max()\n", + " else:\n", + " solver = eval('odespy.' + solver_ID(f)\n", + " solver.set_initial_condition(A)\n", + " ui, ti = solver.solve(time_points)\n", + " # Find error (correct final pos: x=1, y=0)\n", + " orbit_error = np.sqrt(\n", + " (ui[:,0]-u_e[:,0])**2 + (ui[:,2]-u_e[:,1])**2).max()\n", + "\n", + " print ' Orbit no. %d, max error (per cent): %g' % \\\n", + " ((i+1)*orbit_group_size, orbit_error)\n", + "\n", + " E_orbit.append(orbit_error)\n", + "\n", + " # set init. cond. for next time interval\n", + " if solver_ID == 'EC':\n", + " A = [ui[-1,0], ui[-1,1], ui[-1,2], ui[-1,3]]\n", + " elif solver_ID == 'PEFRL':\n", + " I = [ui[-1,0], ui[-1,1]]\n", + " V = [vi[-1,0], vi[-1,1]]\n", + " else: # RK4, adaptive rules, etc.\n", + " A = [ui[-1,0], ui[-1,1], ui[-1,2], ui[-1,3]]\n", + "\n", + "t2 = time.clock()\n", + "CPU_time = (t2 - t1)/(60.0*60.0) # in hours\n", + "return dt, E_orbit, CPU_time\n", + "\n", + "orbit_error_vs_dt(\n", + "f_EC, f_RK4, g, solvers,\n", + "N_orbit_groups=1000,\n", + "orbit_group_size=10):\n", + "'''\n", + "With each solver in list \"solvers\": Simulate\n", + "orbit_group_size*N_orbit_groups orbits with different dt values.\n", + "Collect final 2D position error for each dt and plot all errors.\n", + "'''\n", + "\n", + "for solver_ID in solvers:\n", + " print 'Computing orbit with solver:', solver_ID\n", + " E_values = []\n", + " dt_values = []\n", + " cpu_values = []\n", + " for timesteps_per_period in 200, 400, 800, 1600:\n", + " print '.......time steps per period: ', \\\n", + " timesteps_per_period\n", + " if solver_ID == 'EC':\n", + " dt, E, cpu_time = compute_orbit_and_error(\n", + " f_EC,\n", + " solver_ID,\n", + " timesteps_per_period,\n", + " N_orbit_groups,\n", + " orbit_group_size)\n", + " elif solver_ID == 'PEFRL':\n", + " dt, E, cpu_time = compute_orbit_and_error(\n", + " g,\n", + " solver_ID,\n", + " timesteps_per_period,\n", + " N_orbit_groups,\n", + " orbit_group_size)\n", + " else:\n", + " dt, E, cpu_time = compute_orbit_and_error(\n", + " f_RK4,\n", + " solver_ID,\n", + " timesteps_per_period,\n", + " N_orbit_groups,\n", + " orbit_group_size)\n", + "\n", + " dt_values.append(dt)\n", + " E_values.append(np.array(E).max())\n", + " cpu_values.append(cpu_time)\n", + " print 'dt_values:', dt_values\n", + " print 'E max with dt...:', E_values\n", + " print 'cpu_values with dt...:', cpu_values\n", + "\n", + "\n", + "orbit_error_vs_years(\n", + "f_EC, f_RK4, g, solvers,\n", + "N_orbit_groups=1000,\n", + "orbit_group_size=100,\n", + "N_time_steps = 1000):\n", + "'''\n", + "For each solver in the list solvers:\n", + "simulate orbit_group_size*N_orbit_groups orbits with a fixed\n", + "dt corresponding to N_time_steps steps per year.\n", + "Collect max 2D position errors for each N_time_steps'th run,\n", + "plot these errors and CPU. Finally, make an empirical\n", + "formula for error and CPU as functions of a number\n", + "of cycles.\n", + "'''\n", + "timesteps_per_period = N_time_steps # fixed for all runs\n", + "\n", + "for solver_ID in solvers:\n", + " print 'Computing orbit with solver:', solver_ID\n", + " if solver_ID == 'EC':\n", + " dt, E, cpu_time = compute_orbit_and_error(\n", + " f_EC,\n", + " solver_ID,\n", + " timesteps_per_period,\n", + " N_orbit_groups,\n", + " orbit_group_size)\n", + " elif solver_ID == 'PEFRL':\n", + " dt, E, cpu_time = compute_orbit_and_error(\n", + " g,\n", + " solver_ID,\n", + " timesteps_per_period,\n", + " N_orbit_groups,\n", + " orbit_group_size)\n", + " else:\n", + " dt, E, cpu_time = compute_orbit_and_error(\n", + " f_RK4,\n", + " solver_ID,\n", + " timesteps_per_period,\n", + " N_orbit_groups,\n", + " orbit_group_size)\n", + "\n", + " # E and cpu_time are for every N_orbit_groups cycle\n", + " print 'E_values (fixed dt, changing no of years):', E\n", + " print 'CPU (hours):', cpu_time\n", + " years = np.arange(\n", + " 0,\n", + " N_orbit_groups*orbit_group_size,\n", + " orbit_group_size)\n", + "\n", + " # Now make empirical formula\n", + "\n", + " def E_of_years(x, *coeff):\n", + " return sum(coeff[i]*x**float((len(coeff)-1)-i) \\\n", + " for i in range(len(coeff)))\n", + " E = np.array(E)\n", + " degree = 4\n", + " # note index: polyfit finds p[0]*x**4 + p[1]*x**3 ...etc.\n", + " p = np.polyfit(years, E, degree)\n", + " p_str = map(str, p)\n", + " formula = ' + '.join([p_str[i] + '*x**' + \\\n", + " str(degree-i) for i in range(degree+1)])\n", + "\n", + " print 'Empirical formula (error with years): ', formula\n", + " plt.figure()\n", + " plt.plot(years,\n", + " E, 'b-',\n", + " years,\n", + " E_of_years(years, *p), 'r--')\n", + " plt.xlabel('Number of years')\n", + " plt.ylabel('Orbit error')\n", + " plt.title(solver_ID)\n", + " filename = solver_ID + 'tmp_E_with_years'\n", + " plt.savefig(filename + '.png')\n", + " plt.savefig(filename + '.pdf')\n", + " plt.show()\n", + "\n", + " print 'Predicted CPU time in hours (1 billion years):', \\\n", + " cpu_time*10000\n", + " print 'Predicted max error (1 billion years):', \\\n", + " E_of_years(1E9, *p)\n", + "\n", + "compute_orbit_error_and_CPU():\n", + "'''\n", + "Orbit error and associated CPU times are computed with\n", + "solvers: RK4, Euler-Cromer, PEFRL.'''\n", + "\n", + "def f_EC(u, t):\n", + " '''\n", + " Return derivatives for the 1st order system as\n", + " required by Euler-Cromer.\n", + " '''\n", + " vx, x, vy, y = u # u: array holding vx, x, vy, y\n", + " d = -(x**2 + y**2)**(-3.0/2)\n", + " return [d*x, vx, d*y, vy ]\n", + "\n", + "def f_RK4(u, t):\n", + " '''\n", + " Return derivatives for the 1st order system as\n", + " required by ordinary solvers in Odespy.\n", + " '''\n", + " x, vx, y, vy = u # u: array holding x, vx, y, vy\n", + " d = -(x**2 + y**2)**(-3.0/2)\n", + " return [vx, d*x, vy, d*y ]\n", + "\n", + "def g(u, v):\n", + " '''\n", + " Return derivatives for the 1st order system as\n", + " required by PEFRL.\n", + " '''\n", + " d = -(u[0]**2 + u[1]**2)**(-3.0/2)\n", + " return np.array([d*u[0], d*u[1]])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The standard way of solving the ODE by Odespy is then" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "def u_exact(t):\n", + " \"\"\"Return exact solution at time t.\"\"\"\n", + " return np.array([np.cos(t), np.sin(t)])\n", + "\n", + "u_e = u_exact(time_points).transpose()\n", + "\n", + "solver = odespy.RK4(f_RK4)\n", + "solver.set_initial_condition(A)\n", + "ui, ti = solver.solve(time_points)\n", + "\n", + "# Find error (correct final pos: x=1, y=0)\n", + "orbit_error = np.sqrt(\n", + " (ui[:,0]-u_e[:,0])**2 + (ui[:,2]-u_e[:,1])**2).max()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We develop functions for computing errors and plotting results where we\n", + "can compare different methods. These functions are shown in the solution to\n", + "item d).\n", + "\n", + "Running the code, the time step sizes become" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " dt_values: [0.031415926535897934, 0.015707963267948967,\n", + " 0.007853981633974483, 0.003926990816987242]\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Corresponding maximum errors (per cent) and CPU values (hours) are for the 4th-order Runge-Kutta given in the table below.\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
Quantity $\\Delta t_1$ $\\Delta t_2$ $\\Delta t_3$ $\\Delta t_4$
$\\Delta t$ 0.03 0.02 0.008 0.004
Error 1.9039 0.0787 0.0025 7.7e-05
CPU (h) 0.03 0.06 0.12 0.23
\n", + "For Euler-Cromer we these results:\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
Quantity $\\Delta t_1$ $\\Delta t_2$ $\\Delta t_3$ $\\Delta t_4$
$\\Delta t$ 0.03 0.02 0.008 0.004
Error 2.0162 2.0078 1.9634 0.6730
CPU (h) 0.01 0.02 0.05 0.09
\n", + "\n", + "These results are as expected. The Runge-Kutta implementation is much more accurate than Euler-Cromer, but since it requires more computations, more CPU time is needed. For both methods, accuracy and CPU time both increase as\n", + "the step size is reduced, but the increase is much more pronounced for\n", + "the Euler-Cromer method.\n", + "\n", + "\n", + "\n", + "**d)**\n", + "Implement a solver based on the PEFRL method from\n", + "the section [vib:ode2:PEFRL](#vib:ode2:PEFRL). Verify its 4th-order convergence\n", + "using an equation $u'' + u = 0$.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "Here is a solver function:" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import time\n", + "\n", + "def solver_PEFRL(I, V, g, dt, T):\n", + " \"\"\"\n", + " Solve v' = - g(u,v), u'=v for t in (0,T], u(0)=I and v(0)=V,\n", + " by the PEFRL method.\n", + " \"\"\"\n", + " dt = float(dt)\n", + " Nt = int(round(T/dt))\n", + " u = np.zeros((Nt+1, len(I)))\n", + " v = np.zeros((Nt+1, len(I)))\n", + " t = np.linspace(0, Nt*dt, Nt+1)\n", + "\n", + " # these values are from eq (20), ref to paper below\n", + " xi = 0.1786178958448091\n", + " lambda_ = -0.2123418310626054\n", + " chi = -0.06626458266981849\n", + "\n", + " v[0] = V\n", + " u[0] = I\n", + " # Compare with eq 22 in http://arxiv.org/pdf/cond-mat/0110585.pdf\n", + " for n in range(0, Nt):\n", + " u_ = u[n] + xi*dt*v[n]\n", + " v_ = v[n] + 0.5*(1-2*lambda_)*dt*g(u_, v[n])\n", + " u_ = u_ + chi*dt*v_\n", + " v_ = v_ + lambda_*dt*g(u_, v_)\n", + " u_ = u_ + (1-2*(chi+xi))*dt*v_\n", + " v_ = v_ + lambda_*dt*g(u_, v_)\n", + " u_ = u_ + chi*dt*v_\n", + " v[n+1] = v_ + 0.5*(1-2*lambda_)*dt*g(u_, v_)\n", + " u[n+1] = u_ + xi*dt*v[n+1]\n", + " #print 'v[%d]=%g, u[%d]=%g' % (n+1,v[n+1],n+1,u[n+1])\n", + " return u, v, t" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A proper test function for verification reads" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [], + "source": [ + "def test_solver_PEFRL():\n", + " \"\"\"Check 4th order convergence rate, using u'' + u = 0,\n", + " I = 3.0, V = 0, which has the exact solution u_e = 3*cos(t)\"\"\"\n", + " def g(u, v):\n", + " return np.array([-u])\n", + " def u_exact(t):\n", + " return np.array([3*np.cos(t)]).transpose()\n", + " I = u_exact(0)\n", + " V = np.array([0])\n", + " print 'V:', V, 'I:', I\n", + "\n", + " # Numerical parameters\n", + " w = 1\n", + " P = 2*np.pi/w\n", + " dt_values = [P/20, P/40, P/80, P/160, P/320]\n", + " T = 8*P\n", + " error_vs_dt = []\n", + " for n, dt in enumerate(dt_values):\n", + " u, v, t = solver_PEFRL(I, V, g, dt, T)\n", + " error = np.abs(u - u_exact(t)).max()\n", + " print 'error:', error\n", + " if n > 0:\n", + " error_vs_dt.append(error/dt**4)\n", + " for i in range(1, len(error_vs_dt)):\n", + " #print abs(error_vs_dt[i]- error_vs_dt[0])\n", + " assert abs(error_vs_dt[i]-\n", + " error_vs_dt[0]) < 0.1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "**e)**\n", + "The simulations done previously with the 4th-order Runge-Kutta and\n", + "Euler-Cromer are now to be repeated with the PEFRL solver, so the\n", + "code must be extended accordingly. Then run the simulations and comment\n", + "on the performance of PEFRL compared to the other two.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "With the PEFRL algorithm, we get" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " E max with dt...: [0.0010452575786173163, 6.5310955829464402e-05,\n", + " 4.0475768394248492e-06, 2.9391302503251016e-07]\n", + " cpu_values with dt...: [0.01873611111111106, 0.037422222222222294,\n", + " 0.07511666666666655, 0.14985]\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
Qantity $\\Delta t_1$ $\\Delta t_2$ $\\Delta t_3$ $\\Delta t_4$
$\\Delta t$ 0.03 0.02 0.008 0.004
Error 1.04E-3 6.53E-05 4.05E-6 2.94E-7
CPU (h) 0.02 0.04 0.08 0.15
\n", + "\n", + "The accuracy is now dramatically improved compared to 4th-order Runge-Kutta (and Euler-Cromer).\n", + "With 1600 steps per orbit, the PEFRL maximum error is just below $3.0e-07$ per cent, while\n", + "the corresponding error with Runge-Kutta was about $7.7e-05$ per cent! This is striking,\n", + "considering the fact that the 4th-order Runge-Kutta and the PEFRL schemes are both 4th-order accurate.\n", + "\n", + "\n", + "\n", + "**f)**\n", + "Use the PEFRL solver to simulate 100000 orbits with a fixed time step\n", + "corresponding to 1600 steps per period. Record the maximum error\n", + "within each subsequent group of 1000 orbits. Plot these errors and fit\n", + "(least squares) a mathematical function to the data. Print also the\n", + "total CPU time spent for all 100000 orbits.\n", + "\n", + "Now, predict the error and required CPU time for a simulation of 1\n", + "billion years (orbits). Is it feasible on today's computers to\n", + "simulate the planetary motion for one billion years?\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "The complete code (which also produces the printouts given previously) reads:" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "import scitools.std as plt\n", + "import sys\n", + "import odespy\n", + "import numpy as np\n", + "import time\n", + "\n", + "def solver_PEFRL(I, V, g, dt, T):\n", + " \"\"\"\n", + " Solve v' = - g(u,v), u'=v for t in (0,T], u(0)=I and v(0)=V,\n", + " by the PEFRL method.\n", + " \"\"\"\n", + " dt = float(dt)\n", + " Nt = int(round(T/dt))\n", + " u = np.zeros((Nt+1, len(I)))\n", + " v = np.zeros((Nt+1, len(I)))\n", + " t = np.linspace(0, Nt*dt, Nt+1)\n", + "\n", + " # these values are from eq (20), ref to paper below\n", + " xi = 0.1786178958448091\n", + " lambda_ = -0.2123418310626054\n", + " chi = -0.06626458266981849\n", + "\n", + " v[0] = V\n", + " u[0] = I\n", + " # Compare with eq 22 in http://arxiv.org/pdf/cond-mat/0110585.pdf\n", + " for n in range(0, Nt):\n", + " u_ = u[n] + xi*dt*v[n]\n", + " v_ = v[n] + 0.5*(1-2*lambda_)*dt*g(u_, v[n])\n", + " u_ = u_ + chi*dt*v_\n", + " v_ = v_ + lambda_*dt*g(u_, v_)\n", + " u_ = u_ + (1-2*(chi+xi))*dt*v_\n", + " v_ = v_ + lambda_*dt*g(u_, v_)\n", + " u_ = u_ + chi*dt*v_\n", + " v[n+1] = v_ + 0.5*(1-2*lambda_)*dt*g(u_, v_)\n", + " u[n+1] = u_ + xi*dt*v[n+1]\n", + " #print 'v[%d]=%g, u[%d]=%g' % (n+1,v[n+1],n+1,u[n+1])\n", + " return u, v, t\n", + "\n", + "def test_solver_PEFRL():\n", + " \"\"\"Check 4th order convergence rate, using u'' + u = 0,\n", + " I = 3.0, V = 0, which has the exact solution u_e = 3*cos(t)\"\"\"\n", + " def g(u, v):\n", + " return np.array([-u])\n", + " def u_exact(t):\n", + " return np.array([3*np.cos(t)]).transpose()\n", + " I = u_exact(0)\n", + " V = np.array([0])\n", + " print 'V:', V, 'I:', I\n", + "\n", + " # Numerical parameters\n", + " w = 1\n", + " P = 2*np.pi/w\n", + " dt_values = [P/20, P/40, P/80, P/160, P/320]\n", + " T = 8*P\n", + " error_vs_dt = []\n", + " for n, dt in enumerate(dt_values):\n", + " u, v, t = solver_PEFRL(I, V, g, dt, T)\n", + " error = np.abs(u - u_exact(t)).max()\n", + " print 'error:', error\n", + " if n > 0:\n", + " error_vs_dt.append(error/dt**4)\n", + " for i in range(1, len(error_vs_dt)):\n", + " #print abs(error_vs_dt[i]- error_vs_dt[0])\n", + " assert abs(error_vs_dt[i]-\n", + " error_vs_dt[0]) < 0.1\n", + "\n", + "\n", + "class PEFRL(odespy.Solver):\n", + " \"\"\"Class wrapper for Odespy.\"\"\" # Not used!\n", + " quick_desctiption = \"Explicit 4th-order method for v'=-f, u=v.\"\"\"\n", + "\n", + " def advance(self):\n", + " u, f, n, t = self.u, self.f, self.n, self.t\n", + " dt = t[n+1] - t[n]\n", + " I = np.array([u[1], u[3]])\n", + " V = np.array([u[0], u[2]])\n", + " u, v, t = solver_PFFRL(I, V, f, dt, t+dt)\n", + " return np.array([v[-1], u[-1]])\n", + "\n", + "def compute_orbit_and_error(\n", + " f,\n", + " solver_ID,\n", + " timesteps_per_period=20,\n", + " N_orbit_groups=1000,\n", + " orbit_group_size=10):\n", + " '''\n", + " For one particular solver:\n", + " Calculte the orbits for a multiple of grouped orbits, i.e.\n", + " number of orbits = orbit_group_size*N_orbit_groups.\n", + " Returns: time step dt, and, for each N_orbit_groups cycle,\n", + " the 2D position error and cpu time (as lists).\n", + " '''\n", + " def u_exact(t):\n", + " return np.array([np.cos(t), np.sin(t)])\n", + "\n", + " w = 1\n", + " P = 2*np.pi/w # scaled period (1 year becomes 2*pi)\n", + " dt = P/timesteps_per_period\n", + " Nt = orbit_group_size*N_orbit_groups*timesteps_per_period\n", + " T = Nt*dt\n", + " t_mesh = np.linspace(0, T, Nt+1)\n", + " E_orbit = []\n", + "\n", + " #print ' dt:', dt\n", + " T_interval = P*orbit_group_size\n", + " N = int(round(T_interval/dt))\n", + "\n", + " # set initial conditions\n", + " if solver_ID == 'EC':\n", + " A = [0,1,1,0]\n", + " elif solver_ID == 'PEFRL':\n", + " I = np.array([1, 0])\n", + " V = np.array([0, 1])\n", + " else:\n", + " A = [1,0,0,1]\n", + "\n", + " t1 = time.clock()\n", + " for i in range(N_orbit_groups):\n", + " time_points = np.linspace(i*T_interval, (i+1)*T_interval,N+1)\n", + " u_e = u_exact(time_points).transpose()\n", + " if solver_ID == 'EC':\n", + " solver = odespy.EulerCromer(f)\n", + " solver.set_initial_condition(A)\n", + " ui, ti = solver.solve(time_points)\n", + " # Find error (correct final pos: x=1, y=0)\n", + " orbit_error = np.sqrt(\n", + " (ui[:,1]-u_e[:,0])**2 + (ui[:,3]-u_e[:,1])**2).max()\n", + " elif solver_ID == 'PEFRL':\n", + " # Note: every T_inverval is here counted from time 0\n", + " ui, vi, ti = solver_PEFRL(I, V, f, dt, T_interval)\n", + " # Find error (correct final pos: x=1, y=0)\n", + " orbit_error = np.sqrt(\n", + " (ui[:,0]-u_e[:,0])**2 + (ui[:,1]-u_e[:,1])**2).max()\n", + " else:\n", + " solver = eval('odespy.' + solver_ID(f)\n", + " solver.set_initial_condition(A)\n", + " ui, ti = solver.solve(time_points)\n", + " # Find error (correct final pos: x=1, y=0)\n", + " orbit_error = np.sqrt(\n", + " (ui[:,0]-u_e[:,0])**2 + (ui[:,2]-u_e[:,1])**2).max()\n", + "\n", + " print ' Orbit no. %d, max error (per cent): %g' % \\\n", + " ((i+1)*orbit_group_size, orbit_error)\n", + "\n", + " E_orbit.append(orbit_error)\n", + "\n", + " # set init. cond. for next time interval\n", + " if solver_ID == 'EC':\n", + " A = [ui[-1,0], ui[-1,1], ui[-1,2], ui[-1,3]]\n", + " elif solver_ID == 'PEFRL':\n", + " I = [ui[-1,0], ui[-1,1]]\n", + " V = [vi[-1,0], vi[-1,1]]\n", + " else: # RK4, adaptive rules, etc.\n", + " A = [ui[-1,0], ui[-1,1], ui[-1,2], ui[-1,3]]\n", + "\n", + " t2 = time.clock()\n", + " CPU_time = (t2 - t1)/(60.0*60.0) # in hours\n", + " return dt, E_orbit, CPU_time\n", + "\n", + "def orbit_error_vs_dt(\n", + " f_EC, f_RK4, g, solvers,\n", + " N_orbit_groups=1000,\n", + " orbit_group_size=10):\n", + " '''\n", + " With each solver in list \"solvers\": Simulate\n", + " orbit_group_size*N_orbit_groups orbits with different dt values.\n", + " Collect final 2D position error for each dt and plot all errors.\n", + " '''\n", + "\n", + " for solver_ID in solvers:\n", + " print 'Computing orbit with solver:', solver_ID\n", + " E_values = []\n", + " dt_values = []\n", + " cpu_values = []\n", + " for timesteps_per_period in 200, 400, 800, 1600:\n", + " print '.......time steps per period: ', \\\n", + " timesteps_per_period\n", + " if solver_ID == 'EC':\n", + " dt, E, cpu_time = compute_orbit_and_error(\n", + " f_EC,\n", + " solver_ID,\n", + " timesteps_per_period,\n", + " N_orbit_groups,\n", + " orbit_group_size)\n", + " elif solver_ID == 'PEFRL':\n", + " dt, E, cpu_time = compute_orbit_and_error(\n", + " g,\n", + " solver_ID,\n", + " timesteps_per_period,\n", + " N_orbit_groups,\n", + " orbit_group_size)\n", + " else:\n", + " dt, E, cpu_time = compute_orbit_and_error(\n", + " f_RK4,\n", + " solver_ID,\n", + " timesteps_per_period,\n", + " N_orbit_groups,\n", + " orbit_group_size)\n", + "\n", + " dt_values.append(dt)\n", + " E_values.append(np.array(E).max())\n", + " cpu_values.append(cpu_time)\n", + " print 'dt_values:', dt_values\n", + " print 'E max with dt...:', E_values\n", + " print 'cpu_values with dt...:', cpu_values\n", + "\n", + "\n", + "def orbit_error_vs_years(\n", + " f_EC, f_RK4, g, solvers,\n", + " N_orbit_groups=1000,\n", + " orbit_group_size=100,\n", + " N_time_steps = 1000):\n", + " '''\n", + " For each solver in the list solvers:\n", + " simulate orbit_group_size*N_orbit_groups orbits with a fixed\n", + " dt corresponding to N_time_steps steps per year.\n", + " Collect max 2D position errors for each N_time_steps'th run,\n", + " plot these errors and CPU. Finally, make an empirical\n", + " formula for error and CPU as functions of a number\n", + " of cycles.\n", + " '''\n", + " timesteps_per_period = N_time_steps # fixed for all runs\n", + "\n", + " for solver_ID in solvers:\n", + " print 'Computing orbit with solver:', solver_ID\n", + " if solver_ID == 'EC':\n", + " dt, E, cpu_time = compute_orbit_and_error(\n", + " f_EC,\n", + " solver_ID,\n", + " timesteps_per_period,\n", + " N_orbit_groups,\n", + " orbit_group_size)\n", + " elif solver_ID == 'PEFRL':\n", + " dt, E, cpu_time = compute_orbit_and_error(\n", + " g,\n", + " solver_ID,\n", + " timesteps_per_period,\n", + " N_orbit_groups,\n", + " orbit_group_size)\n", + " else:\n", + " dt, E, cpu_time = compute_orbit_and_error(\n", + " f_RK4,\n", + " solver_ID,\n", + " timesteps_per_period,\n", + " N_orbit_groups,\n", + " orbit_group_size)\n", + "\n", + " # E and cpu_time are for every N_orbit_groups cycle\n", + " print 'E_values (fixed dt, changing no of years):', E\n", + " print 'CPU (hours):', cpu_time\n", + " years = np.arange(\n", + " 0,\n", + " N_orbit_groups*orbit_group_size,\n", + " orbit_group_size)\n", + "\n", + " # Now make empirical formula\n", + "\n", + " def E_of_years(x, *coeff):\n", + " return sum(coeff[i]*x**float((len(coeff)-1)-i) \\\n", + " for i in range(len(coeff)))\n", + " E = np.array(E)\n", + " degree = 4\n", + " # note index: polyfit finds p[0]*x**4 + p[1]*x**3 ...etc.\n", + " p = np.polyfit(years, E, degree)\n", + " p_str = map(str, p)\n", + " formula = ' + '.join([p_str[i] + '*x**' + \\\n", + " str(degree-i) for i in range(degree+1)])\n", + "\n", + " print 'Empirical formula (error with years): ', formula\n", + " plt.figure()\n", + " plt.plot(years,\n", + " E, 'b-',\n", + " years,\n", + " E_of_years(years, *p), 'r--')\n", + " plt.xlabel('Number of years')\n", + " plt.ylabel('Orbit error')\n", + " plt.title(solver_ID)\n", + " filename = solver_ID + 'tmp_E_with_years'\n", + " plt.savefig(filename + '.png')\n", + " plt.savefig(filename + '.pdf')\n", + " plt.show()\n", + "\n", + " print 'Predicted CPU time in hours (1 billion years):', \\\n", + " cpu_time*10000\n", + " print 'Predicted max error (1 billion years):', \\\n", + " E_of_years(1E9, *p)\n", + "\n", + "def compute_orbit_error_and_CPU():\n", + " '''\n", + " Orbit error and associated CPU times are computed with\n", + " solvers: RK4, Euler-Cromer, PEFRL.'''\n", + "\n", + " def f_EC(u, t):\n", + " '''\n", + " Return derivatives for the 1st order system as\n", + " required by Euler-Cromer.\n", + " '''\n", + " vx, x, vy, y = u # u: array holding vx, x, vy, y\n", + " d = -(x**2 + y**2)**(-3.0/2)\n", + " return [d*x, vx, d*y, vy ]\n", + "\n", + " def f_RK4(u, t):\n", + " '''\n", + " Return derivatives for the 1st order system as\n", + " required by ordinary solvers in Odespy.\n", + " '''\n", + " x, vx, y, vy = u # u: array holding x, vx, y, vy\n", + " d = -(x**2 + y**2)**(-3.0/2)\n", + " return [vx, d*x, vy, d*y ]\n", + "\n", + " def g(u, v):\n", + " '''\n", + " Return derivatives for the 1st order system as\n", + " required by PEFRL.\n", + " '''\n", + " d = -(u[0]**2 + u[1]**2)**(-3.0/2)\n", + " return np.array([d*u[0], d*u[1]])\n", + "\n", + " print 'Find orbit error as fu. of dt...(10000 orbits)'\n", + " solvers = ['RK4', 'EC', 'PEFRL']\n", + " N_orbit_groups=1\n", + " orbit_group_size=10000\n", + " orbit_error_vs_dt(\n", + " f_EC, f_RK4, g, solvers,\n", + " N_orbit_groups=N_orbit_groups,\n", + " orbit_group_size=orbit_group_size)\n", + "\n", + " print 'Compute orbit error as fu. of no of years (fixed dt)...'\n", + " solvers = ['PEFRL']\n", + " N_orbit_groups=100\n", + " orbit_group_size=1000\n", + " N_time_steps = 1600 # no of steps per orbit cycle\n", + " orbit_error_vs_years(\n", + " f_EC, f_RK4, g, solvers,\n", + " N_orbit_groups=N_orbit_groups,\n", + " orbit_group_size=orbit_group_size,\n", + " N_time_steps = N_time_steps)\n", + "\n", + "if __name__ == '__main__':\n", + " test_solver_PEFRL()\n", + " compute_orbit_error_and_CPU()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The maximum error develops with number of orbits as seen in the following plot,\n", + "where the red dashed curve is from the mathematical model:\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "We note that the maximum error achieved during the first 100000 orbits is only\n", + "about $1.2e-06$ per cent. Not bad!\n", + "\n", + "For the printed CPU and empirical formula, we get:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " CPU (hours): 1.51591388889\n", + " Empirical formula (E with years):\n", + " 3.15992325978e-26*x**4 + -6.1772567063e-21*x**3 +\n", + " 1.87983349496e-16*x**2 + 2.32924158693e-11*x**1 +\n", + " 5.46989368301e-08*x**0\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since the CPU develops linearly, the CPU time for 100000 orbits can just be multiplied by 10000 to get the\n", + "estimated CPU time required for 1 billion years. This gives 15159 CPU hours (631 days), which is also printed.\n", + "\n", + "With the derived empirical formula, the estimated orbit error after 1 billion years becomes 31593055529 per cent.\n", + "\n", + "[sl 2: Can we really use the plot and the function to predict max E during 1 billion years? Seems hard.]\n", + "\n", + "\n", + "\n", + "Filename: `vib_PEFRL`.\n", + "\n", + "\n", + "\n", + "### Remarks\n", + "\n", + "This exercise investigates whether it is feasible to predict\n", + "planetary motion for the life time of a solar system.\n", + "[hpl 3: Is it???]\n", + "\n", + "\n", + "" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.10.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/fdm-jupyter-book/notebooks/01_vib/vib_gen.ipynb b/fdm-jupyter-book/notebooks/01_vib/vib_gen.ipynb new file mode 100644 index 00000000..539cefb9 --- /dev/null +++ b/fdm-jupyter-book/notebooks/01_vib/vib_gen.ipynb @@ -0,0 +1,2250 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We shall now generalize the simple model problem from Section 1.1 to include a possibly nonlinear damping term $f(u')$, a possibly nonlinear spring (or restoring) force $s(u)$, and some external excitation $F(t)$:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "mu^{\\prime\\prime} + f(u^{\\prime}) + s(u) = F(t),\\quad u(0)=I,\\ u^{\\prime}(0)=V,\\ t\\in (0,T]\n", + "\\thinspace .\n", + "\\label{vib:ode2} \\tag{1}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have also included a possibly nonzero initial value for $u^{\\prime}(0)$.\n", + "The parameters $m$, $f(u^{\\prime})$, $s(u)$, $F(t)$, $I$, $V$, and $T$ are\n", + "input data.\n", + "\n", + "There are two main types of damping (friction) forces: linear $f(u^{\\prime})=bu$, or\n", + "quadratic $f(u^{\\prime})=bu^{\\prime}|u^{\\prime}|$. Spring systems often feature linear\n", + "damping, while air resistance usually gives rise to quadratic damping.\n", + "Spring forces are often linear: $s(u)=cu$, but nonlinear versions\n", + "are also common, the most famous is the gravity force on a pendulum\n", + "that acts as a spring with $s(u)\\sim \\sin(u)$.\n", + "\n", + "\n", + "## A centered scheme for linear damping\n", + "
\n", + "\n", + "Sampling ([1](#vib:ode2)) at a mesh point $t_n$, replacing\n", + "$u^{\\prime\\prime}(t_n)$ by $[D_tD_tu]^n$, and $u^{\\prime}(t_n)$ by $[D_{2t}u]^n$ results\n", + "in the discretization" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "[mD_tD_t u + f(D_{2t}u) + s(u) = F]^n,\n", + "\\label{_auto1} \\tag{2}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "which written out means" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "m\\frac{u^{n+1}-2u^n + u^{n-1}}{\\Delta t^2}\n", + "+ f(\\frac{u^{n+1}-u^{n-1}}{2\\Delta t}) + s(u^n) = F^n,\n", + "\\label{vib:ode2:step3b} \\tag{3}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $F^n$ as usual means $F(t)$ evaluated at $t=t_n$.\n", + "Solving ([3](#vib:ode2:step3b)) with respect to the unknown\n", + "$u^{n+1}$ gives a problem: the $u^{n+1}$ inside the $f$ function\n", + "makes the equation *nonlinear* unless $f(u^{\\prime})$ is a linear function,\n", + "$f(u^{\\prime})=bu^{\\prime}$. For now we shall assume that $f$ is linear in $u^{\\prime}$.\n", + "Then" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "m\\frac{u^{n+1}-2u^n + u^{n-1}}{\\Delta t^2}\n", + "+ b\\frac{u^{n+1}-u^{n-1}}{2\\Delta t} + s(u^n) = F^n,\n", + "\\label{vib:ode2:step3b2} \\tag{4}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "which gives an explicit formula for $u$ at each\n", + "new time level:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "u^{n+1} = (2mu^n + (\\frac{b}{2}\\Delta t - m)u^{n-1} +\n", + "\\Delta t^2(F^n - s(u^n)))(m + \\frac{b}{2}\\Delta t)^{-1}\n", + "\\label{vib:ode2:step4} \\tag{5}\n", + "\\thinspace .\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For the first time step we need to discretize $u^{\\prime}(0)=V$\n", + "as $[D_{2t}u = V]^0$ and combine\n", + "with ([5](#vib:ode2:step4)) for $n=0$. The discretized initial condition\n", + "leads to" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "u^{-1} = u^{1} - 2\\Delta t V,\n", + "\\label{vib:ode2:ic:du} \\tag{6}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "which inserted in ([5](#vib:ode2:step4)) for $n=0$ gives an equation\n", + "that can be solved for\n", + "$u^1$:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "u^1 = u^0 + \\Delta t\\, V\n", + "+ \\frac{\\Delta t^2}{2m}(-bV - s(u^0) + F^0)\n", + "\\thinspace .\n", + "\\label{vib:ode2:step4b} \\tag{7}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## A centered scheme for quadratic damping\n", + "
\n", + "\n", + "When $f(u^{\\prime})=bu^{\\prime}|u^{\\prime}|$, we get a quadratic equation for $u^{n+1}$\n", + "in ([3](#vib:ode2:step3b)). This equation can be straightforwardly\n", + "solved by the well-known formula for the roots of a quadratic equation.\n", + "However, we can also avoid the nonlinearity by introducing\n", + "an approximation with an error of order no higher than what we\n", + "already have from replacing derivatives with finite differences.\n", + "\n", + "\n", + "We start with ([1](#vib:ode2)) and only replace\n", + "$u^{\\prime\\prime}$ by $D_tD_tu$, resulting in" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "[mD_tD_t u + bu^{\\prime}|u^{\\prime}| + s(u) = F]^n\\thinspace .\n", + "\\label{vib:ode2:quad:idea1} \\tag{8}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here, $u^{\\prime}|u^{\\prime}|$ is to be computed at time $t_n$. The idea\n", + "is now to introduce\n", + "a *geometric mean*, defined by" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "(w^2)^n \\approx w^{n-\\frac{1}{2}}w^{n+\\frac{1}{2}},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "for some quantity $w$ depending on time. The error in the geometric mean\n", + "approximation is $\\mathcal{O}({\\Delta t^2})$, the same as in the\n", + "approximation $u^{\\prime\\prime}\\approx D_tD_tu$. With $w=u^{\\prime}$ it follows\n", + "that" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "[u^{\\prime}|u^{\\prime}|]^n \\approx u^{\\prime}(t_{n+\\frac{1}{2}})|u^{\\prime}(t_{n-\\frac{1}{2}})|\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The next step is to approximate\n", + "$u^{\\prime}$ at $t_{n\\pm 1/2}$, and fortunately a centered difference\n", + "fits perfectly into the formulas since it involves $u$ values at\n", + "the mesh points only. With the approximations" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "u^{\\prime}(t_{n+1/2})\\approx [D_t u]^{n+\\frac{1}{2}},\\quad u^{\\prime}(t_{n-1/2})\\approx [D_t u]^{n-\\frac{1}{2}},\n", + "\\label{vib:ode2:quad:idea2} \\tag{9}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "we get" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "[u^{\\prime}|u^{\\prime}|]^n\n", + "\\approx [D_tu]^{n+\\frac{1}{2}}|[D_tu]^{n-\\frac{1}{2}}| = \\frac{u^{n+1}-u^n}{\\Delta t}\n", + "\\frac{|u^n-u^{n-1}|}{\\Delta t}\n", + "\\thinspace .\n", + "\\label{_auto2} \\tag{10}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The counterpart to ([3](#vib:ode2:step3b)) is then" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "m\\frac{u^{n+1}-2u^n + u^{n-1}}{\\Delta t^2}\n", + "+ b\\frac{u^{n+1}-u^n}{\\Delta t}\\frac{|u^n-u^{n-1}|}{\\Delta t}\n", + "+ s(u^n) = F^n,\n", + "\\label{vib:ode2:step3b:quad} \\tag{11}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "which is linear in the unknown $u^{n+1}$. Therefore, we can easily solve\n", + "([11](#vib:ode2:step3b:quad))\n", + "with respect to $u^{n+1}$ and achieve the explicit updating formula" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "u^{n+1} = \\left( m + b|u^n-u^{n-1}|\\right)^{-1}\\times \\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + " \\qquad \\left(2m u^n - mu^{n-1} + bu^n|u^n-u^{n-1}| + \\Delta t^2 (F^n - s(u^n))\n", + "\\right)\n", + "\\thinspace .\n", + "\\label{vib:ode2:step4:quad} \\tag{12}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "In the derivation of a special equation for the first\n", + "time step we run into some trouble: inserting ([6](#vib:ode2:ic:du))\n", + "in ([12](#vib:ode2:step4:quad)) for $n=0$ results in a complicated nonlinear\n", + "equation for $u^1$. By thinking differently about the problem we can\n", + "easily get away with the nonlinearity again. We have for $n=0$ that\n", + "$b[u^{\\prime}|u^{\\prime}|]^0 = bV|V|$. Using this value in ([8](#vib:ode2:quad:idea1))\n", + "gives" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "[mD_tD_t u + bV|V| + s(u) = F]^0\n", + "\\thinspace .\n", + "\\label{_auto3} \\tag{13}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Writing this equation out and using ([6](#vib:ode2:ic:du)) results in the\n", + "special equation for the first time step:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "u^1 = u^0 + \\Delta t V + \\frac{\\Delta t^2}{2m}\\left(-bV|V| - s(u^0) + F^0\\right)\n", + "\\thinspace .\n", + "\\label{vib:ode2:step4b:quad} \\tag{14}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## A forward-backward discretization of the quadratic damping term\n", + "\n", + "The previous section first proposed to discretize the quadratic\n", + "damping term $|u^{\\prime}|u^{\\prime}$ using centered differences:\n", + "$[|D_{2t}|D_{2t}u]^n$. As this gives rise to a nonlinearity in\n", + "$u^{n+1}$, it was instead proposed to use a geometric mean combined\n", + "with centered differences. But there are other alternatives. To get\n", + "rid of the nonlinearity in $[|D_{2t}|D_{2t}u]^n$, one can think\n", + "differently: apply a backward difference to $|u^{\\prime}|$, such that\n", + "the term involves known values, and apply a forward difference to\n", + "$u^{\\prime}$ to make the term linear in the unknown $u^{n+1}$. With\n", + "mathematics," + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "[\\beta |u^{\\prime}|u^{\\prime}]^n \\approx \\beta |[D_t^-u]^n|[D_t^+ u]^n =\n", + "\\beta\\left\\vert\\frac{u^n-u^{n-1}}{\\Delta t}\\right\\vert\n", + "\\frac{u^{n+1}-u^n}{\\Delta t}\\thinspace .\n", + "\\label{vib:ode2:nonlin:fbdiff} \\tag{15}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The forward and backward differences both have an error proportional\n", + "to $\\Delta t$ so one may think the discretization above leads to\n", + "a first-order scheme.\n", + "However, by looking at the formulas, we realize that the forward-backward\n", + "differences in ([15](#vib:ode2:nonlin:fbdiff))\n", + "result in exactly the same scheme as in\n", + "([11](#vib:ode2:step3b:quad)) where we\n", + "used a geometric mean and centered differences and committed errors\n", + "of size $\\mathcal{O}({\\Delta t^2})$. Therefore, the forward-backward\n", + "differences in ([15](#vib:ode2:nonlin:fbdiff))\n", + "act in a symmetric way and actually produce a second-order\n", + "accurate discretization of the quadratic damping term.\n", + "\n", + "\n", + "## Implementation\n", + "
\n", + "\n", + "\n", + "The algorithm arising from the methods in the sections [A centered scheme for linear damping](#vib:ode2:fdm:flin)\n", + "and [A centered scheme for quadratic damping](#vib:ode2:fdm:fquad) is very similar to the undamped case in\n", + "the section [vib:ode1:fdm](#vib:ode1:fdm). The difference is\n", + "basically a question of different formulas for $u^1$ and\n", + "$u^{n+1}$. This is actually quite remarkable. The equation\n", + "([1](#vib:ode2)) is normally impossible to solve by pen and paper, but\n", + "possible for some special choices of $F$, $s$, and $f$. On the\n", + "contrary, the complexity of the\n", + "nonlinear generalized model ([1](#vib:ode2)) versus the\n", + "simple undamped model is not a big deal when we solve the\n", + "problem numerically!\n", + "\n", + "The computational algorithm takes the form\n", + "\n", + "1. $u^0=I$\n", + "\n", + "2. compute $u^1$ from ([7](#vib:ode2:step4b)) if linear\n", + " damping or ([14](#vib:ode2:step4b:quad)) if quadratic damping\n", + "\n", + "3. for $n=1,2,\\ldots,N_t-1$:\n", + "\n", + "a. compute $u^{n+1}$ from ([5](#vib:ode2:step4)) if linear\n", + " damping or ([12](#vib:ode2:step4:quad)) if quadratic damping\n", + "\n", + "Instead of coding the explicit formulas of $u^1$ and $u^{n+1}$ by hand, the code below demonstrates how we can implement this in Devito\n", + "\n", + "Modifying the `solver` function for the undamped case is fairly\n", + "easy, the big difference being many more terms and if tests on\n", + "the type of damping:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from devito import *" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# NOTE:\n", + "# - testConstant gives same errors\n", + "# - testQuadratic and testSinusoidal error increase from 1E-15 to 1E-7\n", + "\n", + "def solver(I, V, m, b, s, F, dt, T, damping='linear'):\n", + " \"\"\"\n", + " Solve m*u'' + f(u') + s(u) = F(t) for t in (0,T],\n", + " u(0)=I and u'(0)=V,\n", + " by a central finite difference method with time step dt.\n", + " If damping is 'linear', f(u')=b*u, while if damping is\n", + " 'quadratic', f(u')=b*u'*abs(u').\n", + " F(t) and s(u) are Python functions.\n", + " \"\"\"\n", + " dt = float(dt); b = float(b); m = float(m) # avoid integer div.\n", + " Nt = int(round(T/dt))\n", + " t = Dimension('t', spacing=Constant('h_t'))\n", + " \n", + " u = TimeFunction(name='u', dimensions=(t,), shape=(Nt+1,), space_order=2, dtype=np.float64)\n", + " u.data[0] = I\n", + " \n", + " # Initialize external function F\n", + " f = TimeFunction(name='f', dimensions=(t,), shape=(Nt+1,), dtype=np.float64)\n", + " for time in range(0, Nt+1):\n", + " f.data[time] = F(time*dt)\n", + " \n", + " # Initial Condition\n", + " if damping == 'linear':\n", + " u.data[1] = u.data[0] + dt*V + dt**2/(2*m)*(-b*V - s(u.data[0]) + f.data[0])\n", + " elif damping == 'quadratic':\n", + " u.data[1] = u.data[0] + dt*V + dt**2/(2*m)*(-b*V*abs(V) - s(u.data[0]) + f.data[0])\n", + " \n", + " # Define equation\n", + " if damping == 'linear':\n", + " eq_u = Eq(m*u.dt2 + b*(u.forward - u.backward)/(2 * t.spacing) + s(u), f)\n", + " elif damping == 'quadratic':\n", + " eq_u = Eq(m*u.dt2 + b*(u.forward - u)*abs(u - u.backward)/t.spacing**2 + s(u), f)\n", + " \n", + " # Solving the equation\n", + " stencil_u = solve(eq_u.evaluate, u.forward)\n", + " update_u = Eq(u.forward, stencil_u)\n", + " \n", + " op = Operator([update_u])\n", + " op.apply(h_t=dt, t_M=Nt-1)\n", + " return u.data, np.linspace(0, Nt*dt, Nt+1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The complete code (`solver` and relevant tests) resides in the file [`vib.py`](${src_vib}/vib.py).\n", + "\n", + "## Verification\n", + "
\n", + "\n", + "### Constant solution\n", + "\n", + "For debugging and initial verification, a constant solution is often\n", + "very useful. We choose $u_e(t)=I$, which implies $V=0$.\n", + "Inserted in the ODE, we get\n", + "$F(t)=s(I)$ for any choice of $f$. Since the discrete derivative\n", + "of a constant vanishes (in particular, $[D_{2t}I]^n=0$,\n", + "$[D_tI]^n=0$, and $[D_tD_t I]^n=0$), the constant solution also fulfills\n", + "the discrete equations. The constant should therefore be reproduced\n", + "to machine precision. The function `test_constant` in `vib.py`\n", + "implements this test.\n", + "\n", + "### Linear solution\n", + "\n", + "Now we choose a linear solution: $u_e(t) = ct + d$. The initial condition\n", + "$u(0)=I$ implies $d=I$, and $u^{\\prime}(0)=V$ forces $c$ to be $V$.\n", + "Inserting $u_e(t)=Vt+I$ in the ODE with linear damping results in" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "0 + bV + s(Vt+I) = F(t),\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "while quadratic damping requires the source term" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "0 + b|V|V + s(Vt+I) = F(t)\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since the finite difference approximations used to compute $u^{\\prime}$ all\n", + "are exact for a linear function, it turns out that the linear $u_e$\n", + "is also a solution of the discrete equations.\n", + "[vib:exer:verify:gen:linear](#vib:exer:verify:gen:linear) asks you to carry out\n", + "all the details.\n", + "\n", + "### Quadratic solution\n", + "\n", + "Choosing $u_e = bt^2 + Vt + I$, with $b$ arbitrary,\n", + "fulfills the initial conditions and\n", + "fits the ODE if $F$ is adjusted properly. The solution also solves\n", + "the discrete equations with linear damping. However, this quadratic\n", + "polynomial in $t$ does not fulfill the discrete equations in case\n", + "of quadratic damping, because the geometric mean used in the approximation\n", + "of this term introduces an error.\n", + "Doing [vib:exer:verify:gen:linear](#vib:exer:verify:gen:linear) will reveal\n", + "the details. One can fit $F^n$ in the discrete equations such that\n", + "the quadratic polynomial is reproduced by the numerical method (to\n", + "machine precision).\n", + "\n", + "# DO WE DELETE VISUALIZATION & USER INTERFACE SECTIONS?\n", + "\n", + "## Visualization\n", + "
\n", + "\n", + "The functions for visualizations differ significantly from\n", + "those in the undamped case in the `vib_undamped.py` program because,\n", + "in the present general case, we do not have an exact solution to\n", + "include in the plots. Moreover, we have no good estimate of\n", + "the periods of the oscillations as there will be one period\n", + "determined by the system parameters, essentially the\n", + "approximate frequency $\\sqrt{s'(0)/m}$ for linear $s$ and small damping,\n", + "and one period dictated by $F(t)$ in case the excitation is periodic.\n", + "This is, however,\n", + "nothing that the program can depend on or make use of.\n", + "Therefore, the user has to specify $T$ and the window width\n", + "to get a plot that moves with the graph and shows\n", + "the most recent parts of it in long time simulations.\n", + "\n", + "The `vib.py` code\n", + "contains several functions for analyzing the time series signal\n", + "and for visualizing the solutions.\n", + "\n", + "## User interface\n", + "
\n", + "\n", + "\n", + "The `main` function is changed substantially from\n", + "the `vib_undamped.py` code, since we need to\n", + "specify the new data $c$, $s(u)$, and $F(t)$. In addition, we must\n", + "set $T$ and the plot window width (instead of the number of periods we\n", + "want to simulate as in `vib_undamped.py`). To figure out whether we\n", + "can use one plot for the whole time series or if we should follow the\n", + "most recent part of $u$, we can use the `plot_empricial_freq_and_amplitude`\n", + "function's estimate of the number of local maxima. This number is now\n", + "returned from the function and used in `main` to decide on the\n", + "visualization technique." + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "def main():\n", + " import argparse\n", + " parser = argparse.ArgumentParser()\n", + " parser.add_argument('--I', type=float, default=1.0)\n", + " parser.add_argument('--V', type=float, default=0.0)\n", + " parser.add_argument('--m', type=float, default=1.0)\n", + " parser.add_argument('--c', type=float, default=0.0)\n", + " parser.add_argument('--s', type=str, default='u')\n", + " parser.add_argument('--F', type=str, default='0')\n", + " parser.add_argument('--dt', type=float, default=0.05)\n", + " parser.add_argument('--T', type=float, default=140)\n", + " parser.add_argument('--damping', type=str, default='linear')\n", + " parser.add_argument('--window_width', type=float, default=30)\n", + " parser.add_argument('--savefig', action='store_true')\n", + " a = parser.parse_args()\n", + " from scitools.std import StringFunction\n", + " s = StringFunction(a.s, independent_variable='u')\n", + " F = StringFunction(a.F, independent_variable='t')\n", + " I, V, m, c, dt, T, window_width, savefig, damping = \\\n", + " a.I, a.V, a.m, a.c, a.dt, a.T, a.window_width, a.savefig, \\\n", + " a.damping\n", + "\n", + " u, t = solver(I, V, m, c, s, F, dt, T)\n", + " num_periods = empirical_freq_and_amplitude(u, t)\n", + " if num_periods <= 15:\n", + " figure()\n", + " visualize(u, t)\n", + " else:\n", + " visualize_front(u, t, window_width, savefig)\n", + " show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The program `vib.py` contains\n", + "the above code snippets and can solve the model problem\n", + "([1](#vib:ode2)). As a demo of `vib.py`, we consider the case\n", + "$I=1$, $V=0$, $m=1$, $c=0.03$, $s(u)=\\sin(u)$, $F(t)=3\\cos(4t)$,\n", + "$\\Delta t = 0.05$, and $T=140$. The relevant command to run is" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " Terminal> python vib.py --s 'sin(u)' --F '3*cos(4*t)' --c 0.03\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This results in a [moving window following the function](${doc_notes}/vib/html//mov-vib/vib_generalized_dt0.05/index.html) on the screen.\n", + "[Figure](#vib:ode2:fig:demo) shows a part of the time series.\n", + "\n", + "\n", + "\n", + "
\n", + "\n", + "

Damped oscillator excited by a sinusoidal function.

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "## The Euler-Cromer scheme for the generalized model\n", + "
\n", + "\n", + "The ideas of the Euler-Cromer method from the section [vib:model2x2:EulerCromer](#vib:model2x2:EulerCromer)\n", + "carry over to the generalized model. We write ([1](#vib:ode2))\n", + "as two equations for $u$ and $v=u^{\\prime}$. The first equation is taken as the\n", + "one with $v^{\\prime}$ on the left-hand side:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "v^{\\prime} = \\frac{1}{m}(F(t)-s(u)-f(v)),\n", + "\\label{vib:ode2:EulerCromer:veq} \\tag{16}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "u^{\\prime} = v\\thinspace .\n", + "\\label{vib:ode2:EulerCromer:ueq} \\tag{17}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Again, the idea is to step ([16](#vib:ode2:EulerCromer:veq)) forward using\n", + "a standard Forward Euler method, while we update $u$ from\n", + "([17](#vib:ode2:EulerCromer:ueq)) with a Backward Euler method,\n", + "utilizing the recent, computed $v^{n+1}$ value. In detail," + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\frac{v^{n+1}-v^n}{\\Delta t} = \\frac{1}{m}(F(t_n)-s(u^n)-f(v^n)),\n", + "\\label{vib:ode2:EulerCromer:dveq0a} \\tag{18}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\frac{u^{n+1}-u^n}{\\Delta t} = v^{n+1},\n", + "\\label{vib:ode2:EulerCromer:dueq0a} \\tag{19}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "resulting in the explicit scheme" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "v^{n+1} = v^n + \\Delta t\\frac{1}{m}(F(t_n)-s(u^n)-f(v^n)),\n", + "\\label{vib:ode2:EulerCromer:dveq} \\tag{20}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "u^{n+1} = u^n + \\Delta t\\,v^{n+1}\\thinspace .\n", + "\\label{vib:ode2:EulerCromer:dueq0} \\tag{21}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We immediately note one very favorable feature of this scheme: all the\n", + "nonlinearities in $s(u)$ and $f(v)$ are evaluated at a previous time\n", + "level. This makes the Euler-Cromer method easier to apply and\n", + "hence much more convenient than the centered scheme for the second-order\n", + "ODE ([1](#vib:ode2)).\n", + "\n", + "The initial conditions are trivially set as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "v^0 = V,\n", + "\\label{_auto4} \\tag{22}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "u^0 = I\\thinspace .\n", + "\\label{_auto5} \\tag{23}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[hpl 1: odespy for the generalized problem]\n", + "\n", + "\n", + "\n", + "\n", + "## The Stoermer-Verlet algorithm for the generalized model\n", + "
\n", + "\n", + "We can easily apply the ideas from the section [vib:model2x2:StormerVerlet](#vib:model2x2:StormerVerlet) to\n", + "extend that method to the generalized model" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{align*}\n", + "v^{\\prime} &= \\frac{1}{m}(F(t)-s(u)-f(v)),\\\\ \n", + "u^{\\prime} &= v\\thinspace .\n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "However, since the scheme is essentially centered differences for\n", + "the ODE system on a staggered mesh, we do not go into detail here,\n", + "but refer to the section [A staggered Euler-Cromer scheme for a generalized model](#vib:ode2:staggered).\n", + "\n", + "## A staggered Euler-Cromer scheme for a generalized model\n", + "
\n", + "\n", + "The more general model for vibration problems," + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "mu'' + f(u') + s(u) = F(t),\\quad u(0)=I,\\ u'(0)=V,\\ t\\in (0,T],\n", + "\\label{_auto6} \\tag{24}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "can be rewritten as a first-order ODE system" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "v' = m^{-1}\\left(F(t) - f(v) - s(u)\\right),\n", + "\\label{vib:ode2:staggered:veq} \\tag{25}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "u' = v\\thinspace .\n", + "\\label{vib:ode2:staggered:ueq} \\tag{26}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It is natural to introduce a staggered mesh (see the section [vib:model2x2:staggered](#vib:model2x2:staggered)) and seek $u$ at mesh points $t_n$ (the numerical value is\n", + "denoted by $u^n$) and $v$ between mesh points at $t_{n+1/2}$ (the numerical\n", + "value is denoted by $v^{n+\\frac{1}{2}}$).\n", + "A centered difference approximation to ([26](#vib:ode2:staggered:ueq))-([25](#vib:ode2:staggered:veq)) can then be written in operator notation as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\lbrack D_tv = m^{-1}\\left(F(t) - f(v) - s(u)\\right)\\rbrack^n,\n", + "\\label{vib:ode2:staggered:dveq} \\tag{27}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\lbrack D_t u = v\\rbrack^{n+\\frac{1}{2}}\\thinspace .\n", + "\\label{vib:ode2:staggered:dueq} \\tag{28}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Written out," + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\frac{v^{n+\\frac{1}{2}} - v^{n-\\frac{1}{2}}}{\\Delta t}\n", + "= m^{-1}\\left(F^n - f(v^n) - s(u^n)\\right),\n", + "\\label{vib:ode2:staggered:dveq2} \\tag{29}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\frac{u^n - u^{n-1}}{\\Delta t} = v^{n+\\frac{1}{2}}\\thinspace .\n", + "\\label{vib:ode2:staggered:dueq2} \\tag{30}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With linear damping, $f(v)=bv$, we can use an arithmetic mean\n", + "for $f(v^n)$: $f(v^n)\\approx = \\frac{1}{2}(f(v^{n-\\frac{1}{2}}) +\n", + "f(v^{n+\\frac{1}{2}}))$. The system\n", + "([29](#vib:ode2:staggered:dveq2))-([30](#vib:ode2:staggered:dueq2))\n", + "can then be solved with respect to the unknowns $u^n$ and $v^{n+\\frac{1}{2}}$:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "v^{n+\\frac{1}{2}} = \\left(1 + \\frac{b}{2m}\\Delta t\\right)^{-1}\\left(\n", + "v^{n-\\frac{1}{2}} + {\\Delta t}\n", + "m^{-1}\\left(F^n - {\\frac{1}{2}}f(v^{n-\\frac{1}{2}}) - s(u^n)\\right)\\right),\n", + "\\label{vib:ode2:staggered:v:scheme:lin} \\tag{31}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "u^n = u^{n-1} + {\\Delta t}v^{n-\\frac{1}{2}}\\thinspace .\n", + "\\label{vib:ode2:staggered:u:scheme:lin} \\tag{32}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In case of quadratic damping, $f(v)=b|v|v$, we can use a geometric mean:\n", + "$f(v^n)\\approx b|v^{n-\\frac{1}{2}}|v^{n+\\frac{1}{2}}$. Inserting this approximation\n", + "in ([29](#vib:ode2:staggered:dveq2))-([30](#vib:ode2:staggered:dueq2)) and\n", + "solving for the unknowns $u^n$ and $v^{n+\\frac{1}{2}}$ results in" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "v^{n+\\frac{1}{2}} = (1 + \\frac{b}{m}|v^{n-\\frac{1}{2}}|\\Delta t)^{-1}\\left(\n", + "v^{n-\\frac{1}{2}} + {\\Delta t}\n", + "m^{-1}\\left(F^n - s(u^n)\\right)\\right),\n", + "\\label{vib:ode2:staggered:v:scheme:quad} \\tag{33}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "u^n = u^{n-1} + {\\Delta t}v^{n-\\frac{1}{2}}\\thinspace .\n", + "\\label{vib:ode2:staggered:u:scheme:quad} \\tag{34}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The initial conditions are derived at the end of\n", + "the section [vib:model2x2:staggered](#vib:model2x2:staggered):" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "u^0 = I,\n", + "\\label{vib:ode2:staggered:u02} \\tag{35}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "v^\\frac{1}{2} = V - \\frac{1}{2}\\Delta t\\omega^2I\n", + "\\label{vib:ode2:staggered:v02} \\tag{36}\\thinspace .\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The PEFRL 4th-order accurate algorithm\n", + "
\n", + "\n", + "A variant of the Euler-Cromer type of algorithm, which provides an\n", + "error $\\mathcal{O}({\\Delta t^4})$ if $f(v)=0$, is called PEFRL\n", + "[[PEFRL_2002]](#PEFRL_2002). This algorithm is very well suited for integrating\n", + "dynamic systems (especially those without damping) over very long time\n", + "periods. Define" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "g(u,v) = \\frac{1}{m}(F(t)-s(u)-f(v))\\thinspace .\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The algorithm is explicit and features these steps:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "u^{n+1,1} = u^n + \\xi\\Delta t v^n,\n", + "\\label{_auto7} \\tag{37}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "v^{n+1,1} = v^n + \\frac{1}{2}(1-2\\lambda)\\Delta t g(u^{n+1,1},v^n),\n", + "\\label{_auto8} \\tag{38}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "u^{n+1,2} = u^{n+1,1} + \\chi\\Delta t v^{n+1,1},\n", + "\\label{_auto9} \\tag{39}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "v^{n+1,2} = v^{n+1,1} + \\lambda\\Delta t g(u^{n+1,2}, v^{n+1,1}),\n", + "\\label{_auto10} \\tag{40}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "u^{n+1,3} = u^{n+1,2} + (1-2(\\chi + \\xi))\\Delta t v^{n+1,2},\n", + "\\label{_auto11} \\tag{41}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "v^{n+1,3} = v^{n+1,2} + \\lambda\\Delta t g(u^{n+1,3}, v^{n+1,2}),\n", + "\\label{_auto12} \\tag{42}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "u^{n+1,4} = u^{n+1,3} + \\chi\\Delta t v^{n+1,3},\n", + "\\label{_auto13} \\tag{43}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "v^{n+1} = v^{n+1,3} + \\frac{1}{2}(1-2\\lambda)\\Delta t g(u^{n+1,4},v^{n+1,3}),\n", + "\\label{_auto14} \\tag{44}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "u^{n+1} = u^{n+1,4} + \\xi\\Delta t v^{n+1}\n", + "\\label{_auto15} \\tag{45}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "The parameters $\\xi$, $\\lambda$, and $\\xi$ have the values" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\xi = 0.1786178958448091,\n", + "\\label{_auto16} \\tag{46}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\lambda = -0.2123418310626054,\n", + "\\label{_auto17} \\tag{47}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\chi = -0.06626458266981849\n", + "\\label{_auto18} \\tag{48}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exercises and Problems\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "## Exercise 1: Implement the solver via classes (DELETE)?\n", + "
\n", + "\n", + "Reimplement the `vib.py` program using a class `Problem` to hold all\n", + "the physical parameters of the problem, a class `Solver` to hold the\n", + "numerical parameters and compute the solution, and a class\n", + "`Visualizer` to display the solution.\n", + "\n", + "\n", + "\n", + "**Hint.**\n", + "Use the ideas and examples from Sections 5.5.1 and 5.5.2 \n", + "in [[Langtangen_decay]](#Langtangen_decay). More specifically, make a superclass\n", + "`Problem` for holding the scalar physical parameters of a problem and\n", + "let subclasses implement the $s(u)$ and $F(t)$ functions as methods.\n", + "Try to call up as much existing functionality in `vib.py` as possible.\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "The complete code looks like this." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Reimplementation of vib.py using classes\n", + "\n", + "import numpy as np\n", + "import scitools.std as plt\n", + "import sympy as sym\n", + "from vib import solver as vib_solver\n", + "from vib import visualize as vib_visualize\n", + "from vib import visualize_front as vib_visualize_front\n", + "from vib import visualize_front_ascii as vib_visualize_front_ascii\n", + "from vib import plot_empirical_freq_and_amplitude as \\\n", + " vib_plot_empirical_freq_and_amplitude\n", + "\n", + "class Vibration:\n", + " '''\n", + " Problem: m*u'' + f(u') + s(u) = F(t) for t in (0,T],\n", + " u(0)=I and u'(0)=V. The problem is solved\n", + " by a central finite difference method with time step dt.\n", + " If damping is 'linear', f(u')=b*u, while if damping is\n", + " 'quadratic', f(u')=b*u'*abs(u'). Zero damping is achieved\n", + " with b=0. F(t) and s(u) are Python functions.\n", + " '''\n", + " def __init__(self, I=1, V=0, m=1, b=0, damping='linear'):\n", + " self.I = I; self.V = V; self.m = m; self.b=b;\n", + " self.damping = damping\n", + " def s(self, u):\n", + " return u\n", + " def F(self, t):\n", + " '''Driving force. Zero implies free oscillations'''\n", + " return 0\n", + "\n", + "class Free_vibrations(Vibration):\n", + " '''F(t) = 0'''\n", + " def __init__(self, s=None, I=1, V=0, m=1, b=0, damping='linear'):\n", + " Vibration.__init__(self, I=I, V=V, m=m, b=b, damping=damping)\n", + " if s != None:\n", + " self.s = s\n", + "\n", + "class Forced_vibrations(Vibration):\n", + " '''F(t)! = 0'''\n", + " def __init__(self, F, s=None, I=1, V=0, m=1, b=0,\n", + " damping='linear'):\n", + " Vibration.__init__(self, I=I, V=V, m=m, b=b,\n", + " damping=damping)\n", + " if s != None:\n", + " self.s = s\n", + " self.F = F\n", + "\n", + "class Solver:\n", + " def __init__(self, dt=0.05, T=20):\n", + " self.dt = dt; self.T = T\n", + "\n", + " def solve(self, problem):\n", + " self.u, self.t = vib_solver(\n", + " problem.I, problem.V,\n", + " problem.m, problem.b,\n", + " problem.s, problem.F,\n", + " self.dt, self.T, problem.damping)\n", + "\n", + "class Visualizer:\n", + " def __init__(self, problem, solver, window_width, savefig):\n", + " self.problem = problem; self.solver = solver\n", + " self.window_width = window_width; self.savefig = savefig\n", + " def visualize(self):\n", + " u = self.solver.u; t = self.solver.t # short forms\n", + " num_periods = vib_plot_empirical_freq_and_amplitude(u, t)\n", + " if num_periods <= 40:\n", + " plt.figure()\n", + " vib_visualize(u, t)\n", + " else:\n", + " vib_visualize_front(u, t, self.window_width, self.savefig)\n", + " vib_visualize_front_ascii(u, t)\n", + " plt.show()\n", + "\n", + "def main():\n", + " # Note: the reading of parameter values would better be done\n", + " # from each relevant class, i.e. class Problem should read I, V,\n", + " # etc., while class Solver should read dt and T, and so on.\n", + " # Consult, e.g., Langtangen: \"A Primer on Scientific Programming\",\n", + " # App E.\n", + " import argparse\n", + " parser = argparse.ArgumentParser()\n", + " parser.add_argument('--I', type=float, default=1.0)\n", + " parser.add_argument('--V', type=float, default=0.0)\n", + " parser.add_argument('--m', type=float, default=1.0)\n", + " parser.add_argument('--b', type=float, default=0.0)\n", + " parser.add_argument('--s', type=str, default=None)\n", + " parser.add_argument('--F', type=str, default='0')\n", + " parser.add_argument('--dt', type=float, default=0.05)\n", + " parser.add_argument('--T', type=float, default=20)\n", + " parser.add_argument('--window_width', type=float, default=30.,\n", + " help='Number of periods in a window')\n", + " parser.add_argument('--damping', type=str, default='linear')\n", + " parser.add_argument('--savefig', action='store_true')\n", + " # Hack to allow --SCITOOLS options\n", + " # (scitools.std reads this argument at import)\n", + " parser.add_argument('--SCITOOLS_easyviz_backend',\n", + " default='matplotlib')\n", + " a = parser.parse_args()\n", + "\n", + " from scitools.std import StringFunction\n", + " if a.s != None:\n", + " s = StringFunction(a.s, independent_variable='u')\n", + " else:\n", + " s = None\n", + " F = StringFunction(a.F, independent_variable='t')\n", + "\n", + " if a.F == '0': # free vibrations\n", + " problem = Free_vibrations(s=s, I=a.I, V=a.V, m=a.m, b=a.b,\n", + " damping=a.damping)\n", + " else: # forced vibrations\n", + " problem = Forced_vibrations(lambda t: np.sin(t),\n", + " s=s, I=a.I, V=a.V,\n", + " m=a.m, b=a.b, damping=a.damping)\n", + "\n", + " solver = Solver(dt=a.dt, T=a.T)\n", + " solver.solve(problem)\n", + "\n", + " visualizer = Visualizer(problem, solver,\n", + " a.window_width, a.savefig)\n", + " visualizer.visualize()\n", + "\n", + "if __name__ == '__main__':\n", + " main()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "Filename: `vib_class`.\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "## Problem 2: Use a backward difference for the damping term\n", + "
\n", + "\n", + "As an alternative to discretizing the damping terms $\\beta u^{\\prime}$ and\n", + "$\\beta |u^{\\prime}|u^{\\prime}$ by centered differences, we may apply\n", + "backward differences:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{align*}\n", + "[u^{\\prime}]^n &\\approx [D_t^-u]^n,\\\\ \n", + "[|u^{\\prime}|u^{\\prime}]^n &\\approx [|D_t^-u|D_t^-u]^n\\\\ \n", + "&= |[D_t^-u]^n|[D_t^-u]^n\\thinspace .\n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The advantage of the backward difference is that the damping term is\n", + "evaluated using known values $u^n$ and $u^{n-1}$ only.\n", + "Extend the [`vib.py`](${src_vib}/vib.py) code with a scheme based\n", + "on using backward differences in the damping terms. Add statements\n", + "to compare the original approach with centered difference and the\n", + "new idea launched in this exercise. Perform numerical experiments\n", + "to investigate how much accuracy that is lost by using the backward\n", + "differences.\n", + "\n", + "\n", + "\n", + "**Solution.**\n", + "The new discretization approach of the linear and quadratic damping\n", + "terms calls for new derivations of the updating formulas (for $u$) in\n", + "the solver. Since backward difference approximations will be used for\n", + "the damping term, we may also use this approximation for the initial\n", + "condition on $u^{\\prime}(0)$ without deteriorating the convergence\n", + "rate any further. Note that introducing backward difference\n", + "approximations for the damping term make our computational schemes\n", + "first order, as opposed to the original second order schemes which\n", + "used central difference approximations also for the damping terms. The\n", + "motivation for also using a backward difference approximation for the\n", + "initial condition on $u^{\\prime}(0)$, is simply that the computational\n", + "schemes get much simpler.\n", + "\n", + "With linear damping, the new discretized form of the equation reads" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "m\\frac{u^{n+1}-2u^n+u^{n-1}}{\\Delta t^2} + b\\frac{u^n - u^{n-1}}{\\Delta t} + s(u^n) = F^n,\n", + "\\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "which gives us" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "u^{n+1} = \\left(2-\\frac{\\Delta t b}{m}\\right)u^n + \\frac{\\Delta t^2}{m}\\left(F^n - s(u^n)\\right) + \\left(\\frac{\\Delta t b}{m} - 1\\right)u^{n-1}.\n", + "\\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With $n=0$, the updating formula becomes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "u^1 = \\left(2-\\frac{\\Delta t b}{m}\\right)u^0 + \\frac{\\Delta t^2}{m}\\left(F^0 - s(u^0)\\right) + \\left(\\frac{\\Delta t b}{m} - 1\\right)u^{-1},\n", + "\\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "which requires some further elaboration because of the unknown $u^{-1}$. We handle this by discretizing the initial condition $u^{\\prime}(0) = V$ by a backward difference approximation as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\frac{u^0 - u^{-1}}{\\Delta t} = V,\n", + "\\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "which implies that" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "u^{-1} = u^0 - \\Delta t V.\n", + "\\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Inserting this expression for $u^{-1}$ in the updating formula for $u^{n+1}$, and simplifying,\n", + "gives us the following special formula for the first time step:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "u^1 = u^0 + \\Delta t V + \\frac{\\Delta t^2}{m}\\left(-bV - s(u^0) + F^0\\right).\n", + "\\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Switching to quadratic damping, the new discretized form of the equations becomes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "m\\frac{u^{n+1}-2u^n+u^{n-1}}{\\Delta t^2} + b|\\frac{u^n - u^{n-1}}{\\Delta t}|\\frac{u^n - u^{n-1}}{\\Delta t} + s(u^n) = F^n,\n", + "\\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "which leads to" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "u^{n+1} = 2u^n - u^{n-1} - \\frac{b}{m}|u^n - u^{n-1}|(u^n - u^{n-1}) + \\frac{\\Delta t^2}{m}\\left(F^n - s(u^n)\\right).\n", + "\\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With $n=0$, this updating formula becomes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "u^1 = 2u^0 - u^{-1} - \\frac{b}{m}|u^0 - u^{-1}|(u^0 - u^{-1}) + \\frac{\\Delta t^2}{m}\\left(F^0 - s(u^0)\\right).\n", + "\\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Again, we handle the unknown $u^{-1}$ via the same expression as above, which be derived from a backward difference approximation to the initial condition on the derivative. Inserting this expression for $u^{-1}$ and simplifying, gives the special updating formula for $u^1$ as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "u^1 = u^0 + \\Delta t V + \\frac{\\Delta t^2}{m}\\left(-b|V|V - s(u^0) + F^0\\right).\n", + "\\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We implement these new computational schemes in a new solver function\n", + "`solver_bwdamping`, so that the discrete solution for $u$ can be found\n", + "by both the original and the new solver. The difference between the\n", + "two different solutions is then visualized in the same way as the\n", + "original solution in `main`.\n", + "\n", + "The convergence rates computed in `test_mms` demonstrates that our\n", + "scheme now is a first order scheme, as $r$ is seen to approach 1.0\n", + "with decreasing $\\Delta t$.\n", + "\n", + "Both solvers reproduce a constant solution exactly (within machine\n", + "precision), whereas sinusoidal and quadratic solutions differ, as\n", + "should be expected after comparing the schemes. Pointing out the\n", + "\"best\" approach is difficult: the backward differences yield a much\n", + "simpler mathematical problem to be solved, while the more complicated\n", + "method converges faster and gives more accuracy for the same cost. On\n", + "the other hand, the backward differences can yield any reasonable\n", + "accuracy by lowering $\\Delta t$, and the results are obtained within a\n", + "few seconds on a laptop.\n", + "\n", + "Here is the complete computer code, arising from copying `vib.py` and changing\n", + "the functions that have to be changed:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# THIS CODE BELOW ?DOESN'T WORK IN THE FIRST PLACE" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import scitools as plt\n", + "\n", + "import os, sys\n", + "sys.path.insert(0, os.path.join(os.path.abspath(os.getcwd()), 'src-vib'))\n", + "from vib import (plot_empirical_freq_and_amplitude,\n", + " visualize_front, visualize_front_ascii,\n", + " minmax, periods, amplitudes,\n", + " solver as solver2)" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.18200549280924477\n", + "LINEAR DIFFERENCE ALL: \n", + "[0. 0.054 0.10314 0.1429974 0.16998503 0.18167401\n", + " 0.17701233 0.15641954 0.12174899 0.07612104 0.02364218 0.03096447\n", + " 0.08278431 0.12715357 0.16007901 0.17859733 0.1810419 0.1671927\n", + " 0.13829615 0.09695295 0.04688398 0.00740454 0.06102666 0.10915637\n", + " 0.14746202 0.17249608 0.18200549 0.17513441 0.15250124 0.11614295\n", + " 0.06933179]\n", + "0.18200549280924297\n" + ] + } + ], + "source": [ + "# both parts of \"testSinusoidal\" doesn't pass using original code (significantly larger error, 1E-14 required, 1E-1 observed)\n", + "# Same problem with the \"linaer\" section of \"testQuadratic\"\n", + "\n", + "def solver_bwdamping(I, V, m, b, s, F, dt, T, damping='linear'):\n", + " \"\"\"\n", + " Solve m*u'' + f(u') + s(u) = F(t) for t in (0,T],\n", + " u(0)=I and u'(0)=V. All terms except damping is discretized\n", + " by a central finite difference method with time step dt.\n", + " The damping term is discretized by a backward diff. approx.,\n", + " as is the init.cond. u'(0). If damping is 'linear', f(u')=b*u,\n", + " while if damping is 'quadratic', f(u')=b*u'*abs(u').\n", + " F(t) and s(u) are Python functions.\n", + " \"\"\"\n", + " dt = float(dt); b = float(b); m = float(m) # avoid integer div.\n", + " Nt = int(round(T/dt))\n", + " u = np.zeros(Nt+1)\n", + " t = np.linspace(0, Nt*dt, Nt+1)\n", + "\n", + " u_original = np.zeros(Nt+1); u_original[0] = I # for testing\n", + "\n", + " u[0] = I\n", + " if damping == 'linear':\n", + " u[1] = u[0] + dt*V + dt**2/m*(-b*V - s(u[0]) + F(t[0]))\n", + " elif damping == 'quadratic':\n", + " u[1] = u[0] + dt*V + \\\n", + " dt**2/m*(-b*V*abs(V) - s(u[0]) + F(t[0]))\n", + " for n in range(1, Nt):\n", + " if damping == 'linear':\n", + " u[n+1] = (2 - dt*b/m)*u[n] + dt**2/m*(- s(u[n]) + \\\n", + " F(t[n])) + (dt*b/m - 1)*u[n-1]\n", + " elif damping == 'quadratic':\n", + " u[n+1] = 2*u[n] - u[n-1] - b/m*abs(u[n] - \\\n", + " u[n-1])*(u[n] - u[n-1]) + dt**2/m*(-s(u[n]) + F(t[n]))\n", + " return u, t\n", + "\n", + "\n", + "import sympy as sym\n", + "\n", + "def test_constant():\n", + " \"\"\"Verify a constant solution.\"\"\"\n", + " u_exact = lambda t: I\n", + " I = 1.2; V = 0; m = 2; b = 0.9\n", + " w = 1.5\n", + " s = lambda u: w**2*u\n", + " F = lambda t: w**2*u_exact(t)\n", + " dt = 0.2\n", + " T = 2\n", + " #u, t = solver(I, V, m, b, s, F, dt, T, 'linear')\n", + " u, t = solver_bwdamping(I, V, m, b, s, F, dt, T, 'linear')\n", + " difference = np.abs(u_exact(t) - u).max()\n", + " print(difference)\n", + " tol = 1E-13\n", + " assert difference < tol\n", + "\n", + " #u, t = solver(I, V, m, b, s, F, dt, T, 'quadratic')\n", + " u, t = solver_bwdamping(I, V, m, b, s, F, dt, T, 'quadratic')\n", + " difference = np.abs(u_exact(t) - u).max()\n", + " print (difference)\n", + " assert difference < tol\n", + "\n", + "def lhs_eq(t, m, b, s, u, damping='linear'):\n", + " \"\"\"Return lhs of differential equation as sympy expression.\"\"\"\n", + " v = sym.diff(u, t)\n", + " if damping == 'linear':\n", + " return m*sym.diff(u, t, t) + b*v + s(u)\n", + " else:\n", + " return m*sym.diff(u, t, t) + b*v*sym.Abs(v) + s(u)\n", + "\n", + "def test_quadratic():\n", + " \"\"\"Verify a quadratic solution.\"\"\"\n", + " I = 1.2; V = 3; m = 2; b = 0.9\n", + " s = lambda u: 4*u\n", + " t = sym.Symbol('t')\n", + " dt = 0.2\n", + " T = 2\n", + "\n", + " q = 2 # arbitrary constant\n", + " u_exact = I + V*t + q*t**2\n", + " F = sym.lambdify(t, lhs_eq(t, m, b, s, u_exact, 'linear'))\n", + " u_exact = sym.lambdify(t, u_exact, modules='numpy')\n", + " #u1, t1 = solver(I, V, m, b, s, F, dt, T, 'linear')\n", + " u1, t1 = solver_bwdamping(I, V, m, b, s, F, dt, T, 'linear')\n", + " diff = np.abs(u_exact(t1) - u1).max()\n", + " print (diff)\n", + " tol = 1E-13\n", + " #assert diff < tol\n", + "\n", + " # In the quadratic damping case, u_exact must be linear\n", + " # in order to exactly recover this solution\n", + " u_exact = I + V*t\n", + " F = sym.lambdify(t, lhs_eq(t, m, b, s, u_exact, 'quadratic'))\n", + " u_exact = sym.lambdify(t, u_exact, modules='numpy')\n", + " #u2, t2 = solver(I, V, m, b, s, F, dt, T, 'quadratic')\n", + " u2, t2 = solver_bwdamping(I, V, m, b, s, F, dt, T, 'quadratic')\n", + " diff = np.abs(u_exact(t2) - u2).max()\n", + " print (diff)\n", + " #assert diff < tol\n", + "\n", + "def test_sinusoidal():\n", + " \"\"\"Verify a numerically exact sinusoidal solution when b=F=0.\"\"\"\n", + " from math import asin\n", + "\n", + " def u_exact(t):\n", + " w_numerical = 2/dt*np.arcsin(w*dt/2)\n", + " return I*np.cos(w_numerical*t)\n", + "\n", + " I = 1.2; V = 0; m = 2; b = 0\n", + " w = 1.5 # fix the frequency\n", + " s = lambda u: m*w**2*u\n", + " F = lambda t: 0\n", + " dt = 0.2\n", + " T = 6\n", + " #u, t = solver(I, V, m, b, s, F, dt, T, 'linear')\n", + " u, t = solver_bwdamping(I, V, m, b, s, F, dt, T, 'linear')\n", + " diff = np.abs(u_exact(t) - u).max()\n", + " print (diff)\n", + " print(\"LINEAR DIFFERENCE ALL: \")\n", + " print(np.abs(u_exact(t) - u))\n", + " tol = 1E-14\n", + " #assert diff < tol\n", + "\n", + " #u, t = solver(I, V, m, b, s, F, dt, T, 'quadratic')\n", + " u, t = solver_bwdamping(I, V, m, b, s, F, dt, T, 'quadratic')\n", + " diff = np.abs(u_exact(t) - u).max()\n", + " print(diff)\n", + " #assert diff < tol\n", + "\n", + "def test_mms():\n", + " \"\"\"Use method of manufactured solutions.\"\"\"\n", + " m = 4.; b = 1\n", + " w = 1.5\n", + " t = sym.Symbol('t')\n", + " u_exact = 3*sym.exp(-0.2*t)*sym.cos(1.2*t)\n", + " I = u_exact.subs(t, 0).evalf()\n", + " V = sym.diff(u_exact, t).subs(t, 0).evalf()\n", + " u_exact_py = sym.lambdify(t, u_exact, modules='numpy')\n", + " s = lambda u: u**3\n", + " dt = 0.2\n", + " T = 6\n", + " errors_linear = []\n", + " errors_quadratic = []\n", + " # Run grid refinements and compute exact error\n", + " for i in range(5):\n", + " F_formula = lhs_eq(t, m, b, s, u_exact, 'linear')\n", + " F = sym.lambdify(t, F_formula)\n", + " #u1, t1 = solver(I, V, m, b, s, F, dt, T, 'linear')\n", + " u1, t1 = solver_bwdamping(I, V, m, b, s,\n", + " F, dt, T, 'linear')\n", + " error = np.sqrt(np.sum((u_exact_py(t1) - u1)**2)*dt)\n", + " errors_linear.append((dt, error))\n", + "\n", + " F_formula = lhs_eq(t, m, b, s, u_exact, 'quadratic')\n", + " #print sym.latex(F_formula, mode='plain')\n", + " F = sym.lambdify(t, F_formula)\n", + " #u2, t2 = solver(I, V, m, b, s, F, dt, T, 'quadratic')\n", + " u2, t2 = solver_bwdamping(I, V, m, b, s,\n", + " F, dt, T, 'quadratic')\n", + " error = np.sqrt(np.sum((u_exact_py(t2) - u2)**2)*dt)\n", + " errors_quadratic.append((dt, error))\n", + " dt /= 2\n", + " # Estimate convergence rates\n", + " tol = 0.05\n", + " for errors in errors_linear, errors_quadratic:\n", + " for i in range(1, len(errors)):\n", + " dt, error = errors[i]\n", + " dt_1, error_1 = errors[i-1]\n", + " r = np.log(error/error_1)/np.log(dt/dt_1)\n", + " # check r for final simulation with (final and) smallest dt\n", + " # note that the method now is 1st order, i.e. r should\n", + " # approach 1.0\n", + " print(abs(r-1.0))\n", + " assert abs(r - 1.0) < tol\n", + "\n", + "def visualize(list_of_curves, legends, title='', filename='tmp'):\n", + " \"\"\"Plot list of curves: (u, t).\"\"\"\n", + " for u, t in list_of_curves:\n", + " plt.plot(t, u)\n", + " plt.hold('on')\n", + " plt.legend(legends)\n", + " plt.xlabel('t')\n", + " plt.ylabel('u')\n", + " plt.title(title)\n", + " plt.savefig(filename + '.png')\n", + " plt.savefig(filename + '.pdf')\n", + " plt.show()\n", + "\n", + "def main():\n", + " import argparse\n", + " parser = argparse.ArgumentParser()\n", + " parser.add_argument('--I', type=float, default=1.0)\n", + " parser.add_argument('--V', type=float, default=0.0)\n", + " parser.add_argument('--m', type=float, default=1.0)\n", + " parser.add_argument('--b', type=float, default=0.0)\n", + " parser.add_argument('--s', type=str, default='4*pi**2*u')\n", + " parser.add_argument('--F', type=str, default='0')\n", + " parser.add_argument('--dt', type=float, default=0.05)\n", + " parser.add_argument('--T', type=float, default=20)\n", + " parser.add_argument('--window_width', type=float, default=30.,\n", + " help='Number of periods in a window')\n", + " parser.add_argument('--damping', type=str, default='linear')\n", + " parser.add_argument('--savefig', action='store_true')\n", + " # Hack to allow --SCITOOLS options\n", + " # (scitools.std reads this argument at import)\n", + " parser.add_argument('--SCITOOLS_easyviz_backend',\n", + " default='matplotlib')\n", + " a = parser.parse_args()\n", + " from scitools import StringFunction\n", + " s = StringFunction(a.s, independent_variable='u')\n", + " F = StringFunction(a.F, independent_variable='t')\n", + " I, V, m, b, dt, T, window_width, savefig, damping = \\\n", + " a.I, a.V, a.m, a.b, a.dt, a.T, a.window_width, a.savefig, \\\n", + " a.damping\n", + "\n", + " # compute u by both methods and then visualize the difference\n", + " u, t = solver2(I, V, m, b, s, F, dt, T, damping)\n", + " u_bw, _ = solver_bwdamping(I, V, m, b, s, F, dt, T, damping)\n", + " u_diff = u - u_bw\n", + "\n", + " num_periods = plot_empirical_freq_and_amplitude(u_diff, t)\n", + " if num_periods <= 40:\n", + " plt.figure()\n", + " legends = ['1st-2nd order method',\n", + " '2nd order method',\n", + " '1st order method']\n", + " visualize([(u_diff, t), (u, t), (u_bw, t)], legends)\n", + " else:\n", + " visualize_front(u_diff, t, window_width, savefig)\n", + " #visualize_front_ascii(u_diff, t)\n", + " plt.show()\n", + "\n", + "if __name__ == '__main__':\n", + " #main()\n", + " #test_constant()\n", + " test_sinusoidal()\n", + " #test_mms()\n", + " #test_quadratic()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here is a comparison of standard method (2nd order) and backward differences for damping (1st order) for 10 (left) and 40 (right) time steps per period:\n", + "\n", + "\n", + "\n", + "\n", + "

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Filename: `vib_gen_bwdamping`.\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "## Exercise 3: Use the forward-backward scheme with quadratic damping\n", + "
\n", + "\n", + "We consider the generalized model with quadratic damping, expressed\n", + "as a system of two first-order equations as in the section [A staggered Euler-Cromer scheme for a generalized model](#vib:ode2:staggered):" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{align*}\n", + "u^{\\prime} &= v,\\\\ \n", + "v' &= \\frac{1}{m}\\left( F(t) - \\beta |v|v - s(u)\\right)\\thinspace .\n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "However, contrary to what is done in the section [A staggered Euler-Cromer scheme for a generalized model](#vib:ode2:staggered),\n", + "we want to apply the idea of a forward-backward discretization:\n", + "$u$ is marched forward by an one-sided Forward Euler scheme applied\n", + "to the first equation, and\n", + "thereafter $v$ can be marched forward by a Backward Euler scheme in the\n", + "second equation, see in the section [vib:model2x2:EulerCromer](#vib:model2x2:EulerCromer).\n", + "\n", + "Express the idea in operator notation and write out the\n", + "scheme. Unfortunately, the backward difference for the $v$ equation\n", + "creates a nonlinearity $|v^{n+1}|v^{n+1}$. To linearize this\n", + "nonlinearity, use the known value $v^n$ inside the absolute value\n", + "factor, i.e., $|v^{n+1}|v^{n+1}\\approx |v^n|v^{n+1}$. Show that the\n", + "resulting scheme is equivalent to the one in the section [A staggered Euler-Cromer scheme for a generalized model](#vib:ode2:staggered) for some time level $n\\geq 1$.\n", + "\n", + "What we learn from this exercise is that the first-order differences\n", + "and the linearization trick play together in \"the right way\" such that\n", + "the scheme is as good as when we (in the section [A staggered Euler-Cromer scheme for a generalized model](#vib:ode2:staggered))\n", + "carefully apply centered differences and a geometric mean on a\n", + "staggered mesh to achieve second-order accuracy.\n", + "There is a\n", + "difference in the handling of the initial conditions, though, as\n", + "explained at the end of the section [vib:model2x2:EulerCromer](#vib:model2x2:EulerCromer).\n", + "Filename: `vib_gen_bwdamping`.\n", + "\n", + "" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Devito", + "language": "python", + "name": "devito" + }, + "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.10.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}