{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "In this blog post, we will look at the vibrating string equation and solve it using two different mathematical methods.\n", "First, we will use a time-domain finite difference method. In a second step, we will use a frequency domain approach (i.e. a Fourier method) to solve the same problem. \n", "\n", "I wrote this post while going through the excellent MOOC Fundamentals of Waves and Vibrations. Feel free to check out what I wrote about the experience [here](http://flothesof.github.io/anki-MOOC-fundamentals-waves-vibrations.html)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# The vibrating string equation " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's start by introducing the [vibrating string](https://en.wikipedia.org/wiki/String_vibration) equation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The equation of the transverse displacement of a vibrating string $y(x,t)$ is, $\\forall x \\in [0, L]$, $\\forall t \\geq 0$:\n", "\n", "$$\n", "m \\dfrac{\\partial^2 y}{\\partial t^2} - T \\dfrac{\\partial^2 y}{\\partial x^2} = 0\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will solve this with so-called Dirichlet boundary conditions, meaning that the string has a zero displacement at its ends, i.e. $\\forall t \\geq 0 \\; y(0, t) = y(L, t) = 0$. Initial conditions will be as follows:\n", "\n", "- $\\forall x \\in [0, L], \\quad y(x, 0) = y_0(x)$\n", "- $\\forall x \\in [0, L], \\quad \\dot{y}(x, 0) = 0$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Time-domain solution " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's implement a time-domain solution known as finite differences. The partial differential equation above, once discretized on a 1D grid, becomes:\n", "\n", "$$\n", "m \\dfrac{y_i^{n+1} - 2 y_i^{n} + y_i^{n-1}}{\\Delta t^2} - T \\dfrac{y_{i+1}^{n} - 2 y_{i}^{n} + y_{i-1}^{n}}{\\Delta x^2} = 0\n", "$$\n", "\n", "After reordering the terms, we obtain the following *propagation step*:\n", "\n", "$$\n", "\\forall i \\in [1, N-1] \\quad y_i^{n+1} = 2 y_i^{n} - y_i^{n-1} + \\dfrac{T}{m} \\dfrac{\\Delta t^2}{\\Delta x^2} \\left ( y_{i+1}^{n} - 2 y_{i}^{n} + y_{i-1}^{n}\\right )\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can setup a simple implementation of the solution as follows:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "class TimeDomainSolver:\n", " \n", " def __init__(self, L, nx, T, m, y0):\n", " \"\"\"Setup of the time domain solver.\"\"\"\n", " self.x = np.linspace(0, L, num=nx)\n", " self.dx = self.x[1] - self.x[0]\n", " c = np.sqrt(T/m)\n", " self.dt = self.dx / c\n", " self.uprev = y0(self.x)\n", " self.ucurr = y0(self.x)\n", " self.alpha = T / m * self.dt ** 2 / self.dx **2\n", " \n", " \n", " def step(self):\n", " \"\"\"Steps the system forward in time.\"\"\"\n", " unext = np.zeros_like(self.ucurr)\n", " unext[1: -1] = 2 * self.ucurr[1:-1] - self.uprev[1:-1] + \\\n", " self.alpha * (self.ucurr[2:] - 2 * self.ucurr[1:-1] + self.ucurr[:-2])\n", " self.ucurr, self.uprev = unext, self.ucurr\n", " \n", " def steps(self, n):\n", " for _ in range(n):\n", " self.step()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To test our implementation, we will consider a triangular initial string with the properties of a low E string of a guitar (such as those that you can find on the back of electric guitar strings, as [here](http://www.daddario.com/DAstringtensionguide.Page?sid=91901025-4136-404a-8c26-e020667e8b5e))." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "L = 0.6477 # m\n", "T = 7.93 * 9.81 # kg * m.s^-2\n", "f0 = 82.4 # Hz\n", "c = 2 * L * f0\n", "m = T / c**2\n", "\n", "triangle = lambda x: (x * 2 / L) * (x < L/2) + (-x * 2 / L + 2) * (x >= L/2) \n", "\n", "td_solver = TimeDomainSolver(L, 200, T, m, triangle)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's visualize the initial position of the string:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "application/javascript": [ "\n", "(function(root) {\n", " function now() {\n", " return new Date();\n", " }\n", "\n", " var force = true;\n", "\n", " if (typeof (root._bokeh_onload_callbacks) === \"undefined\" || force === true) {\n", " root._bokeh_onload_callbacks = [];\n", " root._bokeh_is_loading = undefined;\n", " }\n", "\n", " var JS_MIME_TYPE = 'application/javascript';\n", " var HTML_MIME_TYPE = 'text/html';\n", " var EXEC_MIME_TYPE = 'application/vnd.bokehjs_exec.v0+json';\n", " var CLASS_NAME = 'output_bokeh rendered_html';\n", "\n", " /**\n", " * Render data to the DOM node\n", " */\n", " function render(props, node) {\n", " var script = document.createElement(\"script\");\n", " node.appendChild(script);\n", " }\n", "\n", " /**\n", " * Handle when an output is cleared or removed\n", " */\n", " function handleClearOutput(event, handle) {\n", " var cell = handle.cell;\n", "\n", " var id = cell.output_area._bokeh_element_id;\n", " var server_id = cell.output_area._bokeh_server_id;\n", " // Clean up Bokeh references\n", " if (id != null && id in Bokeh.index) {\n", " Bokeh.index[id].model.document.clear();\n", " delete Bokeh.index[id];\n", " }\n", "\n", " if (server_id !== undefined) {\n", " // Clean up Bokeh references\n", " var cmd = \"from bokeh.io.state import curstate; print(curstate().uuid_to_server['\" + server_id + \"'].get_sessions()[0].document.roots[0]._id)\";\n", " cell.notebook.kernel.execute(cmd, {\n", " iopub: {\n", " output: function(msg) {\n", " var id = msg.content.text.trim();\n", " if (id in Bokeh.index) {\n", " Bokeh.index[id].model.document.clear();\n", " delete Bokeh.index[id];\n", " }\n", " }\n", " }\n", " });\n", " // Destroy server and session\n", " var cmd = \"import bokeh.io.notebook as ion; ion.destroy_server('\" + server_id + \"')\";\n", " cell.notebook.kernel.execute(cmd);\n", " }\n", " }\n", "\n", " /**\n", " * Handle when a new output is added\n", " */\n", " function handleAddOutput(event, handle) {\n", " var output_area = handle.output_area;\n", " var output = handle.output;\n", "\n", " // limit handleAddOutput to display_data with EXEC_MIME_TYPE content only\n", " if ((output.output_type != \"display_data\") || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n", " return\n", " }\n", "\n", " var toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n", "\n", " if (output.metadata[EXEC_MIME_TYPE][\"id\"] !== undefined) {\n", " toinsert[toinsert.length - 1].firstChild.textContent = output.data[JS_MIME_TYPE];\n", " // store reference to embed id on output_area\n", " output_area._bokeh_element_id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n", " }\n", " if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n", " var bk_div = document.createElement(\"div\");\n", " bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n", " var script_attrs = bk_div.children[0].attributes;\n", " for (var i = 0; i < script_attrs.length; i++) {\n", " toinsert[toinsert.length - 1].firstChild.setAttribute(script_attrs[i].name, script_attrs[i].value);\n", " }\n", " // store reference to server id on output_area\n", " output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n", " }\n", " }\n", "\n", " function register_renderer(events, OutputArea) {\n", "\n", " function append_mime(data, metadata, element) {\n", " // create a DOM node to render to\n", " var toinsert = this.create_output_subarea(\n", " metadata,\n", " CLASS_NAME,\n", " EXEC_MIME_TYPE\n", " );\n", " this.keyboard_manager.register_events(toinsert);\n", " // Render to node\n", " var props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n", " render(props, toinsert[toinsert.length - 1]);\n", " element.append(toinsert);\n", " return toinsert\n", " }\n", "\n", " /* Handle when an output is cleared or removed */\n", " events.on('clear_output.CodeCell', handleClearOutput);\n", " events.on('delete.Cell', handleClearOutput);\n", "\n", " /* Handle when a new output is added */\n", " events.on('output_added.OutputArea', handleAddOutput);\n", "\n", " /**\n", " * Register the mime type and append_mime function with output_area\n", " */\n", " OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n", " /* Is output safe? */\n", " safe: true,\n", " /* Index of renderer in `output_area.display_order` */\n", " index: 0\n", " });\n", " }\n", "\n", " // register the mime type if in Jupyter Notebook environment and previously unregistered\n", " if (root.Jupyter !== undefined) {\n", " var events = require('base/js/events');\n", " var OutputArea = require('notebook/js/outputarea').OutputArea;\n", "\n", " if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n", " register_renderer(events, OutputArea);\n", " }\n", " }\n", "\n", " \n", " if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n", " root._bokeh_timeout = Date.now() + 5000;\n", " root._bokeh_failed_load = false;\n", " }\n", "\n", " var NB_LOAD_WARNING = {'data': {'text/html':\n", " \"
\\n\"+\n", " \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n", " \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n", " \"
\\n\"+\n", " \"\\n\"+\n",
" \"from bokeh.resources import INLINE\\n\"+\n",
" \"output_notebook(resources=INLINE)\\n\"+\n",
" \"
\\n\"+\n",
" \"\"),e=0;e<7;e++)n.push(' | '+p(t,e,!0)+\" | \");return\"
---|