diff --git a/Duffing.ipynb b/Duffing.ipynb new file mode 100644 index 0000000..facce78 --- /dev/null +++ b/Duffing.ipynb @@ -0,0 +1,1046 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib notebook\n", + "import numpy as np\n", + "from scipy.integrate import odeint, quad\n", + "from scipy.optimize import brentq\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib import animation, rc\n", + "import seaborn as sbs\n", + "#rc('font', **{'family': 'serif', 'serif': ['Computer Modern'], 'size': 20})\n", + "#rc('text', usetex=True)\n", + "#rc('animation', html='jshtml')\n", + "#plt.rcParams[\"animation.html\"] = \"jshtml\"" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "def deriv(X, t, omega, gamma, alpha, delta, beta):\n", + " \"\"\"Return the derivatives dx/dt and d2x/dt2.\"\"\"\n", + " \n", + " #dVdx = lambda x: x**3 - x\n", + " x, xdot = X\n", + " xdotdot = gamma*np.cos(omega*t) - delta*xdot + alpha*x - beta*x*x*x\n", + " return xdot, xdotdot\n", + "\n", + "def solve_duffing(tmax, dt_per_period, t_trans, x0, v0, params):\n", + " \"\"\"Solve the Duffing equation for parameters gamma, delta, omega.\n", + " \n", + " Find the numerical solution to the Duffing equation using a suitable\n", + " time grid: tmax is the maximum time (s) to integrate to; t_trans is\n", + " the initial time period of transient behaviour until the solution\n", + " settles down (if it does) to some kind of periodic motion (these data\n", + " points are dropped) and dt_per_period is the number of time samples\n", + " (of duration dt) to include per period of the driving motion (frequency\n", + " omega).\n", + " \n", + " Returns the time grid, t (after t_trans), position, x, and velocity,\n", + " xdot, dt, and step, the number of array points per period of the driving\n", + " motion.\n", + " \n", + " \"\"\"\n", + " # Time point spacings and the time grid\n", + "\n", + " omega, gamma, alpha, delta, beta = params\n", + " period = 2*np.pi/omega\n", + " dt = 2*np.pi/omega / dt_per_period\n", + " step = int(period / dt)\n", + " t = np.arange(0, tmax, dt)\n", + " # Initial conditions: x, xdot\n", + " X0 = [x0, v0]\n", + " X = odeint(deriv, X0, t, args=params)\n", + " idx = int(t_trans / dt)\n", + " return t[idx:], X[idx:], dt, step" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the motion for a oscillator with initial position\n", + "# x0 and initially at rest.\n", + "x0, v0 = 1, 0.1\n", + "tmax, t_trans = 10000, 100\n", + "dt_per_period = 500\n", + "\n", + "omega = 1.0\n", + "gamma = 0.4\n", + "alpha = 1.0\n", + "delta = 0.1\n", + "beta = 1.0\n", + "\n", + "omega = 1.4\n", + "gamma, delta = 0.4, 0.1\n", + "\n", + "\n", + "params = (omega, gamma, alpha, delta, beta)\n", + "xs = []" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(1102945,)\n", + "(1102945,)\n" + ] + } + ], + "source": [ + "# Solve the equation of motion.\n", + "params = (omega, gamma, alpha, delta, 0.01)\n", + "t, X, dt, pstep = solve_duffing(tmax, dt_per_period, t_trans, x0, v0, params)\n", + "x, xdot = X.T\n", + "xs.append(x)\n", + "print(t.shape)\n", + "print(x.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# Solve the equation of motion.\n", + "#params = (omega, gamma, 1.2, delta, 0.1)\n", + "#t, X, dt, pstep = solve_duffing(tmax, dt_per_period, t_trans, x0, v0, params)\n", + "#x, xdot = X.T\n", + "#xs.append(x)\n", + "#print(t.shape)\n", + "#print(x.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(1102945,)\n", + "(1102945,)\n" + ] + } + ], + "source": [ + "# Solve the equation of motion.\n", + "params = (omega, gamma, alpha, delta, 1.0)\n", + "t, X, dt, pstep = solve_duffing(tmax, dt_per_period, t_trans, x0, v0, params)\n", + "x, xdot = X.T\n", + "xs.append(x)\n", + "print(t.shape)\n", + "print(x.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "def tsanim(tss):\n", + " \"Plot an array of time series data t, x = [ti], [xi].\"\n", + " nrows = len(tss)\n", + " fig, axes = plt.subplots(nrows=nrows, ncols=1)\n", + " tsmax = lambda ts: max(len(ts[0]), len(ts[1]))\n", + " ts_maxlen = max(map(tsmax, tss))\n", + " tinit = int(0.1*ts_maxlen)\n", + " print(\"Animate time series from t=%d with %d frames\" % (tinit, ts_maxlen))\n", + "\n", + " def tsplot(ax, ts):\n", + " t, x = ts\n", + " ax.set_ylabel(r'$x [\\mathrm{m}]$')\n", + " ax.set_ylim(np.min(x), np.max(x))\n", + " line, = ax.plot(t[:tinit], x[:tinit])\n", + " return line\n", + " \n", + " lines = list(map(lambda axts: tsplot(*axts), zip(axes, tss)))\n", + " axes[-1].set_xlabel(r'$t [\\mathrm{s}]$')\n", + "\n", + " def animate(i):\n", + " \"\"\"Update the image for iteration i of the Matplotlib animation.\"\"\"\n", + " for ax, line, ts in zip(axes, lines, tss):\n", + " t, x = ts\n", + " line.set_data(t[:i+1], x[:i+1])\n", + " ax.set_xlim(t_trans, t[i])\n", + " return\n", + "\n", + " #plt.tight_layout()\n", + " \n", + " #tmp\n", + " fig.suptitle(\"Duffing oscillator, beta = 0.01 (top) and 1.0 (bottom)\")\n", + " \n", + " anim = animation.FuncAnimation(fig, animate, frames=20000, interval=1)\n", + " return anim" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " fig.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", - " button.click(method_name, toolbar_event);\n", - " button.mouseover(tooltip, toolbar_mouse_event);\n", - " nav_element.append(button);\n", - " }\n", - "\n", - " // Add the status bar.\n", - " var status_bar = $('');\n", - " nav_element.append(status_bar);\n", - " this.message = status_bar[0];\n", - "\n", - " // Add the close button to the window.\n", - " var buttongrp = $('
');\n", - " var button = $('');\n", - " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", - " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", - " buttongrp.append(button);\n", - " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", - " titlebar.prepend(buttongrp);\n", - "}\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(el){\n", - " var fig = this\n", - " el.on(\"remove\", function(){\n", - "\tfig.close_ws(fig, {});\n", - " });\n", - "}\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(el){\n", - " // this is important to make the div 'focusable\n", - " el.attr('tabindex', 0)\n", - " // reach out to IPython and tell the keyboard manager to turn it's self\n", - " // off when our div gets focus\n", - "\n", - " // location in version 3\n", - " if (IPython.notebook.keyboard_manager) {\n", - " IPython.notebook.keyboard_manager.register_events(el);\n", - " }\n", - " else {\n", - " // location in version 2\n", - " IPython.keyboard_manager.register_events(el);\n", - " }\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._key_event_extra = function(event, name) {\n", - " var manager = IPython.notebook.keyboard_manager;\n", - " if (!manager)\n", - " manager = IPython.keyboard_manager;\n", - "\n", - " // Check for shift+enter\n", - " if (event.shiftKey && event.which == 13) {\n", - " this.canvas_div.blur();\n", - " event.shiftKey = false;\n", - " // Send a \"J\" for go to next cell\n", - " event.which = 74;\n", - " event.keyCode = 74;\n", - " manager.command_mode();\n", - " manager.handle_keydown(event);\n", - " }\n", - "}\n", - "\n", - "mpl.figure.prototype.handle_save = function(fig, msg) {\n", - " fig.ondownload(fig, null);\n", - "}\n", - "\n", - "\n", - "mpl.find_output_cell = function(html_output) {\n", - " // Return the cell and output element which can be found *uniquely* in the notebook.\n", - " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", - " // IPython event is triggered only after the cells have been serialised, which for\n", - " // our purposes (turning an active figure into a static one), is too late.\n", - " var cells = IPython.notebook.get_cells();\n", - " var ncells = cells.length;\n", - " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", - " data = data.data;\n", - " }\n", - " if (data['text/html'] == html_output) {\n", - " return [cell, data, j];\n", - " }\n", - " }\n", - " }\n", - " }\n", - "}\n", - "\n", - "// Register the function which deals with the matplotlib target/channel.\n", - "// The kernel may be null if the page has been refreshed.\n", - "if (IPython.notebook.kernel != null) {\n", - " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", - "}\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], "source": [ "gamma = 1\n", "alpha = 1\n", @@ -847,812 +57,53 @@ " np.linspace(0, 2, 500)\n", ")\n", "\n", - "plt.figure()\n", - "for beta in (-0.3, 0.0, 1, 4):\n", + "betas = (-0.3, 0.0, 10)\n", + "for i, beta in enumerate(betas):\n", + " fig = plt.figure(i)\n", " G = right(gamma, alpha, beta*beta_mul, delta, z, omega)\n", " plt.contour(omega, z, (z-G), [0])\n", " plt.xlabel(\"$\\omega$\")\n", - " plt.ylabel(\"$z/\\gamma$\")" + " plt.ylabel(\"$z/\\gamma$\")\n", + " fig.suptitle(\"Frequency response, beta=\" + str(beta))\n", + " fig.savefig('frequency-response-%d.png' % i)\n", + "\n", + "fig = plt.figure(len(betas))\n", + "for beta in betas:\n", + " G = right(gamma, alpha, beta*beta_mul, delta, z, omega)\n", + " plt.contour(omega, z, (z-G), [0])\n", + " plt.xlabel(\"$\\omega$\")\n", + " plt.ylabel(\"$z/\\gamma$\")\n", + "fig.suptitle(\"Frequency response, beta=\" + ', '.join(map(str, betas)))\n", + "fig.savefig('frequency-response-%d.png' % len(betas))" ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "application/javascript": [ - "/* Put everything inside the global mpl namespace */\n", - "window.mpl = {};\n", - "\n", - "\n", - "mpl.get_websocket_type = function() {\n", - " if (typeof(WebSocket) !== 'undefined') {\n", - " return WebSocket;\n", - " } else if (typeof(MozWebSocket) !== 'undefined') {\n", - " return MozWebSocket;\n", - " } else {\n", - " alert('Your browser does not have WebSocket support.' +\n", - " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", - " 'Firefox 4 and 5 are also supported but you ' +\n", - " 'have to enable WebSockets in about:config.');\n", - " };\n", - "}\n", - "\n", - "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", - " this.id = figure_id;\n", - "\n", - " this.ws = websocket;\n", - "\n", - " this.supports_binary = (this.ws.binaryType != undefined);\n", - "\n", - " if (!this.supports_binary) {\n", - " var warnings = document.getElementById(\"mpl-warnings\");\n", - " if (warnings) {\n", - " warnings.style.display = 'block';\n", - " warnings.textContent = (\n", - " \"This browser does not support binary websocket messages. \" +\n", - " \"Performance may be slow.\");\n", - " }\n", - " }\n", - "\n", - " this.imageObj = new Image();\n", - "\n", - " this.context = undefined;\n", - " this.message = undefined;\n", - " this.canvas = undefined;\n", - " this.rubberband_canvas = undefined;\n", - " this.rubberband_context = undefined;\n", - " this.format_dropdown = undefined;\n", - "\n", - " this.image_mode = 'full';\n", - "\n", - " this.root = $('
');\n", - " this._root_extra_style(this.root)\n", - " this.root.attr('style', 'display: inline-block');\n", - "\n", - " $(parent_element).append(this.root);\n", - "\n", - " this._init_header(this);\n", - " this._init_canvas(this);\n", - " this._init_toolbar(this);\n", - "\n", - " var fig = this;\n", - "\n", - " this.waiting = false;\n", - "\n", - " this.ws.onopen = function () {\n", - " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", - " fig.send_message(\"send_image_mode\", {});\n", - " if (mpl.ratio != 1) {\n", - " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", - " }\n", - " fig.send_message(\"refresh\", {});\n", - " }\n", - "\n", - " this.imageObj.onload = function() {\n", - " if (fig.image_mode == 'full') {\n", - " // Full images could contain transparency (where diff images\n", - " // almost always do), so we need to clear the canvas so that\n", - " // there is no ghosting.\n", - " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", - " }\n", - " fig.context.drawImage(fig.imageObj, 0, 0);\n", - " };\n", - "\n", - " this.imageObj.onunload = function() {\n", - " fig.ws.close();\n", - " }\n", - "\n", - " this.ws.onmessage = this._make_on_message_function(this);\n", - "\n", - " this.ondownload = ondownload;\n", - "}\n", - "\n", - "mpl.figure.prototype._init_header = function() {\n", - " var titlebar = $(\n", - " '
');\n", - " var titletext = $(\n", - " '
');\n", - " titlebar.append(titletext)\n", - " this.root.append(titlebar);\n", - " this.header = titletext[0];\n", - "}\n", - "\n", - "\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._init_canvas = function() {\n", - " var fig = this;\n", - "\n", - " var canvas_div = $('
');\n", - "\n", - " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", - "\n", - " function canvas_keyboard_event(event) {\n", - " return fig.key_event(event, event['data']);\n", - " }\n", - "\n", - " canvas_div.keydown('key_press', canvas_keyboard_event);\n", - " canvas_div.keyup('key_release', canvas_keyboard_event);\n", - " this.canvas_div = canvas_div\n", - " this._canvas_extra_style(canvas_div)\n", - " this.root.append(canvas_div);\n", - "\n", - " var canvas = $('');\n", - " canvas.addClass('mpl-canvas');\n", - " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", - "\n", - " this.canvas = canvas[0];\n", - " this.context = canvas[0].getContext(\"2d\");\n", - "\n", - " var backingStore = this.context.backingStorePixelRatio ||\n", - "\tthis.context.webkitBackingStorePixelRatio ||\n", - "\tthis.context.mozBackingStorePixelRatio ||\n", - "\tthis.context.msBackingStorePixelRatio ||\n", - "\tthis.context.oBackingStorePixelRatio ||\n", - "\tthis.context.backingStorePixelRatio || 1;\n", - "\n", - " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", - "\n", - " var rubberband = $('');\n", - " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", - "\n", - " var pass_mouse_events = true;\n", - "\n", - " canvas_div.resizable({\n", - " start: function(event, ui) {\n", - " pass_mouse_events = false;\n", - " },\n", - " resize: function(event, ui) {\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " stop: function(event, ui) {\n", - " pass_mouse_events = true;\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " });\n", - "\n", - " function mouse_event_fn(event) {\n", - " if (pass_mouse_events)\n", - " return fig.mouse_event(event, event['data']);\n", - " }\n", - "\n", - " rubberband.mousedown('button_press', mouse_event_fn);\n", - " rubberband.mouseup('button_release', mouse_event_fn);\n", - " // Throttle sequential mouse events to 1 every 20ms.\n", - " rubberband.mousemove('motion_notify', mouse_event_fn);\n", - "\n", - " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", - " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", - "\n", - " canvas_div.on(\"wheel\", function (event) {\n", - " event = event.originalEvent;\n", - " event['data'] = 'scroll'\n", - " if (event.deltaY < 0) {\n", - " event.step = 1;\n", - " } else {\n", - " event.step = -1;\n", - " }\n", - " mouse_event_fn(event);\n", - " });\n", - "\n", - " canvas_div.append(canvas);\n", - " canvas_div.append(rubberband);\n", - "\n", - " this.rubberband = rubberband;\n", - " this.rubberband_canvas = rubberband[0];\n", - " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", - " this.rubberband_context.strokeStyle = \"#000000\";\n", - "\n", - " this._resize_canvas = function(width, height) {\n", - " // Keep the size of the canvas, canvas container, and rubber band\n", - " // canvas in synch.\n", - " canvas_div.css('width', width)\n", - " canvas_div.css('height', height)\n", - "\n", - " canvas.attr('width', width * mpl.ratio);\n", - " canvas.attr('height', height * mpl.ratio);\n", - " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", - "\n", - " rubberband.attr('width', width);\n", - " rubberband.attr('height', height);\n", - " }\n", - "\n", - " // Set the figure to an initial 600x600px, this will subsequently be updated\n", - " // upon first draw.\n", - " this._resize_canvas(600, 600);\n", - "\n", - " // Disable right mouse context menu.\n", - " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", - " return false;\n", - " });\n", - "\n", - " function set_focus () {\n", - " canvas.focus();\n", - " canvas_div.focus();\n", - " }\n", - "\n", - " window.setTimeout(set_focus, 100);\n", - "}\n", - "\n", - "mpl.figure.prototype._init_toolbar = function() {\n", - " var fig = this;\n", - "\n", - " var nav_element = $('
')\n", - " nav_element.attr('style', 'width: 100%');\n", - " this.root.append(nav_element);\n", - "\n", - " // Define a callback function for later on.\n", - " function toolbar_event(event) {\n", - " return fig.toolbar_button_onclick(event['data']);\n", - " }\n", - " function toolbar_mouse_event(event) {\n", - " return fig.toolbar_button_onmouseover(event['data']);\n", - " }\n", - "\n", - " for(var toolbar_ind in mpl.toolbar_items) {\n", - " var name = mpl.toolbar_items[toolbar_ind][0];\n", - " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", - " var image = mpl.toolbar_items[toolbar_ind][2];\n", - " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", - "\n", - " if (!name) {\n", - " // put a spacer in here.\n", - " continue;\n", - " }\n", - " var button = $('');\n", - " button.click(method_name, toolbar_event);\n", - " button.mouseover(tooltip, toolbar_mouse_event);\n", - " nav_element.append(button);\n", - " }\n", - "\n", - " // Add the status bar.\n", - " var status_bar = $('');\n", - " nav_element.append(status_bar);\n", - " this.message = status_bar[0];\n", - "\n", - " // Add the close button to the window.\n", - " var buttongrp = $('
');\n", - " var button = $('');\n", - " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", - " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", - " buttongrp.append(button);\n", - " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", - " titlebar.prepend(buttongrp);\n", - "}\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(el){\n", - " var fig = this\n", - " el.on(\"remove\", function(){\n", - "\tfig.close_ws(fig, {});\n", - " });\n", - "}\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(el){\n", - " // this is important to make the div 'focusable\n", - " el.attr('tabindex', 0)\n", - " // reach out to IPython and tell the keyboard manager to turn it's self\n", - " // off when our div gets focus\n", - "\n", - " // location in version 3\n", - " if (IPython.notebook.keyboard_manager) {\n", - " IPython.notebook.keyboard_manager.register_events(el);\n", - " }\n", - " else {\n", - " // location in version 2\n", - " IPython.keyboard_manager.register_events(el);\n", - " }\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._key_event_extra = function(event, name) {\n", - " var manager = IPython.notebook.keyboard_manager;\n", - " if (!manager)\n", - " manager = IPython.keyboard_manager;\n", - "\n", - " // Check for shift+enter\n", - " if (event.shiftKey && event.which == 13) {\n", - " this.canvas_div.blur();\n", - " event.shiftKey = false;\n", - " // Send a \"J\" for go to next cell\n", - " event.which = 74;\n", - " event.keyCode = 74;\n", - " manager.command_mode();\n", - " manager.handle_keydown(event);\n", - " }\n", - "}\n", - "\n", - "mpl.figure.prototype.handle_save = function(fig, msg) {\n", - " fig.ondownload(fig, null);\n", - "}\n", - "\n", - "\n", - "mpl.find_output_cell = function(html_output) {\n", - " // Return the cell and output element which can be found *uniquely* in the notebook.\n", - " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", - " // IPython event is triggered only after the cells have been serialised, which for\n", - " // our purposes (turning an active figure into a static one), is too late.\n", - " var cells = IPython.notebook.get_cells();\n", - " var ncells = cells.length;\n", - " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", - " data = data.data;\n", - " }\n", - " if (data['text/html'] == html_output) {\n", - " return [cell, data, j];\n", - " }\n", - " }\n", - " }\n", - " }\n", - "}\n", - "\n", - "// Register the function which deals with the matplotlib target/channel.\n", - "// The kernel may be null if the page has been refreshed.\n", - "if (IPython.notebook.kernel != null) {\n", - " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", - "}\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Clipped\n", - "917\n", - "917\n" - ] - }, - { - "data": { - "application/javascript": [ - "/* Put everything inside the global mpl namespace */\n", - "window.mpl = {};\n", - "\n", - "\n", - "mpl.get_websocket_type = function() {\n", - " if (typeof(WebSocket) !== 'undefined') {\n", - " return WebSocket;\n", - " } else if (typeof(MozWebSocket) !== 'undefined') {\n", - " return MozWebSocket;\n", - " } else {\n", - " alert('Your browser does not have WebSocket support.' +\n", - " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", - " 'Firefox 4 and 5 are also supported but you ' +\n", - " 'have to enable WebSockets in about:config.');\n", - " };\n", - "}\n", - "\n", - "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", - " this.id = figure_id;\n", - "\n", - " this.ws = websocket;\n", - "\n", - " this.supports_binary = (this.ws.binaryType != undefined);\n", - "\n", - " if (!this.supports_binary) {\n", - " var warnings = document.getElementById(\"mpl-warnings\");\n", - " if (warnings) {\n", - " warnings.style.display = 'block';\n", - " warnings.textContent = (\n", - " \"This browser does not support binary websocket messages. \" +\n", - " \"Performance may be slow.\");\n", - " }\n", - " }\n", - "\n", - " this.imageObj = new Image();\n", - "\n", - " this.context = undefined;\n", - " this.message = undefined;\n", - " this.canvas = undefined;\n", - " this.rubberband_canvas = undefined;\n", - " this.rubberband_context = undefined;\n", - " this.format_dropdown = undefined;\n", - "\n", - " this.image_mode = 'full';\n", - "\n", - " this.root = $('
');\n", - " this._root_extra_style(this.root)\n", - " this.root.attr('style', 'display: inline-block');\n", - "\n", - " $(parent_element).append(this.root);\n", - "\n", - " this._init_header(this);\n", - " this._init_canvas(this);\n", - " this._init_toolbar(this);\n", - "\n", - " var fig = this;\n", - "\n", - " this.waiting = false;\n", - "\n", - " this.ws.onopen = function () {\n", - " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", - " fig.send_message(\"send_image_mode\", {});\n", - " if (mpl.ratio != 1) {\n", - " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", - " }\n", - " fig.send_message(\"refresh\", {});\n", - " }\n", - "\n", - " this.imageObj.onload = function() {\n", - " if (fig.image_mode == 'full') {\n", - " // Full images could contain transparency (where diff images\n", - " // almost always do), so we need to clear the canvas so that\n", - " // there is no ghosting.\n", - " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", - " }\n", - " fig.context.drawImage(fig.imageObj, 0, 0);\n", - " };\n", - "\n", - " this.imageObj.onunload = function() {\n", - " fig.ws.close();\n", - " }\n", - "\n", - " this.ws.onmessage = this._make_on_message_function(this);\n", - "\n", - " this.ondownload = ondownload;\n", - "}\n", - "\n", - "mpl.figure.prototype._init_header = function() {\n", - " var titlebar = $(\n", - " '
');\n", - " var titletext = $(\n", - " '
');\n", - " titlebar.append(titletext)\n", - " this.root.append(titlebar);\n", - " this.header = titletext[0];\n", - "}\n", - "\n", - "\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._init_canvas = function() {\n", - " var fig = this;\n", - "\n", - " var canvas_div = $('
');\n", - "\n", - " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", - "\n", - " function canvas_keyboard_event(event) {\n", - " return fig.key_event(event, event['data']);\n", - " }\n", - "\n", - " canvas_div.keydown('key_press', canvas_keyboard_event);\n", - " canvas_div.keyup('key_release', canvas_keyboard_event);\n", - " this.canvas_div = canvas_div\n", - " this._canvas_extra_style(canvas_div)\n", - " this.root.append(canvas_div);\n", - "\n", - " var canvas = $('');\n", - " canvas.addClass('mpl-canvas');\n", - " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", - "\n", - " this.canvas = canvas[0];\n", - " this.context = canvas[0].getContext(\"2d\");\n", - "\n", - " var backingStore = this.context.backingStorePixelRatio ||\n", - "\tthis.context.webkitBackingStorePixelRatio ||\n", - "\tthis.context.mozBackingStorePixelRatio ||\n", - "\tthis.context.msBackingStorePixelRatio ||\n", - "\tthis.context.oBackingStorePixelRatio ||\n", - "\tthis.context.backingStorePixelRatio || 1;\n", - "\n", - " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", - "\n", - " var rubberband = $('');\n", - " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", - "\n", - " var pass_mouse_events = true;\n", - "\n", - " canvas_div.resizable({\n", - " start: function(event, ui) {\n", - " pass_mouse_events = false;\n", - " },\n", - " resize: function(event, ui) {\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " stop: function(event, ui) {\n", - " pass_mouse_events = true;\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " });\n", - "\n", - " function mouse_event_fn(event) {\n", - " if (pass_mouse_events)\n", - " return fig.mouse_event(event, event['data']);\n", - " }\n", - "\n", - " rubberband.mousedown('button_press', mouse_event_fn);\n", - " rubberband.mouseup('button_release', mouse_event_fn);\n", - " // Throttle sequential mouse events to 1 every 20ms.\n", - " rubberband.mousemove('motion_notify', mouse_event_fn);\n", - "\n", - " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", - " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", - "\n", - " canvas_div.on(\"wheel\", function (event) {\n", - " event = event.originalEvent;\n", - " event['data'] = 'scroll'\n", - " if (event.deltaY < 0) {\n", - " event.step = 1;\n", - " } else {\n", - " event.step = -1;\n", - " }\n", - " mouse_event_fn(event);\n", - " });\n", - "\n", - " canvas_div.append(canvas);\n", - " canvas_div.append(rubberband);\n", - "\n", - " this.rubberband = rubberband;\n", - " this.rubberband_canvas = rubberband[0];\n", - " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", - " this.rubberband_context.strokeStyle = \"#000000\";\n", - "\n", - " this._resize_canvas = function(width, height) {\n", - " // Keep the size of the canvas, canvas container, and rubber band\n", - " // canvas in synch.\n", - " canvas_div.css('width', width)\n", - " canvas_div.css('height', height)\n", - "\n", - " canvas.attr('width', width * mpl.ratio);\n", - " canvas.attr('height', height * mpl.ratio);\n", - " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", - "\n", - " rubberband.attr('width', width);\n", - " rubberband.attr('height', height);\n", - " }\n", - "\n", - " // Set the figure to an initial 600x600px, this will subsequently be updated\n", - " // upon first draw.\n", - " this._resize_canvas(600, 600);\n", - "\n", - " // Disable right mouse context menu.\n", - " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", - " return false;\n", - " });\n", - "\n", - " function set_focus () {\n", - " canvas.focus();\n", - " canvas_div.focus();\n", - " }\n", - "\n", - " window.setTimeout(set_focus, 100);\n", - "}\n", - "\n", - "mpl.figure.prototype._init_toolbar = function() {\n", - " var fig = this;\n", - "\n", - " var nav_element = $('
')\n", - " nav_element.attr('style', 'width: 100%');\n", - " this.root.append(nav_element);\n", - "\n", - " // Define a callback function for later on.\n", - " function toolbar_event(event) {\n", - " return fig.toolbar_button_onclick(event['data']);\n", - " }\n", - " function toolbar_mouse_event(event) {\n", - " return fig.toolbar_button_onmouseover(event['data']);\n", - " }\n", - "\n", - " for(var toolbar_ind in mpl.toolbar_items) {\n", - " var name = mpl.toolbar_items[toolbar_ind][0];\n", - " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", - " var image = mpl.toolbar_items[toolbar_ind][2];\n", - " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", - "\n", - " if (!name) {\n", - " // put a spacer in here.\n", - " continue;\n", - " }\n", - " var button = $('');\n", - " button.click(method_name, toolbar_event);\n", - " button.mouseover(tooltip, toolbar_mouse_event);\n", - " nav_element.append(button);\n", - " }\n", - "\n", - " // Add the status bar.\n", - " var status_bar = $('');\n", - " nav_element.append(status_bar);\n", - " this.message = status_bar[0];\n", - "\n", - " // Add the close button to the window.\n", - " var buttongrp = $('
');\n", - " var button = $('');\n", - " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", - " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", - " buttongrp.append(button);\n", - " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", - " titlebar.prepend(buttongrp);\n", - "}\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(el){\n", - " var fig = this\n", - " el.on(\"remove\", function(){\n", - "\tfig.close_ws(fig, {});\n", - " });\n", - "}\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(el){\n", - " // this is important to make the div 'focusable\n", - " el.attr('tabindex', 0)\n", - " // reach out to IPython and tell the keyboard manager to turn it's self\n", - " // off when our div gets focus\n", - "\n", - " // location in version 3\n", - " if (IPython.notebook.keyboard_manager) {\n", - " IPython.notebook.keyboard_manager.register_events(el);\n", - " }\n", - " else {\n", - " // location in version 2\n", - " IPython.keyboard_manager.register_events(el);\n", - " }\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._key_event_extra = function(event, name) {\n", - " var manager = IPython.notebook.keyboard_manager;\n", - " if (!manager)\n", - " manager = IPython.keyboard_manager;\n", - "\n", - " // Check for shift+enter\n", - " if (event.shiftKey && event.which == 13) {\n", - " this.canvas_div.blur();\n", - " event.shiftKey = false;\n", - " // Send a \"J\" for go to next cell\n", - " event.which = 74;\n", - " event.keyCode = 74;\n", - " manager.command_mode();\n", - " manager.handle_keydown(event);\n", - " }\n", - "}\n", - "\n", - "mpl.figure.prototype.handle_save = function(fig, msg) {\n", - " fig.ondownload(fig, null);\n", - "}\n", - "\n", - "\n", - "mpl.find_output_cell = function(html_output) {\n", - " // Return the cell and output element which can be found *uniquely* in the notebook.\n", - " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", - " // IPython event is triggered only after the cells have been serialised, which for\n", - " // our purposes (turning an active figure into a static one), is too late.\n", - " var cells = IPython.notebook.get_cells();\n", - " var ncells = cells.length;\n", - " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", - " data = data.data;\n", - " }\n", - " if (data['text/html'] == html_output) {\n", - " return [cell, data, j];\n", - " }\n", - " }\n", - " }\n", - " }\n", - "}\n", - "\n", - "// Register the function which deals with the matplotlib target/channel.\n", - "// The kernel may be null if the page has been refreshed.\n", - "if (IPython.notebook.kernel != null) {\n", - " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", - "}\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/plain": [ - "[]" - ] - }, - "execution_count": 100, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "f_resonance = 8e3\n", "omega_resonance = 2*np.pi*f_resonance\n", diff --git a/duffing-stuffing.ipynb b/duffing-stuffing.ipynb index 2658d15..010ea70 100644 --- a/duffing-stuffing.ipynb +++ b/duffing-stuffing.ipynb @@ -26,9 +26,9 @@ "c1 = 0.1\n", "c2 = 0.1\n", "c3 = 0.1\n", - "a1 = 0.1\n", - "a2 = 0.1\n", - "a3 = 0.1\n", + "a1 = 0.12\n", + "a2 = 0.05\n", + "a3 = 0.08\n", "k1 = 0.05\n", "k2 = 0.21\n", "k3 = 0.11\n", @@ -85,8 +85,8 @@ "outputs": [], "source": [ "# Set up the motion for a oscillator with initial positions and velocities\n", - "x0, v0 = (0.5, -0.5), (0.5, -0.1)\n", - "tmax, t_trans = 150, 0\n", + "x0, v0 = (-1, -2.5), (1.5, -0.1)\n", + "tmax, t_trans = 200, 0\n", "dt_per_period = 100" ] }, @@ -132,49 +132,67 @@ "outputs": [], "source": [ "%%capture\n", - "fig2, ((ax11, ax12), (ax21, ax22)) = plt.subplots(nrows=2, ncols=2)\n", + "fig2, ((ax11, ax12, ax13), (ax21, ax22, ax23)) = plt.subplots(nrows=2, ncols=3)\n", + "\n", + "#fig2.suptitle(\"Duffing\")\n", "\n", "# Positions\n", "ax11.set_xlabel(r'$x_1$')\n", "ax11.set_ylabel(r'$x_2$')\n", "ln11, = ax11.plot([], [])\n", - "ax11.set_xlim([np.min(x1), np.max(x1)])\n", - "ax11.set_ylim([np.min(x2), np.max(x2)])\n", + "ax11.set_xlim([1.2*np.min(x1), 1.2*np.max(x1)])\n", + "ax11.set_ylim([1.2*np.min(x2), 1.2*np.max(x2)])\n", "\n", "# Velocities\n", - "ax12.set_xlabel(r'$v_1$')\n", - "ax12.set_ylabel(r'$v_2$')\n", + "ax21.set_xlabel(r'$v_1$')\n", + "ax21.set_ylabel(r'$v_2$')\n", + "ln21, = ax21.plot([], [])\n", + "ax21.set_xlim([1.2*np.min(xd1), 1.2*np.max(xd1)])\n", + "ax21.set_ylim([1.2*np.min(xd2), 1.2*np.max(xd2)])\n", + "\n", + "# x1(t)\n", + "ax12.set_xlabel(r'$t$')\n", + "ax12.set_ylabel(r'$x_1$')\n", "ln12, = ax12.plot([], [])\n", - "ax12.set_xlim([np.min(xd1), np.max(xd1)])\n", - "ax12.set_ylim([np.min(xd2), np.max(xd2)])\n", + "ax12.set_xlim([t[0], t[-1]])\n", + "ax12.set_ylim([1.2*np.min(x1), 1.2*np.max(x1)])\n", + "\n", + "# x2(t)\n", + "ax22.set_xlabel(r'$t$')\n", + "ax22.set_ylabel(r'$x_2$')\n", + "ln22, = ax22.plot([], [])\n", + "ax22.set_xlim([t[0], t[-1]])\n", + "ax22.set_ylim([1.2*np.min(x2), 1.2*np.max(x2)])\n", "\n", "# Phase spaces\n", - "ax21.set_xlabel(r'$\\mathbf{x}$')\n", - "ax21.set_ylabel(r'$\\mathbf{v}$')\n", - "ln21a, = ax21.plot([], [])\n", - "ln21b, = ax21.plot([], [])\n", - "ax21.set_xlim([min(np.min(x1), np.min(x2)), max(np.max(x1), np.max(x2))])\n", - "ax21.set_ylim([min(np.min(xd1), np.min(xd2)), max(np.max(xd1), np.max(xd2))])\n", + "ax13.set_xlabel(r'$\\mathbf{x}$')\n", + "ax13.set_ylabel(r'$\\mathbf{v}$')\n", + "ln13a, = ax13.plot([], [])\n", + "ln13b, = ax13.plot([], [])\n", + "ax13.set_xlim([1.2*min(np.min(x1), np.min(x2)), 1.2*max(np.max(x1), np.max(x2))])\n", + "ax13.set_ylim([1.2*min(np.min(xd1), np.min(xd2)), 1.2*max(np.max(xd1), np.max(xd2))])\n", "\n", "# F(t)\n", "F1 = f1(t)\n", "F2 = f2(t)\n", - "ax22.set_xlabel(r'$t$')\n", - "ax22.set_ylabel(r'$\\mathbf{F}(t)$')\n", - "ln22a, = ax22.plot([], [])\n", - "ln22b, = ax22.plot([], [])\n", - "ax22.set_xlim([t[0], t[-1]])\n", - "ax21.set_ylim([min(np.min(F1), np.min(F2)), max(np.max(F1), np.max(F2))])\n", + "ax23.set_xlabel(r'$t$')\n", + "ax23.set_ylabel(r'$\\mathbf{F}(t)$')\n", + "ln23a, = ax23.plot([], [])\n", + "ln23b, = ax23.plot([], [])\n", + "ax23.set_xlim([t[0], t[-1]])\n", + "ax23.set_ylim([1.2*min(np.min(F1), np.min(F2)), 1.2*max(np.max(F1), np.max(F2))])\n", "\n", "plt.tight_layout()\n", "\n", "def animate(i):\n", " ln11.set_data(x1[:i], x2[:i])\n", - " ln12.set_data(xd1[:i], xd2[:i])\n", - " ln21a.set_data(x1[:i], xd1[:i])\n", - " ln21b.set_data(x2[:i], xd2[:i])\n", - " ln22a.set_data(t[:i], F1[:i])\n", - " ln22b.set_data(t[:i], F2[:i])\n", + " ln21.set_data(xd1[:i], xd2[:i])\n", + " ln12.set_data(t[:i], x1[:i])\n", + " ln22.set_data(t[:i], x2[:i])\n", + " ln13a.set_data(x1[:i], xd1[:i])\n", + " ln13b.set_data(x2[:i], xd2[:i])\n", + " ln23a.set_data(t[:i], F1[:i])\n", + " ln23b.set_data(t[:i], F2[:i])\n", " \n", "\n", "anim = animation.FuncAnimation(fig2, animate, frames=len(t), interval=100, repeat_delay=1000)" @@ -189,6 +207,15 @@ "anim" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "anim.save('animations/coupled-duffing.mp4', writer='ffmpeg')" + ] + }, { "cell_type": "code", "execution_count": null, diff --git a/duffing-stuffing2.ipynb b/duffing-stuffing2.ipynb new file mode 100644 index 0000000..d1337ae --- /dev/null +++ b/duffing-stuffing2.ipynb @@ -0,0 +1,344 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib notebook\n", + "\n", + "import numpy as np\n", + "from scipy.integrate import odeint\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib import animation\n", + "plt.rcParams[\"animation.html\"] = \"jshtml\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def deriv(X, t, gamma, delta, omega):\n", + " \"\"\"Return the derivatives dx/dt and d2x/dt2.\"\"\"\n", + " \n", + " x, xdot = X\n", + " xdotdot = -dVdx(x) -delta * xdot + gamma * np.cos(omega*t)\n", + " return xdot, xdotdot\n", + "\n", + "def solve_duffing(tmax, dt_per_period, t_trans, x0, v0, gamma, delta, omega):\n", + " \"\"\"Solve the Duffing equation for parameters gamma, delta, omega.\n", + " \n", + " Find the numerical solution to the Duffing equation using a suitable\n", + " time grid: tmax is the maximum time (s) to integrate to; t_trans is\n", + " the initial time period of transient behaviour until the solution\n", + " settles down (if it does) to some kind of periodic motion (these data\n", + " points are dropped) and dt_per_period is the number of time samples\n", + " (of duration dt) to include per period of the driving motion (frequency\n", + " omega).\n", + " \n", + " Returns the time grid, t (after t_trans), position, x, and velocity,\n", + " xdot, dt, and step, the number of array points per period of the driving\n", + " motion.\n", + " \n", + " \"\"\"\n", + " # Time point spacings and the time grid\n", + "\n", + " period = 2*np.pi/omega\n", + " dt = 2*np.pi/omega / dt_per_period\n", + " step = int(period / dt)\n", + " t = np.arange(0, tmax, dt)\n", + " # Initial conditions: x, xdot\n", + " X0 = [x0, v0]\n", + " X = odeint(deriv, X0, t, args=(gamma, delta, omega))\n", + " idx = int(t_trans / dt)\n", + " return t[idx:], X[idx:], dt, step" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the motion for a oscillator with initial position\n", + "# x0 and initially at rest.\n", + "x0, v0 = 0, 0\n", + "tmax, t_trans = 18000, 300\n", + "omega = 1.4\n", + "gamma, delta = 0.4, 0.1\n", + "dt_per_period = 100" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Solve the equation of motion.\n", + "delta = 0.0\n", + "t, X, dt, pstep = solve_duffing(tmax, dt_per_period, t_trans, x0, v0, gamma, delta, omega)\n", + "x1, x1dot = X.T" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Solve the equation of motion.\n", + "delta = 0.1\n", + "t, X, dt, pstep = solve_duffing(tmax, dt_per_period, t_trans, x0, v0, gamma, delta, omega)\n", + "x2, x2dot = X.T" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m1 = 1\n", + "m2 = 1\n", + "c1 = 0.1\n", + "c2 = 0.1\n", + "c3 = 0.1\n", + "a1 = 0.12\n", + "a2 = 0.05\n", + "a3 = 0.08\n", + "k1 = 0.05\n", + "k2 = 0.21\n", + "k3 = 0.11\n", + "C1 = 0.15\n", + "C2 = 0.30\n", + "omega = 0.1\n", + "\n", + "f1 = lambda t: C1*np.cos(omega*t)\n", + "f2 = lambda t: C2*np.cos(omega*t)\n", + "\n", + "def deriv(X, t):\n", + " \"\"\"Return the derivatives dx/dt and d2x/dt2.\"\"\"\n", + " x1, x2, xd1, xd2 = X\n", + " F1 = f1(t)\n", + " F2 = f2(t)\n", + " xdd1 = F1 - xd1*(c1 + c2) + xd2*c2 - x1*(k1 + k2) + x2*k2 - a1*x1**3 + a2*(x2 - x1)**3\n", + " xdd2 = F2 - xd2*(c2 + c3) + xd1*c1 - x2*(k2 + k3) + x1*k2 - a3*x2**3 + a2*(x2 - x1)**3\n", + " return xd1, xd2, xdd1/m1, xdd2/m2\n", + "\n", + "\n", + "def solve_duffing(tmax, dt_per_period, t_trans, x0, v0):\n", + " \"\"\"Solve the Duffing equation with the standard odeint solver.\n", + " \n", + " Find the numerical solution to the Duffing equation using a suitable\n", + " time grid: tmax is the maximum time (s) to integrate to; t_trans is\n", + " the initial time period of transient behaviour until the solution\n", + " settles down (if it does) to some kind of periodic motion (these data\n", + " points are dropped) and dt_per_period is the number of time samples\n", + " (of duration dt) to include per period of the driving motion (frequency\n", + " omega).\n", + " \n", + " Returns the time grid, t (after t_trans), position, x, and velocity,\n", + " xdot, dt, and step, the number of array points per period of the driving\n", + " motion.\n", + " \n", + " \"\"\"\n", + " # Time point spacings and the time grid\n", + "\n", + " period = 2*np.pi/omega\n", + " dt = 2*np.pi/omega / dt_per_period\n", + " step = int(period / dt)\n", + " t = np.arange(0, tmax, dt)\n", + " # Initial conditions: x, xdot\n", + " X0 = x0 + v0\n", + " X = odeint(deriv, X0, t)\n", + " idx = int(t_trans / dt)\n", + " return t[idx:], X[idx:], dt, step" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the motion for a oscillator with initial positions and velocities\n", + "x0, v0 = (-1, -2.5), (1.5, -0.1)\n", + "tmax, t_trans = 200, 0\n", + "dt_per_period = 100" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Solve the equation of motion.\n", + "t, X, dt, pstep = solve_duffing(tmax, dt_per_period, t_trans, x0, v0)\n", + "print(\"# samples: %d\" % len(t))\n", + "print(\"dt: %.4f\" % dt)\n", + "print(\"Steps per period: %d\" % pstep)\n", + "x1, x2, xd1, xd2 = X.T" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig1, ((ax11, ax12), (ax21, ax22)) = plt.subplots(nrows=2, ncols=2)\n", + "\n", + "ax11.plot(t, x1)\n", + "ax11.set_xlabel(r'$x_1$')\n", + "ax12.plot(t, x2)\n", + "ax12.set_xlabel(r'$x_2$')\n", + "ax21.plot(t, xd1)\n", + "ax21.set_xlabel(r'$v_1$')\n", + "ax22.plot(t, xd2)\n", + "ax22.set_xlabel(r'$v_2$')\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "%%capture\n", + "fig2, ((ax11, ax12, ax13), (ax21, ax22, ax23)) = plt.subplots(nrows=2, ncols=3)\n", + "\n", + "#fig2.suptitle(\"Duffing\")\n", + "\n", + "# Positions\n", + "ax11.set_xlabel(r'$x_1$')\n", + "ax11.set_ylabel(r'$x_2$')\n", + "ln11, = ax11.plot([], [])\n", + "ax11.set_xlim([1.2*np.min(x1), 1.2*np.max(x1)])\n", + "ax11.set_ylim([1.2*np.min(x2), 1.2*np.max(x2)])\n", + "\n", + "# Velocities\n", + "ax21.set_xlabel(r'$v_1$')\n", + "ax21.set_ylabel(r'$v_2$')\n", + "ln21, = ax21.plot([], [])\n", + "ax21.set_xlim([1.2*np.min(xd1), 1.2*np.max(xd1)])\n", + "ax21.set_ylim([1.2*np.min(xd2), 1.2*np.max(xd2)])\n", + "\n", + "# x1(t)\n", + "ax12.set_xlabel(r'$t$')\n", + "ax12.set_ylabel(r'$x_1$')\n", + "ln12, = ax12.plot([], [])\n", + "ax12.set_xlim([t[0], t[-1]])\n", + "ax12.set_ylim([1.2*np.min(x1), 1.2*np.max(x1)])\n", + "\n", + "# x2(t)\n", + "ax22.set_xlabel(r'$t$')\n", + "ax22.set_ylabel(r'$x_2$')\n", + "ln22, = ax22.plot([], [])\n", + "ax22.set_xlim([t[0], t[-1]])\n", + "ax22.set_ylim([1.2*np.min(x2), 1.2*np.max(x2)])\n", + "\n", + "# Phase spaces\n", + "ax13.set_xlabel(r'$\\mathbf{x}$')\n", + "ax13.set_ylabel(r'$\\mathbf{v}$')\n", + "ln13a, = ax13.plot([], [])\n", + "ln13b, = ax13.plot([], [])\n", + "ax13.set_xlim([1.2*min(np.min(x1), np.min(x2)), 1.2*max(np.max(x1), np.max(x2))])\n", + "ax13.set_ylim([1.2*min(np.min(xd1), np.min(xd2)), 1.2*max(np.max(xd1), np.max(xd2))])\n", + "\n", + "# F(t)\n", + "F1 = f1(t)\n", + "F2 = f2(t)\n", + "ax23.set_xlabel(r'$t$')\n", + "ax23.set_ylabel(r'$\\mathbf{F}(t)$')\n", + "ln23a, = ax23.plot([], [])\n", + "ln23b, = ax23.plot([], [])\n", + "ax23.set_xlim([t[0], t[-1]])\n", + "ax23.set_ylim([1.2*min(np.min(F1), np.min(F2)), 1.2*max(np.max(F1), np.max(F2))])\n", + "\n", + "plt.tight_layout()\n", + "\n", + "def animate(i):\n", + " ln11.set_data(x1[:i], x2[:i])\n", + " ln21.set_data(xd1[:i], xd2[:i])\n", + " ln12.set_data(t[:i], x1[:i])\n", + " ln22.set_data(t[:i], x2[:i])\n", + " ln13a.set_data(x1[:i], xd1[:i])\n", + " ln13b.set_data(x2[:i], xd2[:i])\n", + " ln23a.set_data(t[:i], F1[:i])\n", + " ln23b.set_data(t[:i], F2[:i])\n", + " \n", + "\n", + "anim = animation.FuncAnimation(fig2, animate, frames=len(t), interval=100, repeat_delay=1000)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "anim" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "anim.save('animations/coupled-duffing.mp4', writer='ffmpeg')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/frequency-response/frequency-response-0.png b/frequency-response/frequency-response-0.png new file mode 100644 index 0000000..d5771bd Binary files /dev/null and b/frequency-response/frequency-response-0.png differ diff --git a/frequency-response/frequency-response-1.png b/frequency-response/frequency-response-1.png new file mode 100644 index 0000000..8d45c15 Binary files /dev/null and b/frequency-response/frequency-response-1.png differ diff --git a/frequency-response/frequency-response-2.png b/frequency-response/frequency-response-2.png new file mode 100644 index 0000000..2830501 Binary files /dev/null and b/frequency-response/frequency-response-2.png differ diff --git a/frequency-response/frequency-response-3.png b/frequency-response/frequency-response-3.png new file mode 100644 index 0000000..2fc4d18 Binary files /dev/null and b/frequency-response/frequency-response-3.png differ diff --git a/tsanim.py b/tsanim.py new file mode 100644 index 0000000..14f8a89 --- /dev/null +++ b/tsanim.py @@ -0,0 +1,42 @@ +import numpy as np +from scipy.integrate import odeint, quad +from scipy.optimize import brentq +import matplotlib.pyplot as plt +from matplotlib import animation, rc +import seaborn as sbs +#rc('font', **{'family': 'serif', 'serif': ['Computer Modern'], 'size': 20}) +#rc('text', usetex=True) +#rc('animation', html='jshtml') +#plt.rcParams["animation.html"] = "jshtml" + + +def tsanim(tss): + "Plot an array of time series data t, x = [ti], [xi]." + nrows = len(tss) + fig, axes = plt.subplots(nrows=nrows, ncols=1) + tsmax = lambda ts: max(len(ts[0]), len(ts[1])) + ts_maxlen = max(map(tsmax, tss)) + tinit = int(0.1*ts_maxlen) + print("Animate time series from t=%d with %d frames" % (tinit, ts_maxlen)) + + def tsplot(ax, ts): + t, x = ts + ax.set_ylabel(r'$x$') + ax.set_ylim(np.min(x), np.max(x)) + line, = ax.plot(t[:tinit], x[:tinit]) + return line + + lines = list(map(lambda axts: tsplot(*axts), zip(axes, tss))) + axes[-1].set_xlabel(r'$t$') + + def animate(i): + """Update the image for iteration i of the Matplotlib animation.""" + for ax, line, ts in zip(axes, lines, tss): + t, x = ts + line.set_data(t[:i+1], x[:i+1]) + ax.set_xlim(0, t[i]) + return + + plt.tight_layout() + anim = animation.FuncAnimation(fig, animate, frames=ts_maxlen, interval=10) + return anim