{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Recently, I've come across the term Verlet integration several times: an [interesting numerical exploration of a moon explosion by Jason Cole](https://jasmcole.com/2017/09/20/the-moon-blew-up-without-warning-and-for-no-apparent-reason/) and the [wonderful Lennard-Jones molecular dynamics simulator by Daniel Schroeder](http://physics.weber.edu/schroeder/md/). I've been interested in numerical simulation for a long time (and wrote about it [here](http://flothesof.github.io/harmonic-oscillator-three-methods-solution.html) regarding the harmonic oscillator as well as [here](http://flothesof.github.io/charged-particle-trajectories-E-and-B-fields.html) regarding the movement of charged particles in E and B fields) so I think it's time to explore what exactly this numerical ODE integration technique is." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Verlet integration " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "According to [Wikipedia](https://en.wikipedia.org/wiki/Verlet_integration), Verlet integration is named after [Loup Verlet](https://en.wikipedia.org/wiki/Loup_Verlet), a French physicist who published [a paper in 1967 about a model of the Argon gas using 864 particles](https://journals.aps.org/pr/abstract/10.1103/PhysRev.159.98) following forces obeying the Lennard-Jones potential. The integration method that he used in the original paper is written as follows:\n", "\n", "$$\n", "\\vec{r}_i (t+h) = -\\vec{r}_i(t-h) + 2 \\vec{r}_i (t) + \\sum_{j \\neq i} \\vec{f}(r_{ij}(t)) h^2\n", "$$\n", "\n", "For our purpose, we will use the formulation found in Wikipedia:\n", "$$\n", "\\vec{r}_i (t+h) = -\\vec{r}_i(t-h) + 2 \\vec{r}_i (t) + \\vec{a}(r_{i}(t)) h^2\n", "$$\n", "\n", "In this equation, $\\vec{r}$ is the position vector for the particle we are modeling and $\\vec{a}$ is its acceleartion.\n", "\n", "As stated in the Wikipedia article, the interesting thing about this integration scheme is that it is an order more accurate than integration by simple Taylor expansion alone and that is does not use the velocity vector at all. For me, this is quite surprising since most numerical integration techniques I have used before integrate acceleration to get velocity and integrate velocity to get positions.\n", "\n", "This is all the theory we need. As we will see in the remainder of this post, this formulation is very general and will allow us to play with many examples: oscillators, particle systems, pendulums and many more." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Applications" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Harmonic oscillator" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's start with the harmonic oscillator since it is a good scalar example. The harmonic oscillator's equation is: \n", "$$\n", "\\ddot{x} + \\omega_0^2 x = 0\n", "$$\n", "\n", "Therefore, the scalar acceleration $a$ is computed as $-\\omega_0^2 x$." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "class VerletSolver:\n", " def __init__(self, dt, x0, xd0, acc_func):\n", " \"Inits the solver.\"\n", " self.dt = dt\n", " self.dt_squared = dt**2\n", " self.t = dt\n", " self.acc_func = acc_func\n", " self.x0 = x0\n", " self.xd0 = xd0\n", " # the second position vector is obtained by double integration of acceleration\n", " # using the initial velocity xd0\n", " x1 = acc_func(x0) * dt**2 / 2. + xd0 * dt + x0\n", " self.x = [x1, x0]\n", " \n", " def step(self):\n", " \"Steps the solver.\"\n", " xt, xtm1 = self.x\n", " xtp1 = - xtm1 + 2 * xt + self.acc_func(xt) * self.dt_squared\n", " self.x = (xtp1, xt)\n", " self.t += self.dt\n", " \n", " def step_until(self, tmax, snapshot_dt):\n", " \"Steps the solver until a given time, returns snapshots.\"\n", " ts = [self.t]\n", " vals = [self.x[0]]\n", " niter = max(1, int(snapshot_dt // self.dt))\n", " while self.t < tmax:\n", " for _ in range(niter):\n", " self.step()\n", " vals.append(self.x[0])\n", " ts.append(self.t)\n", " return np.array(ts), np.array(vals)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def acc_func_harmonic(x):\n", " \"\"\"Acceleration for harmonic oscilator.\"\"\"\n", " omega_0 = 1.\n", " return -omega_0 * x" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "solver = VerletSolver(0.01, 2, -0.01, acc_func_harmonic)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ts, vals = solver.step_until(12, snapshot_dt=0.1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(ts, vals, '-o')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's see how \"stable\" this is if we vary the integration step size." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for step_size in [0.001, 0.01, 0.05, 0.25]:\n", " solver = VerletSolver(step_size, 2, -0.01, acc_func=acc_func_harmonic)\n", " ts, vals = solver.step_until(12, snapshot_dt=0.5)\n", " plt.plot(ts, vals, \"-o\", label='step_size: {}'.format(step_size))\n", "plt.legend()\n", "plt.ylim(-2.1, 2.1);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As the above graph shows, the Verlet integration works very well in this case." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A two-dimensional spring " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's move on to a two-dimensional spring. The difference with the previous example is that the spring is described by a vector instead of a scalar: its position depends on two coordinates. This integrates easily with the vector Verlet update equation that was explained in the theory section.\n", "\n", "This time, we thus use a 2d NumPy array to store initial positions and velocities `x0`, `xd0`.\n", "\n", "One more thing we need is the description of the behavior of a spring. To quote H-P Langtangen:\n", "\n", "> Stretching the elastic wire a distance ΔL gives rise to a spring force kΔL in the opposite direction of the stretching.\n", "\n", "This is what we implement below:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def acc_func_spring(x):\n", " \"\"\"Acceleration for a 2D spring.\"\"\"\n", " k = 10\n", " stretch = np.linalg.norm(x) \n", " n = x / np.linalg.norm(x)\n", " acc = - k * stretch * n\n", " return acc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's define initial conditions with a position and a speed and then run the simulation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0 = np.array([0, -1], dtype=np.float)\n", "xd0 = np.array([0.0, 0.1], dtype=np.float)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "spring = VerletSolver(0.01, x0, xd0, acc_func_spring)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ts, vals = spring.step_until(12, snapshot_dt=0.1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(vals[:, 0], vals[:, 1])\n", "plt.axis('equal');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's plot the trajectory according to time:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(ts, vals[:, 0])\n", "plt.plot(ts, vals[:, 1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or as a frame-by-frame animation:" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "from IPython.display import HTML\n", "import matplotlib.animation as mpl_anim" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 50\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " ax.plot(vals[index, 0], vals[index, 1], 'o')\n", " ax.plot(vals[index-5:index+1, 0], vals[index-5:index+1, 1], '-')\n", " ax.set_xlim(-1, 1)\n", " ax.set_ylim(-1, 1)\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we give the spring a slight nudge to the right, we can get elliptical oscillations:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xd0 = np.array([-0.5, .1], dtype=np.float)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "spring = VerletSolver(0.1, x0, xd0, acc_func_spring)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ts, vals = spring.step_until(12, snapshot_dt=0.1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(vals[:, 0], vals[:, 1])\n", "plt.axis('equal');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's plot the trajectory according to time:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(ts, vals[:, 0])\n", "plt.plot(ts, vals[:, 1])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 50\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " ax.plot(vals[index, 0], vals[index, 1], 'o')\n", " ax.plot(vals[index-5:index+1, 0], vals[index-5:index+1, 1], '-')\n", " ax.set_xlim(-1, 1)\n", " ax.set_ylim(-1, 1)\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Shooting a cannonball " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def acc_func_cannonball(x):\n", " \"\"\"Force function for a cannonball.\"\"\"\n", " g = 10.\n", " acc = g * np.array([0., -1.])\n", " return acc" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0 = np.array([0, 0], dtype=np.float)\n", "xd0 = np.array([10.0, 10.0], dtype=np.float)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cannonball = VerletSolver(0.001, x0, xd0, acc_func_cannonball)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ts, vals = cannonball.step_until(2, snapshot_dt=0.1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(vals[:, 0], vals[:, 1])\n", "plt.axis('equal');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's plot the trajectory according to time:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(ts, vals[:, 0], label='x')\n", "plt.plot(ts, vals[:, 1], label='y')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's now fire a lot of cannonbals." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for angle in np.deg2rad(np.linspace(0, 90, num=20)):\n", " x0 = np.array([0, 0], dtype=np.float)\n", " xd0 = 15 * np.array([np.cos(angle), np.sin(angle)], dtype=np.float)\n", " cannonball = VerletSolver(0.001, x0, xd0, acc_func_cannonball)\n", " ts, vals = cannonball.step_until(2, snapshot_dt=0.1)\n", " plt.plot(vals[:, 0], vals[:, 1], '-ok')\n", "plt.axis('equal')\n", "plt.axis([0, 20, 0, 7])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's do an animation of cannonball isochrones. We first precompute the trajectories:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t_max = 3.5\n", "trajectories = {}\n", "angles = np.deg2rad(np.linspace(25, 85, num=20))\n", "for angle in angles:\n", " x0 = np.array([0, 0], dtype=np.float)\n", " xd0 = 15 * np.array([np.cos(angle), np.sin(angle)], dtype=np.float)\n", " cannonball = VerletSolver(0.001, x0, xd0, acc_func_cannonball)\n", " ts, vals = cannonball.step_until(t_max, snapshot_dt=0.01)\n", " trajectories[angle] = (ts, vals)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And then animate the plot of the simultaneous trajectories as a function of time." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 50\n", "def animate(i):\n", " ax.clear() \n", " for angle in angles:\n", " ts, vals = trajectories[angle]\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.plot(vals[:index, 0], vals[:index, 1], '-k')\n", " ax.plot(vals[index, 0], vals[index, 1], 'ok')\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", " ax.set_xlim(-1, 19)\n", " ax.set_ylim(0, 20)\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Asteroids around the earth " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also simulate some asteroids orbiting around the earth (considered fixed).\n", "\n", "First, we generate a set of asteroids as well as their velocities." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def make_asteroids(n, rmin, rmax):\n", " \"\"\"Creates a bunch of asteroids.\"\"\"\n", " rs = np.random.uniform(rmin, rmax, size=n) \n", " thetas = np.random.uniform(0, 360, size=n)\n", " x = np.stack([rs * np.cos(thetas), rs * np.sin(thetas)], axis=1)\n", " n = np.stack([rs * np.sin(thetas), -rs * np.cos(thetas)], axis=1)\n", " xd = n \n", " return x, xd" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0, xd0 = make_asteroids(200, 0.5, 1.1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "ax.quiver(x0[:, 0], x0[:, 1], xd0[:, 0], xd0[:, 1])\n", "ax.plot([0], [0], 'ro', ms=10, label='earth')\n", "ax.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We then create a plotting function:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_asteroids(x, ax, **plot_kwargs):\n", " xy = x.reshape(-1, 2)\n", " ax.plot(xy[:, 0], xy[:, 1], 'o', **plot_kwargs)\n", " ax.plot([0], [0], 'ro', ms=10, label='earth')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "plot_asteroids(x0, ax, ms=4, alpha=0.5)\n", "ax.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now define the gravitational force, which always pulls the asteroids towards the earth." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def acc_func_asteroids(x):\n", " \"\"\"Acceleration for asteroids orbiting the earth.\"\"\"\n", " xy = x.reshape(-1, 2)\n", " r = np.linalg.norm(xy - np.array([0., 0.]), axis=1)\n", " n = xy / r[:, np.newaxis]\n", " acc = - 0.0000001 * 1/r[:, np.newaxis]**2 * n\n", " return acc.ravel()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "asteroid_field = VerletSolver(1., \n", " x0.ravel(), \n", " (xd0.ravel() + np.random.randn(xd0.size) / 20) / 2000, \n", " acc_func_asteroids)\n", "ts, vals = asteroid_field.step_until(40000, 100)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 100\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " plot_asteroids(vals[index], ax)\n", " if index>4:\n", " plot_asteroids(vals[index-1], ax, alpha=0.5, ms=5)\n", " plot_asteroids(vals[index-2], ax, alpha=0.4, ms=4)\n", " plot_asteroids(vals[index-3], ax, alpha=0.3, ms=3)\n", " plot_asteroids(vals[index-4], ax, alpha=0.2, ms=2)\n", " ax.axis([-2.5, 2.5, -2.5, 2.5])\n", " ax.legend(loc='upper right')\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=200, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's also render a close-up view:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 100\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " plot_asteroids(vals[index], ax)\n", " if index>4:\n", " plot_asteroids(vals[index-1], ax, alpha=0.5, ms=5)\n", " plot_asteroids(vals[index-2], ax, alpha=0.4, ms=4)\n", " plot_asteroids(vals[index-3], ax, alpha=0.3, ms=3)\n", " plot_asteroids(vals[index-4], ax, alpha=0.2, ms=2)\n", " ax.axis([-1.25, 1.25, -1.25, 1.25])\n", " ax.legend(loc='upper right')\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=200, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As one can see, it seems that particles closer to the earth have faster orbits, while particles more farther away go more slowly." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A pendulum " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's now turn to pendulums. A pendulum is basically a spring with an additional force: the pull of gravity along the vertical axis. Let's write that as a function:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def acc_func_pendulum(x):\n", " \"\"\"Force function with spring of fixed size.\"\"\"\n", " g = 1.\n", " l0 = 1.\n", " elongation = np.linalg.norm(x) - l0\n", " k = 100.\n", " acc = - k * elongation * x / np.linalg.norm(x) + g * np.array([0, -1])\n", " return acc" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0 = np.array([1, 0], dtype=np.float)\n", "xd0 = np.array([0.0, 0.0], dtype=np.float)\n", "pendulum = VerletSolver(0.01, x0, xd0, acc_func_pendulum)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ts, vals = pendulum.step_until(15, snapshot_dt=0.05)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(vals[:, 0], vals[:, 1])\n", "plt.axis('equal');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's plot the trajectory according to time:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(ts, vals[:, 0], label='x')\n", "plt.plot(ts, vals[:, 1], label='y')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or as an animation:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 50\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " ax.plot([0, vals[index, 0]], [0, vals[index, 1]])\n", " ax.plot(vals[index, 0], vals[index, 1], 'o')\n", " ax.plot(vals[index-10:index+1, 0], vals[index-10:index+1, 1], '-')\n", " ax.axis('equal')\n", " ax.axis([-1.5, 1.5, -2, 1])\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One interesting thing about pendulums is that when you make the spring less stiff, you can obtain more varied behaviours:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def acc_func_pendulum_elastic(x):\n", " \"\"\"Force function with spring of fixed size.\"\"\"\n", " g = 1.\n", " l0 = 1.\n", " elongation = np.linalg.norm(x) - l0\n", " k = 5.\n", " acc = - k * elongation * x / np.linalg.norm(x) + g * np.array([0, -1])\n", " return acc" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0 = np.array([1, 0], dtype=np.float)\n", "xd0 = np.array([0.0, 0.0], dtype=np.float)\n", "\n", "pendulum = VerletSolver(0.01, x0, xd0, acc_func_pendulum_elastic)\n", "\n", "ts, vals = pendulum.step_until(30, snapshot_dt=0.1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 50\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " ax.plot([0, vals[index, 0]], [0, vals[index, 1]])\n", " ax.plot(vals[index, 0], vals[index, 1], 'o')\n", " ax.plot(vals[index-10:index+1, 0], vals[index-10:index+1, 1], '-')\n", " ax.axis('equal')\n", " ax.axis([-1.5, 1.5, -2, 1])\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Two particle systems" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another interesting example consists of two particles attached by a spring." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def acc_func_two_spring(x):\n", " \"\"\"Force function with spring of fixed size.\"\"\"\n", " g = 1.\n", " l0 = 1.\n", " x1 = x[:2]\n", " x2 = x[2:]\n", " elongation1 = np.linalg.norm(x1 - x2) - l0\n", " elongation2 = np.linalg.norm(x2 - x1) - l0\n", " k = 10.\n", " acc = np.zeros_like(x)\n", " acc[:2] = - k * elongation1 * (x1 - x2) / np.linalg.norm(x1 - x2)\n", " acc[2:] = - k * elongation2 * (x2 - x1) / np.linalg.norm(x2 - x1)\n", " return acc" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0 = np.array([0.75, 0, -0.75, 0], dtype=np.float)\n", "xd0 = np.array([0.0, 0.0, 0.0, 0.0], dtype=np.float)\n", "twosprings = VerletSolver(0.01, x0, xd0, acc_func_two_spring)\n", "ts, vals = twosprings.step_until(2.8, snapshot_dt=0.05)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's plot the trajectory according to time:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(ts, vals[:, 0], label='x')\n", "plt.plot(ts, vals[:, 1], label='y')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's animate this." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 50\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " ax.plot(vals[index, [0, 2]], vals[index, [1, 3]], 'o', ms=50)\n", " ax.axis('equal')\n", " ax.axis([-1.5, 1.5, -1, 1])\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What happens if the particles start with a little spin?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0 = np.array([0.75, 0, -0.75, 0], dtype=np.float)\n", "xd0 = np.array([0.0, 0.1, 0.0, -.1], dtype=np.float)\n", "twosprings = VerletSolver(0.01, x0, xd0, acc_func_two_spring)\n", "ts, vals = twosprings.step_until(8, snapshot_dt=0.05)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 50\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " ax.plot(vals[index, [0, 2]], vals[index, [1, 3]], 'o')\n", " # tail 1\n", " ax.plot(vals[0:index+1, 0], vals[0:index+1, 1], '-')\n", " # tail 2\n", " ax.plot(vals[0:index+1, 2], vals[0:index+1, 3], '-')\n", " ax.axis('equal')\n", " ax.axis([-1.5, 1.5, -1, 1])\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "They still pull on each other, however the two particles now also \"dance in a circle\" following their initial momentum." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A chain of oscillators strung together " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's build a chain of oscillators that are fixed on both ends." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def build_chain(n_free):\n", " \"\"\"Returns x0 and xd0 for a chain of n_free horizontal particles.\"\"\"\n", " l0 = 1\n", " xs = np.arange(n_free + 2, dtype=np.float) * l0\n", " ys = np.zeros_like(xs)\n", " x0 = np.c_[xs.reshape(-1, 1), ys.reshape(-1, 1)].ravel()\n", " xd0 = np.zeros_like(x0) \n", " return x0, xd0" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib.collections import LineCollection\n", "\n", "def plot_chain(ax, x):\n", " xy = x.reshape(-1, 2)\n", " line = LineCollection([[start, stop] for start, stop in zip(xy[:-1,:], xy[1:,:])], linestyles='-.')\n", " ax.add_artist(line)\n", " ax.plot(xy[:, 0], xy[:, 1], 'o')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0, xd0 = build_chain(3)\n", "x0[4] += 50 / 100" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "plot_chain(ax, x0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's write the acceleration function." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def acc_func_chain(x):\n", " \"\"\"Acceleration of a chain of oscillators strung together. \n", " Node 0 and -1 are assumed fixed.\"\"\"\n", " k = 10.\n", " l0 = 1\n", " xy = x.reshape(-1, 2)\n", " acc = np.zeros_like(xy)\n", " for index in range(1, xy.shape[0] - 1):\n", " x0 = xy[index-1, :]\n", " x1 = xy[index, :]\n", " x2 = xy[index+1, :]\n", " \n", " elongation1 = (np.linalg.norm(x1 - x0) - l0) / l0\n", " elongation2 = (np.linalg.norm(x2 - x1) - l0) / l0\n", " acc[index, :] = - k * elongation1 * (x1 - x0) / np.linalg.norm(x1 - x0) \\\n", " - k * elongation2 * (x1 - x2) / np.linalg.norm(x2 - x1)\n", " \n", " return acc.ravel()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now let's test this:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "chain = VerletSolver(0.01, x0, xd0, acc_func_chain)\n", "ts, vals = chain.step_until(5, 0.1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's see what the trajectories of the individual particles are like:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "reshaped_vals = vals.reshape(-1, 5, 2)\n", "for particle_index in range(reshaped_vals.shape[1]):\n", " traj = reshaped_vals[:, particle_index, :]\n", " ax.plot(traj[:, 0])\n", " ax.plot(traj[:, 1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's animate this:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 50\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " plot_chain(ax, vals[index])\n", " ax.axis('equal')\n", " ax.axis([-0.5, 4.5, -0.1, 0.1])\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What if one of the extremities is loose? This is just a matter of updating the acceleration of the last particle in the chain with the pulling of its neighbor." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def acc_func_chain2(x):\n", " \"\"\"Acceleration of a chain of oscillators strung together. \n", " Node 0 and -1 are assumed fixed.\"\"\"\n", " k = 10.\n", " l0 = 1\n", " xy = x.reshape(-1, 2)\n", " acc = np.zeros_like(xy)\n", " for index in range(1, xy.shape[0] - 1):\n", " x0 = xy[index-1, :]\n", " x1 = xy[index, :]\n", " x2 = xy[index+1, :]\n", " \n", " elongation1 = (np.linalg.norm(x1 - x0) - l0) / l0\n", " elongation2 = (np.linalg.norm(x2 - x1) - l0) / l0\n", " acc[index, :] = - k * elongation1 * (x1 - x0) / np.linalg.norm(x1 - x0) \\\n", " - k * elongation2 * (x1 - x2) / np.linalg.norm(x2 - x1)\n", " index += 1\n", " x0 = xy[index-1, :]\n", " x1 = xy[index, :]\n", " elongation1 = (np.linalg.norm(x1 - x0) - l0) / l0\n", " acc[index, :] = - k * elongation1 * (x1 - x0) / np.linalg.norm(x1 - x0) \\\n", " \n", " return acc.ravel()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0, xd0 = build_chain(3)\n", "x0[4] += 50 / 100\n", "chain = VerletSolver(0.01, x0, xd0, acc_func_chain2)\n", "ts, vals = chain.step_until(5, 0.1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 50\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " plot_chain(ax, vals[index])\n", " ax.axis('equal')\n", " ax.axis([-0.5, 4.5, -0.1, 0.1])\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What if we add a little bit of noise to this chain of oscilators?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xd1 = xd0.copy()\n", "xd1[2:] += np.random.randn(8) / 100\n", "chain = VerletSolver(0.01, x0, xd1, acc_func_chain2)\n", "ts, vals = chain.step_until(40, 0.1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 50\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " plot_chain(ax, vals[index])\n", " ax.axis('equal')\n", " ax.axis([-4.5, 4.5, -0.1, 0.1])\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We get a behaviour that resembles that of a N-link pendulum!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What if we use more than 5 particles?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0, xd0 = build_chain(28)\n", "x0[-2] += .5\n", "chain = VerletSolver(0.01, x0, xd0, acc_func_chain2)\n", "ts, vals = chain.step_until(40, 0.1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(10, 5))\n", "\n", "FRAMES = 50\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " plot_chain(ax, vals[index])\n", " ax.axis([-0.5, 30.5, -0.1, 0.1])\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compression waves appear! Let's plot this data in the \"earth science style\" to see the wave advancing across time and space at the same time." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 8))\n", "reshaped_vals = vals.reshape(-1, 30, 2)\n", "for time_index in range(0, reshaped_vals.shape[0], 10):\n", " delta_traj = reshaped_vals[time_index, :, 0] - np.arange(30)\n", " ax.plot(delta_traj - time_index / 40, 'k')\n", "ax.axis('off');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What we see in the above graph is a double reflexion on the left side of the particle chain, as well as the physical phenomenon of dispersion: the wave started as a sharp peak but ended up as a diffuse velocity field." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2D waves in a lattice " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can of course extend this sort of chaining to a grid. Let's first write a function that creates a grid of particles." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "def build_grid(nx, ny):\n", " \"\"\"Builds a grid of nx times ny nodes.\"\"\"\n", " x = np.arange(nx, dtype=np.float)\n", " y = np.arange(ny, dtype=np.float)\n", " X, Y = np.meshgrid(x, y, indexing='ij')\n", " return X, Y" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "from numba import njit" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "@njit\n", "def grid2linear(X, Y):\n", " \"\"\"Returns a linear vector from a X,Y displacement representation.\"\"\"\n", " return np.concatenate((X.ravel(), Y.ravel()))" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "@njit\n", "def linear2grid(linear_repr, nx, ny):\n", " \"\"\"Returns a grid X,Y displacement from a linear vector.\"\"\"\n", " n = nx * ny\n", " X, Y = linear_repr[:n], linear_repr[n:2*n]\n", " X = X.reshape((nx, ny))\n", " Y = Y.reshape((nx, ny))\n", " return X, Y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's test these functions:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "nx, ny = 3, 2\n", "X, Y = build_grid(nx, ny)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 0. 0.]\n", " [ 1. 1.]\n", " [ 2. 2.]]\n", "[[ 0. 1.]\n", " [ 0. 1.]\n", " [ 0. 1.]]\n" ] } ], "source": [ "print(X)\n", "print(Y)" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 0., 0., 1., 1., 2., 2., 0., 1., 0., 1., 0., 1.])" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "linear_repr = grid2linear(X, Y)\n", "linear_repr" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([[ 0., 0.],\n", " [ 1., 1.],\n", " [ 2., 2.]]), array([[ 0., 1.],\n", " [ 0., 1.],\n", " [ 0., 1.]]))" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "linear2grid(linear_repr, nx, ny)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's now code the acceleration function, which will connect each node with its neighbors." ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "@njit\n", "def acc_func_grid_template(x, nx, ny, k):\n", " \"\"\"Acceleration of a grid of oscillators strung together to its 4 neighbors.\"\"\"\n", " l0 = 1.\n", " l1 = np.sqrt(2)\n", " X, Y = linear2grid(x, nx, ny)\n", " acc_X, acc_Y = np.zeros_like(X), np.zeros_like(Y)\n", " for i in range(nx):\n", " for j in range(ny):\n", " x_center = np.array([X[i, j], Y[i, j]])\n", " for neighbor in [(i+1, j, l0), (i-1, j, l0),\n", " (i, j+1, l0), (i, j-1, l0),\n", " (i+1, j+1, l1), (i+1, j-1, l1),\n", " (i-1, j+1, l1), (i-1, j-1, l1)]:\n", " ii, jj, l = neighbor\n", " if (ii >= 0 and ii < nx) and (jj >=0 and jj < ny):\n", " x_neighb = np.array([X[ii, jj], Y[ii, jj]])\n", " stretch = np.linalg.norm(x_center - x_neighb) - l\n", " acc = - k * l0 / l * stretch * (x_center - x_neighb) / np.linalg.norm(x_center - x_neighb)\n", " acc_X[i, j] += acc[0]\n", " acc_Y[i, j] += acc[1]\n", " return grid2linear(acc_X, acc_Y)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "from functools import partial" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/plain": [ "array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "acc_func_grid = partial(acc_func_grid_template, nx=nx, ny=ny, k=20.)\n", "acc_func_grid(linear_repr)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [], "source": [ "linear_repr[1] += 0.1" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 0.00992562, -2.66129819, 0.65137257, 2. , 0. ,\n", " 0. , 0.0992562 , 0.62449111, -0.7237473 , 0. ,\n", " 0. , 0. ])" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "acc_func_grid(linear_repr)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [], "source": [ "def plot_grid(x, nx, ny, ax=None):\n", " X, Y = linear2grid(x, nx, ny)\n", " if ax is None:\n", " ax = plt.gca()\n", " ax.scatter(X.ravel(), Y.ravel())" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAEWZJREFUeJzt3X2MZXddx/H3x90WiiIFd4ywD2yJ\nS7UYTWFSEYwWwXRbQxcjylaJgNUNalGjadKmppoa40P/8CHWh0oIlGhLqbWuZMmKUIIBFzul0KWt\ni8vy0OkSu9a2SljpQ77+cc+W29uZnTO792H66/uVTPac3/ndc777u7/5zJlz7p2bqkKS1JZvmnUB\nkqTxM9wlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDVo/qwNv2LChtm7dOqvDS9LT\n0u233/5fVTW3Ur+ZhfvWrVtZWFiY1eEl6WkpyZf69POyjCQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3\nSWqQ4S5JDVox3JO8K8n9ST67zPYk+dMkB5PcmeTl4y9TkrQafd7E9G7gz4Drltl+PrCt+/p+4C+6\nf9ekW+64j6v3HuDwQ0d50emncel5Z/KGszfOuiytAc4NTcos5taK4V5VH0uy9ThddgDX1eCTtvcl\nOT3JC6vqK2OqcWxuueM+Lr95P0cffRyA+x46yuU37wfwm/gZzrmhSZnV3BrHNfeNwL1D64td25pz\n9d4DTwzwMUcffZyr9x6YUUVaK5wbmpRZza1xhHuWaKslOya7kiwkWThy5MgYDr06hx86uqp2PXM4\nNzQps5pb4wj3RWDz0Pom4PBSHavq2qqar6r5ubkV/6jZ2L3o9NNW1a5nDueGJmVWc2sc4b4b+Nnu\nVTOvBB5ei9fbAS4970xOO2Xdk9pOO2Udl5535owq0lrh3NCkzGpurXhDNcn1wLnAhiSLwG8BpwBU\n1V8Ce4ALgIPA14C3TarYk3Xs5oWviNAo54YmZVZzK4MXuUzf/Px8+ffcJWl1ktxeVfMr9fMdqpLU\nIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y\n3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNd\nkhpkuEtSgwx3SWqQ4S5JDTLcJalBvcI9yfYkB5IcTHLZEtu3JLk1yR1J7kxywfhLlST1tWK4J1kH\nXAOcD5wFXJTkrJFuvwncWFVnAzuBPx93oZKk/vqcuZ8DHKyqQ1X1CHADsGOkTwHf2i0/Dzg8vhIl\nSavVJ9w3AvcOrS92bcN+G3hzkkVgD/COpXaUZFeShSQLR44cOYFyJUl99An3LNFWI+sXAe+uqk3A\nBcB7kzxl31V1bVXNV9X83Nzc6quVJPXSJ9wXgc1D65t46mWXi4EbAarqX4FnAxvGUaAkafX6hPtt\nwLYkZyQ5lcEN090jfb4MvBYgyXczCHevu0jSjKwY7lX1GHAJsBe4h8GrYu5KclWSC7tuvwH8QpLP\nANcDb62q0Us3kqQpWd+nU1XtYXCjdLjtyqHlu4FXj7c0SdKJ8h2qktQgw12SGmS4S1KDDHdJapDh\nLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S\n1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkN\n6hXuSbYnOZDkYJLLlunzU0nuTnJXkr8db5mSpNVYv1KHJOuAa4AfBRaB25Lsrqq7h/psAy4HXl1V\nDyb59kkVLElaWZ8z93OAg1V1qKoeAW4Adoz0+QXgmqp6EKCq7h9vmZKk1egT7huBe4fWF7u2YS8F\nXprk40n2Jdk+rgIlSau34mUZIEu01RL72QacC2wC/iXJ91TVQ0/aUbIL2AWwZcuWVRcrSeqnz5n7\nIrB5aH0TcHiJPv9QVY9W1ReAAwzC/kmq6tqqmq+q+bm5uROtWZK0gj7hfhuwLckZSU4FdgK7R/rc\nArwGIMkGBpdpDo2zUElSfyuGe1U9BlwC7AXuAW6sqruSXJXkwq7bXuCBJHcDtwKXVtUDkypaknR8\nqRq9fD4d8/PztbCwMJNjS9LTVZLbq2p+pX6+Q1WSGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCX\npAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lq\nkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1qFe4J9me\n5ECSg0kuO06/NyapJPPjK1GStForhnuSdcA1wPnAWcBFSc5aot9zgV8BPjnuIiVJq9PnzP0c4GBV\nHaqqR4AbgB1L9Psd4A+B/xtjfZKkE9An3DcC9w6tL3ZtT0hyNrC5qj5wvB0l2ZVkIcnCkSNHVl2s\nJKmfPuGeJdrqiY3JNwF/BPzGSjuqqmurar6q5ufm5vpXKUlalT7hvghsHlrfBBweWn8u8D3AR5N8\nEXglsNubqpI0O33C/TZgW5IzkpwK7AR2H9tYVQ9X1Yaq2lpVW4F9wIVVtTCRiiVJK1ox3KvqMeAS\nYC9wD3BjVd2V5KokF066QEnS6q3v06mq9gB7RtquXKbvuSdfliTpZPgOVUlqkOEuSQ0y3CWpQYa7\nJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtS\ngwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXI\ncJekBvUK9yTbkxxIcjDJZUts//Ukdye5M8mHk7x4/KVKkvpaMdyTrAOuAc4HzgIuSnLWSLc7gPmq\n+l7gJuAPx12oJKm/Pmfu5wAHq+pQVT0C3ADsGO5QVbdW1de61X3ApvGWKUlajT7hvhG4d2h9sWtb\nzsXAB0+mKEnSyVnfo0+WaKslOyZvBuaBH15m+y5gF8CWLVt6lihJWq0+Z+6LwOah9U3A4dFOSV4H\nXAFcWFVfX2pHVXVtVc1X1fzc3NyJ1CtJ6qFPuN8GbEtyRpJTgZ3A7uEOSc4G/opBsN8//jIlSaux\nYrhX1WPAJcBe4B7gxqq6K8lVSS7sul0NfAvw/iSfTrJ7md1JkqagzzV3qmoPsGek7cqh5deNuS5J\n0knwHaqS1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJ\napDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QG\nGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQev7dEqyHfgTYB3wzqr6/ZHtzwKuA14BPAC8qaq+\nON5Sv+GWO+7j6r0HOPzQUV50+mlcet6ZvOHsjZM6nJ4hnFealFnMrRXDPck64BrgR4FF4LYku6vq\n7qFuFwMPVtV3JtkJ/AHwpkkUfMsd93H5zfs5+ujjANz30FEuv3k/gN+IOmHOK03KrOZWn8sy5wAH\nq+pQVT0C3ADsGOmzA3hPt3wT8NokGV+Z33D13gNPDNIxRx99nKv3HpjE4fQM4bzSpMxqbvUJ943A\nvUPri13bkn2q6jHgYeDbRneUZFeShSQLR44cOaGCDz90dFXtUh/OK03KrOZWn3Bf6gy8TqAPVXVt\nVc1X1fzc3Fyf+p7iRaeftqp2qQ/nlSZlVnOrT7gvApuH1jcBh5frk2Q98Dzgv8dR4KhLzzuT005Z\n96S2005Zx6XnnTmJw+kZwnmlSZnV3OrzapnbgG1JzgDuA3YCPz3SZzfwFuBfgTcCH6mqp5y5j8Ox\nGxC+qkHj5LzSpMxqbqVPBie5APhjBi+FfFdV/W6Sq4CFqtqd5NnAe4GzGZyx76yqQ8fb5/z8fC0s\nLJz0f0CSnkmS3F5V8yv16/U696raA+wZabtyaPn/gJ9cbZGSpMnwHaqS1CDDXZIaZLhLUoMMd0lq\nkOEuSQ0y3CWpQYa7JDWo15uYJnLg5AjwpZPczQbgv8ZQzjitxZrAulZjLdYE1rVaa7GucdT04qpa\n8Y9zzSzcxyHJQp93ak3TWqwJrGs11mJNYF2rtRbrmmZNXpaRpAYZ7pLUoKd7uF876wKWsBZrAuta\njbVYE1jXaq3FuqZW09P6mrskaWlP9zN3SdIS1mS4J9me5ECSg0kuW2L7s5K8r9v+ySRbh7Zd3rUf\nSHLelOv69SR3J7kzyYeTvHho2+NJPt197Z5yXW9NcmTo+D8/tO0tSf6j+3rLFGv6o6F6PpfkoaFt\nkxyrdyW5P8lnl9meJH/a1X1nkpcPbZvUWK1U0890tdyZ5BNJvm9o2xeT7O/GaqwfkNCjrnOTPDz0\nXF05tO24z/+E67p0qKbPdvPpBd22iYxXks1Jbk1yT5K7kvzqEn2mO7eqak19MfhAkM8DLwFOBT4D\nnDXS55eAv+yWdwLv65bP6vo/Czij28+6Kdb1GuA53fIvHqurW//qDMfrrcCfLfHYFwCHun+f3y0/\nfxo1jfR/B4MPgZnoWHX7/iHg5cBnl9l+AfBBBp8L/Ergk5Mcq541verYsYDzj9XUrX8R2DCjsToX\n+MDJPv/jrmuk7+sZfDLcRMcLeCHw8m75ucDnlvg+nOrcWotn7ucAB6vqUFU9AtwA7BjpswN4T7d8\nE/DaJOnab6iqr1fVF4CD3f6mUldV3VpVX+tW9zH4vNlJ6zNeyzkP+FBV/XdVPQh8CNg+g5ouAq4f\nw3FXVFUf4/if77sDuK4G9gGnJ3khkxurFWuqqk90x4Tpzas+Y7Wck5mT465rKnOrqr5SVZ/qlv8X\nuAcY/Ry9qc6ttRjuG4F7h9YXeeogPdGnqh4DHga+redjJ1nXsIsZ/JQ+5tlJFpLsS/KGMdW0mrp+\novtV8KYkxz7wfFLj1Xu/3aWrM4CPDDVPaqz6WK72Sc6t1RidVwX8U5Lbk+yaQT0/kOQzST6Y5GVd\n25oYqyTPYRCSfzfUPPHxyuAy8dnAJ0c2TXVu9fqYvSnLEm2jL+lZrk+fx56o3vtO8mZgHvjhoeYt\nVXU4yUuAjyTZX1Wfn1Jd/whcX1VfT/J2Br/1/EjPx06qpmN2AjdV1eNDbZMaqz5mMbd6SfIaBuH+\ng0PNr+7G6tuBDyX59+7Mdho+xeCt8F/N4HOWbwG2sQbGqvN64ONVNXyWP9HxSvItDH6Y/FpV/c/o\n5iUeMrG5tRbP3BeBzUPrm4DDy/VJsh54HoNf0/o8dpJ1keR1wBXAhVX19WPtVXW4+/cQ8FEGP9mn\nUldVPTBUy18Dr+j72EnVNGQnI782T3Cs+liu9knOrRUl+V7gncCOqnrgWPvQWN0P/D3juwy5oqr6\nn6r6are8BzglyQZmPFZDjje3xj5eSU5hEOx/U1U3L9FlunNr3DcWxnBjYj2DGwpn8I2bMS8b6fPL\nPPmG6o3d8st48g3VQ4zvhmqfus5mcCNp20j784FndcsbgP9gTDeYetb1wqHlHwf21Tdu5Hyhq+/5\n3fILplFT1+9MBje4Mo2xGjrGVpa/SfhjPPmm179Ncqx61rSFwf2jV420fzPw3KHlTwDbpzhW33Hs\nuWMQkl/uxq3X8z+purrtx074vnka49X9v68D/vg4faY6t8Y22GN+4i5gcLf588AVXdtVDM6GAZ4N\nvL+b8P8GvGTosVd0jzsAnD/luv4Z+E/g093X7q79VcD+bpLvBy6ecl2/B9zVHf9W4LuGHvtz3Tge\nBN42rZq69d8Gfn/kcZMeq+uBrwCPMjhjuhh4O/D2bnuAa7q69wPzUxirlWp6J/Dg0Lxa6Npf0o3T\nZ7rn94opj9UlQ/NqH0M/fJZ6/qdVV9fnrQxeXDH8uImNF4NLZQXcOfQ8XTDLueU7VCWpQWvxmrsk\n6SQZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNej/Ac4Q57g9RYGnAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_grid(linear_repr, nx=nx, ny=ny)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's try this with some random noise." ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [], "source": [ "nx, ny = 3, 3\n", "X, Y = build_grid(nx, ny)\n", "x0 = grid2linear(X, Y)\n", "x0[1] += 0.3\n", "xd0 = np.zeros_like(x0)\n", "acc_func_grid = partial(acc_func_grid_template, nx=nx, ny=ny, k=10.)\n", "grid = VerletSolver(0.01, x0, xd0, acc_func_grid)\n", "ts, vals = grid.step_until(2.5, snapshot_dt=0.1)" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 0.12652114, -4.82279424, 0.12652114, 0.78487598, 3. ,\n", " 0.78487598, 0. , 0. , 0. , 0.42173715,\n", " 0. , -0.42173715, -1.12125139, 0. , 1.12125139,\n", " 0. , 0. , 0. ])" ] }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ "acc_func_grid(x0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's now plot the grid of particles, between its initial and final state in the simulation. " ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAEJ1JREFUeJzt3X+oXGedx/H3d5NbueCSC+aCbZIa\nd7eE9UekdbZWXJZCkFTZNkVLzP6xWlkpuHa7ukvA+kcs+cddArrWiiWuRbOINdQQUmkJWBd0Fyy9\naWrSGrIbZSX3puC13dwqDppkv/vHTMzNZG5mJpmZc+6T9wsud85zntzz7TPPfHrm/JiJzESSVJY/\nqLoASdLwGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAq2sasOrV6/O9evXV7V5\nSVqWDh069MvMnO7Vr7JwX79+PTMzM1VtXpKWpYj4eT/9PCwjSQUy3CWpQIa7JBXIcJekAhnuklQg\nw12SCmS4S1KBDHdJKpDhLkmjdGQvfOFt8NBU6/eRvWPZbGV3qEpS8Y7shScfgDPN1vLCydYywMat\nI920e+6SNCrP7LwQ7OedabbaR8xwl6RRWZgdrH2Ieh6WiYh1wB7gjcD/Absz84sdfQL4IvB+4DfA\nvZn5/PDLbdl/eI5dB49z6nSTG6Ym2b55A3ffvGZUm1PBnEsalq5zadXa1qGYTqvWjryefvbczwL/\nmJl/CtwGfCIi3tLR533ATe2f+4CvDLXKRfYfnuPBfUeZO90kgbnTTR7cd5T9h+dGtUkVyrmkYVlq\nLj33x38HE5MXd56YhE07Rl5Tz3DPzJfP74Vn5q+AY0Dnrs0WYE+2/AiYiojrh14tsOvgcZpnzl3U\n1jxzjl0Hj49icyqYc0nDstRc+uRPboI7H4ZV64Bo/b7z4ZGfTIUBr5aJiPXAzcCzHavWAIvfe8y2\n217u+Pf30dqz58Ybbxys0rZTp5sDtUtLcS5pWC47lzZuHUuYd+r7hGpEvB74DvDJzHytc3WXf5KX\nNGTuzsxGZjamp3t+kUhXN0xNDtQuLcW5pGGp41zqK9wjYoJWsH8zM/d16TILrFu0vBY4dfXlXWr7\n5g1MTqy4qG1yYgXbN2+o7GYBLU+XnUvSAOo4l/q5WiaArwHHMvPzS3Q7ANwfEY8D7wIWMvPlJfpe\nlfNXMlxyVnrFf1Z2s4CWpyXnklfLaEB1nEuRecnRk4s7RPw58EPgKK1LIQE+A9wIkJmPtv8H8Ahw\nB61LIT+amZf9gtRGo5FD/Q7VL7xtiUuO1sGnXhzediSpQhFxKDMbvfr13HPPzP+g+zH1xX0S+ET/\n5Y1AhTcLqDBH9rbuIFyYbV2PvGmH7/607JRzh+pSNwWM4WYBFeT8Z4EsnATywuE9z99omSkn3Dft\nqOxmARWkws8CkYapnHDfuLWymwVUEA/vqRBlfeRvRTcLqCAVfhaINEzl7LlLw+DhPRXCcJcW8/Ce\nClHWYRlpGDy8pwK45y5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWp\nQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpk\nuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIK1DPcI+KxiPhFRLy4xPrbI2IhIl5o\n/+wYfpmSpEGs7KPP14FHgD2X6fPDzPzLoVQkSbpqPffcM/MHwKtjqEWSNCTDOub+7oj4cUQ8HRFv\nXapTRNwXETMRMTM/Pz+kTUuSOg0j3J8H3pSZ7wC+BOxfqmNm7s7MRmY2pqenh7BpSVI3Vx3umfla\nZv66/fgpYCIiVl91ZZKkK3bV4R4Rb4yIaD++tf03X7navytJunI9r5aJiG8BtwOrI2IW+CwwAZCZ\njwL3AB+PiLNAE9iWmTmyiiVJPfUM98z8qx7rH6F1qaQkqSa8Q1WSCmS4S1KBDHdJKpDhLkkFMtwl\nqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIK\nZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCG\nuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBeoZ7hHxWET8IiJeXGJ9RMTDEXEiIo5ExC3D\nL1OSNIiVffT5OvAIsGeJ9e8Dbmr/vAv4Svt3PRzZC8/shIVZWLUWNu2AjVvHWsL+w3PsOnicU6eb\n3DA1yfbNG7j75jVjrUHd9XxuajB/VAN9zIO6vc57hntm/iAi1l+myxZgT2Ym8KOImIqI6zPz5SHV\neOWO7IUnH4AzzdbywsnWMoztBbr/8BwP7jtK88w5AOZON3lw31EAA75iPZ+bGswf1UAf86COr/Nh\nHHNfA5xctDzbbqveMzsvPCHnnWm22sdk18Hjv3/Cz2ueOceug8fHVoO66/nc1GD+qAb6mAd1fJ0P\nI9yjS1t27RhxX0TMRMTM/Pz8EDbdw8LsYO0jcOp0c6B2jU/P56YG80c10Mc8qOPrfBjhPgusW7S8\nFjjVrWNm7s7MRmY2pqenh7DpHlatHax9BG6YmhyoXePT87mpwfxRDfQxD+r4Oh9GuB8APty+auY2\nYKEWx9uhddJjomNwJyZb7WOyffMGJidWXNQ2ObGC7Zs3jK0GddfzuanB/FEN9DEP6vg673lCNSK+\nBdwOrI6IWeCzwARAZj4KPAW8HzgB/Ab46KiKHdj5k14VXu1w/mRKnc6iq6Xnc1OD+aMa6GMe1PF1\nHq2LXMav0WjkzMxMJduWpOUqIg5lZqNXP+9QlaQCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWp\nQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpk\nuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7\nJBXIcJekAhnuklQgw12SCmS4S1KBDHdJKlBf4R4Rd0TE8Yg4ERGf7rL+3oiYj4gX2j8fG36pkqR+\nrezVISJWAF8G3gvMAs9FxIHM/ElH129n5v0jqFGSNKB+9txvBU5k5s8y83fA48CW0ZYlSboa/YT7\nGuDkouXZdlunD0bEkYh4IiLWDaU6SdIV6Sfco0tbdiw/CazPzI3A94BvdP1DEfdFxExEzMzPzw9W\nqSSpb/2E+yyweE98LXBqcYfMfCUzf9te/Crwzm5/KDN3Z2YjMxvT09NXUq8kqQ/9hPtzwE0R8eaI\nuA7YBhxY3CEirl+0eBdwbHglSpIG1fNqmcw8GxH3AweBFcBjmflSROwEZjLzAPBARNwFnAVeBe4d\nYc2SpB4is/Pw+Xg0Go2cmZmpZNuStFxFxKHMbPTq5x2qklQgw12SCmS4S1KBDHdJKpDhLkkFMtwl\nqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIK\nZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCG\nu9TpyF74wtvgoanW7yN7q65IGtjKqguQauXIXnjyATjTbC0vnGwtA2zcWl1d0oDcc5cWe2bnhWA/\n70yz1S4tI2WFu2+ndbUWZgdrl2qqnHA//3Z64SSQF95OG/AaxKq1g7VLNVVOuPt2WsOwaQdMTF7c\nNjHZapeWkXLC3bfTGoaNW+HOh2HVOiBav+982JOpWnb6ulomIu4AvgisAP41M/+pY/3rgD3AO4FX\ngA9l5v8Mt9QeVq1tH5Lp0i4NYuNWw1zLXs9wj4gVwJeB9wKzwHMRcSAzf7Ko298A/5uZfxIR24B/\nBj40ioIB9h+eY9fB45w63eSGqUm2b97A3Zt2XHwJG/h2Wj11nUs3r6m6LC1DdZtL/RyWuRU4kZk/\ny8zfAY8DWzr6bAG+0X78BLApImJ4ZV6w//AcD+47ytzpJgnMnW7y4L6j7D/3Ht9OayBLzqXDc1WX\npmWmjnOpn3BfAyw+3jHbbuvaJzPPAgvAG4ZRYKddB4/TPHPuorbmmXPsOni8FeSfehEeOt36bbDr\nMi47l6QB1HEu9RPu3fbA8wr6EBH3RcRMRMzMz8/3U98lTp1uDtQuLcW5pGGp41zqJ9xngXWLltcC\np5bqExErgVXAq51/KDN3Z2YjMxvT09NXVPANU5MDtUtLcS5pWC47lyq6ubKfcH8OuCki3hwR1wHb\ngAMdfQ4AH2k/vgf4fmZesuc+DNs3b2ByYsVFbZMTK9i+ecMoNqeCOZc0LEvNpX95y39XdnNlz3Bv\nH0O/HzgIHAP2ZuZLEbEzIu5qd/sa8IaIOAH8A/DpURV8981r+NwH3s6aqUkCWDM1yec+8HavcNDA\nnEsalqXm0p/99EuV3VwZI9rB7qnRaOTMzEwl25aksXhoii6nH4FoXfhxBSLiUGY2evUr5w5VSaqb\nCj+ryHCXpFGp8LOKDHdJGpUKP6vIb2KSpFGq6LOK3HOXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5J\nBTLcJalAhrskFaiyDw6LiHng51f5Z1YDvxxCOSVzjPrjOPXHcept1GP0pszs+YUYlYX7METETD+f\njnYtc4z64zj1x3HqrS5j5GEZSSqQ4S5JBVru4b676gKWAceoP45Tfxyn3moxRsv6mLskqbvlvucu\nSeqi9uEeEXdExPGIOBERl3zxdkS8LiK+3V7/bESsH3+V1etjnO6NiPmIeKH987Eq6qxSRDwWEb+I\niBeXWB8R8XB7DI9ExC3jrrEO+hin2yNiYdFcGv3XCtVMRKyLiH+PiGMR8VJE/H2XPtXOp8ys7Q+w\nAvgp8EfAdcCPgbd09Plb4NH2423At6uuu6bjdC/wSNW1VjxOfwHcAry4xPr3A08DAdwGPFt1zTUd\np9uB71ZdZ8VjdD1wS/vxHwL/1eU1V+l8qvue+63Aicz8WWb+Dngc2NLRZwvwjfbjJ4BNERFjrLEO\n+hmna15m/gB49TJdtgB7suVHwFREXD+e6uqjj3G65mXmy5n5fPvxr4BjwJqObpXOp7qH+xrg5KLl\nWS4dwN/3ycyzwALwhrFUVx/9jBPAB9tvD5+IiHXjKW1Z6XccBe+OiB9HxNMR8daqi6lS+1DwzcCz\nHasqnU91D/due+Cdl/f006d0/YzBk8D6zNwIfI8L73Z0gXOpP8/TugX+HcCXgP0V11OZiHg98B3g\nk5n5WufqLv9kbPOp7uE+Cyzew1wLnFqqT0SsBFZx7b2l7DlOmflKZv62vfhV4J1jqm056We+XfMy\n87XM/HX78VPARESsrrissYuICVrB/s3M3NelS6Xzqe7h/hxwU0S8OSKuo3XC9EBHnwPAR9qP7wG+\nn+2zGdeQnuPUcazvLlrHCHWxA8CH21c53AYsZObLVRdVNxHxxvPntSLiVlo58kq1VY1X+7//a8Cx\nzPz8Et0qnU8rx7WhK5GZZyPifuAgrStCHsvMlyJiJzCTmQdoDfC/RcQJWnvs26qruBp9jtMDEXEX\ncJbWON1bWcEViYhv0brSY3VEzAKfBSYAMvNR4ClaVzicAH4DfLSaSqvVxzjdA3w8Is4CTWDbNbhD\n9R7gr4GjEfFCu+0zwI1Qj/nkHaqSVKC6H5aRJF0Bw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQV\nyHCXpAL9P9RRuB58sM8iAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_grid(vals[0], nx, ny)\n", "plot_grid(vals[-1], nx, ny)" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " Once \n", " Loop \n", " Reflect \n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 50\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " plot_grid(vals[index], nx, ny, ax=ax)\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", " ax.axis('off')\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What if we give some of the particles a specific directional nudge? " ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [], "source": [ "nx, ny = 11, 12\n", "X, Y = build_grid(nx, ny)\n", "x0 = grid2linear(X, Y)\n", "x0[:ny] += 0.6\n", "x0[nx*ny:nx*ny+ny] += 0.2\n", "xd0 = np.zeros_like(x0)\n", "acc_func_grid = partial(acc_func_grid_template, nx=nx, ny=ny, k=1.)\n", "grid = VerletSolver(0.05, x0, xd0, acc_func_grid)\n", "ts, vals = grid.step_until(20, snapshot_dt=0.05)" ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJztnX9wVOd97p+vfqzRxtfGQXKEJO4l\n6WSc2oQYLDJurcQdUwb7kgDpxArNcMfTKSEzdiB4Oig4E8sZ0l4IdBriTOi1rpw2GTIB1bFBDnUd\njXOHRG4Tg4CLsR1PEucHSBCEcwVttUQSeu8fuwd2V+fsz/N9392j5/OPYBH7vOfdc77nPd8953nE\nGANCCCHVT43rARBCCAkHFnRCCIkILOiEEBIRWNAJISQisKATQkhEYEEnhJCIwIJOCCERgQWdEEIi\nAgs6IYREhDqbYo2NjWbhwoU2JQkhpOoZGhq6aIxpyvd7Vgv6woULcezYMZuShBBS9YjIrwv5PbZc\nCCEkIrCgE0JIRGBBJ4SQiMCCTgghEYEFnRBCIgILOiGERAQWdEIIiQgs6IQQEhGsPlhUyRw8MYzd\nL76JkbEEWuY2YOvK27B2SavrYZFTfcBL24FLZ4Gb24Dl3cDiTtejIrOdCt0vo1fQS5jogyeG8diz\nryIxeRUAMDyWwGPPvgoAlVXUXe1ELnWf3wxMJpJ/v3Qm+XegIg4eMkup4P0yWi0Xb6IvnQFgrk/0\nqb6c/233i29eK+Yeicmr2P3im/n1vrII+OLc5M88OmVR4rZVrS6QPIl4B43HZCL5Oql8bB4fNqng\n/TJaK/RcE53jzDkylijq9YMnhnHycA+6JvciLhPJF5XP0uMvdCNewraVw8ETw7j70OfRDHu66a2v\nX8w567/iuHRWTfOhG19BV/0BxBPn1a9GrLb5UldZ5tJZ/BaN2DHxII7dtEJP02cVO/7dR7Cr/zXc\nuWqj+pWv6twG7H/m0ll07PyB07ZttFboQQd6ngLQMreh4Ne99syGiX3Xi7mH0ln64IlhzBk/7/+P\nIRe3dM3Hnn0Vt5pRa7qe5vBYAgbAyPQ8/1+8uU1F86M1g8mTdOIctK9GsrfVa/MdPDEculb6VZbA\noBmj2FHfi7suD+hp+iyu4jKBDRP79DQB4FQfxr/8Pqw+eAcOjH8KH60ZDH9uA/a/ETPPzueZg2gV\n9KADPU8B2LryNjTU12a81lBfi60rb5vxu157pkUu+r+ZQqHb/eKbGDH6xS1bMzF5FSOm0Zpudutr\n11Qnxk0s85fqG5KrZgXNrro+ayfpk4d7MCCP4K0bPonB2GasrhksrM1XCgHFtauuT08z4Dhokbf1\nNFMnrnjiHGrEoK3mInbW94Y/t8u7k/thGgncgC9PZl7JqW1nDqJV0H0mupACsHZJK3b82fvROrcB\nAqB1bgN2/Nn7fS+XvDaMzUI3MpbwLW7jJhZqccvWBOwU1WxNj/7pDmyb3ICz040ABLh5AfDRJ0Nt\ngaRrWjtJn+pD1+RetNVcRI0go/AEtfnKIkdxBYJbi2WRYxWrppnjxBWq5uLO5H548wJ4++W2ib9E\n/3THjF9V2c4cRKuH7h3oJdyRsXZJa0H9rpa5DRhOFdid9b2ZKzqlQtcytwH9Yx3AZHIV2SJvY8TM\nQ29sPb6o1N/1trN/OlP3gjSi+aP/U6Wv7Gmm0z/dgaH4Cry87b7Q9bI1R0wj2vyKetgn6Ze2z7gS\n8ArPUHxFuFpAcvyXzsx42SuuQS3HsljendlDR3IBsmuqU08zz4krVM3FnRnHwLGdPwB8irfKduYg\nWit0IDnJj54GvjiW/Bly4fHaM+mrx2kjGG+YH/rq0U+zY+JJvOf338YK83XcuWpj6FrZmgCu6d4x\nvR8/XnNE7UvCYlpfGprWrkZyFB6VbfW5cvWKq9r8plax4w3zMW0EZ6cbsW1yA/qnO/Q0c1wV2NyP\nPLQ1/YjWCt0C3ip+94tv4vmx5OpR+9vsdE1b36DPRs3nxzrwzvqY/l0uASvmK/FmnW1Nu3K9dpfL\n5IMYumkFdmjO7+JOxBd3Ztxx0qr5mQZcFfTG1mPHKv8Wali42Hf9EGOMNbH29nbDCDoy68m+pQ9I\nrqCVrvBmFRX6BGe5iMiQMaY93+9xhU6Ibcr4rofkIau3PdtgQSfEBbO88BAdovelKCGEzFJY0Akh\nJCKwoBNCSERgQSeEkIjAL0VTMOCCEFLtcIUOy853YWHba9qWXlQ9tEl0qaB9Nm9BF5FviMgFETmd\n9to7RWRARH6W+nmL7jALpMSJLTngwhW2Qyds6bnYrgo5EEmV4jIAxoe8T4qKyIcB/AeAbxljFqVe\n2wXgd8aYnSKyDcAtxpjP5RNTfVK0jKfv3r3tMPxmQQD8cueqGa+7CERI1/y3OZ9FM3x8ym9ekPSv\nCVnzwPin0FbjZ1oVrh6+ssj3kfjQdQDf/WXcxLCr/mErAQzA9fltvzyAx2L/hHfhIsRSqIY1zVN9\nGH+hG3MS5zEynTSUU59fm5oB++x5NOHHa46Eplnok6J5V+jGmB8C+F3Wy2sAfDP1528CWFv0CMOm\njFioUgIubAYiZLeEbIROpGvaspU1OZJgQsdVAEMKb37vujyAHfW9aMYoxFKohjXNU32YOrQp6U+O\npD951+ReDD63VzXgwqpmwL55q7lYVQEX7zLGnAOA1M9bwxtSiZSYVgSUFnAB2AtEyG4J2fBiT9e0\n5f3+W/jrBL1eFi4CGNLw5tdmqIZ1zZe2o+7qlYyX4jKBLdivN7+2NXM4PEYy4EJENorIMRE5Njoa\nsLIMgxLTioDSAi4Ae4EI2Sb5Nmxe0zVt2crumHjQN8Rjx8SDoeoAcBPAkK6Ten+byVfWNXOcNNXm\n17ZmDmtiwH7ARakF/bciMh8AUj8vBP2iMabHGNNujGlvamoqUa4ASkwr8li7pBUvb7sPv9y5Ci9v\nuy+w95XehrG1cs1u/Xhe7OfRBK0kn3TNDO93JT0AOHbTigyPec9D+9hNCqEPeQ5E7WAC7/1tJl9Z\n18xx0lSbX9uaKd/382ia4fsOVE/ART+Ah1J/fgjAoXCGUwY+sVAaRcdFIIJfS2ig9t5k2IRykIdH\n/3QHVpivo3/Nayp6nuZA7b3XQjw6Jp7EQO290QlgSMObX5sRf9Y1l3djqnZOxkvjJoY9WKc3vy40\nF3fix2uO4I7p/eiYePJaMXcRcFHIXS7fAfAnABoB/BbAEwAOAugD8F8B/AbAg8aY7C9OZxAVP3TX\nd7nYevCJmrrwLhclXGhCdz8q9C4XBlwQQkiFE9pti4QQQqoDFnRCCIkILOiEEBIRWNAJISQisKAT\nQkhEYEEnhJCIwICLFAy4IIRUO9Vf0E/1JU2FLp0t+cEez4XOM6PyAi4AsKiTaBHC8UIql+puuYRk\nLl9WwEUlhCRELb2o2t+/UsdRYWEMkcTxvlXdK/RcHuhFrDqCHNGCXk9/fHpn7Gk04PfJf/AOECD0\nVU9gSyg7qCHEMfhq1r6spgcAR/ufwqLjj+vNqc98TR3ahL/ufw3f/I8P2mu3pXy7r1m9psZRB+it\nmEM6Xoql4uwclK5S1PfdAqjuFXoZHujplBpwsbWu7/qH56HgK50z87SMYI9SNMdf6FbR8zRbhnbp\nzqnPfNVdvYINE/us5smOv9A9w7e77uqV5PwqYTVAJEV2IMyB8U9h9cE7MP7l96mtXnMeL0pXKVb2\n3QKo7oJehgd6OqUGXNjylc7ZEgrppFao5pzEeRU9T3M+lOc0h1+2h41ggqB5DJzfELAaIJLC249W\n1wxiZ30v2mouokZMMuVLqd2T83hRWgBZ2XcLoLoLepke6B6lBlzY8pXO2RIK6aRWsOb0PBU9T1N9\nTvMEW6SPRZOgeQyc3xCwGiCSwptHm8lMOY8XpQWQlX23AKq7oIfogV5KwIUtX+mcLaGQTmqFavbG\n1qvoeZp+c5rADeHNaZ5gi/SxaNIbW+9bXHtj69U0rQaIpPDm0WYyU87jRWkBZGXfLYDqLuhAsng/\nelot6CGb9PaMl+QzbBphLIVqeFxrCSkFewRp3rlqo1qQiBdwkV50hk0jTi/9UnhzmjVf4w3z0W02\nXgslAOwEE9y5aiO6zcaM4tptNibnVwmrASJpmg31tVZXrzmPF6UFkJV9twDoh14CFfetPTWrTtOV\nrivNk4d70DW5N7PtUt+gsgDyNG3f5cKAC0LI7IEPNZVMoQW9uu9DJ4RUD4s7WcCVqf4eOiGEEAAs\n6IQQEhlY0AkhJCKwoBNCSERgQSdkNlEpTpNEBd7lkoIBFyTyKDpzksqAK3TkcWcjJCooGVORyqGs\ngi4ij4rIayJyWkS+IyJzwhpYwYRwCcmAC+IMm5+dkjEVqRxKbrmISCuAzQBuN8YkRKQPwDoA/xjS\n2PIT0iUkAy6K1ASq8vHpStIEHAQi3NyW8gH3eV2LU30Yf6EbcxLnMTI9D72x9bhz1Ubrn+me23+G\nZb/4mupTqumaD934CrrqDyCeOG/1qdhyWy51ABpEpA5AHMBI+UMqgpAuIRlwUaSmYkiAn+bR/qfU\nVrEHTwxj8Lm9ODD+Kfzihk/iwPinMPjcXvV2m4tAhKN/sAmJbDdAE8PRP9ikouelMsUT51ADg7aa\ni+ia3Ks+v9n70V2XB7Bo6Auq0XvZQR5dk3uTnu+Wo/5KLujGmGEAfwvgNwDOAbhkjPl+WAMriJAu\nIRlwUaSmYkhAtuaKq0ew6PjjagfjycM92C49qeAFoK3mIrZLD04e7gnl/YNwEYiw5fX34nNZ9rmf\nm9yALa+/V0UPL22fkcoUlwlswX7VAJHs/airrg8Nyl7s6Zo2vd+zKbmgi8gtANYAeDeAFgDvEJEZ\nZs4islFEjonIsdHR0dJH6kdI3sYMuChSUzEkIJsu5augDRP7Zhx8cZnAhol9obx/EC4CEUbGEuif\n7siwz+2f7tAL88iRDqUZIJL93jYWXumaNr3fsymn5fKnAH5pjBk1xkwCeBbAH2f/kjGmxxjTboxp\nb2pqKkPOhxC9jRlwUYSmYkjAjNeUD46WmreLej0sXAQiFNNaDIUc6VCaASLZ723jxJmu6TK5qJyC\n/hsAd4tIXEQEwHIAb4QzrAJRCnfIxWwPuNAOCcjWPBeUdxnSwXGlobmo18PCRSBCMa3FUFjejana\nzBvfxk0Me7DOSqiGx66pzhnfHYS98ErXtLXQ86Pku1yMMT8RkWcAHAcwBeAEAN3Gox+WLTm9lbv3\nbfbQTStwdOVnVL+1z9accSeGwhzk1kxphXyXi5/myO1daH31icyefYgHR/yB7Zg6tCmj1ztVOwfx\nB3T7nde3NYYPjXVYubsm734UNos7UQfMuMulQ/kuF79j9PTtC1XvcknXfH6sA++sjzm5y4UBF6Ty\n0Q5GYPACqXCYWEQIIRGh0ILOR/8JISQisKATQkhEYEEnhJCIwIJOCCERgX7oKeiHTgipdljQcd1Y\nx/Ni8EyhALCoE0KqBrZcUKYfOiHVCD30I0l1F/SQdspi/dDD1g+FShoLqWyU7I8rDhfHhOPjsHpb\nLiEGO7TMbcCwT/EOMhA6eGIYJw/3JD2PPac+5XCCvGETCiEXvpq1L6s+VZkeHvJY7J/wLlyEKD+9\n6er7E1dhHncf+jyaEWB/XO4c+zx1e/DqPda303p4iCPNbKp3hR6iJ3cxpkVev93PdlXL8zhv5qmC\nP7mf5uBzezF1aJPays7TvOvyAHbU96IZoxDlFaSrPNmCdENe7Xmat5oAG+tyHSx9Vv5ThzZh8Lm9\nVufXRXiIC00/qregh+jJXYwfutdvt+l5nLfHr+BP7qe5BftnBBaEucN6mjYDArK3c3XNIAbkEaw+\ndIfqJXPez1ShLeJpqtm7+iws6q5ewRbsz3hN+/spF+EhLjT9qN6WS8j5iGuXtBZ0Gej11UdMI9r8\nirqC53HeHr9CVqSfpvZJzNO0ebJM387VNYPYWd9rpY2W9zPNddVV4li899411Zm5nUA4DpY5Ai2C\nxqLByFgCIzF7x6crTT+qd4Wu5MmdD6+vbtPzOG8wgcJc+GlqG/d7mjYDAtK30+aVQd7PVOGqy3tv\nz8ff82E/j6ZwfPxzBFoEjUUDV+EhtjX9qN6C7iDcArjeb88+KMYb5qvp5+3xK8yFn+YerJsRWBDm\nSczTtHmyTN9Om1cGeT9ThVSo7HCWjoknccf0fvx4zZFw9lufhcVU7RzswbqM11RDNeAuPMS2ph/V\n23IBrIdbADON7IfiKyojmCDkufDT7Fj5MOpqP6B2l0t66MNjl2HlLpeM7Ry3d8mc9zNd3p155xJQ\n9klNPeBi8czgk7rl3ei4eg/+zeJdLm7DQ+xp+kE/dEI8sm//BJJF1MKVX+B4GLxBULgfenWv0AkJ\nE58VptMi6uAKlFQ3LOiEpMMiSqqY6v1SlBBCSAYs6IQQEhFY0AkhJCKwh56CAReEkGqHBR0MuCCE\nRAO2XBCBgAt6oRNCUGZBF5G5IvKMiPxURN4QkT8Ka2CBKBSvqg64mC1hBYSUio3jtBJqAcpvuXwV\nwL8YYz4uIjEA8RDGFIxSkENVB1wouPLl1VR8gpEBF/oBF7MptKT98gB2xp7WDZ3wqUuJZz+D07/6\nf1i2+tPhaBRIyY/+i8hNAP4vgPeYAt+k7Ef/v7IowCZ2AfDo6ZLfNruHDiQNhPw80b3fHZBH0Fbj\n5/tR3liKHt+hOwD4Tb8AXxwLXfNby36NZa8+ofJ4vKe54uoRf3tXhUfw07dzdc0guur60CJv40q8\nGfEHQkjwKUDXI2OfUzhp+ml+PPav2Fnfm+lzH+Jcu/5MB2Ob9Y/TgLo0bBpxdO0PQzlJF/rofzkt\nl/cAGAXwDyJyQkR6ReQdZbxffhQsRYEqD7hQcOXLpbng+O7Q05GyNV0EXHhe6G01F1EjBvHEOdXW\nVc7PVKmNNhtDS6wcpwHvNR9vW/8erpyCXgdgKYC/N8YsAfCfALZl/5KIbBSRYyJybHQ0IPqqUJSK\nF5As6i9vuw+/3LkKL2+7L/Csmh5woTWWIE3f15V84YM01eLL4DbgwmbBSdf1fV0hUjBIM+qhJVaO\n0xw+8JpBHn6UU9DPAjhrjPlJ6u/PIFngMzDG9Bhj2o0x7U1NTWXIwVmoRToVF3Ch5AsfpHlBAj7D\nEA4QlwEXNgtOuq7v60pXorMxtMTKcbq8Oxlkkca4iWHXVKdqkIcfJRd0Y8x5AGdExHOqXw7g9VBG\nFYSjUIt0Kjbg4tHTyZ75o6dDGUOQ5pmlW9VOqi4DLmwWnHTddK59pkpXorMxtMQ7TodNI4xWzVjc\nidNLv4Rhk6wFZ6cbsW1yAwZq71UN8vCj3LtcNgH4duoOl7cA/EX5Q8qDYze8ig24sKS5bMn9wMJb\nVO6KcBlw0Xt4feZdS4Dq1V/Oz7Q2/HCLIM3Ih5aMJTB00wocXfkZ1eNl2epP4+CCj2TM7Q4GXBDi\nkEoKlKiksRDnFHqXCws6IYRUODZuWySEEFJBsKATQkhEYEEnhJCIwIJOCCERgX7oKRhwQQipdljQ\nwYALQkg0qO6WS0gexCUHXFSIBzIhhADVvEIP0Ru9pIALJW/2onHxAAofetHH1hzzs4wU1VvQQwx2\nKCXg4u5Dn0czdIIlgjR9gwkUTyouNV2FIdjSTNf1DRBRmuN0zYdufAVfMP/run2u0qKkYoI8NDVP\n9WH8hW7MSZzHyPQ89MbW485VG623bKu35RKiI11e86s0vH67po1skObwWAIG13v84y90q3mTu9S8\n6/IAdtT3ohmjEOVYPRea6brZ83vwxLCafW625oaJfape6Omad10ewI9im/GjxMew7OCHcbT/qdA0\ngjR951aDU32YOrQJ8cQ51MCgreYiuib3YvC5vXqaAVRvQQ/Rka6UgAub7nxBPf45ifP+/yGEk4pL\nTRdhCLb90HN+b6Nkn5utacMyePeLb15LK0qGhwCtchGLjj/uJjxEg5e2zzgxxmUCW7C/qgIu3BKy\nN3qxARc27UADe/zT8/z/QwgnFZeaLsIQbPuh5/zeRsk+N1vTxqJkZCzhe7JswO+thYesrhnEYOrq\nQOXmhYB9pEXerqqAC7c48kb3+urZfujn0aSmH9TL742tV/Mmd6npIgzBth96zoALpSCXbE0bi5KW\nuQ1Ow0MyowWh00rLkVhUNQEXFYFCsEM+ss3zOyaexB3T+/HjNUfU9IN6/Heu2qh2UnOp6SIMwaZm\num461763UVqsZGv2T3eg22xMhrMoLYq2rrwN5+AuPMRKK21594yQkHETwx6ssx5wQfvcEpgV39o7\n1pzVd7ko4kLzaP9TWHT88WSbxaO+QfWK2tvOHyU+llyZz0CSC8GwUL7LhX7ohJDKwdX97l9ZlGyz\nZHPzguRVfZVQaEGv3vvQCSHVg6voyOU6cX6VSnX30AkhJBcVECxvE67QCSHRxnGwvE24QieEkIjA\ngk4IIRGBLZcUDLgghFQ7LOhgwAUhJBpUb8slxHAJ62Y+YcOgDUIIQijoIlIrIidE5HthDKggPL/o\nS2eAEKxOSw64qIQiGvJc5NWysc0u57ZSPtdKGwupCsJouXwWwBsAbgrhvQojxHALoPiAixmPMltI\nKwrs8Yc8F0GatsIQXM5t++UB7Iw9bU075/c2SiEXR/ufwoLju3GrGcUFacKZpVuxbPWny9yS3Mwm\niwPbc5tNWSt0EWkDsApAbzjDKZCQ/aKLDbhoGdqV6UsBqHpn5zTsV/LOdhWG4HJut9b1WdPOG8Kg\nEHJxtP8pLBr6ApoxihoBmjGKRUNfCC9swueKwnrYBBwEXMDC3BZIuS2XPQC6AEyHMJbCCdkvutiA\ni/mwaweas8ev5J3tKgzB5dzatHnN+72Nwol6wfHdaMj2JZcJLDi+u+T3vEZA6+/k4R7r30+5+E5M\ndW6LoOSCLiIfAXDBGDOU5/c2isgxETk2OhoQ21YsCn7RxQRc2PbOztnjV/LOdhWG4HJubWrn/d5G\n4UQdFJt4qwk4kRVDwBXFhol9vr+uGfwwMpa4Fmrx1g2fxGBsM1bXDKpqqs5tEZSzQr8HwGoR+RWA\n/QDuE5EZn54xpscY026MaW9qaipDLg2H/gwtcxt8vbMTuEHN8CdnGILSXLgKQ3A5tzY90XN+poDK\nifqC+B9/FyTgRFYMQak9NW/7v64Y/PDQja9khFq01VzEzvpePHTjK2qaqnNbBCUXdGPMY8aYNmPM\nQgDrAPzAGLM+tJHlw0G4BZDstw/U3puRVjRsGnF66ZesB1xc6/ErzIWrMASXc+ulUA2bRhjlhUJB\nn2nIJ+ozS7cikX2yNDGcWbq15Pe8RsCVw5WG5oK/nwqLrvoDM0It4jKBrvoDapqqc1sEfLCoSLxW\nzO4XY/jQWIeVb9Cva9r71t5Ps2Plw4gv+RsLmu7mduimFTi68jPqd0QU9JmGbCq1bPWncRRI3Ylx\nERekEWfuCulOjACb2vgD27Hj6vut7rvxgCDzoNfDQHVui4ABF4SQcHAVYpFNREIt0mHABSHELpVi\nUzvLQi3Sqd5H/wkhxI9ZFmqRDlfohJDoUSlXC5bhCp0QQiICCzohhEQEtlxSMOCCEFLtsKCDAReE\nkGhQXS0XJX/oqg+4iDL0BCekYKpnha7kDw0UH3BRUZ7Sig9z+GrWvqz68Ih1D/ZTfRh/oRtzEucx\nMj0PvbH1uHPVRitXZrPFJ3y2aFaCH3r1PCmq+PTXPTt/4Btw0Tq3AS9vuy/jNc/3ON0qM2FiOH3X\nX6t9eNktISDph/GtZb/GslefmPkARQj33Pppfjz2r9hZ35vpix6Snp/mYGwz2mp83OrCeuLvVB+m\nDm3K2J5xE0O32YiOjz2sWgCCPtMg2+aq00wtNMylsxgx8/DlyU70T3foaqZwMbfadaHQJ0Wrp+Wi\nFOQAFBdw4cL3OKgltOD47tBDEHJpbsF+1ZAL6x7sL22fsT1xmcAW7Fdvt3nbmm7zOiCP4OThHjXN\nk4d7MCCPZFjKqrQW07zRBQatknQ7XF0zCMCuH7o3v6/VfAJ3H7pXrWVX9X7o1lEKcgCKC7hw4Xsc\n1PoJGksYBc9PU7vAWvdgD7J8lbdVvbOB657d2TavXZN71fJguyb3zrCUVfEJ9/FGj8sEuuqub5e2\nHzqAGfPbjFG1vN0o+KHbRSnIwaPQgAsXvsdB3tFBYwmj4PlpahdY6x7sAeMeMfNU/bqB5LZ21fX5\n2ryqxO29tN3fUrauL/xtzXGivPZnxfn13ttvfrUiBaveD906FeLP4ML3OKgldGbpVrWTnJ/mHqzD\nVO0cFT0/TXUP9uXdM7Zn3MSwB+tU/bqB5LamF7gMNOL2chTZ0Lc1x4kS0PdD9/Yjm5GC9EMvhQrw\nZ3Dhexzknb1syf3AwltU7joJ8kOvq/2A2l0u1j3YF3eiDphxl0uHhbtc1i5pxfj3mxFPnJv5jxpx\neze3+d5UcCXeHP62+rgdJnADdk91otWix/2FQ03JNks2CvNLP3RCZjvZt+ICod415EzL03PtjW57\nmxWhHzohlY5XVGwUPptanp7roml7mysArtAJIaTCid596IQQQnLCgk4IIRGBBZ0QQiICCzohhEQE\n3uWSggEXhJBqhwUdDLgghEQDtlxQRsCF6/AF1/qEkIqi5BW6iCwA8C0AzQCmAfQYY74a1sCKosyn\n0ooNuDh4YhgnD/ega3LvdfMfjfCFLM30ltCe23+W6YWuoO8i4KKiwkOUcRX8cPJwDzZM7ENLzdu4\n0tCM+APbVR+2caU5G0I1simn5TIF4K+MMcdF5L8AGBKRAWPM6yGNrTBCSDJqmdvgG3Dh5wjntWcG\nZB/iNQFObiHvqH4toZahXYAEeKErhE0MjyUw+NxefCQ94CLkk0hGSEDK7vTmoS/gKBB+UU8LYFhm\n5uGuyU4Mo8Nau81Fm+/giWEMPrcX26Xn2r4bT5xLhnwAagsRF5q25/Zo/1NYNrQLP8JFjMQasety\nJx57dkJV04+SWy7GmHPGmOOpP/87gDcA2G84+3gvF2uRWUzAhdeesenk5tcSmg9dfRcBF9ZCAhwH\nMABuAi52v/gmtmD/DEvZuqtXdCx7HWpm77srrh5JBlxotCdP9WHR8cfRKple8yuuHrGeSxxKD11E\nFgJYAuAnPv+2UUSOicix0dGoP0tdAAALfElEQVSAQIZyCCHJqJiAC68Nox6+4KOZ8ZqyvouAC2sh\nAY4DGLz3txpwkdK0uRBxqZmON89J50Vz/coyrHl+aTsa8PuMl7z9SXs/yqbsgi4iNwL4LoAtxpjL\n2f9ujOkxxrQbY9qbmgICGcohpCSjQgMuvDaMeviCj2Y6u6Y6kcANavouAi6shQQ4DmDw3t9qwEVK\n0+ZCxKVmOupBFzn2J+39KJuyCrqI1CNZzL9tjHk2nCEViXKSUTZee6Z/ugPbJjfg7HQjpo0kQxiU\nbDn9WkIDtffi9NIvqQV+uAi4sBYS4DiAAXAQcJHS3IN1MxYiU7VzVI8XF5rp+676FULA/nQO89T3\no2zKuctFADwN4A1jzN+FN6QisWyRmR7C8PxYB4biK6wZ9vsGXEDnDhAXARfWQgIcBzAADgIu4H2m\nD2PX4Tprd5y407y+714Q5aCLgP1p5K4u63e5lGyfKyIdAH4E4FUkb1sEgM8bY/456P/QPpdUDAxg\nmD3YmGfl/alQ+1z6oRPikko4scwGqnyeWdAJISQiMOCCEEJmGSzohBASEVjQCSEkIrCgE0JIRKAf\neopKcEojhJByiE5BL+O2JAZcEEKiQDRaLmnueaWY75QccOFpM2SCVBvcbyNJNFbouSx0C1illxJw\nsfvFN9F+eQA7Y09fd1pTDLmYFYb9p/ow/kI35iTOY2R6Hnpj63Hnqo3R204f3YdufAVd9QcQT5xX\nffClEsJZIrnvOtLMJhoFvUwL3VICLhKTV3Eg1jfDNlMj5CJdc3XNILrG+9By8O2kF4iSJ4anueLq\nERyI9aElcRHnDjbi6JkunQShU32YOrQJ8ZTfumcl2/3cFICH1UMfDmA/Wm64iJHxRux5bp2qpqeb\n8ZlO9iI+pVtgXYSzuAh+cBUeUglt22i0XMq00C0l4AKw4OKWpZnpn22Sxk5h+jpnaa64eiTDr7tV\nLmLR8cd1Ls9f2j4jPCMuE9iC/aohAScP92C79GR4km+XHtWgCSBzP1K3d83StOZP7iL44VQf7j50\nL16r+QQGY5vthJa40AwgGgW9TAvdUgIuAHshF56mrQPf0/TTa8Dvdfy6c3hKa4YEbJjY5+tJvmFi\nn5omkLkf2Sqw1sNZbAc/pL5La8ZoxgnEK7CR0cxBNAr64s6kc1oZ3uDFBlwA9kIuPE2byS8tcxvs\nJs3k8CjXDAloqfH3JA96PTTdtG2yVWCth7PYDn7Ik0QVGc0cRKOgA8ni/ehp4ItjyZ9KTmrp7Rkv\n5GLYNMIohExka9pMftm68jacg8WkmeXdM8Izxk0Me7BONSTgSkNzUa+HRfp+ZKvAWg9nsR38kOME\nohZa4kIzB9Ep6JbIbs8M3bQCR9f+EKJ4IvE0e2PrrcXerV3SipG7ulRj7jJY3Im6NV/DeMN8TENw\ndroRu+ofRsfHdL+cjD+wfcaJZKp2TvLLZkXS96Pnpzuwq/7hZGFVXBhka34i/r/Rv/Y1xD/3U50F\nkE8rVDX4IeAEckEaA1uoVamZA9rnVhu2fZ2r3Ee6IGbDNrrC5ty6CAyxpEk/dELI7MPFydmCJgs6\nIYREBAZcEELILIMFnRBCIgILOiGERAQWdEIIiQjRMOcKgUpwSiOEkHKITkFnwAUhZJYTjZaLy4AL\nEk1sBUAwaIKESFkrdBG5H8BXAdQC6DXG7AxlVMXiKOAi6ub5s1Vzz+0/w7JXn7i+Tyn5kx/tfwqL\njj9uJSAFwKwKEJktmtmUXNBFpBbA1wGsAHAWwFER6TfGvB7W4ArGcsCF7UAEV5qVEhLQeuZ7WPaL\nr6k8ieen2TK0C5DSFwiF6i4b2oUG0Q9IAeAmQCR1Alk9fh7tZh521XSif6zDyn40G44XP8ppuXwQ\nwM+NMW8ZYyYA7AewJpxhFYnFgAsXgQguNNPbUKtrBjEY24zXaj6Buw/dq9YW8Gt9rbh6JBmqUWI7\nrRTN+dC3Dd794ptWdK5hO0Ak1QaNJ86hRkyGT7h2O9P18eJRbQEXrQDOpP39bOq1DERko4gcE5Fj\no6OjZcjlwGLAhYtABBeaXrspMyUJaMaoWkqSX4urqy5HzJ+Spg2b4pGxhFU7ZOsBInl8wqMWWpJ+\nvAzGNuOtGz6JwdhmtF8eUNP0o5yCLj6vzTCGMcb0GGPajTHtTU1NZcjlwGbAhYNABCeaqXaTzZQk\nvxaXdsiGn+auqU512+CWuQ2+PugJ3KBjT2w7QCTHCQTQDX5wdbxkL37aai5iZ+xpq190l1PQzwJY\nkPb3NgAj5Q2nDCwFXLgIRHCh6bWhbKYW+bW+tEM2/DQHau/F6aVfKmuBUIjuQO29GUETw6Yxqavk\nTW41QCTHCUQ7+MHV8fK5eouRjQGUc5fLUQDvFZF3AxgGsA7AJ0MZVQUTf2A7pg5tyuhHagciuND0\nrlAuHGpKtlmyUWgLeJrpdwqM3N6F1vQ7ToBQV8t+mltX3oZlS+4H8OlQNHLrxvChsQ79uyIWd6IO\nmHGXS4fWXS7Lu2f4hI+bGHpj67FjlW7wg6vjxRwKuALQ+E4kgLLsc0XkvwPYg+Rti98wxvxNrt+P\njH1uRD2XA3Vthwb4jYEBFNWHy8/NhfZXFqW+vM/i5gXJrkEZ0A+dhAcLKiH5UVz8FFrQo/PoP9Fj\ncScLOCH58I4Rh4sfFnRCCAkLx4ufaHi5EEIIYUEnhJCowIJOCCERgQWdEEIiAgs6IYREBBZ0QgiJ\nCCzohBASEVjQCSEkIlh99F9ERgH8WuntG4GgxIBIwe2MHrNlW7mdpfPfjDF5/cetFnRNRORYIV4H\n1Q63M3rMlm3ldurDlgshhEQEFnRCCIkIUSroegmwlQW3M3rMlm3ldioTmR46IYTMdqK0QieEkFlN\nJAq6iNwvIm+KyM9FZJvr8WggIgtE5P+IyBsi8pqIfNb1mDQRkVoROSEi33M9Fi1EZK6IPCMiP019\nrn/kekwaiMijqX32tIh8R0Tm5P9f1YGIfENELojI6bTX3ikiAyLys9TPW2yNp+oLuojUAvg6gAcA\n3A7gz0XkdrejUmEKwF8ZY/4QwN0AHonodnp8FsAbrgehzFcB/Isx5n0APoAIbq+ItALYDKDdGLMI\nyfzhdW5HFSr/COD+rNe2AXjJGPNeAC+l/m6Fqi/oAD4I4OfGmLeMMRMA9gNY43hMoWOMOWeMOZ76\n878jefDrRac7RETaAKwC0Ot6LFqIyE0APgzgaQAwxkwYY8bcjkqNOgANIlIHIA5gxPF4QsMY80MA\nv8t6eQ2Ab6b+/E0Aa22NJwoFvRVAetT2WUS00HmIyEIASwD8xO1I1NgDoAvAtOuBKPIeAKMA/iHV\nWuoVkXe4HlTYGGOGAfwtgN8AOAfgkjHm+25Hpc67jDHngORCDMCttoSjUNDF57XI3rojIjcC+C6A\nLcaYy67HEzYi8hEAF4wxQ67HokwdgKUA/t4YswTAf8LipbktUv3jNQDeDaAFwDtEZL3bUUWXKBT0\nswAWpP29DRG6pEtHROqRLObfNsY863o8StwDYLWI/ArJ9tl9IrLP7ZBUOAvgrDHGu8p6BskCHzX+\nFMAvjTGjxphJAM8C+GPHY9LmtyIyHwBSPy/YEo5CQT8K4L0i8m4RiSH5hUu/4zGFjogIkv3WN4wx\nf+d6PFoYYx4zxrQZYxYi+Vn+wBgTuRWdMeY8gDMiclvqpeUAXnc4JC1+A+BuEYmn9uHliOCXv1n0\nA3go9eeHAByyJVxnS0gLY8yUiHwGwItIfoP+DWPMa46HpcE9AP4HgFdF5GTqtc8bY/7Z4ZhIeWwC\n8O3UQuQtAH/heDyhY4z5iYg8A+A4kndqnUCEnhgVke8A+BMAjSJyFsATAHYC6BORv0TyhPagtfHw\nSVFCCIkGUWi5EEIIAQs6IYREBhZ0QgiJCCzohBASEVjQCSEkIrCgE0JIRGBBJ4SQiMCCTgghEeH/\nA+oppgxJcWhFAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_grid(vals[0], nx, ny)\n", "plot_grid(vals[-1], nx, ny)" ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " Once \n", " Loop \n", " Reflect \n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 64, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 40\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " plot_grid(vals[index], nx, ny, ax=ax)\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", " ax.axis([-1, nx+1, -1, ny+1])\n", " ax.axis('off')\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another version might be to add random movements to the particles." ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [], "source": [ "nx, ny = 11, 12\n", "X, Y = build_grid(nx, ny)\n", "x0 = grid2linear(X, Y)\n", "x0 += np.random.randn(*x0.shape) / 10\n", "xd0 = np.zeros_like(x0)\n", "acc_func_grid = partial(acc_func_grid_template, nx=nx, ny=ny, k=1.)\n", "grid = VerletSolver(0.05, x0, xd0, acc_func_grid)\n", "ts, vals = grid.step_until(20, snapshot_dt=0.05)" ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " Once \n", " Loop \n", " Reflect \n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 68, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 40\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " plot_grid(vals[index], nx, ny, ax=ax)\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", " ax.axis([-1, nx+1, -1, ny+1])\n", " ax.axis('off')\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Lennard Jones particles " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Trying to mess with the png encoding in to_jshtml " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Simulating the spring." ] }, { "cell_type": "code", "execution_count": 320, "metadata": {}, "outputs": [], "source": [ "x0 = np.array([0, -1], dtype=np.float)\n", "xd0 = np.array([0.0, 0.1], dtype=np.float)\n", "spring = VerletSolver(0.01, x0, xd0, acc_func_spring)\n", "ts, vals = spring.step_until(12, snapshot_dt=0.1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Making the animation." ] }, { "cell_type": "code", "execution_count": 321, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "FRAMES = 50\n", "def animate(i):\n", " index = np.argmin((np.linspace(0, 1, vals.shape[0]) - i / FRAMES)**2)\n", " ax.clear()\n", " ax.plot(vals[index, 0], vals[index, 1], 'o')\n", " ax.plot(vals[index-5:index+1, 0], vals[index-5:index+1, 1], '-')\n", " ax.set_xlim(-1, 1)\n", " ax.set_ylim(-1, 1)\n", " ax.set_title(f\"t = {ts[index]:.2f}\")\n", "\n", "animation = mpl_anim.FuncAnimation(fig, animate, init_func=lambda: None,\n", " frames=FRAMES, interval=100, blit=False)\n", "plt.close(fig)" ] }, { "cell_type": "code", "execution_count": 322, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 322, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plt.rcParams['animation.html_args']" ] }, { "cell_type": "code", "execution_count": 323, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'h264'" ] }, "execution_count": 323, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plt.rcParams['animation.codec']" ] }, { "cell_type": "code", "execution_count": 324, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 324, "metadata": {}, "output_type": "execute_result" } ], "source": [ "HTML(animation.to_html5_video())" ] }, { "cell_type": "code", "execution_count": 325, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 325, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mpl_anim.HTMLWriter(fps=10, extra_args={'frame_format':'jpeg'})" ] }, { "cell_type": "code", "execution_count": 326, "metadata": {}, "outputs": [], "source": [ "mpl_anim.Animation.to_jshtml = mpl_anim" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import tempfile\n", "HTMLWriter = mpl_anim.HTMLWriter\n", "import os\n", "\n", "def to_jshtml_jpeg(self, fps=None, embed_frames=True, default_mode=None):\n", " \"\"\"Generate HTML representation of the animation\"\"\"\n", " if fps is None and hasattr(self, '_interval'):\n", " # Convert interval in ms to frames per second\n", " fps = 1000 / self._interval\n", "\n", " # If we're not given a default mode, choose one base on the value of\n", " # the repeat attribute\n", " if default_mode is None:\n", " default_mode = 'loop' if self.repeat else 'once'\n", "\n", " if hasattr(self, \"_html_representation\"):\n", " return self._html_representation\n", " else:\n", " # Can't open a second time while opened on windows. So we avoid\n", " # deleting when closed, and delete manually later.\n", " with tempfile.NamedTemporaryFile(suffix='.html',\n", " delete=False) as f:\n", " self.save(f.name, writer=HTMLWriter(fps=fps,\n", " embed_frames=embed_frames,\n", " default_mode=default_mode,\n", " extra_args={'frame_format':'jpeg'}))\n", " # Re-open and get content\n", " with open(f.name) as fobj:\n", " html = fobj.read()\n", "\n", " # Now we can delete\n", " os.remove(f.name)\n", "\n", " self._html_representation = html\n", " return html" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "toc": { "colors": { "hover_highlight": "#DAA520", "navigate_num": "#000000", "navigate_text": "#333333", "running_highlight": "#FF0000", "selected_highlight": "#FFD700", "sidebar_border": "#EEEEEE", "wrapper_background": "#FFFFFF" }, "moveMenuLeft": true, "nav_menu": {}, "navigate_menu": true, "number_sections": true, "sideBar": true, "skip_h1_title": false, "threshold": 4, "toc_cell": false, "toc_position": {}, "toc_section_display": "block", "toc_window_display": true, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 2 }