paqueteviajero-diagonalizacion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discretizando la ecuación de Schrödinger 2. Operador evolucion.\n",
    "**Autor: Juan Mauricio Matera (Facultad de Ciencias Exactas, UNLP)** \n",
    "\n",
    "23 de agosto de 2020\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En esta segunda versión, en lugar de resolver la ecuación de Schrödinger discretizada en x y t, sólo discretizaremos en x. Para esto, notamos que la ecuación de Schrödinger para Hamiltonianos independientes de t puede resolverse por separación de variables:\n",
    "$$\n",
    "\\psi(\\vec{x},t)=\\sum_{E} a(E) e^{-{\\bf i}E t/\\hbar}\\varphi_E(\\vec{x})\n",
    "$$\n",
    "con \n",
    "$$\n",
    "E \\varphi_i(\\vec{x})=-\\frac{\\hbar^2}{2 m}\\nabla^2 \\varphi_i(\\vec{x}) +V(x)\\varphi_i(\\vec{x})\n",
    "$$\n",
    "que es una ecuación de autovalores para el operador ${\\bf H}=-\\frac{\\hbar^2}{2 m}\\nabla^2+V(x)$. \n",
    "Como ${\\bf H}$ es un operador autoadjunto, sus autofunciones $\\varphi_E(\\vec{x})$ expanden una base ortogonal,\n",
    "lo que permite fijar los coeficientes $a(E)$ a partir de la condición inicial usando la ortogonalidad:\n",
    "$$\n",
    "a(E)=\\frac{\\int \\varphi_E(\\vec{x})^*\\psi(\\vec{x},0) dx}{\\int \\varphi_E(\\vec{x})^*\\varphi_E(\\vec{x}) dx}\n",
    "$$\n",
    "\n",
    "Podemos también expresar en forma compacta la solución, introduciendo el operador evolución\n",
    "$$\n",
    "{\\bf U}(t)=\\exp(-{\\bf i}{\\bf H}t/\\hbar)\n",
    "$$\n",
    "de manera que\n",
    "$$\n",
    "\\psi(\\vec{x},t+ \\Delta t)={\\bf U}(\\Delta t)\\psi(\\vec{x},t)=\n",
    "\\sum_{E} a(E) e^{-{\\bf i}E \\Delta t/\\hbar}\\varphi_E(\\vec{x})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Discretización\n",
    "============\n",
    "El problema es entonces construir ${\\bf U}(t)$. Para esto, vamos a construir una versión discretizada del operador ${\\bf H}$, sobre un conjunto de puntos equiespaciados en el intervalo $(-L,L)$. Para esto, introducimos una base finita de funciones en ese intervalo, definidas por\n",
    "$$\n",
    "\\varphi_n(x)=\\exp({\\bf i} \\frac{\\pi n}{L} x)/\\sqrt{2L}\n",
    "$$\n",
    "con $2N+1$ el número de puntos en que dividimos el intervalo, y $-N\\leq n\\leq N$. Estas funciones son autofunciones de $\\nabla^2$  por lo que el término cinético es diagonal en esta base.\n",
    "\n",
    "Por otro lado, podemos combinar estas funciones para obtener una base complementaria\n",
    "$$\n",
    "\\xi_{m}(x)=\\sum_{n=-N}^{N}  \\frac{e^{-{\\bf i}\\frac{\\pi n m}{N+1/2}}}{\\sqrt{(2N+1)}}\\varphi_n(x)=\n",
    "\\frac{1}{\\sqrt{2L(2N+1)}}\\frac{\\sin\\left(\\pi(N+1/2)(x/L-m/(N+1/2))\\right)}\n",
    "{\\sin\\left(\\pi/2(x/L-m/(N+1/2))\\right)}\n",
    "$$\n",
    "que cumplen las condiciones\n",
    "\\begin{eqnarray}\n",
    "\\int_{-L}^L \\xi_{m}^*(x)\\xi_{m'}(x) dx&=&\\frac{L}{N+1/2}\\sum_{n=-N}^{N}\\xi_{m}^*(x_n)\\xi_{m'}(x_n)=\\delta_{m,m'}\\\\\n",
    "\\int_{-L}^{L} \\xi_{m}^*(x)\\phi(x)dx&\\approx&\\frac{(\\phi(x_n) -  \\left(\\frac{L}{\\pi(N+1/2)}\\right)^2\\phi''(x_n)+\\ldots)}{\\sqrt{(N+1/2)/L}}\\\\\n",
    "\\frac{\\xi_{m}(x_n)}{\\sqrt{(N+1/2)/L}}&\\approx& \\delta_{m,n}\n",
    "\\end{eqnarray}\n",
    "con $x_n=\\frac{n L}{N+1/2}$. De estas propiedades, se sigue que \n",
    "$$\n",
    "\\psi(x)\\approx \\sum_{m=-N}^{N}  \\frac{\\xi_m(x)}{\\sqrt{(N+1/2)/L}}\\psi(x_m)\n",
    "$$\n",
    "y\n",
    "$$\n",
    "V(x)\\psi(x)\\approx \\sum_{m=-N}^{N}  \\frac{\\xi_m(x)}{\\sqrt{(N+1/2)/L}}V(x_m)\\psi(x_m)\n",
    "$$\n",
    "de lo que el potencial es (para $N$ suficientemente grande) diagonal en la base $\\xi_m(x)$. \n",
    "\n",
    "\n",
    "De la misma manera, podemos determinar el efecto del operador energía cinética ($T=-\\frac{\\hbar^2}{2M}\\nabla^2$) en esta base, a partir de su acción sobre la base $\\varphi_n(x)$:\n",
    "$$\n",
    "T \\xi_m(x)=\\frac{\\hbar^2}{2M}\\sum_{n=-N}^{N} \\frac{\\pi^2 n^2}{L^2} \\frac{e^{-{\\bf i}\\frac{\\pi n m}{N+1/2}}}{\\sqrt{(2N+1)}}\\varphi_n(x)=\n",
    "\\frac{\\hbar^2}{2M}\\sum_{n=-N}^{N}\\sum_{m'=-N}^{N} \\frac{\\pi^2 n^2}{L^2} \\frac{e^{{\\bf i}\\frac{\\pi n (m-m')}{N+1/2}}}{{(2N+1)}}\\xi_{m'}(x)=\\sum_{m'}T_{mm'}\\xi_{m'}(x)\n",
    "$$\n",
    "con\n",
    "$$\n",
    "T_{mm'}=\\frac{\\hbar^2}{2M}\\sum_{n=-N}^{N}\\frac{\\pi^2 n^2}{L^2} \\frac{e^{{\\bf i}\\frac{\\pi n (m-m')}{N+1/2}}}{{(2N+1)}}\n",
    "$$\n",
    "\n",
    "Restringidos a esta base,  cualquier solución de la ecuación de Schrödinger para una condición inicial suficientemente suave se expresa como\n",
    "\n",
    "$$\n",
    "\\psi(x,t)\\approx \\sum_{m,m'}    U_{m',m}(t-t_0)  \\psi(x_{m},t_0)\\frac{\\xi_{m}(x)}{\\sqrt{(N+1/2)/L}}\n",
    "$$\n",
    "\n",
    "con \n",
    "$$\n",
    "U_{m,m'}=[\\exp(-{\\bf i}H(t-t_0)/\\hbar]_{m,m'}\n",
    "$$\n",
    "el operador de evolución en el subespacio de $\\xi_{m}(x)$, con\n",
    "$$\n",
    "H_{m,m'}=\\sum_{k=-N}^{N} \\frac{\\hbar^2 \\pi^2 n^2}{2 M L^2}\\frac{e^{{\\bf i}(m-m')\\frac{\\pi n}{N+1/2}}}{2N+1}+U(m)\\delta_{m,m'}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "De esta manera, el problema se reduce a \n",
    "1.  Construir la matriz H. \n",
    "2. Diagonalizar la matriz $H$ como $H={\\cal V}{\\Omega}{\\cal V}^\\dagger$ con $\\Omega$ una matriz diagonal y \n",
    "${\\cal V}^\\dagger$ construir el operador de evolución con un paso constante $\\Delta t$: \n",
    "$$\n",
    "U(\\Delta t)={\\cal V}\\exp(-{\\bf i}{\\Omega}\\Delta t/\\hbar){\\cal V}^\\dagger\n",
    "$$\n",
    "3. Construir el estado inicial como un vector $\\Psi_{t_0}$ con los valores de $\\psi(x,t_0)$ evaluados en los puntos $x_n$\n",
    "4. Construir el estado en el instante $t+n\\Delta t$ mediante la recurrencia\n",
    "$$\n",
    "\\Psi_{t+\\Delta t}=U(\\Delta t)\\Psi_{t}\n",
    "$$\n",
    "\n",
    "La ventaja de este método respecto a la integración discreta en tiempos es que preserva la norma de la función de onda. Es importante sin embargo que el estado inicial esté localizado lejos de los bordes, y que durante la evolución, nunca alcance un valor apreciable en estos. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementación\n",
    "\n",
    "En esta notebook, implementamos el algoritmo anterior de manera más o menos directa. El único detalle es que explotamos la función fft de la librería numpy para evaluar eficientemente los elementos de matriz de $T$. Esto sin embargo implica tener que convertir el índice $n$ correspondiente a los estados de $k$ del intervalo (0,2N+1) al intervalo (-N,N).\n",
    "Otro detalle de la implementación es que en lugar de fijar N, se fijó el espaciado $dx=L/N$. Esto tiene consecuencias menores en los desarrollos previos, que involucran correcciones ${\\cal O}(1/N)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulación"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-08-23T04:14:08.427491Z",
     "start_time": "2020-08-23T04:14:08.276944Z"
    },
    "hide_input": true,
    "init_cell": true,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d76d820473004d548e3bc92749b03f34",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(VBox(children=(Label(value='Paquete'), FloatSlider(value=-150.0, description='centro del paquet…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "from IPython.display import HTML\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "\n",
    "L = 200\n",
    "dx = .5\n",
    "U = lambda x:(0.*x)\n",
    "\n",
    "buttonrunsim = widgets.Button(description=\"Simulate\")\n",
    "buttongenerateu = widgets.Button(description=\"Generate U\")\n",
    "outwdt = widgets.Output()\n",
    "\n",
    "centropaquete = widgets.FloatSlider(value=-3*L/4,min=-L,max=L,step=1,\n",
    "    description='centro del paquete',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "anchopaquete = widgets.FloatSlider(value=L/8,min=1,max=L,step=1,\n",
    "    description='ancho del paquete',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "\n",
    "valuek0 = widgets.FloatSlider(value=0,min=-np.pi/dx,max=np.pi/dx,step=.1,\n",
    "    description='k0',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "valueV0 = widgets.FloatSlider(value=0,min=-4,max=4,step=.1,\n",
    "    description='V0',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "potencial = widgets.Dropdown(options=[\"Lineal\",\"Cuadrado\",\"Rampa\",\"Cuadrático\",\"Gaussiano\"],\n",
    "                            description=\"Forma del Potencial\",\n",
    "                             ensure_option=True,\n",
    "                            disabled=False)\n",
    "\n",
    "valueW = widgets.FloatSlider(value=1.,min=1.,max=300.,step=.5,\n",
    "    description='Alcance potencial',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "\n",
    "items = [[widgets.Label(\"Paquete\"),centropaquete,anchopaquete,valuek0,buttonrunsim,\n",
    "          widgets.Label(\"Potencial\"),valueW,valueV0,potencial,buttongenerateu],[outwdt]]\n",
    "\n",
    "widgets.HBox([widgets.VBox(it) for it in items])\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false
   },
   "source": [
    "# Código que controla las simulaciones\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-08-23T04:14:08.436163Z",
     "start_time": "2020-08-23T04:14:08.430651Z"
    },
    "hide_input": false,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "# Para fijar parámetros a mano\n",
    "\n",
    "\n",
    "#### Parámetros de la discretización\n",
    "dx = .5\n",
    "dt = 10.\n",
    "L = 200\n",
    "\n",
    "####  Estado inicial\n",
    "a = 20.\n",
    "x0 = -L/3\n",
    "k0 = 10./a\n",
    "\n",
    "#### Potencial\n",
    "U0= -1.\n",
    "w = 10.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-08-12T16:07:09.782844Z",
     "start_time": "2020-08-12T16:07:09.778530Z"
    }
   },
   "source": [
    "## Rutinas que generan inicializan y generan la simulación"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-08-23T04:16:05.371884Z",
     "start_time": "2020-08-23T04:16:04.889087Z"
    },
    "hide_input": false,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "import numpy.fft as fft\n",
    "import numpy.linalg as la\n",
    "############    Inicializar  #############\n",
    "\n",
    "\n",
    "def init_evolop(L=L,dx=dx,dt=dt,pot=None):\n",
    "    global rho0\n",
    "    global psi\n",
    "    global psi0\n",
    "    global xs\n",
    "    global evol_op\n",
    "    global P, ham\n",
    "    xs = np.linspace(-L,L,int(2*L/dx))\n",
    "    ## Construyo el operador evolución\n",
    "    nhalf = int(L/dx)\n",
    "    # Energía cinética\n",
    "    ks = np.pi/L*np.array([(1.*q)  if q<nhalf else  (q-2*nhalf) for q in range(2*nhalf)])\n",
    "    Uft = np.array([fft.ifft(v)*np.sqrt(2*nhalf) for v in  np.eye(2*nhalf)])\n",
    "    P = (ks*(Uft).transpose()).dot(Uft.conj())    \n",
    "    ham = .5*P.dot(P)\n",
    "    # Energía potencial\n",
    "    if pot is None:\n",
    "        U = 0.*xs\n",
    "    else:\n",
    "        U = pot(xs)\n",
    "    ham = ham + np.diag(U)\n",
    "    ens,states = np.linalg.eigh(ham)\n",
    "    evol_op = (states.conj()).dot((np.exp(-1j*dt*ens)*states).transpose())\n",
    "    return evol_op\n",
    "\n",
    "init_evolop()\n",
    "\n",
    "def init_state(L=L,dx=dx,dt=dt,a=a,x0=x0,k0=k0):\n",
    "    global rho0\n",
    "    global psi\n",
    "    global psi0\n",
    "    global xs\n",
    "    global evol_op\n",
    "    global P\n",
    "    xs = np.linspace(-L,L,int(2*L/dx))\n",
    "    rho0 = np.array([np.exp(-(x-x0)**2/(2*a**2))  for x in xs]) \n",
    "    rho0[0] = 0\n",
    "    rho0[-1] = .5*rho0[-2]\n",
    "    rho0[1] = .5*rho0[2]\n",
    "    rho0 = rho0/sum(rho0)/dx\n",
    "    psi0 = rho0**.5 * np.exp(1j*k0*(xs-x0))\n",
    "    psi = psi0\n",
    "    return psi\n",
    "\n",
    "\n",
    "def Ulineal(xs):\n",
    "    return -U0*xs/w\n",
    "\n",
    "def Ugaussiano(xs):\n",
    "    return U0*np.exp(-xs**2/(.5*w**2))\n",
    "\n",
    "\n",
    "def Uescalon(xs):\n",
    "    return np.array([U0 if abs(x)<.5*w else 0 for x in xs])\n",
    "\n",
    "\n",
    "def Uarmonico(xs):\n",
    "    return np.array([U0*(1-4*(x/w)**2) if abs(x)<.5 * w else 0 for x in xs])\n",
    "\n",
    "\n",
    "def Urampadoble(xs):\n",
    "    return np.array([U0*(1-2*abs(x)/w) if abs(x)<.5 * w else 0 for x in xs])\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def make_animation_new(ts=50, progress=None):\n",
    "    global psi\n",
    "    global xs\n",
    "    global U0,k0,w,x0,a,U, evol_op\n",
    "    \n",
    "    psi = init_state(a=a,dt=dt,k0=k0,x0=x0)\n",
    "    rho0 = abs(psi)**2\n",
    "    fig1 = plt.figure()\n",
    "    plt.xlim(-L,L)\n",
    "    plt.ylim(-abs(U0)-k0**2-.5,.5+abs(U0)+k0**2)\n",
    "    plt.plot(xs, U(xs),ls=\"-.\")\n",
    "    density,= plt.plot(xs,50*(U0+.1)*rho0+.5*k0**2)\n",
    "\n",
    "    def update_graphic(t,density):\n",
    "        global psi, evol_op\n",
    "        if progress is not None:\n",
    "            progress.value = t\n",
    "\n",
    "        if t>=1:\n",
    "            psi = evol_op.dot(psi)\n",
    "        rho = abs(psi)**2\n",
    "        density.set_data(xs,50*(U0+.1)*rho+.5*k0**2)\n",
    "        \n",
    "        \n",
    "        return density,\n",
    "\n",
    "    line_ani = animation.FuncAnimation(fig1, update_graphic, ts, fargs=(density,),\n",
    "                                       interval=100, blit=True)\n",
    "\n",
    "    plt.close()\n",
    "    return HTML(line_ani.to_jshtml())\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rutinas que conectan los controles de la primera celda con la simulacion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-08-23T04:14:09.104226Z",
     "start_time": "2020-08-23T04:14:09.090501Z"
    },
    "hide_input": false,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "def on_buttongenerateu(b):\n",
    "    global U, evol_op, w, U0,ham\n",
    "    U0 = valueV0.value\n",
    "    w = valueW.value\n",
    "    if potencial.value == \"Lineal\":\n",
    "        U = Ulineal\n",
    "    elif potencial.value == \"Gaussiano\":\n",
    "        U = Ugaussiano\n",
    "    elif potencial.value == \"Cuadrático\":\n",
    "        U = Uarmonico\n",
    "    elif potencial.value == \"Cuadrado\":\n",
    "        U = Uescalon\n",
    "    elif potencial.value == \"Rampa\":\n",
    "        U = Urampadoble\n",
    "    else:\n",
    "        print(\"No encontrado\")\n",
    "        return\n",
    "    outwdt.clear_output()\n",
    "    with outwdt:\n",
    "        print(\"U=\",U)\n",
    "        print(\"U0=\",U0)\n",
    "        print(\"w=\",w)\n",
    "        evol_op = init_evolop(pot=U)\n",
    "    outwdt.clear_output()\n",
    "    with outwdt:\n",
    "        plt.plot(xs, U(xs))\n",
    "        display(plt.show())\n",
    "        plt.ylim(-5,5)\n",
    "        plt.xlim(0,100)\n",
    "        plt.scatter(range(len(ham)),sorted(la.eigvalsh(ham)))\n",
    "        plt.title(\"Hamiltonian Spectrum\")\n",
    "        display(plt.show())\n",
    "        \n",
    "\n",
    "def on_button_start_sim(b):\n",
    "    global U0, w, k0, a, x0, U, U0, w\n",
    "    k0 = valuek0.value\n",
    "    a = anchopaquete.value\n",
    "    x0 = centropaquete.value\n",
    "    \n",
    "    outwdt.clear_output()\n",
    "    progress = widgets.FloatProgress(value=0,min=0,max=50,step=1,description='Simulating:',\n",
    "    bar_style='info',orientation='horizontal')\n",
    "    \n",
    "    with outwdt:\n",
    "        print([a,x0,k0,U0,w])\n",
    "        display(progress)\n",
    "        result = make_animation_new(progress=progress)     \n",
    "        outwdt.clear_output()\n",
    "        display(result)\n",
    "        \n",
    "\n",
    "buttonrunsim.on_click(on_button_start_sim)\n",
    "buttongenerateu.on_click(on_buttongenerateu)\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Initialization Cell",
  "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.6"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": false,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}