04/009x_43.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Stochastic Processes: <br>Data Analysis and Computer Simulation\n",
    "<br>\n",
    "\n",
    "\n",
    "# Brownian motion 2: computer simulation\n",
    "<br>\n",
    "\n",
    "\n",
    "# 3. Simulations with on-the-fly animation\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.1. Simulation code with on-the-fly animation\n",
    "\n",
    "## Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np # import numpy library as np\n",
    "import matplotlib.pyplot as plt # import pyplot library as plt \n",
    "import matplotlib.mlab as mlab # import mlab module to use MATLAB commands with the same names\n",
    "import matplotlib.animation as animation # import animation modules from matplotlib\n",
    "from mpl_toolkits.mplot3d import Axes3D # import Axes3D from mpl_toolkits.mplot3d\n",
    "from IPython.display import HTML\n",
    "plt.style.use('ggplot') # use \"ggplot\" style for graphs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Define `init` function for `FuncAnimation`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def init():\n",
    "    global R,V,W,Rs,Vs,Ws,time\n",
    "    R[:,:] = 0.0 # initialize all the variables to zero\n",
    "    V[:,:] = 0.0 # initialize all the variables to zero\n",
    "    W[:,:] = 0.0 # initialize all the variables to zero\n",
    "    Rs[:,:,:] = 0.0 # initialize all the variables to zero\n",
    "    Vs[:,:,:] = 0.0 # initialize all the variables to zero\n",
    "    Ws[:,:,:] = 0.0 # initialize all the variables to zero\n",
    "    time[:]   = 0.0 # initialize all the variables to zero\n",
    "    title.set_text(r'') # empty title\n",
    "    line.set_data_3d([],[],[]) # set line data to show the trajectory of particle n in 3d (x,y,z)\n",
    "    particles.set_data_3d([],[],[]) # set position current (x,y,z) position data for all particles\n",
    "    return particles,title,line # return listed objects that will be drawn by FuncAnimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Define `animate` function for `FuncAnimation`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def animate(i):\n",
    "    global R,V,W,Rs,Vs,Ws,time # define global variables\n",
    "    time[i]=i*dt # store time in each step in an array time\n",
    "    W = std*np.random.randn(nump,dim) # generate an array of random forces accordingly to Eqs.(F10) and (F11)\n",
    "    R, V = R + V*dt, V*(1-zeta/m*dt)+W/m # update R & V via Eqs.(F5)&(F9)\n",
    "    Rs[i,:,:]=R # accumulate particle positions at each step in an array Rs\n",
    "    Vs[i,:,:]=V # accumulate particle velocitys at each step in an array Vs\n",
    "    Ws[i,:,:]=W # accumulate random forces at each step in an array Ws\n",
    "    title.set_text(r\"t = \"+str(time[i])) # set the title to display the current time\n",
    "    line.set_data(Rs[:i+1,n,0],Rs[:i+1,n,1]) # set the line in 2D (x,y)\n",
    "    line.set_3d_properties(Rs[:i+1,n,2]) # add z axis to set the line in 3D\n",
    "    particles.set_data(R[:,0],R[:,1]) # set the current position of all the particles in 2d (x,y)\n",
    "    particles.set_3d_properties(R[:,2]) # add z axis to set the particle in 3D\n",
    "    return particles,title,line # return listed objects that will be drawn by FuncAnimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Set parameters and initialize variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "dim  = 3 # system dimension (x,y,z)\n",
    "nump = 1000 # number of independent Brownian particles to simulate \n",
    "nums = 1024 # number of simulation steps\n",
    "dt   = 0.05 # set time increment, \\Delta t\n",
    "zeta = 1.0 # set friction constant, \\zeta\n",
    "m    = 1.0 # set particle mass, m\n",
    "kBT  = 1.0 # set temperatute, k_B T\n",
    "std  = np.sqrt(2*kBT*zeta*dt) # calculate std for \\Delta W via Eq.(F11)\n",
    "np.random.seed(0) # initialize random number generator with a seed=0\n",
    "R = np.zeros([nump,dim]) # array to store current positions and set initial condition Eq.(F12)\n",
    "V = np.zeros([nump,dim]) # array to store current velocities and set initial condition Eq.(F12)\n",
    "W = np.zeros([nump,dim]) # array to store current random forcces\n",
    "Rs = np.zeros([nums,nump,dim]) # array to store positions at all steps\n",
    "Vs = np.zeros([nums,nump,dim]) # array to store velocities at all steps\n",
    "Ws = np.zeros([nums,nump,dim]) # array to store random forces at all steps\n",
    "time = np.zeros([nums]) # an array to store time at all steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Perform and animate the simulation using `FuncAnimation` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,10)) # set fig with its size 10 x 10 inch\n",
    "ax = fig.add_subplot(111,projection='3d') # creates an additional axis to the standard 2D axes\n",
    "box = 40 # set draw area as box^3\n",
    "ax.set_xlim(-box/2,box/2) # set x-range\n",
    "ax.set_ylim(-box/2,box/2) # set y-range\n",
    "ax.set_zlim(-box/2,box/2) # set z-range\n",
    "ax.set_xlabel(r\"x\",fontsize=20) # set x-lavel\n",
    "ax.set_ylabel(r\"y\",fontsize=20) # set y-lavel\n",
    "ax.set_zlabel(r\"z\",fontsize=20) # set z-lavel\n",
    "ax.view_init(elev=12,azim=120) # set view point\n",
    "particles, = ax.plot([],[],[],'ro',ms=8,alpha=0.5) # define object particles\n",
    "title = ax.text(-180.,0.,250.,r'',transform=ax.transAxes,va='center') # define object title\n",
    "line, = ax.plot([],[],[],'b',lw=1,alpha=0.8) # define object line\n",
    "n = 0 # trajectry line is plotted for the n-th particle\n",
    "anim = animation.FuncAnimation(fig,func=animate,init_func=init,frames=nums,interval=5,blit=True,repeat=False)\n",
    "## If you have ffmpeg installed on your machine \n",
    "## you can save the animation by uncomment the last line\n",
    "## You may install ffmpeg by typing the following command in command prompt\n",
    "## conda install -c menpo ffmpeg\n",
    "## \n",
    "# anim.save('movie.mp4',fps=50,dpi=100)\n",
    "\n",
    "HTML(anim.to_html5_video())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary of simulation methods\n",
    "\n",
    "### Original differential equation\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{d\\mathbf{R}(t)}{dt}=\\mathbf{V}(t)\\tag{F1}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "m\\frac{d\\mathbf{V}(t)}{dt}=\\color{black}{-\\zeta\\mathbf{V}(t)}+\\color{black}{\\mathbf{F}(t)} \n",
    "\\tag{F2}\n",
    "\\end{equation}\n",
    "$\\hspace{80mm}$with\n",
    "\\begin{equation}\n",
    "\\langle \\mathbf{F}(t)\\rangle=\\mathbf{0} \n",
    "\\tag{F3}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\langle \\mathbf{F}(t)\\mathbf{F}(0)\\rangle = {2k_B T\\zeta}\\mathbf{I}\\delta(t)  \n",
    "\\tag{F4}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Euler method\n",
    "\n",
    "$$\n",
    "\\mathbf{V}_{i+1}\n",
    "=\\left(1-\\frac{\\zeta}{m}\\Delta t\\right)\\mathbf{V}_i + \\frac{1}{m} {\\Delta \\mathbf{W}_i} \n",
    "\\tag{F9}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{R}_{i+1}=\\mathbf{R}_i+\\mathbf{V}_i \\Delta t \\hspace{15mm}\\tag{B3}\n",
    "$$\n",
    "$\\hspace{80mm}$with\n",
    "\\begin{equation}\n",
    "\\langle \\Delta \\mathbf{W}_i\\rangle=\\mathbf{0} \n",
    "\\tag{F10}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\langle \\Delta \\mathbf{W}_i\\Delta \\mathbf{W}_j\\rangle = {2k_B T\\zeta}\\Delta t\\mathbf{I}\\delta_{ij}  \n",
    "\\tag{F11}\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 2nd order Runge-Kutta method\n",
    "\n",
    "$$\n",
    "\\mathbf{V}'_{i+\\frac{1}{2}}\n",
    "=\\mathbf{V}_i-\\frac{\\zeta}{m}\\frac{\\Delta t}{2}\\mathbf{V}_{i}  \n",
    "=\\left(1-\\frac{\\zeta}{m}\\frac{\\Delta t}{2}\\right)\\mathbf{V}_{i}  \n",
    "\\tag{F12}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{V}_{i+1}\n",
    "=\\mathbf{V}_i-\\frac{\\zeta}{m}\\Delta t\\mathbf{V}'_{i+\\frac{1}{2}} + \\frac{1}{m} {\\Delta \\mathbf{W}_i} \n",
    "\\tag{F13}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{R}_{i+1}=\\mathbf{R}_i+\\mathbf{V}'_{i+\\frac{1}{2}} \\Delta t \\hspace{15mm}\n",
    "\\tag{F14}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 4th order Runge-Kutta method\n",
    "\n",
    "$$\n",
    "\\mathbf{V}'_{i+\\frac{1}{2}}\n",
    "=\\mathbf{V}_i-\\frac{\\zeta}{m}\\frac{\\Delta t}{2}\\mathbf{V}_{i}  \n",
    "\\tag{F15}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{V}''_{i+\\frac{1}{2}}\n",
    "=\\mathbf{V}_i-\\frac{\\zeta}{m}\\frac{\\Delta t}{2}\\mathbf{V}'_{i+\\frac{1}{2}}  \n",
    "\\tag{F16}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{V}'''_{i+1}\n",
    "=\\mathbf{V}_i-\\frac{\\zeta}{m}{\\Delta t}\\mathbf{V}''_{i+\\frac{1}{2}}  \n",
    "\\tag{F17}\n",
    "$$\n",
    "\n",
    "\n",
    "$$\n",
    "\\mathbf{V}_{i+1}\n",
    "=\\mathbf{V}_i-\\frac{\\zeta}{m}\\frac{\\Delta t}{6}\\left(\\mathbf{V}+\\mathbf{V}'_{i+\\frac{1}{2}}+\\mathbf{V}''_{i+\\frac{1}{2}}+\\mathbf{V}'''_{i+1}\\right) + \\frac{1}{m} {\\Delta \\mathbf{W}_i} \n",
    "\\tag{F18}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{R}_{i+1}=\\mathbf{R}_i+\n",
    "\\frac{\\Delta t}{6}\\left(\\mathbf{V}+\\mathbf{V}'_{i+\\frac{1}{2}}+\\mathbf{V}''_{i+\\frac{1}{2}}+\\mathbf{V}'''_{i+1}\\right)\n",
    "\\hspace{15mm}\n",
    "\\tag{F19}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}