Resolviendo una PDE: Ecuación de Schrödinger

La Ecuación de Schrödinger fue desarrollada por el célebre físico austríaco Erwin Schrödinger en la década de los 20 y describe la evolución espacial y temporal de una partícula a través de su función de onda $\Psi(x,t)$. En general la ecuación posee la forma: \begin{align} i \hbar \frac{\partial\Psi}{\partial t} &= \hat{H} \Psi \end{align} donde $i$ es la unidad imaginaria y $\hbar$ es la constante reducida de Planck $\displaystyle \frac{h}{2\pi}$. Para una partícula no relativista, el operador $\hat{H}$ conocido como el Hamiltoniano (que es el total de energías cinéticas y potenciales de partículas en un sistema) queda descrito en función de un potencial $V(x,t)$, la masa inercial $m$ y el Laplaciano $\nabla^2$. La ecuación (1) entonces toma la siguiente forma: \begin{align} i \hbar \frac{\partial\Psi}{\partial t} &= \left(\frac{-\hbar^2}{2m}\nabla^2 + V\right) \Psi\\ \end{align} donde $x \in [-L, L]$ y $t > 0 $. La solución a esta ecuación entonces es la descripción de un sistema cuántico, la función de onda de una partícula. En general, esta onda es compleja, posee una parte real y una imaginaria. Este código desarrolla la ecuación en ausencia de potencial, es decir $V(x,t)=0$.



In [1]:

    
import numpy as np
from scipy.constants import hbar, electron_mass as me, proton_mass as mp
from scipy.integrate import fixed_quad
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
%matplotlib notebook

1. Gráfico de Condición inicial



In [2]:

    
# Initial conditions, outside for plotting.
f_re = lambda x: amp*np.exp(-(x-xc)**2.0/s)*np.cos(2.0*np.pi*(x-xc)/wl)
f_im = lambda x: amp*np.exp(-(x-xc)**2.0/s)*np.sin(2.0*np.pi*(x-xc)/wl)



In [3]:

    
# Center of wave
xc = 0
# Wavelength
wl = 4
# Exponential decay coeficient
s = 3
# Amplitude
amp = 0.5

xx = np.linspace(-3, 3, 400)
yy_re = f_re(xx)
yy_im = f_im(xx)

plt.figure(figsize=(12,6))
plt.ylim(-1, 1)
plt.plot(xx, yy_re, 'b-', lw=2.0, label=r"$\Re(\Psi_0)$")
plt.plot(xx, yy_im, 'r-', lw=2.0, label=r"$\Im(\Psi_0)$")
plt.grid(True)
plt.legend(loc='best')
plt.draw()

2. Resolución mediante Forward Difference con diferencias de primer orden



In [4]:

    
def schroedinger1D(xl, xr, yb, yt, M, N, D, xc, wl, s, amp):
    """
        Schrödinger Equation Simulation (with no potential, i.e. V(x,t) = 0)
        Args:
            [xl, xr]:    Space interval
            [yb, yt]:    Time interval
            M, N:        Space and time steps
            D:           Custom coefficient (should be hbar/(2m))
            xc:          Center for initial condition
            wl:          Wavelength for initial condition
            s:           Decay for initial condition
            amp:         Amplitude for initial condicion
    """
    # Normalize wave
    f = lambda x: f_re(x)**2 + f_im(x)**2
    area = fixed_quad(f, xl, xr, n=5)[0]
    f_real = lambda x: f_re(x)/area
    f_imag = lambda x: f_im(x)/area
    # Boundary conditions for all t
    l = lambda t: 0*t
    r = lambda t: 0*t
    # Effective mass
    #m = mp*me/(mp+me)
    # "Diffusion coefficient"
    #D = hbar / (2*m)
    # Step sizes and sigma constant
    h, k = (xr-xl)/M, (yt-yb)/N
    m, n = M-1, N
    sigma = D*k/(h**2)
    print("k=%f" % k)
    print("h=%f" % h)
    print("Sigma=%f" % sigma)
    # Finite differences matrix
    A_real = np.diag(2*sigma*np.ones(m)) + np.diag(-sigma*np.ones(m-1),1) + np.diag(-sigma*np.ones(m-1),-1)
    A_imag = -A_real
    # Left boundary condition u(xl,t) from time yb
    lside = l(yb+np.arange(0,n+1)*k)
    # Right boundary condition u(xr,t) from time tb
    rside = r(yb+np.arange(0,n+1)*k)
    # Initial conditions
    W_real = np.zeros((m, n+1))
    W_imag = np.zeros((m, n+1))
    W_real[:,0] = f_real(xl + np.arange(1,m+1)*h)
    W_imag[:,0] = f_imag(xl + np.arange(1,m+1)*h)
    # Solving for imaginary and real part
    for j in range(0,n-1):
        b_cond = np.concatenate(([lside[j]], np.zeros(m-2), [rside[j]]))
        W_real[:,j+1] =  W_real[:,j] + A_real.dot(W_imag[:,j]) - sigma*b_cond
        W_imag[:,j+1] =  W_imag[:,j] + A_imag.dot(W_real[:,j]) + sigma*b_cond
    print("System solved.")
    return np.vstack([lside, W_real, rside]).T, np.vstack([lside, W_imag, rside]).T



In [5]:

    
xl, xr, yb, yt, M, N, D, xc, wl, s, amp = (-3, 3, 0, 4, 16, 2000, 1.0, 0, 4.0, 3.0, 1.0)
W_real, W_imag = schroedinger1D(xl, xr, yb, yt, M, N, D, xc, wl, s, amp)









    



k=0.002000
h=0.375000
Sigma=0.014222
System solved.

3. Resultados gráficos



In [6]:

    
[X, T] = np.meshgrid(np.linspace(xl, xr, M+1), np.linspace(yb, yt, N+1))
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)
matrix = ax.imshow(W_real, aspect='auto')
ax.set_xlabel("$x$", fontsize=20)
ax.set_ylabel("$t$", fontsize=20)
ax.set_xticks(np.linspace(0,M,5))
ax.set_xticklabels(np.linspace(xl,xr,5))
ax.set_title(r"Real part of $\Psi(x,t)$", fontsize=20)
fig.colorbar(matrix)
plt.draw()



In [7]:

    
# Plot results
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel("$x$", fontsize=20)
ax.set_ylabel("$t$", fontsize=20)
ax.set_zlabel("$\Psi(x,t)$", fontsize=20)
ax.set_xlim(xl, xr)
ax.set_title("Real part of $\Psi(x,t)$")
surface = ax.plot_surface(X, T, W_real, cmap=cm.jet, linewidth=0, antialiased=True, rstride=1, cstride=1)
fig.colorbar(surface)
plt.tight_layout()
plt.draw()



In [8]:

    
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)
matrix = ax.imshow(W_imag, aspect='auto')
ax.set_xlabel("$x$", fontsize=20)
ax.set_ylabel("$t$", fontsize=20)
ax.set_xticks(np.linspace(0,M,5))
ax.set_xticklabels(np.linspace(xl,xr,5))
ax.set_title(r"Imaginary part of $\Psi(x,t)$", fontsize=20)
fig.colorbar(matrix)
plt.draw()



In [9]:

    
# Plot results
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel("$x$", fontsize=20)
ax.set_ylabel("$t$", fontsize=20)
ax.set_zlabel("$\Psi(x,t)$", fontsize=20)
ax.set_xlim(xl, xr)
ax.set_title("Imaginary part of $\Psi(x,t)$")
surface = ax.plot_surface(X, T, W_imag, cmap=cm.jet, linewidth=0, antialiased=True, rstride=1, cstride=1)
fig.colorbar(surface)
plt.tight_layout()
plt.draw()

Una idea interesante es discutir la estabilidad del método, puesto que ajustar parámetros para lograr estabilidad en este caso parece dificil. El análisis nos dirá si es posible obtener combinaciones de $k$ y $h$ tales que nuestros errores no se amplifiquen.

Disclaimer

El presente notebook ha sido creado para el curso ILI286 - Computación Científica 2, del Departamento de Informática, Universidad Técnica Federico Santa María. El material ha sido creado por Alejandro Sazo (asazo@alumnos.inf.utfsm.cl). En caso de encontrar un error, por favor no dude en contactar al email especificado. Puede encontrar la última versión del código en https://github.com/asazo/CC2



In [ ]: