Soluci´ on on Num´ erica erica al Oscilador Oscilado r Arm´ onico onico Cu´ antico Unidimensional antico Neri Ne ria a P´erez er ez Jos´ Jos´e Ange An gell Centro de Investigaci´ on on en Matem´ aticas aticas
15 de diciembre de 2014
Resumen
1 ψ
En este trabajo se presenta la implementaci´on on realizada para la soluci´on on al problema del oscilador arm´onico onico cu´antianti2 d co unidimensional ψ (x) + V + V (x (x)ψ (x) = E ψ (x) con 2m dx2 m V (x (x) = ω 2 x2 para un electr´on on con ω = 1 10 15 rad/s. rad/s. 2 Dicha ecuaci´on on se discretiz´o mediante diferencias finitas. Se comparan resultados de la implementaci´on on y los de un software comercial. Se muestran las primeras tres funciones de onda obtenidas num´ ericamente ericamente y se muestra que la diferencia entre estas y las anal´ıticas ıticas es del orden de 1 10 6 .
−
√ −
2
2
2
(5)
d2 ˆ ψ (x) = Eψ + V (x (x) ψ (x) = H E ψ (x) 2m dx2
(6)
−
−
La ecuaci´on on de Schr¨oedinger oedinger dependiente del tiempo es la ecuaci´ on on fundamental de la f´ısica que describe el comportacomp ortamiento miento de los cuerpos peque˜ nos. nos. Tambi´en en se le denota como la ecuaci´on on de onda de Schr¨oedinger oedinger y es una ecuaci´on on diferencial parcial que describe como la funci´on on de onda de un sistema sistem a f´ısico ısico evoluciona en el tiempo. tiemp o. La ecuaci´ ecuaci´ on dependiente del tiempo est´a dada como on ˆΨ − 2 m ∂ ∂xΨ + V (x, (x, t)Ψ(x, )Ψ(x, t) = H
Donde en la segunda igualdad de la ecuaci´on on (6) se redefi2 d ˆ ni´ o al operador hamiltoniano como H = + V (x (x). 2m dx2 Y a las soluciones de dicha ecuaci´on on se les denomina funciones de onda. Como se puede apreciar la ecuaci´on on (6) representa un problema de valores y vectores propios, donde los vectores propios del operador hamiltoniano son las soluciones a la ecuaci´ on on independiente del tiempo. Conviene hacer ´enfasis enfasis en que todos los valores valores propios de ´esta esta ecuaci´ on on son reales . Esto es una propiedad necesaria de cualquier teor´ıa ıa v´alida alida ya que en los experimentos no se miden, por ejemplo, distancias de 1 + 6i 6i metros. Es importante imp ortante mencionar que a pesar que la energ´ energ´ıa potencial no depende del tiempo, hay pocos casos en los que existe soluci´ on on anal´ anal´ıtica y en las aplicaciones se busca el calcular los valores propios m´as as peque˜ nos ya que en las funciones de nos onda asociadas a estos es m´as as plausible pl ausible hallar a la part´ıcula, ıcula, esto debido al principio de m´ınima acci´ on. on. En este trabajo m 2 2 se abordar´a el caso V (x (x) = ω x , que se denomina osci2 lador arm´onico onico cu´antico. antico. Las funciones de onda asociadas a este problema son
Intr Introdu oducc cci´ i´ on on
∂ Ψ i = ∂t
d2 ψ + V (x (x)ψ (x) = E = E 2m dx2
iEt La soluci´on on a la ecuaci´on o n (4) es T ( T (t) = e . Mientras que la ecuaci´on on (5) se reescribe como
×
1.
2
−
−
×
−
(1)
donde i = 1, es la consta constant ntee reduci reducida da de Planc Planck, k, ˆ V (x, (x, t) representa la energ´ energ´ıa potencial del sistema y H se se le denomina operador hamiltoniano, que representa la energ´ energ´ıa total del sistema. En el caso en que V que V (x, (x, t) = V (x (x), se puede escribir Ψ(x, Ψ( x, t) = ψ (x)T ( T (t), por lo que al substituir en la ecuaci´on on (1) se tiene
2 1 ξ 2 dT d2 ψ = T ( T (t) 2 + V (x (x)ψ (x)T ( T (t) (2) mω 4 1 dt 2m dx ψn (x) = H n (ξ ) e 2 , π 2n n! Y al dividir esta ´ultima ultima ecuaci´on on por Ψ(x, Ψ(x, t) = ψ( ψ (x)T ( T (t) se tiene mω 2 2 i dT 1 d ψ ξ = = x, n = 0, 1,... (7) = + V (x (x)ψ (x) (3) 2 T dt ψ 2m dx Dado que el lado derecho de la ecuaci´on (3) depende s´ olo olo de y los valores propios la posici´on on y el lado izquierdo del tiempo, entonces deben ser 1 E n = n + (8) ω constantes para que esto ocurra, esto implica que la ecuaci´on on 2 (3) se convierte en las ecauciones Donde H n (ξ ) son los polinomios asociados de Hermite. i dT = E (4) La forma de onda de las primeras cuatro funciones de onda T dt
i ψ (x)
−
−
√
1
−
del oscilador arm´onico, as´ı como sus cuadrados se muestran en la figura 1.
Figura 2: Cuadrado de una funci´on de onda. Ser´a relativamente m´ as f´ acil de hallar cerca de la posici´on A que cerca de la posici´on B. El ´area sombreada representa la probabilidad de hallar la part´ıcula en el rango dx .
2.
Soluci´ on Num´ erica
En este traba jo se discretiz´o el intervalo 7 7 6 10 , 6 10 en N sub intervalos Por otro lado, se discretiz´ o la ecuaci´o n (6) mediante diferencias centrales, la matriz resultante es una matriz bandada tridiagonal de tama˜ no N N , y el siguiente es el stencil que define a la matriz:
− ×
−
×
−
×
Figura 1: (a) Primeras cuatro funciones de onda del oscilador arm´onico; (b) Cuadrado de sus amplitudes
1.1.
2m∆2
[ 1, 2, 1] + [0, V (xi ), 0]
−
−
(9)
donde ∆ es el incremento en la malla y x i son los puntos en la malla. y los u ´ nicos renglones distintos son el primero y el ´ultimo, que ser´ an: Interpretaci´ on Estad´ıstica de la Solu-
ci´ on de la Funci´ on de Onda
2m∆2 Una pregunta natural que surge es: una vez que se ha obtenido la o las soluciones a las ecuaciones (1) o (6) ¿ como interpretarlas? Ya que una part´ıcula est´a instintivamente localizada en un punto en el espacio, mientras que una onda, como su nombre lo indica, est´a dispersa en el espacio. Entonces ¿como es posible que la funci´on de onda describa el estado de una part´ıcula? La respuesta se obtiene mediante la interpretaci´on estad´ıstica de Born de la funci´on de onda, que dice que Ψ(x, t) 2 = Ψ(x, t) Ψ(x, t) representa la probabilidad de hallar a la part´ıcula en el punto x al tiempo t1 , o m´as precisamente:
|
|
[2, 1] + [V (x0 ), 0] ,
−
[ 1, 2] + [0, V (xN 1 )] 2m∆2 La soluci´ o n al oscilador arm´ o nico cu´antico ser´an los vectores propios asociados a esta matriz y las soluciones de inter´es son las asociadas a los valores propios m´as peque˜ nos de dicha matriz. Donde el vector propio asociado al valor propio m´ as peque˜ n o se le denomina estado base y a los vectores propios subsecuentes se les denomina estados excitados.
∗
2.1.
−
−
Implementaci´ on
2
|Ψ(x, t)| dx = Probabilidad de hallar a la part´ıcula
La implementaci´ on de la soluci´on sigue el siguiente pseudo c´ odigo:
en la posici´on entre x y x + dx al tiempo t;
Proveer un archivo que contenga la matriz formada por 2 en la diagonal principal ,-1 en las diagonales adyacentes y 0 en los otras entradas.
Esta interpretaci´o n es la m´as aceptada. En la figura 2 se muestra el cuadrado de una funci´on de onda y se ilustra la interpretaci´ on estad´ıstica.
Leer el archivo provisto, determinar si el archivo define una matriz sim´etrica y almacenar la informaci´on de los elementos no nulos en un objeto denominado matriz rala.
1 La funci´ o n de onda Ψ(x, t) es compleja, pero |Ψ|2 = Ψ Ψ, donde Ψ denota al complejo conjugado de Ψ, es una funci´o n real y no negativa, como una funci´ on de probabilidad debe ser. ∗
∗
2
Se define el intervalo sobre el que se desea trabajar y se particiona de acuerdo al tama˜no de la matriz.
Formar la matriz hamiltoniana que ser´a la suma del stencil de la segunda derivada con el potencial evaluado en los puntos de la malla como se enuncia en la ecuaci´on (9).
Calcular los primeros valores propios de la matriz mediante deflaci´on2 .
Figura 3: Diferencia entre valores propios num´ericos y anal´ıticos. Azul, implementaci´on. Rojo, software comercial, N=200.
Almacenar los valores y vectores propios para su posterior an´ alisis.
Se reportan los primeros 10 valores propios en entrada est´andar y se comparan con los valores anal´ıticos.
3.
Resultados
Se realiz´o la partici´on del intervalo 7 7 6 10 , 6 10 en 200 y 1000 partes iguales y se calcularon los primeros 200 valores propios m´a s peque˜ nos de las matrices asociadas con la implementaci´ on realizada y con un software comercial.
− ×
−
×
−
En la figura 3 se muestra la diferencia entre los valores propios calculados mediante la implementaci´on y los anal´ıticos en azul; en rojo la diferencia entre los valores propios obtenidos por el software comercial para N=200. Mientras que en la figura 4 se muestra un acercamiento de la misma. Se aprecia que en ambos casos se predicen correctamente los primeros 4 valores propios. Finalmente, se realiza una comparaci´ o n an´ aloga con una matriz de tama˜ no 1000 1000 en la figura 5 y en la figura 6 se presenta un acercamiento a esta. Se observa que el comportamiento de ambas curvas es similar y se predicen correctamente los primeros 12 valores propios.
×
2
Figura 4: Acercamiento de la figura 3 .
Ver el apendice A
3
Figura 5: Diferencia entre valores propios num´ericos y anal´ıticos. Azul, implementaci´on. Rojo, software comercial, N=1000. Figura 7: Diferencia entre Ψ 0 (x) anal´ıtica y Ψ0 (x) numerica .
Figura 6: Acercamiento de la figura 5 .
Finalmente se estudiaron los primeros tres vectores propios m´ as peque˜ n os de la matriz de tama˜ n o 1000. En las figuras 7 - 9 se muestra la diferencia entre los valores anal´ıticos de las primeras tres funciones de onda y los obtenidos. Mientras que en la figura 10 se muestran las primeras tres funciones de onda. Por otro lado, en las figuras 11 - 13 se muestran las diferencias entre Ψi (x) Ψi (x) num´erica y Ψi (x) Ψi (x) anal´ıtica y en la figura 14 se muestran Ψ 0 (x) Ψ0 (x), Ψ1 (x) Ψ1 (x), Ψ2 (x) Ψ2 (x), que son las densidades de probabilidad del estado base, el primer estado excitado y el segundo estado Figura 8: Diferencia entre Ψ 1 (x) anal´ıtica y Ψ1 (x) numerica excitado, respectivamente. . ∗
∗
∗
∗
∗
4
Figura 9: Diferencia entre Ψ 2 (x) anal´ıtica y Ψ 2 (x) numerica .
Figura 11: Diferencia entre Ψ0 (x) Ψ0 (x) anal´ıtica y Ψ0 (x) Ψ0 (x) numerica .
Figura 10: Primeras tres funciones de onda.
Figura 12: Diferencia entre Ψ1 (x) Ψ1 (x) anal´ıtica y Ψ1 (x) Ψ1 (x) numerica .
∗
∗
∗
∗
5
Deflaci´ on para quitar los valores propios ya calculados anteriormente. Con respecto al primer cuello de botella se pueden realizar al menos dos mejoras: Dado que el problema a estudiar es una matriz bandada, que el programa reciba como par´ametros el stencil que define a la misma y el tama˜no de la matriz a estudiar. Que el programa lea la matriz en el formato definido en la direcci´on electr´onica http://math.nist.gov/MatrixMarket/formats.html
Figura 13: Diferencia entre Ψ2 (x) Ψ2 (x) anal´ıtica y Ψ2 (x) Ψ2 (x) numerica . ∗
∗
Mientras que en el caso del segundo cuello de botella se puede paralelizar la secci´ o n de c´odigo que realiza la deflaci´ on y realizar la comparaci´on de tiempos de ejecuci´on. Por otro lado la implementaci´on realizada se puede extender al estudio de otros problemas, por ejemplo el caso de la denominada part´ıcula en una ca ja de potencial infinito, con V (x) de la forma V (x) =
0
∞
si 0 x L en otro caso .
≤ ≤
(10)
Cabe mencionar que se realiz´o la prueba con este potencial, pero los resultados, en el caso de los valores propios no son f´ aciles de interpretar, ya que se encuentran en funci´on de la constante L, sin embargo los obtenidos num´ericamente son independientes de este par´ametro, ya que siempre se calcular´ an los valores propios de la matriz con stencil
Figura 14: Densidades de probabilidad del estado base y los primeros dos estados excitados.
4.
− [−1, 2, 1]. 2m∆ 2
5.
Discusi´ on
Conclusiones
Se implement´ o la soluci´on a la ecuaci´on de Schr¨oedinger De las figuras mostradas en la secci´on de resultados se independiente del tiempo con potencial de oscilador arm´onim aprecia que en los primeros tres vectores propios de la ma- co de la forma V (x) = ω 2 x2 . La ecuaci´on fue discretizada 2 triz estudiada se tiene una discrepancia m´axima de hasta mediante diferencias centrales. 6 10 6 con los valores anal´ıticos, mientras que en los valo- Las soluciones de la ecuaci´on son los vectores prores propios de la matriz, se tiene una concordancia con los pios asociados a la matriz hamiltoniana cuyo stencil es 2 primeros 12 valores propios anal´ıticos para una discretiza [ 1, 2, 1] + [0, V (xi ), 0]. ci´on de tama˜ no 1000 y 4 para una de tama˜no 200. 2m∆2 Se observa que la curva de los valores propios es similar a Los resultados muestran que se predicen de manera correcta nos con una matriz la del software comercial, aunque en ambos casos se tiene los primeros 3 valores propios m´as peque˜ no 200 200, y al crecer el n´umero de valor propio un comportamiento no lineal despu´es de los valores propios de tama˜ erico incrementa de una forma no calculados de forma correcta. Se cree que esto se debe al calculado, su valor num´ stencil que define a la matriz, una modificaci´on que valdr´ıa lineal. Mientras que con la matriz de tama˜no 1000 1000 se la pena realizar es la comparaci´on entre este stencil y uno predicen correctamente los 12 primeros valores propios. En ambos casos se observa que el comportamiento de la implede cinco puntos. Por otro lado se observ´o que durante la ejecuci´on del pro- mentaci´on y el software comercial es similar. Con respecto a los vectores propios se muestran resultados para los tres asograma hay esencialmente dos cuellos de botella: ciados a los tres valores propios m´as peque˜ nos. Se observa Lectura del archivo de la matriz que representa el pro- que la diferencia m´axima entre los obtenidos num´ericamente blema. y los anal´ıticos difieren a los m´as en 6 10 6 .
×
−
−
−
×
×
×
6
−
A.
M´ etodo de la Potencia Inversa
Este m´ etodo consiste en dar una entrada x0 tal que x0 = 1 y posteriormente realizar la iteracion xk+1 A 1 xk , de tal forma que en cada iteracion, el vector x k+1 A 1 xk es normalizado. Se puede mostrar que si la subsucesi´on formada por xk converge, entonces la subsucesi´on µk definida como µ k = < xT , xk > converger´a a el valor propio m´as peque˜ no, < xT , Axk > k+1 denotado como λ 1 .
←
−
−
{ }
{ }
A.1.
Pseudo C´ odigo
Figura 15: Pseudo codigo del gradiente conjugado precondicionado [4].
El pseudo c´odigo para realizar el m´etodo de la potencia inversa es el siguiente:
A.2.
Definir x0 tal que x0 y λ = 0, λ1 = 1, que ser´an los primeros valores aproximados de λ 1 .
Determinaci´ on de los m Valores Propios m´ as Peque˜ nos de una Matriz
Se busca el valor propio λ m al conocer los anteriores mas peque˜ nos λ m 1 con 1 m 1.
≤ −
−
A.2.1.
Para cada i hasta i max .
Pseudo C´ odigo
El pseudo c´odigo para realizar el calculo de los m valores propios mas peque˜ nos de A es el que se enuncia:
• Realizar la asignaci´on x ← A
1
−
1
Para cada i menor que m
x0
• Definir una sugerencia x que ser´a la aproximaci´on 0
x1 . x1
a φ i .
• Realizar la asignaci´on x ← 1
•
• Definir λ = 1000,
λ1 = 400 que ser´an las aproxi-
maciones a λ i .
1 Definir λ = . < x 1 , x0 >
i
• Asignar x ∗ ← x − < x , φ > φ • Realizarx ∗el m´etodo de la potencia inversa con input x ∗ . • Guardar λ y φ para accesar a ellos en la siguiente 0
0
0
k
k
k =1
◦ Si |λ |λ− λ| | < ε, donde ε es una tolerancia de0
1
0
finida por usuario, salir y reportar λ y
0
x1 . x1
0
iteraci´ on.
◦ En caso contrario definir λ = λ y x = xx .
i
1
1
0
1
B.
M´ etodos Implementados
Como la matriz del problema a estudiar es bandada En la implementacion realizada no se calcula la inversa de (rala), entre m´ a s fino sea el particionado del intervalo 1 la matriz rala A , sino que se resuelve el sistema mediante a estudiar, m´ as entradas nulas habr´ a en la matriz que el m´etodo del gradiente conjugado precondicionado, cuyo representar´ a a este problema, por este motivo se decidi´o impseudo c´odigo se muestra en la figura 15 plementar el uso de la plantilla de clase std :: map, en −
7
donde se almacenar´an las entradas no nulas de dicha matriz.
B.3.
B.1.
Una vez definido el tipo Matriz se defini´o la clase M rala que contendr´a los elementos no nulos de la matriz a estudiar con las siguientes variables protegidas:
Plantilla std::map
Clase M rala
bool isSymmetric: Esta variable identificar´a cuando la Los maps son contenedores asociativos que almacenan elematriz rala sea sim´etrica. mentos formados por una combinaci´o n de una clave y un valor mapeado, siguiendo un orden espec´ıfico. int m,n: Con estas variables se identifica a las dimenEn un map, los valores clave son generalmente usados para siones de la matriz, en el caso del problema abordado ordenar e identificar u ´ nivocamente a los elementos, mientras son cuadradas, por lo que basta con un ´ındice, pero se que los valores mapeados almacenan el contenido asociado usan los dos para futuras implementaciones. a esta clave. Los tipos de clave y el valor mapeado pueden Los siguientes son los m´etodos p´ ublicos de la clase: ser diferentes, pero se agrupan juntos en un tipo. Internamente, los elementos en un map son siempre ordenaM rala(int m, int n, bool is): Constructor de la clase dos mediante su clave siguiendo un criterio de ordenamiento. matriz, define las dimensiones y si es sim´etrica. Los valores mapeados en un map pueden ser accesados directamente por su clave correspondiente usando el operador M rala(): Constructor de la clase matriz, define bracket []. Con lo anterior en mente se defini´o un tipo de par´ ametros por default. matriz, a saber: M rala(M rala &mat): Constructor de copia de la clase, copia las variables de la matriz mat a la instanciada. typedf map < int, map < int, double > > M atriz double obtiene(int i, int j): M´ etodo que devuelve la entrada (i,j) de la matriz en base si es sim´etrica o no.
En los contenedores del nuevo tipo se almacenar´an los elementos no nulos de la matriz. La ventaja que presentan sobre los arreglos es que si se accede a ´ındices no definidos, el map por default informa que su contenido es cero, un ejemplo de esto se muestra a continuaci´on:
int size(): M´etodo que devuelve el n´ umero de renglones de la matriz rala; void asigna(int i,int j,double dato): M´ etodo que asigna el dato a la entrada (i,j) de la matriz.
Matriz My mat; My mat[0][1] = 5; cout << My mat[0][1] << endl; cout << My mat[1][1] << endl;
void modifica(int m,int n, bool is): M´ etodo que modifica los par´ametros protegidos de la clase.
bool estado(): M´etodo que devuelve la variable proteLa salida de este programa ser´a 5 en la primera linea y 0 gida isSymetric. en la siguiente. La ventaja que presenta lo anterior es que basta con almavoid imprime(); M´ etodo que imprime la matriz rala a cenar en el map los elementos no nulos de la matriz que se consola. desea estudiar, pero esta ventaja es su debilidad: se pueden hacer operaciones que no son v´alidas, por ejemplo multipliB.4. Operaciones Definidas para los Elecar una matriz de 5 5 con un vector de 6 1, por lo que mentos de la Clase M rala siempre que se defina una funci´on que involucre elementos de tipo Matriz se debe verificar si dicha operaci´on es v´alida Se definieron las siguientes operaciones entre elementos (esto se debe realizar con operaciones matriz-vector) de la clase M rala, elementos de la clase vector y variables doubles:
×
B.2.
×
Plantilla std::vector
vector operator*(M rala &A, vector &b): Sobrecarga del operador para realizar el producto matriz rala con vector, el resultado ser´a un vector. En la definici´on de esta sobrecarga se verifica primero que esta operaci´on sea v´alida mediante las dimensiones del vector y la matriz.
∗
Los vectores son una secuencia de contenedores que representan arreglos que pueden cambiar de tama˜no. De la misma forma que los arreglos, en los vectores se usa un almacenamiento contiguo para sus elementos, pero a diferencia de los arreglos, su tama˜ no puede cambiar din´ amicamente, manejando la informaci´ on contenida en ellos.
double operator*(vector &a,vector &b): Sobrecarga del operador para realizar el producto punto de dos vectores.
∗
8
vector operator/(vector &a,vector &b): Sobrecarga del operador / para realizar el cociente de vectores punto a punto. void imprime(vector &b): Imprime el contenido del vector b. vector construye(int n): Toma un elemento de la clase vector y lo asigna de tama˜ no n, inicializa sus valores a cero. vector diag(M rala &A): Obtiene la diagonal de la matriz rala A. void multiplica mrala(double a,M rala &A): Multiplica a los elementos de la matriz rala con el escalar a. void aniade diagonal(M rala &A,vector &b): Suma los elementos del vector b a los elementos de la diagonal de la matriz rala.
Referencias [1] Griffiths, D. J., and Harris, E. G. (1995). Introduction to quantum mechanics (Vol. 2). New Jersey: Prentice Hall. [2] De La Pe˜ na, L. (1996). Introducci´ o n a la Mec´anica Cu´antica. Fondo de Cultura econ´omica. [3] Asmar, N. H. (2005). Partial differential equations with Fourier series and boundary value problems. Pearson Prentice Hall. [4] Jorge, Nocedal, and J. Wright Stephen. Numerical optimization. Springerverlang, USA (1999).
9
