{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "<center><h2>KdV Solitons</h2></center>\n",
    "<p>\n",
    "Solitons arising from explicit solutions of the  Korteweg and deVries (KdeV)  equation (nonlinear and dispersive):\n",
    "\n",
    "$$\n",
    "\\frac{\\partial u(x,t)}{\\partial t} + \\varepsilon u(x,t)\\,\n",
    "\\frac{\\partial u(x,t)}{\\partial x} + \\mu  \\, \\frac{\\partial ^3\n",
    "u(x,t)}{\\partial x^3} = 0. \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "\"\"\"From \"COMPUTATIONAL PHYSICS\" & \"COMPUTER PROBLEMS in PHYSICS\" \n",
    "by RH Landau, MJ Paez, and CC Bordeianu (deceased). Copyright R Landau, \n",
    "Oregon State Unv, MJ Paez, Univ Antioquia, C Bordeianu (deceased), \n",
    "Univ Bucharest, 2020. Please respect copyright & acknowledge our work.\"\"\"\n",
    "\n",
    "# KdVsolitons.ipynb\n",
    "\n",
    "% matplotlib notebook\n",
    "\n",
    "from numpy import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "\n",
    "ds = 0.4                              # Delta x\n",
    "dt = 0.1                              # Delta t\n",
    "max = 2000                            # Numb t steps\n",
    "mu = 0.1;                             # Mu from KdeV equation\n",
    "eps = 0.2;                            # Epsilon from KdeV eq\n",
    "mx = 131                              # grid in x max\n",
    "fac = mu*dt/(ds**3)                   # combo\n",
    "u = np.zeros( (mx, 3), float)         # Soliton amplitude \n",
    "\n",
    "def init():    # Initial wave\n",
    "   for  i in range(0, mx):      \n",
    "      u[i, 0] = 0.5*(1.-((math.exp(2*(0.2*ds*i - 5.))-1)/(math.exp(2*(0.2*ds*i-5.) )+1)))\n",
    "   u[0, 1] = 1.\n",
    "   u[0, 2] = 1.\n",
    "   u[130, 1] = 0. # End points\n",
    "   u[130, 2] = 0.\n",
    "    \n",
    "init()   \n",
    "k = range(0,mx)   \n",
    "fig = plt.figure()          # figure to plot (a changing line)\n",
    "ax = fig.add_subplot(111, autoscale_on=False, xlim=(0, mx), ylim=(0, 3))\n",
    "ax.grid()                                                               \n",
    "plt.ylabel(\"Height\")        \n",
    "plt.title(\"Soliton (runs slowly)\")\n",
    "line, = ax.plot(k, u[k,0],\"b\", lw=2)    \n",
    "\n",
    "for  i in range (1, mx - 1):                                       # First time step\n",
    "    a1 = eps*dt*(u[i + 1, 0] + u[i, 0] + u[i - 1, 0])/(ds*6.)     \n",
    "    if i>1 and  i < 129: a2 = u[i+2, 0] + 2.*u[i-1, 0] - 2.*u[i+1, 0] - u[i-2, 0]\n",
    "    else:  a2 = u[i - 1, 0] - u[i + 1, 0]\n",
    "    a3 = u[i + 1, 0] - u[i - 1, 0] \n",
    "    u[i, 1] = u[i, 0] - a1*a3 -3*fac*a2/8.\n",
    " \n",
    "def animate(num):  \n",
    "    for i in range(1, mx - 2):\n",
    "          a1 = eps*dt*(u[i + 1, 1] + u[i, 1] + u[i - 1, 1])/(3.*ds)\n",
    "          if i>1 and i < mx - 2: a2 = u[i+2, 1] + 2.*u[i-1, 1] - 2.*u[i+1, 1] - u[i-2, 1]\n",
    "          else:  a2 = u[i-1, 1] - u[i+1, 1]  \n",
    "          a3 = u[i+1, 1] - u[i-1, 1] \n",
    "          u[i, 2] = u[i, 0] - a1*a3 - 3.*fac*a2/4.\n",
    "    line.set_data(k,u[k,2])     # data to plot (x,y)=(position,height)                        \n",
    "    u[k, 0] = u[k, 1]    # Recycle array         \n",
    "    u[k, 1] = u[k, 2]\n",
    "    return line,\n",
    "\n",
    "ani = animation.FuncAnimation(fig, animate,1)                                                                # Finish plot\n",
    "plt.show()\n",
    "print(\"finished\") "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}