<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Discretizando la ecuación de Schrödinger 2. Operador evolucion.
Autor
Step1:  Código que controla las simulaciones
Step2:  Rutinas que generan inicializan y generan la simulación
Step3:  Rutinas que conectan los controles de la primera celda con la simulacion
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
import ipywidgets as widgets
from IPython.display import display
L = 200
dx = .5
U = lambda x:(0.*x)
buttonrunsim = widgets.Button(description="Simulate")
buttongenerateu = widgets.Button(description="Generate U")
outwdt = widgets.Output()
centropaquete = widgets.FloatSlider(value=-3*L/4,min=-L,max=L,step=1,
    description='centro del paquete',
    ensure_option=True,
    disabled=False)
anchopaquete = widgets.FloatSlider(value=L/8,min=1,max=L,step=1,
    description='ancho del paquete',
    ensure_option=True,
    disabled=False)
valuek0 = widgets.FloatSlider(value=0,min=-np.pi/dx,max=np.pi/dx,step=.1,
    description='k0',
    ensure_option=True,
    disabled=False)
valueV0 = widgets.FloatSlider(value=0,min=-4,max=4,step=.1,
    description='V0',
    ensure_option=True,
    disabled=False)
potencial = widgets.Dropdown(options=["Lineal","Cuadrado","Rampa","Cuadrático","Gaussiano"],
                            description="Forma del Potencial",
                             ensure_option=True,
                            disabled=False)
valueW = widgets.FloatSlider(value=1.,min=1.,max=300.,step=.5,
    description='Alcance potencial',
    ensure_option=True,
    disabled=False)
items = [[widgets.Label("Paquete"),centropaquete,anchopaquete,valuek0,buttonrunsim,
          widgets.Label("Potencial"),valueW,valueV0,potencial,buttongenerateu],[outwdt]]
widgets.HBox([widgets.VBox(it) for it in items])
Explanation: Discretizando la ecuación de Schrödinger 2. Operador evolucion.
Autor: Juan Mauricio Matera (Facultad de Ciencias Exactas, UNLP) 
23 de agosto de 2020
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:
$$
\psi(\vec{x},t)=\sum_{E} a(E) e^{-{\bf i}E t/\hbar}\varphi_E(\vec{x})
$$
con 
$$
E \varphi_i(\vec{x})=-\frac{\hbar^2}{2 m}\nabla^2 \varphi_i(\vec{x}) +V(x)\varphi_i(\vec{x})
$$
que es una ecuación de autovalores para el operador ${\bf H}=-\frac{\hbar^2}{2 m}\nabla^2+V(x)$. 
Como ${\bf H}$ es un operador autoadjunto, sus autofunciones $\varphi_E(\vec{x})$ expanden una base ortogonal,
lo que permite fijar los coeficientes $a(E)$ a partir de la condición inicial usando la ortogonalidad:
$$
a(E)=\frac{\int \varphi_E(\vec{x})^\psi(\vec{x},0) dx}{\int \varphi_E(\vec{x})^\varphi_E(\vec{x}) dx}
$$
Podemos también expresar en forma compacta la solución, introduciendo el operador evolución
$$
{\bf U}(t)=\exp(-{\bf i}{\bf H}t/\hbar)
$$
de manera que
$$
\psi(\vec{x},t+ \Delta t)={\bf U}(\Delta t)\psi(\vec{x},t)=
\sum_{E} a(E) e^{-{\bf i}E \Delta t/\hbar}\varphi_E(\vec{x})
$$
Discretizació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
$$
\varphi_n(x)=\exp({\bf i} \frac{\pi n}{L} x)/\sqrt{2L}
$$
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.
Por otro lado, podemos combinar estas funciones para obtener una base complementaria
$$
\xi_{m}(x)=\sum_{n=-N}^{N}  \frac{e^{-{\bf i}\frac{\pi n m}{N+1/2}}}{\sqrt{(2N+1)}}\varphi_n(x)=
\frac{1}{\sqrt{2L(2N+1)}}\frac{\sin\left(\pi(N+1/2)(x/L-m/(N+1/2))\right)}
{\sin\left(\pi/2(x/L-m/(N+1/2))\right)}
$$
que cumplen las condiciones
\begin{eqnarray}
\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'}\
\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}}\
\frac{\xi_{m}(x_n)}{\sqrt{(N+1/2)/L}}&\approx& \delta_{m,n}
\end{eqnarray}
con $x_n=\frac{n L}{N+1/2}$. De estas propiedades, se sigue que 
$$
\psi(x)\approx \sum_{m=-N}^{N}  \frac{\xi_m(x)}{\sqrt{(N+1/2)/L}}\psi(x_m)
$$
y
$$
V(x)\psi(x)\approx \sum_{m=-N}^{N}  \frac{\xi_m(x)}{\sqrt{(N+1/2)/L}}V(x_m)\psi(x_m)
$$
de lo que el potencial es (para $N$ suficientemente grande) diagonal en la base $\xi_m(x)$. 
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)$:
$$
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)=
\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)
$$
con
$$
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)}}
$$
Restringidos a esta base,  cualquier solución de la ecuación de Schrödinger para una condición inicial suficientemente suave se expresa como
$$
\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}}
$$
con 
$$
U_{m,m'}=[\exp(-{\bf i}H(t-t_0)/\hbar]{m,m'}
$$
el operador de evolución en el subespacio de $\xi{m}(x)$, con
$$
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'}
$$
De esta manera, el problema se reduce a 
1.  Construir la matriz H. 
2. Diagonalizar la matriz $H$ como $H={\cal V}{\Omega}{\cal V}^\dagger$ con $\Omega$ una matriz diagonal y 
${\cal V}^\dagger$ construir el operador de evolución con un paso constante $\Delta t$: 
$$
U(\Delta t)={\cal V}\exp(-{\bf i}{\Omega}\Delta t/\hbar){\cal V}^\dagger
$$
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$
4. Construir el estado en el instante $t+n\Delta t$ mediante la recurrencia
$$
\Psi_{t+\Delta t}=U(\Delta t)\Psi_{t}
$$
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. 
Implementació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).
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)$.
Simulación
End of explanation
# Para fijar parámetros a mano
#### Parámetros de la discretización
dx = .5
dt = 10.
L = 200
####  Estado inicial
a = 20.
x0 = -L/3
k0 = 10./a
#### Potencial
U0= -1.
w = 10.
Explanation: Código que controla las simulaciones
End of explanation
import numpy.fft as fft
import numpy.linalg as la
############    Inicializar  #############
def init_evolop(L=L,dx=dx,dt=dt,pot=None):
    global rho0
    global psi
    global psi0
    global xs
    global evol_op
    global P, ham
    xs = np.linspace(-L,L,int(2*L/dx))
    ## Construyo el operador evolución
    nhalf = int(L/dx)
    # Energía cinética
    ks = np.pi/L*np.array([(1.*q)  if q<nhalf else  (q-2*nhalf) for q in range(2*nhalf)])
    Uft = np.array([fft.ifft(v)*np.sqrt(2*nhalf) for v in  np.eye(2*nhalf)])
    P = (ks*(Uft).transpose()).dot(Uft.conj())    
    ham = .5*P.dot(P)
    # Energía potencial
    if pot is None:
        U = 0.*xs
    else:
        U = pot(xs)
    ham = ham + np.diag(U)
    ens,states = np.linalg.eigh(ham)
    evol_op = (states.conj()).dot((np.exp(-1j*dt*ens)*states).transpose())
    return evol_op
init_evolop()
def init_state(L=L,dx=dx,dt=dt,a=a,x0=x0,k0=k0):
    global rho0
    global psi
    global psi0
    global xs
    global evol_op
    global P
    xs = np.linspace(-L,L,int(2*L/dx))
    rho0 = np.array([np.exp(-(x-x0)**2/(2*a**2))  for x in xs]) 
    rho0[0] = 0
    rho0[-1] = .5*rho0[-2]
    rho0[1] = .5*rho0[2]
    rho0 = rho0/sum(rho0)/dx
    psi0 = rho0**.5 * np.exp(1j*k0*(xs-x0))
    psi = psi0
    return psi
def Ulineal(xs):
    return -U0*xs/w
def Ugaussiano(xs):
    return U0*np.exp(-xs**2/(.5*w**2))
def Uescalon(xs):
    return np.array([U0 if abs(x)<.5*w else 0 for x in xs])
def Uarmonico(xs):
    return np.array([U0*(1-4*(x/w)**2) if abs(x)<.5 * w else 0 for x in xs])
def Urampadoble(xs):
    return np.array([U0*(1-2*abs(x)/w) if abs(x)<.5 * w else 0 for x in xs])
def make_animation_new(ts=50, progress=None):
    global psi
    global xs
    global U0,k0,w,x0,a,U, evol_op
    
    psi = init_state(a=a,dt=dt,k0=k0,x0=x0)
    rho0 = abs(psi)**2
    fig1 = plt.figure()
    plt.xlim(-L,L)
    plt.ylim(-abs(U0)-k0**2-.5,.5+abs(U0)+k0**2)
    plt.plot(xs, U(xs),ls="-.")
    density,= plt.plot(xs,50*(U0+.1)*rho0+.5*k0**2)
    def update_graphic(t,density):
        global psi, evol_op
        if progress is not None:
            progress.value = t
        if t>=1:
            psi = evol_op.dot(psi)
        rho = abs(psi)**2
        density.set_data(xs,50*(U0+.1)*rho+.5*k0**2)
        
        
        return density,
    line_ani = animation.FuncAnimation(fig1, update_graphic, ts, fargs=(density,),
                                       interval=100, blit=True)
    plt.close()
    return HTML(line_ani.to_jshtml())
    
Explanation: Rutinas que generan inicializan y generan la simulación
End of explanation
def on_buttongenerateu(b):
    global U, evol_op, w, U0,ham
    U0 = valueV0.value
    w = valueW.value
    if potencial.value == "Lineal":
        U = Ulineal
    elif potencial.value == "Gaussiano":
        U = Ugaussiano
    elif potencial.value == "Cuadrático":
        U = Uarmonico
    elif potencial.value == "Cuadrado":
        U = Uescalon
    elif potencial.value == "Rampa":
        U = Urampadoble
    else:
        print("No encontrado")
        return
    outwdt.clear_output()
    with outwdt:
        print("U=",U)
        print("U0=",U0)
        print("w=",w)
        evol_op = init_evolop(pot=U)
    outwdt.clear_output()
    with outwdt:
        plt.plot(xs, U(xs))
        display(plt.show())
        plt.ylim(-5,5)
        plt.xlim(0,100)
        plt.scatter(range(len(ham)),sorted(la.eigvalsh(ham)))
        plt.title("Hamiltonian Spectrum")
        display(plt.show())
        
def on_button_start_sim(b):
    global U0, w, k0, a, x0, U, U0, w
    k0 = valuek0.value
    a = anchopaquete.value
    x0 = centropaquete.value
    
    outwdt.clear_output()
    progress = widgets.FloatProgress(value=0,min=0,max=50,step=1,description='Simulating:',
    bar_style='info',orientation='horizontal')
    
    with outwdt:
        print([a,x0,k0,U0,w])
        display(progress)
        result = make_animation_new(progress=progress)     
        outwdt.clear_output()
        display(result)
        
buttonrunsim.on_click(on_button_start_sim)
buttongenerateu.on_click(on_buttongenerateu)
    
    
Explanation: Rutinas que conectan los controles de la primera celda con la simulacion
End of explanation
<END_TASK>