Ondas y relación de dispersión.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Resolución de la ecuación de onda en medios dispersivos. Condiciones iniciales."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "La función\n",
    "\\begin{equation}\\label{soluciongeneral}\n",
    "\\Psi(x,t)=\\Re\\int F(k) e^{{\\bf i}(k x-\\omega(k)t)}dk\n",
    "\\end{equation}\n",
    "es una solución de la ecuación de onda modificada, porque es una combinación lineal de soluciones. Su derivada es\n",
    "$$\n",
    "\\dot{\\Psi}(x,t)=\\Re\\frac{\\partial}{\\partial t}{\\Psi}(x,t)=\n",
    "\\Re\\int  -{\\bf i}\\omega(k) F(k) e^{{\\bf i}(k x-\\omega(k)t)}dk\n",
    "$$\n",
    "\n",
    "Para fijar las condiciones iniciales, evaluamos estas integrales en $t=0$, de manera que obtenemos las siguientes condiciones:\n",
    "$$\n",
    "\\Psi(x,0)=\\Re\\int F(k) e^{{\\bf i}k x}dk=\\int \\frac{F^*(k)+F(k)}{2} e^{{\\bf i}k x}dk\n",
    "$$\n",
    "$$\n",
    "\\dot{\\Psi}(x,0)=\\Re\\int  -{\\bf i}\\omega(k) F(k) e^{{\\bf i}(k x)}dk\n",
    "=\\int  \\omega(k){\\bf i} \\frac{F^*(k)-F(k)}{2} e^{{\\bf i}k x}dk\n",
    "$$\n",
    "\n",
    "Para recuperar $F(k)$, transformamos Fourier a ambos lados:\n",
    "\n",
    "$$\n",
    "\\int \\Psi(x,0)e^{-{\\bf i}k x}\\frac{dx}{2\\pi}=\\frac{F^*(k)+F(k)}{2} =\\Re B(k)\n",
    "$$\n",
    "$$\n",
    "\\int \\dot{\\Psi}(x,0)e^{-{\\bf i}k x}\\frac{dx}{2\\pi}=\\omega(k){\\bf i} \\frac{F^*(k)-F(k)}{2}\n",
    "=\\omega(k)\\frac{F(k)-F^*(k)}{2{\\bf i }}=\\omega(k)\\Im F(k)\n",
    "$$\n",
    "Luego,\n",
    "$$\n",
    "F(k)=\\int \\Psi(x,0)e^{-{\\bf i}k x}\\frac{dx}{2\\pi} + \\frac{{\\bf i}}{\\omega(k)}\\int \\dot{\\Psi}(x,0)e^{-{\\bf i}k x}\\frac{dx}{2\\pi}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Podemos escribir directamente \\ref{soluciongeneral} en términos de las condiciones iniciales como\n",
    "\n",
    "$$\n",
    "\\Psi(x,t)=\\Re \\int\\int  (\\Psi(x',0)+\\frac{{\\bf i}}{\\omega(k)} \\dot{\\Psi}(x',0))   e^{{\\bf i}k (x-x')-\\omega(k)t}     \\frac{dk}{2\\pi} dx'  \n",
    "$$\n",
    "ó, definiendo\n",
    "$$\n",
    "G(x,t)=\\Re \\int \\frac{{\\bf i} e^{{\\bf i}(k x-\\omega(k)t)}}{\\omega(k)}\\frac{dk}{2\\pi}=\n",
    "\\int \\frac{\\sin(\\omega(k)t-k x)}{\\omega(k)}\\frac{dk}{2\\pi}\n",
    "$$\n",
    "de forma que\n",
    "$$\n",
    "\\frac{\\partial}{\\partial t}G(x,t)=\\Re \\int e^{{\\bf i}(k x-\\omega(k)t)}\\frac{dk}{2\\pi}=\n",
    "\\int \\cos(\\omega(k)t-k x)\\frac{dk}{2\\pi}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\Psi(x,t)=\\int \\left(\\frac{\\partial G(x-x',t)}{\\partial t} \\Psi(x',0)- G(x-x',t) \\dot{\\Psi}(x',0))\\right) dx' \n",
    "$$\n",
    "\n",
    "Nótese la analogía *formal* respecto a la identidad de Green que nos permitía relacionar la solución\n",
    "del problema de Laplace para ciertas condiciones de borde\n",
    "\n",
    "$$\n",
    "\\Psi(\\vec{r})=\\int_{\\partial {\\cal R}} \\left( (\\check{n}\\cdot \\nabla G(\\vec{r},\\vec{r}') \\Psi(x',0)- G(\\vec{r},\\vec{r}') \\check{n}\\cdot \\nabla\\Psi(\\vec{r'})\\right) d\\vec{r}' \n",
    "$$\n",
    "donde la integral se definía sobre el borde de la región ${\\cal R}$ a resolver. En el caso de la solución para la ecuación de onda, la *derivada normal* se remplaza por la derivada temporal, y el \"borde\" es todo el espacio (a t constante). Por otro lado, hay dos diferencias importantes entre ambas relaciones\n",
    "* En el caso de la ecuación de onda, estamos definiendo una *condición inicial*, mientras que en el problema de Laplace estamos fijando una *condición de borde*.\n",
    "* Para el caso de la ecuación de onda, la condición inicial hay que fijarla tanto para la función como su derivada, mientras que en el problema de Laplace, no se pueden fijar condiciones en la función y su derivada normal simultaneamente.\n",
    "* La derivación fue distinta: una salió como consecuencia del teorema de Stokes, mientras que la otra la construimos \"a mano\" despejando los coeficientes de la solución general.\n",
    "\n",
    "Más adelante, veremos que efectivamente estas dos expresiones están íntimamente relacionadas, si tomamos al tiempo como una \"coordenada espacial\" más. \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Paquetes de onda"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Un paquete de ondas es una solución de la ecuación de onda que en cada instante se encuentra localizada en una región del espacio. Para describirla, conviene trabajar con la \"versión compleja\" \n",
    "$$\n",
    "\\psi(x,t)=\\int F(k) e^{{\\bf i}(k x-\\omega(k)t)}dt\n",
    "$$\n",
    "de manera que si definimos \n",
    "$$\n",
    "\\rho(x,t)=\\frac{|\\psi(x,t)|^2}{2}\n",
    "$$\n",
    "una \"densidad de energía\", la cantidad\n",
    "$$\n",
    "U=\\int \\rho(x,t)dx =\\int \\frac{2\\pi |F(k)|^2}{2} dk \n",
    "$$\n",
    "es independiente de t. Notemos que\n",
    "$$\n",
    "\\rho(x,t)=\\frac{\\Re(\\psi(x,t))^2+\\Im(\\psi(x,t))^2}{2}\n",
    "$$\n",
    "de manera que \n",
    "$$\n",
    "\\rho_p(x,t)=\\frac{\\Re(\\psi(x,t))^2}{2}\n",
    "$$\n",
    "puede interpretarse como una *densidad de energía potencial*, y\n",
    "$$\n",
    "\\rho_k(x,t)=\\frac{\\Im(\\psi(x,t))^2}{2}\n",
    "$$\n",
    "una *densidad de energía cinética*. \n",
    "\n",
    "Para describir el movimiento del paquete de ondas, podemos usar a la densidad de energía de la misma forma que definíamos la densidad de masa en un sistema de partículas para definir un \"centro de masa\"\n",
    "$$\n",
    "x_{CM}(t)=\\frac{\\int \\rho(x,t) x dx}{\\int \\rho(x,t)dx}\n",
    "$$\n",
    "Como en cada instante $\\rho_p(x,t)$ y $\\rho_k(x,t)$ son aproximadamente  proporcionales a $\\rho(x,t)$, podemos\n",
    "usar cualquiera de estas para estimar la posición del \"centro de masa\" del paquete, y por lo tanto, la ubicación de la energía, cantidad de movimiento e información que este transporta.\n",
    "\n",
    "En un medio no dispersivo, la velocidad de fase $v_f=\\omega/k$ es constante, y \n",
    "$$\n",
    "\\psi(x,t)=\\psi(x-v_f t,0)\n",
    "$$\n",
    "de manera que el paquete se mueve sin distorsión, y \n",
    "$$\n",
    "x_{CM}(t)=x_{CM}(0) + v_f t\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Por otro lado, en un medio dispersivo, la velocidad de $x_{cm}$ dependerá de la composición en frecuencias\n",
    "que tenga el paquete. Por ser un paquete, $|F(k)|$ será no nula en un intervalo de números de onda\n",
    "de, al menos, la inversa del ancho espacial del paquete de ondas. Por lo tanto, las fases  varían de manera diferente para las diferentes componentes armónicas que componen al paquete, y este sufre \"dispersión\": el paquete se deforma.\n",
    "\n",
    "Si en el intervalo de números de onda sobre el que $|F(k)|$ tiene soporte, la derivada de $\\omega(k)$ respecto a $k$ es aproximadamente constante,\n",
    "\\begin{equation}\\label{condvg}\n",
    "\\omega(k)\\approx \\omega(k_0) + (k-k_0)v_g  \n",
    "\\end{equation}\n",
    "y\n",
    "$$\n",
    "\\psi(x,t)\\approx \\int F(k)e^{{\\bf i}(k x - ({\\bf i \\omega(k_0)} + (k-k_0)v_g)t)  }dk=\n",
    "\\int F(k)e^{{\\bf i}k (x-v_g t) -  {\\bf i}(\\omega(k_0) -k_0 v_g) t}   dk\n",
    "\\approx e^{{\\bf i}(\\omega(k_0)-k_0 v_g)t} \n",
    "\\psi(x-v_g,0)\n",
    "$$\n",
    "Llamamos entonces a\n",
    "$$\n",
    "v_g =\\frac{\\partial \\omega}{\\partial k}\n",
    "$$\n",
    "la *velocidad de grupo* del paquete de ondas. Cuando la condición \\ref{condvg} se cumple, entonces\n",
    "$$\n",
    "x_{cm}(t)\\approx x_{cm}(0)+ v_g t \n",
    "$$\n",
    "\n",
    "Por otro lado, si el paquete de ondas es muy angosto, o la relación de dispersión presenta un comportamiento singular en el entorno de $k_0$, \\ref{condvg} falla y $v_g$ deja de tener que ver con la velocidad del paquete. En particular, si $v_g$ sufre un cambio de signo en el entorno de $k_0$, lo que obtenemos es algo semejante a una *onda estacionaria*: las componentes armónicas se propagan en sentidos contrarios dando como resultado que el paquete deja de propagarse."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T15:46:28.202134Z",
     "start_time": "2020-05-08T15:46:28.197852Z"
    }
   },
   "source": [
    "# Demostraciones numéricas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T17:28:41.849936Z",
     "start_time": "2020-05-08T17:28:41.825585Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import numpy.fft as fft\n",
    "import matplotlib.animation as animation\n",
    "from matplotlib.animation import FuncAnimation\n",
    "from IPython.display import HTML\n",
    "import ipywidgets as widgets\n",
    "\n",
    "\n",
    "nmax = 200\n",
    "xmax = 60\n",
    "tmax = 10\n",
    "\n",
    "xs = None\n",
    "ks = None\n",
    "aes = None\n",
    "\n",
    "def paquete_inicial_gaussiano(x0, k0, a):\n",
    "    \"\"\"\n",
    "    Dada una función que corresponde a un paquete gaussiano de ancho a,\n",
    "    centrado en x0, y con una fase k0(x-x0)\n",
    "    \"\"\"\n",
    "    return (lambda x: np.exp(-(x-x0)**2/(2*a**2)+1j*k0*(x-x0) ))\n",
    "\n",
    "def paquete_inicial_cuadrado(x0, k0, a):\n",
    "    def pq(x):\n",
    "        if type(x) is np.ndarray:\n",
    "            return np.array([pq(xx) for xx in x])\n",
    "        else:\n",
    "            return np.exp(1j*k0*(x-x0))  if abs(x-x0)< a else 0 \n",
    "    return pq\n",
    "\n",
    "\n",
    "def paquete_inicial_triangular(x0, k0, a):\n",
    "    def pq(x):\n",
    "        if type(x) is np.ndarray:\n",
    "            return np.array([pq(xx) for xx in x])\n",
    "        else:\n",
    "            return np.exp(1j*k0*(x-x0))*(1-abs(x-x0)/a)  if abs(x-x0)< a else 0 \n",
    "    return pq\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def evol_paquete(paquete0, reldisp):\n",
    "    \"\"\"\n",
    "    Dada una función que define al estado inicial del paquete de ondas, y una \n",
    "    relación de dispersión, devuelve una función de x y t que describe la \n",
    "    evolución del paquete dado.\n",
    "    \"\"\"\n",
    "    global nmax, xmax, aes, xs, ks\n",
    "    # nmax = 1000\n",
    "    # xmax = 60\n",
    "    deltax = xmax / (2*nmax+1)\n",
    "    ks = 2*np.pi/xmax * (((np.array(range(0, 2*nmax+1))+nmax) % (2*nmax+1))-nmax)\n",
    "    xs = np.linspace(0,xmax,2*nmax+1)\n",
    "    aes = np.fft.fft(np.array(paquete0(xs)))\n",
    "    def evol(t):\n",
    "        return  fft.ifft(aes * np.exp(-1j*reldisp(ks)*t))\n",
    "    \n",
    "    return evol\n",
    "    \n",
    "    \n",
    "def rel_dispersion_lineal(k):\n",
    "    return k\n",
    "\n",
    "def rel_dispersion_cuadratica(k):\n",
    "    return k**2+1\n",
    "\n",
    "\n",
    "def rel_dispersion_relativista(k):\n",
    "    return np.sqrt(k**2+1)\n",
    "\n",
    "\n",
    "def rel_dispersion_anomala(k):\n",
    "    return -30*(1/(k+.1*1j)).real\n",
    "\n",
    "\n",
    "paquetes={\"gaussiano\":paquete_inicial_gaussiano,\n",
    "          \"cuadrado\": paquete_inicial_cuadrado,\n",
    "          \"triangular\": paquete_inicial_triangular,}\n",
    "\n",
    "reldisps={\"lineal\":rel_dispersion_lineal,\n",
    "        \"cuadratica\": rel_dispersion_cuadratica,\n",
    "         \"anomala\": rel_dispersion_anomala,}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T17:28:47.474251Z",
     "start_time": "2020-05-08T17:28:42.392229Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "73a2669976094ba9a927e07e371f959c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Label(value='listo...')"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "afe170ad2a124019a082a05463d0e44c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='tipo pulso:', options=('gaussiano', 'cuadrado', 'triangular'), val…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "paquete = None\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "tipopulso = widgets.Dropdown(options=['gaussiano','cuadrado',\"triangular\"],\n",
    "    description='tipo pulso:',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "tiporeldisp = widgets.Dropdown(options=['lineal','cuadratica',\"anomala\"],\n",
    "    description='relación de dispersión:',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "\n",
    "klimslider = widgets.FloatLogSlider(base=2.,min=6.,max=13.,step=1.,\n",
    "    description='k maximo',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "k0slider = widgets.FloatSlider(min=0.,max=10.,step=1.,\n",
    "    description='k0',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "x0slider = widgets.FloatSlider(min=2,max=xmax,step=1.,\n",
    "    description='x0',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "aslider = widgets.FloatLogSlider(base=2.,min=-2.,max=np.log(xmax)/np.log(2.),step=.5,\n",
    "    description='semi ancho pulso',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "\n",
    "status = widgets.Label(\"listo...\")\n",
    "display(status)\n",
    "\n",
    "def definir_paquete(pulso,reldisp,x0,k0,a,kmax):\n",
    "    global paquete\n",
    "    global nmax\n",
    "    global xs\n",
    "    reldisp = reldisps[reldisp]\n",
    "    nmax = int(round(kmax * (xmax)/(2*np.pi)))\n",
    "    status.value = \"Espere: fabricando el pulso... \"\n",
    "    paquete = evol_paquete(paquetes[pulso](x0,k0,a), reldisp)\n",
    "    status.value = \"... listo\"\n",
    "    print(\"Paquete \"+ pulso + \" definido.\")\n",
    "    print(\"x0= \"+ str(x0))\n",
    "    print(\"k0= \"+ str(k0))\n",
    "    vg = 50*(reldisp(k0+.01)-reldisp(k0-.01) )\n",
    "    print(\"vel grupo en k0 = \"+ str(vg))\n",
    "    if k0==0:\n",
    "        vf = reldisp(.0001)/.0001\n",
    "    else:\n",
    "        vf = reldisp(k0)/k0\n",
    "        \n",
    "    print(\"vel fase en k0 = \"+ str( vf   ))\n",
    "    fig, axs = plt.subplots(1,2,figsize=(8,4))\n",
    "    axs[0].plot(xs,paquete(0).real)\n",
    "    axs[0].set_title(\"Pulso a t=0\")\n",
    "    axs[1].plot(ks,reldisp(ks),label=\"$\\\\omega(k)$\")\n",
    "    axs[1].plot(ks,abs(aes),label=\"F(k)\")\n",
    "    axs[1].legend()\n",
    "    axs[1].set_title(\"relación de dispersión y amplitudes\")\n",
    "    axs[1].set_xlabel(\"k\")\n",
    "    axs[1].set_xlim(-20/a,20/a)\n",
    "    \n",
    "    plt.show()\n",
    "    \n",
    "    ts = np.linspace(0,tmax,100)\n",
    "    rhos = np.array([abs(paquete(t))**2 for t in ts])\n",
    "    \n",
    "    energias = np.array([sum(rho) for rho in rhos])\n",
    "    print(\"U(0)=\",energias[0])\n",
    "    cms = np.array([np.dot(rho,xs) for rho in rhos])/energias\n",
    "    \n",
    "    fig, axs = plt.subplots(1,2,figsize=(8,4))\n",
    "    axs[0].plot(ts,energias)\n",
    "    axs[0].set_title(\"$\\\\rho_{p}(t)$\")\n",
    "    axs[0].set_ylim(0,2*energias[0])\n",
    "    axs[1].plot(ts,cms,label=\"$x_{cm}$\")\n",
    "    axs[1].plot(ts,ts*vg + cms[0],ls=\"-.\",label=\"$t v_{g}$\")\n",
    "    axs[1].plot(ts,ts*vf + cms[0],ls=\":\",label=\"$t v_{f}$\")\n",
    "    plt.legend()\n",
    "    axs[1].set_title(\"$x_{cm}(t)$\")\n",
    "    plt.show()\n",
    "    # Animación\n",
    "    fig, ax=plt.subplots(1,1)\n",
    "    gpaq, =ax.plot(xs,paquete(0),animated=True)\n",
    "    gcm = ax.scatter(np.array([cms[0]]),np.array([0]),animated=True)\n",
    "    ax.set_ylim(-5,5)\n",
    "\n",
    "    def init():\n",
    "        ax.set_xlim(0, xmax)\n",
    "        ax.set_ylim(-5, 5)\n",
    "        return gpaq,\n",
    "\n",
    "    def update(t):\n",
    "        gpaq.set_data(xs,paquete(t))\n",
    "        gcm.set_offsets([[cms[int(t/tmax*99)],0.]])\n",
    "        return (gcm,gpaq)\n",
    "\n",
    "    ani = FuncAnimation(fig, update, frames=ts,\n",
    "                        init_func=init, blit=True)\n",
    "\n",
    "    ani = ani.to_jshtml()\n",
    "    display(HTML(ani))\n",
    "    plt.close()\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "    \n",
    "interact=widgets.interact(definir_paquete,\n",
    "                          pulso=tipopulso,\n",
    "                          reldisp=tiporeldisp,\n",
    "                          x0=x0slider,\n",
    "                          k0=k0slider,\n",
    "                          a=aslider,\n",
    "                          kmax=klimslider\n",
    "                          )"
   ]
  }
 ],
 "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.6"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "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
}