El oscilador armónico cuántico en Python 3

Hola amigx cuánticx! Seguro que si trabajás o sos un entusiasta de la Física algo de osciladores armónicos viste en tu vida. Y sí… sabemos resolver tan pocas cosas de manera analítica que las exprimimos hasta el fin. Hoy voy a contarte sobre el oscilador armónico cuántico, pero desde un punto de vista completamente numérico. Me parece que es algo que si bien está en muchos textos, puede ser interesante para tener una primera visión de los problemas numéricos que aparecen en las simulaciones del quehacer científico. Pasemos a definir nuestro problema: El Hamiltoniano de un oscilador unidimensional tiene la siguiente pinta:

(1) $\begin{equation*}H = \frac{p^2}{2m}+\frac{1}{2}m\omega^2x^2\end{equation*}$

Si utilizamos unidades atómicas, en la representación de coordenadas la Ecuación de Schrödinger independiente del tiempo será

(2) $\begin{equation*}-\frac{1}{2m}\frac{d^2}{dx^2}\psi\left(x\right)+\frac{1}{2}m\omega^2x^2 \psi\left(x\right) = E \psi\left(x\right), \end{equation*}$

que es un problema de autovalores, y que es bien conocido. Si consideramos que la partícula es un electrón $(m=1)$ , las energías toman los valores $E_n = \omega(n+1/2)$ , y las autofunciones correspondientes son:

(3) $\begin{equation*}\psi_n(x)=\sqrt{\frac{1}{2^nn!}}\left(\frac{\omega}{\pi}\right)^{\frac{1}{4}}e^{-\omega x^2 /2}\end{equation*}$

Las soluciones analíticas se pueden encontrar en la mayoría de los libros de texto [1,2] inclusive en Wikipedia en la versión en inglés donde está desarrollado de manera bastante completa y una versión reducida en español. En este artículo exploraremos un método para resolverlo de manera numérica en Python 3, utilizando la librería NumPy.

Como todos los métodos numéricos, lo que se hace es discretizar todo el problema en una grilla, que llamaremos $x_i$ , con $i={0,1,...N-1}$ . En este caso utilizaremos una grilla equi-espaciada. Esto nos permite utilizar el método de Numerov para realizar la integración. Este método nos dice que en una ecuación con la forma $\psi''(x)=f(x)\psi(x)$ , evaluada sobre una grilla, se puede integrar de manera iterativa de la siguiente manera:

(4) $\begin{equation*}\phi_{i+1} = 2\phi_i-\phi_{i-1}+\Delta^2 f_i \psi_i, \end{equation*}$

donde:

(5) $\begin{equation*}\phi_{i} = \psi_i\left(1-\frac{\Delta^2}{12} f_i \right). \end{equation*}$

En las ecuaciones (4) y (5), $\Delta$ es el espaciado de la grilla y $\psi_i$ es la función evaluada en el punto $x_i$ de la grilla. Esto nos permitirá integrar a partir de dos puntos previos el punto siguiente. Este procedimiento se puede repetir de manera iterativa partiendo de los dos puntos del borde izquierdo (o derecho) de la grilla hasta llegar al otro extremo. En el siguiente snippet les muestro las implementaciones hacia adelante y hacia atrás.

def numerov_fwd(u, f, delta, n_cut):  #u es la función, n_cut es el punto sobre la grilla hasta donde integrar
    #forward
    fi_0 = 0
    fi_1 = u[1]*(1-delta**2/12*f[1])
    for i in range(1,n_cut+1):
        fi_2 = 2*fi_1 - fi_0 + delta**2*f[i]*u[i]
        u[i+1] = fi_2/(1-delta**2/12*f[i+1])
        fi_0 = fi_1
        fi_1 = fi_2
    return u
def numerov_bwd(u, f, delta, n_cut):
    #backward
    fi_0 = u[n_grid - 1]*(1-delta**2/12*f[n_grid-1])
    fi_1 = u[n_grid - 2]*(1-delta**2/12*f[n_grid-2])
    for i in range(n_grid-2,n_cut-1,-1):
        fi_2 = 2*fi_1 - fi_0 + delta**2*f[i]*u[i]
        u[i-1] = fi_2/(1-delta**2/12*f[i-1])
        fi_0 = fi_1
        fi_1 = fi_2
    return u

Es claro que para que la solución encontrada sea la solución correspondiente a una autofunción de nuestro problema, el autovalor provisto en $f(x)$ y las condiciones de contorno deben ser las exactas. Para resolver un problema de autovalores como el que describe la Ecuación (2), donde no se conocen los autovalores, conviene explotar la continuidad de la función de onda y su derivada, utilizando un método de disparo como veremos en el apartado que sigue.

Método de disparo y el problema de autovalores y autovectores

Para resolver este problema, integraremos la Ecuación (2) hasta un punto intermedio, partiendo tanto por la derecha, como por la izquierda. Como las soluciones, si existen son únicas, para encontrarlas seguiremos el siguiente esquema (voy a utilizar la notación de los arrays de NumPy para los índices donde los valores negativos significa contar desde el último):

Dadas dos condiciones de borde por izquierda: $\psi[0], \psi[1]$ , y otras por derecha: $\psi[-1], \psi[-2]$ , integraremos hasta un punto intermedio para un valor propuesto de la energía $E$ , que determina la integración de la ecuación como parte de $f(x) = -2 m (E-\frac{1}{2}m \omega^2 x^2)$
En ese punto intermedio igualamos las soluciones encontradas por izquierda y por derecha, y calculamos la derivada en ese punto mediante alguna aproximación. Si la derivada no es igual por derecha e izquierda dentro de un error determinado, quiere decir que la energía no corresponde a un autovalor. Esto se puede asegurar porque la derivada logarítmica de la función de onda debe ser continua, es decir: $\left. \frac{d (ln\psi)}{dx}\right|_{izq} = \left. \frac{d (ln\psi)}{dx}\right|_{der}$
Este esquema se puede repetir, variando el valor de la energía de prueba hasta conseguir que la diferencia en los valores de la derivada obtenidos por izquierda y derecha difieran en un valor menor a una tolerancia determinada. Una vez que se alcanza este valor, podemos decir que hemos encontrado el autovalor, y la autofunción correspondiente (no olvidar que se deben normalizar).

Primero construyamos nuestra grilla radial

#parametros de grilla
n_grid = 1024     #numero de puntos de grilla
x_max = 15        #distancia del origen donde se considera el infinito
x_grid = np.linspace(-x_max,x_max,n_grid) #grilla espacial
delta = x_grid[1]-x_grid[0] #paso de la grilla

Resolvamos para un sistema de masa $m=1$ y frecuencia $\omega=1$ . Para esto veamos qué sucede para una energía $E = 0.3$ a.u.

E_try=0.3
f = -2*(E_try-0.5*x_grid**2)         #funcional de la energía del integrador de Numerov
n_corte=n_grid//2                       #lugar donde se deben intersectar las funciones
psi_fwd=np.zeros(n_grid)           #defino el array de la función
psi_fwd[0] = 1e-4
psi_fwd[1] = 1 #defino las condiciones de borde
psi_fwd = numerov_fwd(psi_fwd.copy(), f, delta, n_corte+1) #evoluciono un paso mas para derivar por el punto medio
psi_fwd = psi_fwd/psi_fwd[n_corte] #normalizo a 1 el punto donde corto

Figura 1. Integración hacia adelante para E=0.3

En la Figura (1) vemos la integración numérica de izquierda a derecha utilizando el algoritmo de Numerov. Hemos elegido integrar hasta la mitad de los puntos más uno, y forzamos a la solución para que tome el valor 1 en la mitad de la grilla. Es importante aclarar que las condiciones de borde no tienen que ser exactas, ya que el proceso de normalización lleva a un cambio de escala en la totalidad de la función. Esta es una de las grandes ventajas de este método donde no se conocen exactamente las condiciones asintóticas de las soluciones. Por otra parte siempre hay que tener una estimación adecuada debido a que la solución puede tomar valores que escapan de los 64 bits que tenemos para representar números de punto flotante.

Se puede hacer lo mismo partiendo desde el lado derecho, y fijarnos en la diferencia entre las derivadas por derecha y por izquierda para el punto elegido para hacer el encuentro entre las funciones.

psi_bwd=np.zeros(n_grid)
psi_bwd[n_grid-1] = 1e-4
psi_bwd[n_grid-2] = 1
psi_bwd = numerov_bwd(psi_bwd.copy(), f, delta, n_corte-1)
psi_bwd = psi_bwd/psi_bwd[n_corte]

Figura 2. Integración hacia adelante y atrás para una E=0.3

En este caso vemos que si bien la función está forzada a ser continua, porque le hemos exigido que en el punto medio de la grilla tome el valor unidad, las derivadas por derecha e izquierda no son las mismas. Para ver los valores que toma la derivada podemos hacer el cálculo por diferencias finitas en el punto medio (por esto es que integramos un paso extra del punto de encuentro):

der_fwd = (psi_fwd[n_corte+1]-psi_fwd[n_corte-1])/2/delta
der_bwd = -(psi_bwd[n_corte-1]-psi_bwd[n_corte+1])/2/delta
print('deriv por izquierda:',der_fwd)
print('deriv por derecha:',der_bwd)
#out:
deriv por izquierda: 0.3028597909587141
deriv por derecha: -0.32331840338242146

Donde vemos que la diferencia entre los valores de la derivada es notable. Veamos qué sucede cuando calculamos para el valor exacto de la energía del estado fundamental siguiendo el mismo procedimiento. En la Figura (3) vemos la función integrada por derecha e izquierda, donde se aprecia que por derecha e izquierda ambas soluciones se acercan de modo suave. Los valores que se obtienen para la derivada son los siguientes:
derivada por izq: -0.014656459323199651
derivada por der: -0.014656449080109192.

Figura 4. Soluciones numéricas para la energía exacta del estado fundamental E=0.5

¡La diferencia en el error absoluto de las derivadas se encuentra recién en la octava cifra decimal! Lo que nos dice que nos encontramos ante un método sumamente potente para resolver este tipo de problemas. Para completar el procedimiento basta con normalizar la función. Recordemos que si:

(6) $\begin{equation*}\int_{-\infty}^{\infty}\left| \psi \left( x\right) \right|^2 dx = N\end{equation*}$

La función se puede normalizar dividiendo por $\sqrt{N}$ , como vemos a continuación mediante el método del trapecio:

psi = np.hstack([psi_fwd[:n_corte],psi_bwd[n_corte:]]) #unifico las soluciones
psi = psi/np.sqrt(np.sum(psi[1:]**2+psi[:-1]**2)*delta/2) #normalizo integrando por trapecio

Y así hemos encontrado de un modo completamente numérico el estado fundamental del oscilador armónico. No está de más que verifiquen que se corresponde con la solución analítica, y lo recomiendo enfáticamente.

Hasta aquí llega este artículo, en uno próximo veremos qué sucede con los errores como función de la energía, y cómo sistematizar la búsqueda de autovalores y autofunciones de los estados excitados. Dejen comentarios, compartan y esas cosas ¡Hasta pronto!

[1] Eisberg, R. M., & Resnick, R. (1985). Quantum physics of atoms, molecules, solids, nuclei, and particles (2nd ed.). New York: Wiley.
[2] Cohen-Tannoudji, Claude, Bernard Diu, and Franck Laloë. 2005. Quantum mechanics. New York: Wiley

Importante

La guitarra científica

El oscilador armónico cuántico en Python 3

Método de disparo y el problema de autovalores y autovectores

Deja un comentario Cancelar respuesta