duffing-stuffing.ipynb

@@ -26,9 +26,9 @@
     "c1 = 0.1\n",
     "c2 = 0.1\n",
     "c3 = 0.1\n",
-    "a1 = 0.12\n",
+    "a1 = 0.1\n",
-    "a2 = 0.05\n",
+    "a2 = 0.1\n",
-    "a3 = 0.08\n",
+    "a3 = 0.1\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 = (-1, -2.5), (1.5, -0.1)\n",
+    "x0, v0 = (0.5, -0.5), (0.5, -0.1)\n",
-    "tmax, t_trans = 200, 0\n",
+    "tmax, t_trans = 150, 0\n",
     "dt_per_period = 100"
    ]
   },
@@ -132,67 +132,49 @@
    "outputs": [],
    "source": [
     "%%capture\n",
-    "fig2, ((ax11, ax12, ax13), (ax21, ax22, ax23)) = plt.subplots(nrows=2, ncols=3)\n",
+    "fig2, ((ax11, ax12), (ax21, ax22)) = plt.subplots(nrows=2, ncols=2)\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_xlim([np.min(x1), np.max(x1)])\n",
-    "ax11.set_ylim([1.2*np.min(x2), 1.2*np.max(x2)])\n",
+    "ax11.set_ylim([np.min(x2), np.max(x2)])\n",
     "\n",
     "# Velocities\n",
-    "ax21.set_xlabel(r'$v_1$')\n",
+    "ax12.set_xlabel(r'$v_1$')\n",
-    "ax21.set_ylabel(r'$v_2$')\n",
+    "ax12.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_xlim([np.min(xd1), np.max(xd1)])\n",
-    "ax12.set_ylim([1.2*np.min(x1), 1.2*np.max(x1)])\n",
+    "ax12.set_ylim([np.min(xd2), np.max(xd2)])\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",
+    "ax21.set_xlabel(r'$\\mathbf{x}$')\n",
-    "ax13.set_ylabel(r'$\\mathbf{v}$')\n",
+    "ax21.set_ylabel(r'$\\mathbf{v}$')\n",
-    "ln13a, = ax13.plot([], [])\n",
+    "ln21a, = ax21.plot([], [])\n",
-    "ln13b, = ax13.plot([], [])\n",
+    "ln21b, = ax21.plot([], [])\n",
-    "ax13.set_xlim([1.2*min(np.min(x1), np.min(x2)), 1.2*max(np.max(x1), np.max(x2))])\n",
+    "ax21.set_xlim([min(np.min(x1), np.min(x2)), 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",
+    "ax21.set_ylim([min(np.min(xd1), np.min(xd2)), 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",
+    "ax22.set_xlabel(r'$t$')\n",
-    "ax23.set_ylabel(r'$\\mathbf{F}(t)$')\n",
+    "ax22.set_ylabel(r'$\\mathbf{F}(t)$')\n",
-    "ln23a, = ax23.plot([], [])\n",
+    "ln22a, = ax22.plot([], [])\n",
-    "ln23b, = ax23.plot([], [])\n",
+    "ln22b, = ax22.plot([], [])\n",
-    "ax23.set_xlim([t[0], t[-1]])\n",
+    "ax22.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",
+    "ax21.set_ylim([min(np.min(F1), np.min(F2)), 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(xd1[:i], xd2[:i])\n",
-    "    ln12.set_data(t[:i], x1[:i])\n",
+    "    ln21a.set_data(x1[:i], xd1[:i])\n",
-    "    ln22.set_data(t[:i], x2[:i])\n",
+    "    ln21b.set_data(x2[:i], xd2[:i])\n",
-    "    ln13a.set_data(x1[:i], xd1[:i])\n",
+    "    ln22a.set_data(t[:i], F1[:i])\n",
-    "    ln13b.set_data(x2[:i], xd2[:i])\n",
+    "    ln22b.set_data(t[:i], F2[: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)"
@@ -207,15 +189,6 @@
     "anim"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "anim.save('animations/coupled-duffing.mp4', writer='ffmpeg')"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,

duffing-stuffing2.ipynb

@@ -1,344 +0,0 @@
-{
- "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
-}

tsanim.py

@@ -1,42 +0,0 @@
-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