UNA
SOLUCION A LA ECUACION
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
|
= |
é ê ê ê ë
|
- |
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
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
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
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.
- Resolver la ecuación de
Schrödinger unimensional para un potencial cuadrado.
Discretizaremos el espacio y el tiempo de forma que la función de ondas es:
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
o lo que es lo
mismo, dada Fj,n, j = 0,..,N, cj,n es la solución de la ecuación
que en forma
explícita puede escribirse como:
cj+1,n+ |
é ê ê ê ë
|
-2+ |
2i
|
- |
~ V
|
j
|
ù ú ú ú û
|
cj,n+cj-1,n = |
4i
|
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
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:
- (0) Dar los parámetros iniciales: N, nciclos y l. Generar [s], [k]0, [V]j y Fj,0 (incluyendo las condiciones de contorno
F0,0 = FN,0 = 0.
- (1) Calcular a's y b's utilizando la recurrencias (16) y (17)
respectivamente.
- (2) Calcular c's de (14).
- (3) Calcular Fj,n+1 de (9).
- (4) n=n+1, Ir a (2).
©
1985 Javier de Lucas