{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## $N$ interacting masses in $\\cos$ potential\n", "Additional requirements for this example: scipy\n", "\n", "### Sine Gordon equation\n", "\n", "Let $x$ be an array with position of $N$ masses in $\\cos$ potential linearly coupled to nearest neighbours. Then the Newton's equations read:\n", "\n", "$$\\begin{cases}\n", "\\displaystyle\\frac{dx}{dt} = v \\\\\n", "\\displaystyle\\frac{dv}{dt} = - M x + \\sin(x),\\;\\;\\;(1)\n", "\\end{cases}$$\n", "\n", "where M is a discrete Laplacian. \n", "\n", "On the other hand this is an aproximation to a sine-Gordon equation: \n", "\n", "$$x_{tt}- x_{\\xi\\xi} + \\sin x= 0$$\n", "\n", "which has the following soliton solutions:\n", "\n", "$$x_\\text{soliton}(\\xi, t) := 4 \\arctan \\exp\\left\\{ \\frac{ \\xi - v t }{\\sqrt{1 - v^2}}\\right\\}$$\n", "\n", "Below we solve numerically the system (1) and animate results in 3D.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.integrate import odeint\n", "import numpy as np " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 153\n", "M = np.diag((N-1)*[ 1.0],-1)+np.diag((N-1)*[ 1.0],1)+np.diag(N*[-2.0],0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "y0 = np.zeros(2*N)\n", "x = y0[:N] # Fist N variables are positions\n", "v = y0[N:]\n", "ksi = np.linspace(-5,25,N, dtype=np.float32)\n", "h = np.diff(ksi)[0]\n", "\n", "instanton = lambda x,v,t: 4 * np.arctan(np.exp( (x-v*t)/np.sqrt(1-v**2) ))\n", "v1,v2 = 0.4,0.05\n", "x[:] = instanton(ksi,v1,0)\n", "v[:] = -v1/h*np.gradient( instanton(ksi,v1,0) )\n", "\n", "x[:] += instanton(ksi,v2,50)\n", "v[:] += -v2/h*np.gradient( instanton(ksi,v2,50) )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def lhs(y_,t):\n", " y = y_.copy()\n", " x = y_[:N]\n", " v = y_[N:]\n", " y[:N] = v \n", " y[N:] = -np.sin(x) + 1/h**2* np.dot(M,x) \n", " y[0] = 0\n", " y[N-1] = 0 \n", " return y\n", "ts = np.linspace(0,140,50) \n", "%time xt = odeint(lhs,y0, ts).astype(np.float32)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import k3d\n", "import numpy as np\n", "import time \n", "plot = k3d.plot()\n", "r = .6\n", "origins = np.vstack([ksi,np.zeros(N),np.zeros(N)]).T.copy().astype(np.float32)\n", "vectors = np.vstack([np.zeros(N),r*np.sin(xt[0,:N]),r*np.cos(xt[0,:N])] ).T.astype(np.float32)\n", "\n", "vector_plot = k3d.vectors(origins, vectors)\n", "line_plot = k3d.line(vectors + origins,color=0xff0000)\n", "\n", "plot += vector_plot\n", "plot += line_plot\n", "plot.display()\n", "\n", "def update_plot(xx):\n", " vectors = np.vstack([np.zeros_like(xx),r*np.sin(xx),r*np.cos(xx)] ).T\n", " vector_plot.vectors = vectors.copy()\n", " line_plot.positions = vectors+origins" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot.camera_auto_fit = False\n", "plot.grid_auto_fit = False" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time \n", "import time\n", "for xx in xt[:,:N]:\n", " update_plot(xx)\n", " time.sleep(0.05)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from ipywidgets import widgets,interact\n", "\n", "@interact(ith = widgets.IntSlider(min=0,max=(ts.size-1)))\n", "def draw(ith):\n", " update_plot(xt[ith,:N])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }