Física Computacional

Resolución de ecuaciones en derivadas parciales: la ecuación de

Schrödinger
Física Computacional Resolución de ecuaciones en derivadas parciales: la ecuación de Schrödinger 2
1. Fundamento Teórico

En Mecánica cuántica, el estado físico de un sistema unidimensional de una partícula viene descrito
completamente por una función de onda compleja Φ(x, t) que obedece la denominada ecuación de Schrödinger
∂Φ(x, t)  ~2 ∂ 2
 

i~ = − + V (x) Φ(x, t) = HΦ (1)

∂t 2m ∂x2
donde V (x) es el potencial al que está sometida la partícula, que supondremos que tiene una masa m
independiente del tiempo, ~ es la constante de Planck reducida y H es el operador Hamiltoniano, que es hermítico.
La densidad de probabilidad de encontrar a la partícula en un volumen dV alrededor del punto x y en el instante
t es
dP = |Φ(x, t)|2 dV. (2)
Normalización: Para un sistema compuesto únicamente por una partícula, como es el nuestro, la probabilidad
total de encontrar la partícula en algún punto del espacio, en cualquier instante t, es igual la unidad,
Z

V
|Φ(x, t)|2 dV = 1. (3)
Nótese que la integral de probabilidad es constante durante la evolución temporal, lo que habrá que tenerse en
cuenta al plantear la resolución numérica del problema.
Podemos observar que la ecuación de Schrödinger es lineal y que se trata de una ecuación diferencial de primer
orden en el tiempo.
√
Para evitar el uso continuado de m y ~, haremos el cambio de variable t → t~ en el tiempo y x → x~/ 2m en
el espacio, lo que conduce a una ecuación equivalente a (1) pero en la que ~ = 1 y m = 1/2:
∂2
 
∂Φ
− HΦ ≡  − V (x)  Φ(x, t) = −i . (4)
∂x2 ∂t
Física Computacional Resolución de ecuaciones en derivadas parciales: la ecuación de Schrödinger 3
Esta ecuación puede resolverse términos de las autofunciones del hamiltoniano H, que, recordemos, es
independiente del tiempo.
Llamando un a los estados acotados y uk a los estados del continuo, entonces la solución dependiente del tiempo
viene dada de forma rigurosa por
Z
an e−iEn t un (x) + dka(k)e−iE(k)t uk (x),
X
Φ(x, t) = (5)
n
con
an = hun , φ(x, 0)i, a(k) = huk , φ(x, 0)i, (6)
donde Φ(x, 0) es el estado inicial del sistema (y por lo tanto an y a(k) son los coeficientes Φ en las bases {un } y
{uk }, respectivamente).
De seguir este procedimiento, toparíamos con los inconvenientes de que:
• En la mayoría de los casos, las autofunciones y los autovalores del hamiltoniano no se pueden calcular
analíticamente y han de determinarse numéricamente.
• Por motivos numéricos, el estado inicial debe restringirse necesariamente a una combinación lineal finita
de autofunciones y, en cualquier caso, el número de éstas debe ser el menor posible para reducir el trabajo
computacional. Truncar la base tiene como inconveniente añadido que la evolución deja de ser unitaria.
Por todo ello, en lugar de seguir el procedimiento anterior integrararemos directamente la ecuación (4). Pero hay
que ser muy cuidadosos con el método numérico porque, en general, la discretización provoca que no se conserve
la normalización de la función de onda (ecuación 3).
La solución formal de la ecuación de Schrödinger es

Φ(x, t) = e−i(t−t0 )H Φ(x, t0 ) (7)

Física Computacional Resolución de ecuaciones en derivadas parciales: la ecuación de Schrödinger 4
donde el operador e−itH es unitario (su inverso coindice con su adjunto) por ser H hermítico. Esta última
propiedad será de utilidad a la hora de hallar la discretización más adecuada para resolver el problema
numéricamente como veremos a continuación.

2. Método numérico

Discretizaremos el espacio y el tiempo tomando xj = jh y tn = ns, con j = 0, 1, ..., N y n = 0, 1, 2, ..., donde

h es el espaciado en la discretización espacial y s es el espaciado temporal.
La función de onda viene dada en cada punto del retículo espacio-temporal por:

Φ(xj , tn ) → Φ(jh, ns) = Φj,n con j = 0, 1, ..., N y n = 0, 1, 2, .... (8)

Supondremos que las condiciones de contorno para la función de onda en j = 0 y j = N son las correspondientes
a la existencia de un potencial infinito en esos puntos, esto es, la densidad de probabilidad de encontrar la
partícula en dichos puntos es cero. Así, Φ0,n = ΦN,n = 0, en cualquier instante.
Puede definirse de forma inmediata un primer algoritmo a partir de la expresión

Φj,n+1 = e−isH Φj,n . (9)

Si s es muy pequeño, el operador de evolución exp (−isH) puede aproximarse por su desarrollo de Taylor a
primer orden en t. Si además observamos que el operador Hamiltoniano no es más que la aplicación de una
derivada segunda, su discretización espacial es obvia y en total obtenemos el siguiente algoritmo

Φj,n+1 = 1 − isHD + O(sHD )2 Φj,n ,

 
(10)
Física Computacional Resolución de ecuaciones en derivadas parciales: la ecuación de Schrödinger 5
donde
∂ 2Φ 1
Φ00j ≡ = (Φ j+1 − 2Φ j + Φ j−1 ) + O(h2
) (11)
∂x2 h2
y HD viene dado por
1
HD fj = − (fj+1 − 2fj + fj−1 ) + Vj fj (12)
h2
con Vj = V (jh).
Este algoritmo tan simple y natural adolece de un grave problema: el operador (1 − isH) no es unitario. Esto
implica que durante la iteración del algoritmo para encontrar la evolución de la función de onda, la normalización,
2
j h|Φj,n | , irá variando con el tiempo, lo que es incompatible con la restricción de normalización a la unidad
P

expresada en (3) y viola la interpretación de Born.

Así pues, el objetivo será encontrar un operador evolución similar al anterior pero que sea exactamente unitario.
Esto se consigue haciendo uso de la aproximación de Cayley para el operador evolución temporal
1 − isHD /2
e−isH ≈ , (13)
1 + isHD /2
que conduce al algoritmo
1 − isHD /2
Φj,n+1 = Φj,n . (14)
1 + isHD /2
Nótese que además de ser unitario, este operador es exacto hasta orden (sHD )2 .
Para completar el algoritmo, sólo resta diseñar la estrategia para utilizar eficazmente el operador 1/(1 + isHD ).
Para ello reescribimos la anterior ecuación como
 
2
Φj,n+1 =  − 1 Φj,n = χj,n − Φj,n , (15)
1 + isHD /2
Física Computacional Resolución de ecuaciones en derivadas parciales: la ecuación de Schrödinger 6
donde
2
χj,n ≡ Φj,n , (16)
1 + isHD /2
o lo que es lo mismo, dada Φj,n , j = 0, .., N , χj,n es la solución de la ecuación

[1 + isHD /2] χj,n = 2Φj,n , (17)

que en forma explícita puede escribirse como

 
2i 4i
χj+1,n + −2 + − Ṽj  χj,n + χj−1,n = Φj,n (18)
s̃ s̃

donde s̃ = s/h2 y Ṽj = h2 Vj .

De esta forma, el algoritmo de evolución es claro: dada una función de onda en el instante n para toda posición
j se resuelve el conjunto de ecuaciones (18) y se obtiene χj,n (j = 0, .., N ). Con esta solución se recurre a la
ecuación (15) para obtener las nuevas Φj,n+1 (j = 0, ..., N ), iterándose el proceso.
Sólo falta por conocer cómo resolver ecuaciones del tipo (18), es decir, cómo invertir matrices tridiagonales.
En este caso, hemos de resolver un conjunto de ecuaciones del tipo

A− 0 +
j χj−1,n + Aj χj,n + Aj χj+1,n = bj,n j = 1, . . . , N − 1 (19)

donde A− 0 +
j = 1, Aj = −2 + 2i/s̃ − Ṽj , Aj = 1 y bj,n = 4iΦj,n /s̃.

Las condiciones de contorno son χ0 = χN = 0 (notemos que estas condiciones de contorno implican que en (15)
Φ0 y ΦN son nulas en cualquier instante).
Física Computacional Resolución de ecuaciones en derivadas parciales: la ecuación de Schrödinger 7
Para resolver la recurrencia (19) suponemos que su solución es del tipo
χj+1,n = αj χj,n + βj,n j = 0, . . . , N − 1 (20)
donde, para garantizar que se cumple χN,n = 0, tomaremos αN −1 = βN −1,n = 0.
Sustituyendo esta expresión en (19) obtenemos
A−j bj,n − A+j βj,n
χj,n =− 0 χ j−1 + 0 + A+ α . (21)
Aj + A+j α j A j j j

Si identificamos las ecuaciones (20) y (21) podemos definir las recurrencias para los coeficientes α y β así
 
αj−1 = −A−
j γj , βj−1,n = γj bj,n − A+
j βj,n . (22)

donde γj−1 = A0j + A+

j αj .

Estas ecuaciones nos dan la forma de obtener todas las α y β partiendo de j = N − 1 y obteniendo en orden
decreciente αj y βj,n con j = N − 2, . . . , 1, 0.
Obsérvese que α no depende del tiempo y sólo es necesario calcularlas una vez.
Una vez obtenidas las α y β, se usa la recurrencia (20) para hallar las χj en orden de j crecientes.
Conocida χj,n , y con ella Φj,n+1 , queda discutir qué función de onda inicial se usa y con qué potencial. La función
de onda inicial que vamos a usar es una onda plana con una amplitud gaussiana, esto es,
2 2
Φ(x, 0) = eik0 x e−(x−x0 ) /2σ . (23)
Notemos que con esta elección, la densidad de probabilidad de encontrar inicialmente la partícula en un punto
x es una gaussiana centrada en x0 y de anchura σ.
Física Computacional Resolución de ecuaciones en derivadas parciales: la ecuación de Schrödinger 8
El número de oscilaciones completas que la función de onda tiene sobre la red depende de k0 . En particular,
k0 N h = 2πnciclos . En lugar de dar como parámetro inicial k0 , daremos nciclos .
Obviamente nciclos = 0, 1, ...., N . Pero físicamente no tendría mucho sentido que se produjese una oscilación
completa de la función de onda entre dos puntos del retículo, pues querría decir que la discretización no es
suficientemente fina. Así pues, restringiremos el parámetro nciclos a los valores 1, ..., N/4. De esta forma, un ciclo
tendrá 4 puntos como mínimo.
La posición media inicial y la anchura de la gaussiana serán por ejemplo x0 = N h/4 y σ = N h/16, aunque es
recomendable escribir la función de onda inicial de manera genérica y jugar con x0 y σ.
Por último, el potencial que usaremos tendrá una anchura N/5, estará entrado en N/2 y su altura será
proporcional a la energía de la función de onda incidente: λk02 (puede usarse, por ejemplo, λ = 0,3).
En resumen, utilizando las constantes reescaladas, la función de onda en el retículo se tomará
2 2
Φj,0 = eik̃0 j e−8(4j−N ) /N , (24)

donde k̃0 = k0 h = 2πnciclos /N donde nciclos = 1, ..., N/4.

El potencial es 

0 si j ∈
6 [2N/5, 3N/5],
Ṽj = Vj h2 = 

(25)
 λk̃02 si j ∈ [2N/5, 3N/5].

Sólo queda por fijar el parámetro s̃ = s/h2 . Puesto que la energía es proporcional a k02 y el operador dinámico
discreto tiende a ser exacto en potencias de Hs, lo óptimo es elegir ||H||s < 1, esto es, k02 s < 1. Así deducimos
que s̃ < 1/k̃02 . En particular tomaremos s̃ = 1/4k̃02 .
Física Computacional Resolución de ecuaciones en derivadas parciales: la ecuación de Schrödinger 9
Resumiendo, los parámetro que se han de fijar inicialmente son N , nciclos y λ, pues todos los demás se determinan
a partir de ellos.
Por lo tanto, el algoritmo para resolver la ecuación de Schrödinger unidimensional puede esquematizarse del
siguiente modo:
1. Dar los parámetros iniciales: N , nciclos y λ. Generar s̃, k̃0 , Ṽj , Φj,0 (incluyendo las condiciones de contorno
Φ0,0 = ΦN,0 = 0) y α.
2. Calcular β utilizando la recurrencia (22).
3. Calcular χ a partir de (20).
4. Calcular Φj,n+1 de (15).
5. n = n + 1, ir a al paso 2.

3. Problemas

Obligatorio: Resolver la ecuación de Schrödinger unidimensional para un potencial cuadrado. Comprobar que
se conserva la norma.
Voluntario 1: Estudiar el coeficiente de transmisión
El coeficiente de transmisión es la probabilidad de encontrar a la partícula al otro lado del obstáculo para tiempos
largos. En sistemas clásicos este coeficiente es 1 si la energía de la partícula es mayor que la energía del escalón y
0 si es menor. En sistemas cuánticos esto cambia por la acción del efecto túnel. Para determinar si una partícula
es reflejada o transmitida debemos colocar detectores al final del sistema. Como ya hemos visto, la probabilidad
de encontrar a la partícula entre los puntos x1 y x2 viene dada por
Física Computacional Resolución de ecuaciones en derivadas parciales: la ecuación de Schrödinger 10
Z x
|Φ(x)|2 dx
2
P (x1 , x2 ) = x1
(26)

Suponemos que tenemos detectores finitos, actuando a derecha e izquierda de la barrera. Si estos detectores
tienen un ancho de N/5 la probabilidad a tiempo n de detectar la partícula a la derecha vendrá dado por
2 PN/5 2
PT (n) = N i=N/5 |Φi,n | , y la probabilidad de detectarla a la izquierda sería PR (n) = i=0 |Φi,n | . Después de
P

realizar el experimento m veces el coeficiente de transmisión se calcula como K = NT /m, donde NT es el número
de veces que hemos detectado la partícula a la derecha del potencial.
Si la partícula es detectada no hace falta continuar con la simulación, pero si no se detecta hay que proyectarla.
Esto significa que a cada paso que no haya detección hay que hacer los coeficientes Φi = 0 para i ∈ [N/5, N ],
si hemos aplicado el detector derecho, y Φi = 0 para i ∈ [0, N/5] si ha sido el izquierdo. Para garantizar la
normalización tendremos que calcular tras cada paso el valor

N
|Φj |2
X
k= (27)
j=0

y multiplicar todos los elementos de la función de onda de la manera Φj = √1 Φj .

k
3
Para obtener una buena estadística habrá que simular el sistema al menos 10 veces.
Hay que realizar:

• Estudiar la dependencia en N de K realizando simulaciones para N = 500, 1000, 2000.

• Estudiar la dependencia en V (x) de K, usando valores λ = 0, 1; 0,3; 0,5; 1; 5; 10.
• Comparar con resultados teóricos.
Física Computacional Resolución de ecuaciones en derivadas parciales: la ecuación de Schrödinger 11
Hay que tener en cuenta que un detector que funcione a cada paso afecta fuertemente la dinámica del sistema
(comprobadlo). Así habrá que aplicar los detectores sólo tras un intervalo de tiempo. Este intervalo debe ser lo
suficientemente grande para dejar el sistema evoluciona, y lo suficientemente pequeño como para evitar reflexiones
en las paredes. Los resultados también se pueden analizar en función de este parámetro. Discutid este fenónmeno.

Números complejos en Fortran.

Declaración: implicit double complex (h)

Asignación de valores constantes: h = (1,0d0, 0,3d0) (=1+i0.3).
Asignación de variables: h = complex(a, b) (=a+ib).