notebooks/verlet_2body.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Introduction to Verlet integration (exmaple: 2-body system)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />\n",
    "\n",
    "<p style=\"text-align: justify;font-size:15px\"> \n",
    "In the molecular dyanmics simulations, the velocity, position and acceleration are sloved numerically by using the Verlet algorithm.\n",
    "<p>\n",
    "\n",
    "<p style=\"text-align: justify;font-size:20px\"> \n",
    "$\\mathbf{r}(t+\\Delta t) = 2\\mathbf{r}(t)-\\mathbf{r}(t-\\Delta t) + \\mathbf{a}(t)\\Delta t^2 + O(\\Delta t^4)$\n",
    "<p>\n",
    "\n",
    "<p style=\"text-align: justify;font-size:20px\">\n",
    "$\\mathbf{a}(t) = - \\frac{1}{m}\\nabla\\mathbf{v}(\\mathbf{r}(t))$\n",
    "<p>\n",
    "    \n",
    "<p style=\"text-align: justify;font-size:23px\">\n",
    "$\\mathbf{v}(t+ \\Delta t) = \\frac{\\mathbf{r}(t+\\Delta t)-\\mathbf{r}(t-\\Delta t)}{2\\Delta t}$\n",
    "<p>\n",
    "    \n",
    "<p style=\"text-align: justify;font-size:15px\"> \n",
    "Since this is a numerical method, the choosing of the time step $\\Delta t$ is very important. \n",
    "    If one chooses a very small time step to get accurate results, the simulations will take \n",
    "    long time to run. On the other hand, if one chooses a very large time step, the molecular \n",
    "    simulations will go to wrong results.\n",
    "<p>\n",
    "    \n",
    "<p style=\"text-align: justify;font-size:15px\"> \n",
    "Here, we domonstrates a two-body system (sun and earth) with the Verlet algorithm. The sun is \n",
    "    fixed at the centre positon. By tunning the time step with the slider, one can see how the\n",
    "    choosing of the time step $\\Delta t$ influences the trajectory of the earth.\n",
    "<p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.animation import FuncAnimation\n",
    "from ipywidgets import Button, IntSlider\n",
    "from numpy.linalg import norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_force(r1, r2, m1, m2):\n",
    "    x = G*m1*m2/np.linalg.norm(r1-r2)**2\n",
    "    v = (r2-r1)/np.linalg.norm(r2-r1)\n",
    "    return x*v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_acceleration(f, m):\n",
    "    return f/m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verlet_update(r1, r2, v, m1, m2, dt, a):\n",
    "    \n",
    "    rn = 2.0*r1[-1] - r1[-2]  + dt*dt*a\n",
    "    vn = (rn - r1[-2])/(2*dt)\n",
    "    \n",
    "    r1[-2] = r1[-1]\n",
    "    r1[-1] = rn\n",
    "    \n",
    "    f = compute_force(rn, r2, m1, m2)\n",
    "    a = compute_acceleration(f, m1)\n",
    "    \n",
    "    return r1, vn, a\n",
    "    \n",
    "    #return np.append(r1, [rn], axis=0), vn\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "\n",
    "m1 = 1.0\n",
    "m2 = 1000.0\n",
    "G = 0.2\n",
    "v = [0.0, 3.0]\n",
    "dt = 0.01\n",
    "\n",
    "r1 = np.array([[10, 0],[10.0, 3.0*dt]])\n",
    "r2 = np.array([0, 0])\n",
    "\n",
    "f = compute_force(r1[-1], r2, m1, m2)\n",
    "a = compute_acceleration(f, m1)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "fig.canvas.header_visible = False\n",
    "\n",
    "sun = plt.scatter([0.0],[0.0], s=23**2, \n",
    "               marker=\"o\", linewidth=0, color='gold', label = 'Sun')\n",
    "earth, = plt.plot([10.0],[0.0], 'r.', alpha = 1.0, markersize=15, label = 'Earth')\n",
    "origin = plt.plot([10.0],[0.0], 'kx', alpha = 1.0, markersize=15, label = 'Origin point')\n",
    "avector = plt.quiver([10.0], [0], [0], [0], color = ['blue'], units='inches', scale = 10.0, label=\"a\")\n",
    "fvector = plt.quiver([10.0], [0], [0], [0], color = ['gold'], units='inches', scale = 10.0, label=\"v\")\n",
    "trace, = plt.plot([10.0], [0.0], \"k\", linewidth=0.5)\n",
    "title = ax.set_title(\"\")\n",
    "\n",
    "plt.legend(loc=2, fontsize=12)\n",
    "\n",
    "plt.xlim(-11.0, 11.0)\n",
    "plt.ylim(-11.0, 11.0)\n",
    "plt.grid()\n",
    "\n",
    "\n",
    "xt = [];\n",
    "yt = [];\n",
    "\n",
    "def init():\n",
    "    return (earth, avector, fvector, trace)\n",
    "\n",
    "    \n",
    "def update(i):\n",
    "    global r1, v, a, xt, yt\n",
    "    r1, v, a = verlet_update(r1, r2, v, m1, m2, dt, a)\n",
    "\n",
    "    earth.set_data([r1[-1, 0]], [r1[-1, 1]])\n",
    "    xt.append(r1[-1, 0])\n",
    "    yt.append(r1[-1, 1])\n",
    "    trace.set_data(xt, yt)\n",
    "    avector.set_offsets(np.array([r1[-1, 0], r1[-1, 1]]))\n",
    "    fvector.set_offsets(np.array([r1[-1, 0], r1[-1, 1]]))\n",
    "    avector.set_UVC(a[0], a[1])\n",
    "    fvector.set_UVC(v[0], v[1])\n",
    "    \n",
    "    title.set_text(\"v(t) = %-5.3f\" % norm(v))\n",
    "\n",
    "    return (earth, avector, fvector, trace)\n",
    "\n",
    "\n",
    "anim = FuncAnimation(fig, update, frames = 20000, interval = 10,\n",
    "                    init_func=init, blit=True, repeat = True)\n",
    "\n",
    "\n",
    "def onClick(event):\n",
    "    if button_pause.description == \"Pause\":\n",
    "        button_pause.description = \"Play\"\n",
    "        anim.event_source.stop()\n",
    "    else:\n",
    "        button_pause.description = \"Pause\"\n",
    "        anim.event_source.start()\n",
    "        \n",
    "\n",
    "def clear_data(b):\n",
    "    global xt, yt\n",
    "    \n",
    "    xt = [];\n",
    "    yt = [];\n",
    "\n",
    "button_pause = Button(description=\"Pause\");\n",
    "button_pause.on_click(onClick)\n",
    "button_clear = Button(description=\"Clear trace\")\n",
    "button_clear.on_click(clear_data)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "def on_dt_change(c):\n",
    "    global dt, r1, xt, yt\n",
    "    \n",
    "    anim.event_source.stop()\n",
    "    \n",
    "    dt = w_dt.value*0.01\n",
    "    r1 = np.array([[10, 0],[10.0, 3.0*dt]])\n",
    "    \n",
    "    xt = []\n",
    "    yt = []\n",
    "    \n",
    "    anim.event_source.start()\n",
    "\n",
    "w_dt = IntSlider(value = 1, min = 1, max = 20, step=1, description=\"dt: \")\n",
    "w_dt.observe(on_dt_change, names='value')\n",
    "\n",
    "display(w_dt, button_pause, button_clear)"
   ]
  },
  {
   "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}