La ecuación de Schrödinger describe la evolución (determinista) de un campo
F llamado función de ondas. Esta ecuación en una dimensión es de la forma:

i _
h
 
  F(x,t)
t
= é
ê
ê
ê
ë
-
_
h
2

2m
  2
x2
+V(x) ù
ú
ú
ú
û
F(x,t) = HF
(1)

donde `h es la constante de Planck, F(x,t) es el campo que representa a una partícula en el instante t, m es la masa de dicha partícula, V(x) es el potencial al que se ve sometida y H es el operador Hamiltoniano. La función de ondas es la magnitud que contiene toda la información accesible a un observador externo al sistema (notar que es una variable compleja). Para evitar el uso continuado de m y `h, reescalamos el tiempo t® t`h y el espacio x® x`h/Ö[2m] lo que nos da una ecuación de Schrödinger equivalente a (1) pero con `h = 1 y m = 1/2. Esta ecuación reescalada será la que utilicemos a partir de ahora.

El carácter físico de la función de ondas fue dado por Max Born al interpretar que

P(x,t) = |F(x,t)|2
(2)
es proporcional a la probabilidad de encontrar la partícula en el punto x en el instante t. Para que sea una probabilidad, debemos de garantizar que
  ó
õ
¥

-
¥ 
dx P(x,t) = 1
(3)
esto es, las funciones de onda han de ser de cuadrado integrable luego han de pertenecer a un espacio de Hilbert. Notar que la integral de probabilidad es una constante durante la evolución.

Finalmente decir que la solución formal de la ecuación de Schrödinger es

F(x,t) = e-itH F(x,0)
(4)
donde el operador eitH es unitario al ser H hermítico. Esta última propiedad nos será útil a la hora de hallar la discretización más adecuada.

Discretizaremos el espacio y el tiempo de forma que la función de ondas es:

F(x,t)®F(jh,ns) = fj,n
(5)
donde n = 0,1,2,...¥, j = 0,1,...,N, h es el espaciado del retículo y s es el espaciado temporal (ambos supuestamente pequeños. Supondremos que las condiciones de contorno para la función de ondas en j = 0 y j = N son las correspondientes a la existencia de un potencial infinito en esos puntos, esto es, la probabilidad de encontrar la partícula en dichos puntos es cero. Asi, F0,n = FN,n = 0, para todo tiempo. En este contexto, podemos definir de forma inmediata un primer algoritmo discreto a partir de la expresión (4). Según esa expresión, si t es muy pequeña, podemos aproximar el operador de evolución exp(-itH) por su desarrollo de Taylor quedándonos 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
F(j,n+1) = (1-isHD+O((sHD)2)F(j,n)
(6)
donde
HD fj = (- 1
h2
dj2+Vj)fj = - 1
h2
(fj+1-2fj+fj-1)+Vj
(7)
y Vj = V(jh). Este algoritmo tan simple y natural tiene un problema grave: 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 ondas, su normalización, åj h |Fj,n|2, irá variando con el tiempo lo que es incompatible con el carácter de probabilidad que se le tiene que dar. Asi pues, nuestro objetivo será encontrar un operador evolución similar al anterior pero que sea unitario. Esto lo conseguimos con el siguiente algoritmo:
Fj,n+1 = 1-isHD/2
1+isHD/2
Fj,n
(8)
Notar que además de ser unitario, este operador es exácto hasta orden (sHD)3. El único problema que nos resta con este algoritmo es el de diseñar la estrategia utilizar eficazmente el operador 1/(1+isHD). Para ello reescribimos la anterior ecuación como:
Fj,n+1 = é
ê
ë
  2
1+isHD/2
-1 ù
ú
û
Fj,n = cj,n-Fj,n
(9)
donde
  2
1+isHD/2
Fj,n = cj,n
(10)
o lo que es lo mismo, dada Fj,n, j = 0,..,N, cj,n es la solución de la ecuación
[ 1+isHD/2 ]cj,n = 2Fj,n
(11)
que en forma explícita puede escribirse como:
cj+1,n+ é
ê
ê
ê
ë
-2+ 2i
~
s
- ~
V
 

j 
ù
ú
ú
ú
û
cj,n+cj-1,n = 4i
~
s
Fj,n
(12)
donde [s\tilde] = s/h2 y [V\tilde]j = h2Vj. De esta forma el algoritmo básico de evolución queda bastante claro: dada una función de ondas en el instante n para toda posición j se resuelve el conjunto de ecuaciones (12) y se obtiene cj,n j = 0,..,N. Con esta solución vamos a la ecuación (9) y obtenemos las nuevas Fj,n+1 j = 0,...,N, volviendo a iterar el proceso. Lo único que nos queda por conocer es como resolver ecuaciones del tipo (12), osea, como invertir matrices tridiagonales.

En nuestro caso hemos de resolver un conjunto de ecuaciones del tipo:

Aj-cj-1+Aj0cj+Aj+cj+1 = bj     j = 1,¼,N-1
(13)
donde Aj- = 1, Aj0 = -2+2i/[s\tilde]-[V\tilde]j, Aj+ = 1 y bj = 4iFj,n/[s\tilde]. Las condiciones de contorno c0 = cN = 0 (notar que estas condiciones de contorno implican que en (9) F0 y FN son cero para todo tiempo.

Para resolver la recurrencia (13) suponemos que su solución es del tipo:

cj+1 = ajcj+bj    j = 0,¼,N-1
(14)
donde para garantizar que se cumple cN = 0 tomaremos que aN-1 = bN-1 = 0. Sustituyendo esta expresión en (13) obtenemos:
cj = - Aj-
Aj0+Aj+aj
cj-1+ bj-Aj+bj
Aj0+Aj+aj
 
(15)
Si identificamos las ecuaciones (14) y (15) podemos definir las recurrencias para los coeficientes a y b, asi
aj-1 = -Aj -gj
(16)
 
bj-1 = gj(bj-Aj+bj )
(17)
donde gj-1 = Aj0+Aj+aj. Estas ecuaciones nos dan la forma de obtener todas las a y b partiendo de j = N-1 y obteniendo en orden decreciente aj y bj j = N-2, ...,1. Una vez que tenemos todas las a y b, utilizamos la recurrencia (14) para hallar las cj en orden creciente en j.

Una vez ya sabemos cómo obtener c y con ella Fn+1 nos queda discutir qué función de onda inicial introducimos y qué potencial usamos. La función de onda inicial que vamos a usar es una onda plana con una amplitud gausiana, esto es:

F(x,0) = eik0xe-(x-x0)2/2s2
(18)
notar que con esta elección, la probabilidad de encontrar inicialmente a la partícula en un punto x es una gausiana centrada en x0 y de anchura s. El número de oscilaciones completas que la función de ondas tiene sobre la red depende de k0 , así k0Nh = 2pnciclos. En lugar de dar como parámetro inicial k0, darémos nciclos. Obviamente nciclos = 0,1,....,N, pero físicamente no tendría mucho sentido que una oscilación completa de la función de ondas se produjese entre dos puntos del retículo pues querría decir que no estamos haciendo lo suficientemente fina nuestra discretización. Así pues el parámetro nciclos lo restringiremos a valores 1,...,N/4. De esta forma, un ciclo, como mínimo tendrá cuatro puntos. La posición media inicial y la anchura de la gausiana las tomamos: x0 = Nh/4 y s = Nh/8 respectivamente. Por útimo, el potencial que usaremos será de anchura N/5, centrado en N/2 y altura proporcional a la energía de la función de ondas incidente: lk02 donde l = 0.3, por ejemplo.

En resumen, utilizando las constantes reescaladas, la función de ondas en el retículo se tomará:

Fj,0 = ei[k]0 j e-2(4j-N)2/N2
(19)
donde [k]0 = k0 h = 2pnciclos/N donde nciclos = 1,...,N/4. El potencial
  ~
V
 

j 
= Vj h2 = 0    si     j not in [2N/5,3N/5]
(20)
        = l
[k]
 
2
0 
    si     j Î [2N/5,3N/5]
(21)
Lo único que nos queda por fijar es el parámetro [s] = s/h2. Puesto que la energía es proporcional a k02 y nuestro operador dinámico discreto tiende a ser exacto en potencias de Hs, lo óptimo es elegir Hs < 1, esto es, k02s < 1, asi deducimos que [s] < 1/[k\tilde]02. En particular tomaremos [s] = 1/4[k]02. Notar que los parámetro que hemos de dar al sistema inicialmente son N, nciclos y l pues todos los demás se fijan a partir de ellos.

En resumen, el algoritmo para resolver la ecuación de Schrödinger unidimensional es:

                    © 1985 Javier de Lucas