Método Variacional com Monte Carlo aplicado ao oscilador harmônico quântico

The Monte Carlo variational method applied to the quantum harmonic oscillator

Resumos

O método variacional com Monte Carlo é ilustrado em uma aplicação simples da mecânica quântica: o oscilador harmônico unidimensional. Pormenores tanto do método de Monte Carlo como do variacional são discutidos. Além disto, um simulador gráfico do problema está disponível no endereço <A HREF="http://www.ifi.unicamp.br/%7Emreis/oh/oh.html">http://www.ifi.unicamp.br/~mreis/oh/oh.html</A> da rede de computadores de tal forma a permitir ao leitor variar diversos parâmetros do cálculo e acompanhar nossas discussões dos resultados de forma a ganhar uma maior familiaridade com o material apresentado.

método variacional; Monte Carlo

The Monte Carlo variational method is illustrated in the simple quantum-mechanics application: the one-dimensional harmonic oscillator. Details of Monte Carlo and variational methods are explained. Besides, a graphic simulator for the problem is available at <A HREF="http://www.ifi.unicamp.br/%7Emreis/oh/oh.html">http://www.ifi.unicamp.br/~mreis/oh/oh.html</A> to allow the reader changes several calculations parameters and follows our explanation of results to acquire more familiarity with the text presented.

variational method; Monte Carlo

ARTIGOS GERAIS

Método Variacional com Monte Carlo aplicado ao oscilador harmônico quântico

The Monte Carlo variational method applied to the quantum harmonic oscillator

Marcelo A. dos Reis¹; S.A. Vitiello

Instituto de Física, Universidade Estadual de Campinas, Campinas, SP, Brazil

RESUMO

O método variacional com Monte Carlo é ilustrado em uma aplicação simples da mecânica quântica: o oscilador harmônico unidimensional. Pormenores tanto do método de Monte Carlo como do variacional são discutidos. Além disto, um simulador gráfico do problema está disponível no endereço http://www.ifi.unicamp.br/~mreis/oh/oh.html da rede de computadores de tal forma a permitir ao leitor variar diversos parâmetros do cálculo e acompanhar nossas discussões dos resultados de forma a ganhar uma maior familiaridade com o material apresentado.

Palavras-chave: método variacional, Monte Carlo.

ABSTRACT

The Monte Carlo variational method is illustrated in the simple quantum-mechanics application: the one-dimensional harmonic oscillator. Details of Monte Carlo and variational methods are explained. Besides, a graphic simulator for the problem is available at http://www.ifi.unicamp.br/~mreis/oh/oh.html to allow the reader changes several calculations parameters and follows our explanation of results to acquire more familiarity with the text presented.

Keywords: variational method, Monte Carlo.

1. Introdução

Esse trabalho tem por objetivo ilustrar um poderoso método numérico que pode ser implementado na resolução de problemas físicos. O método variacional com Monte Carlo é aplicado com sucesso a diversos problemas de muitos corpos quânticos. Entre eles, podemos mencionar os sistemas formados por átomos de hélio que a baixas temperaturas apresenta o fenômeno da superfluidez, a capacidade de sob certas condições escoar sem viscosidade. Neste artigo com a finalidade de ilustrar o método, vamos aplicá-lo ao importante problema da mecânica quântica: o oscilador harmônico unidimensional. Embora o problema do oscilador possua solução analítica conhecida e desse modo não se justifica qualquer tratamento numérico na busca de uma solução, o método será empregado para esse sistema e seus resultados confrontados com os valores esperados para ilustrar a aplicação do método de Monte Carlo na obtenção de resultados variacionais.

Inicialmente, faremos breve revisão teórica do oscilador harmônico quântico [1] e em seguida apresentamos o método de Monte Carlo e o algoritmo de Metropolis [2]. A partir dessas considerações é apresentada a implementação do método de Monte Carlo ao problema proposto. Na segunda parte desse artigo estão contidos os desenvolvimentos computacionais e na terceira os resultados são abordados. Além disto, o leitor poderá acessar no endereço http://www.ifi.unicamp.br/~mreis/oh/oh.html da rede de computadores um softaware gráfico que realiza esta simulação. Essa simulação foi desenvolvida tendo em vista sua utilização nos diferentes sistemas operacionais, mais especificamente foi construí­do um Applet-Java [3] que pode rodar em praticamente todos os navegadores atuais inclusive naqueles pertencentes à classe "software livre", o que significa uma grande portabilidade do programa desenvolvido. O usuário do programa poderá interagir com parâmetros do método e obter valores da energia e seu erro padrão assim como visualizar a função de onda tentativa escolhida.

2. Oscilador harmônico quântico unidimensional

Nesta seção faremos uma breve descrição teórica do sistema oscilador harmônico quântico unidimensional. Os resultados aqui apresentados serão confrontados com aqueles de nossas simulações numéricas.

O operador Hamiltoniano desse sistema é dado por:

onde P é o operador momento e X, posição da partícula de massa m oscilando com frequência angular w.

Uma vez que H é independente do tempo, para obter os níveis de energia desse sistema basta resolver a equação de autovalor:

Podemos escrever essa equação na representação das coordenadas |x > e obter a equação de Schrodinger independente do tempo:

onde E_n é a energia para um dado número quântico n e y_n é o autoestado correspondente a essa energia:

As autofunções do oscilador harmônico são dadas por:

e para o estado fundamental n = 0 temos:

3. Cálculo variacional

A seguir, apresentamos a idéia do método variacional e mostramos como é possível implementá-lo numericamente através do método de Monte Carlo.

3.1. O método variacional

Em situações onde não é possível determinarmos a solução exata da equação de Schrödinger, o método variacional é de grande utilidade. Ao construirmos uma função de onda tentativa Y_T, é importante que ela incorpore da melhor maneira possível os conhecimentos que temos do sistema.

A energia variacional do sistema pode ser escrita na representação das coordenadas como:

Logo,

É interessante notar que a função tentativa Y_T(x) não precisa estar normalizada, pois essa constante de normalização é cancelada algebricamente em (8).

Neste trabalho onde queremos ilustrar a aplicação da teoria variacional em uma implementação com o método de Monte Carlo, vamos imaginar que as autofunções do oscilador harmônico sejam desconhecidas.

Usaremos a seguinte função tentativa:

que localiza a partícula na região de menor energia potencial e que depende do parâmetro variacional a.

Se substituirmos a função tentativa (9) na expressão da energia variacional em (8), obtemos:

onde adotamos um sistema de unidades onde

²/m = mw² = 1. Determinamos então a energia variacional minimizando a expressão da Eq. (10) com respeito a a, isto é, impondo = 0. Obtemos o mínimo em a = , resultando em:

que é a energia do estado fundamental do oscilador harmônico. Não há nada de surpreendente nesse resultado já que o valor a = 1/2 corresponde àquele da função de onda do estado fundamental do problema. Tudo que foi feito até agora foi apenas supor que a função de onda do estado fundamental desse sistema seja proporcional a uma Gaussiana e calcular analiticamente a energia associada a essa função variacional.

3.2. O método variacional com Monte Carlo [4]

Uma vez construída a função tentativa , a integral expressa em (8) deve ser calculada. Aqui admitimos que esta integral não pode ser tratada analiticamente com a finalidade de ilustrarmos a aplicação do método variacional de Monte Carlo.

Reescrevemos a integral dada em (8) de outra forma:

onde poderíamos ter abandonado o módulo já que nossa função tentativa é real.

Levando em conta a densidade de probabilidade p(x) associada a função tentativa e definindo a energia local :

podemos reescrever

Utilizando o Hamiltoniano da Eq. (3) no sistema de unidades onde

²/m = mw² = 1 e a função tentativa (9), obtemos:

O método variacional com Monte Carlo consiste em amostrar a densidade de probabilidade p(x) de modo a calcular numericamente (13) da seguinte forma:

onde entendemos amostrar p(x) como construir uma sequência {x₁, x₂, ..., x_M} onde a frequência dos valores x é proporcional a p(x) (ver próxima seção).

Os resultados estão sujeitos a incertezas estatísticas, o que torna essencial a estimativa da variância:

O erro padrão é dado por:

3.3. Algoritmo de Metropolis

Para amostrar as M configurações do problema a partir da densidade de probabilidade p(x) usamos o algoritmo de Metropolis [2]. Uma observação aqui é notar que este algoritmo pode ser aplicado a qualquer função tentativa.

Um fluxograma dos passos básicos desse algoritmo é apresentado na Fig. 1 onde z e z₁ são números aleatórios no intervalo [0,1], D é um parâmetro que determina a amplitude do deslocamento da partícula e no qual os resultados são independentes. O número de elementos da amostra é M e p(x) é a densidade de probabilidade. Note que o fator de normalização da função tentativa não precisa ser determinado, já que na razão das probabilidades ele é cancelado.

3.4. Minimização das correlações produzidas pelo algoritmo

Fica claro que o próprio algoritmo da Fig. 1 traz correlações nas amostras. Devemos considerar esse fato e para eliminarmos as correlações utilizaremos o que chamamos de médias por blocos. Além disso, devemos também descartar algumas amostras iniciais para evitarmos qualquer tendência nos resultados devido a escolha arbitrária da posição inicial.

A média por blocos consiste em dividir os M elementos amostrados em N blocos e calcular valores médios de energia para M/N elementos do i-ésimo bloco (onde i = 1, 2,..., N) e tratar cada um destes como valores independentes para o cálculo da variância. O valor médio da energia é dado pela média dos valores e independe do número N, ao contrário da variância que depende diretamente do número de blocos escolhidos. O parâmetro N, portanto, permite que a correlação entre as amostras seja levada em conta e que as incertezas estatísticas sejam corretamente avaliadas.

4. Simulações

Em trabalhos onde a física computacional desempenha um papel relevante, a programação de computadores é em geral um aspecto importante de uma forma ou de outra. Assim, encorajamos vivamente o leitor a programar os algoritmos aqui descritos para ganhar familiaridade com os poderosos métodos de cálculo que eles representam. Infelizmente, mesmo com um programa que esteja livre de erros lógicos e, portanto funcionando corretamente, não é imediata sua utilização de forma ótima. Os cálculos dependem de parâmetros que a princípio não devem influenciar nos resultados. Na prática, entretanto, uma escolha desastrada dos mesmos pode levar a uma convergência muito lenta dos resultados ou até mesmo em casos extremos, a uma resposta errada. Com a finalidade de o leitor poder vivenciar um pouco a situação de utilizar o programa que aqui descrevemos, no endereço http://www.ifi.unicamp.br/~mreis/oh/oh.html oferecemos acesso a um applet-java que realiza a simulação aqui descrita.

O programa que implementa o cálculo variacional que descrevemos para um potencial harmônico foi inicialmente escrito em Fortran 90 e posteriormente, adaptado ao Java. Nesta seção discutiremos alguns aspectos do programa utilizado neste trabalho.

No Applet-Java disponível, além de poder alterar o valor do parâmetro variacional a que determina o valor da energia tentativa, o usuário pode alterar outros três parâmetros do cálculo: o número de amostras M que determina a precisão dos resultados, a amplitude do deslocamento das configurações D que é responsável pela convergência dos resultados e o número de blocos N associado a estimativa da variância. O valor destes últimos parâmetros a princípio não podem influenciar o resultado da energia variacional calculada.

O programa produz M amostras conforme descrito na seção 3.3 e calcula a energia e seu respectivo erro padrão. Além dessa grandeza, o programa também apresenta a fração de amostras propostas x_novo que é efetivamente incorporada à sequência de amostras x_i utilizada no cálculo da energia. O valor deste parâmetro (D) é ajustado de tal forma a permitir uma aceitação em torno de 50%. Outro resultado apresentado é o desvio percentual que a energia calculada pelo método de Monte Carlo possui da energia calculada analiticamente (10). Além da visualização da função de onda, esse desvio ajuda o usuário a perceber mais facilmente como os parâmetros fornecidos ao programa estão influenciando os resultados.

4.1. O parâmetro delta (D)

Para analisar o comportamento da fração de passos aceitos em função do deslocamento D, vamos fixar a = 0.49 de modo que para todas as simulações que realizaremos tenhamos a mesma função de onda tentativa e mesma energia dentro das incertezas de natureza estatística. Os resultados estão ilustrados na Fig. 2. Como esperado, quanto menor o valor de delta maior aceitação se tem. Isso ocorre porque quando permitimos somente deslocamentos pequenos, a probabilidade da candidata à nova configuração tem um valor próximo ao da probabilidade do elemento anterior que já pertence à amostra. Como a condição de aceitação verifica se a razão das probabilidades do candidato à amostra e o seu elemento anterior é maior que um número aleatório em um intervalo de [0,1] é muito provável que haja aceitação dessa nova configuração já que a razão é aproximadamente a unidade. O contrário ocorre para delta muito grande, a razão das probabilidades pode ser muito diferente uma da outra, pois com um delta grande permitimos grandes deslocamentos e assim comparando com o número aleatório no intervalo [0,1] pode ocorrer muita rejeição, como de fato ocorre. Uma regra empírica mostra que um bom valor para o delta é manter uma faixa de aceitação em torno de 50%, o que corresponde a um D = 4 em nosso caso.

Na Fig. 2 apresentamos a fração de aceitação para 10² e 10⁴ passos Monte Carlo, podemos ver também que para um maior número de passos, mais suave é a variação da fração de aceitação em função do delta. Isso porque temos um maior número de amostras e, portanto menores flutuações estatísticas.

4.2. O número de blocos (N)

Utilizando nosso programa escrito em Fortran, fizemos diversos testes para determinar o número de blocos adequado para estimarmos um erro padrão confiável. Após essa análise disponibilizamos alguns valores de blocos para o usuário interagir.

A metodologia para essa análise foi a seguinte: primeiro estipulamos um valor para a próximo de 0.5 e um determinado número de passos. Na Fig. 3 apresentamos o erro padrão em função do número de blocos. Para obter os resultados da Fig. 3(a) escolhemos um valor de D inferior a 4 (o valor adequado para esse parâmetro) para aumentar a aceitação e por em evidência o comportamento do erro padrão. Apresentamos na Fig. 3(b) resultados com o parâmetro D = 4. A partir de uma mesma amostragem em cada caso, as energias são calculadas com diferentes números de blocos e a energia média estimada será a mesma independentemente do número de blocos enquanto que o comportamento do erro padrão depende deste número.

O comportamento de DE em função do número de blocos varia com o a valor inicial dos parâmetros. As Figs. 3(a) e 3(b) ilustram tal situação. O número de bloco adequado para cada conjunto de simulações é aquele onde o erro padrão é máximo. Nos conjuntos de parâmetros ilustrados na Fig. 3, uma boa estimativa do erro seria usar um número de blocos entre 10 e 100, portanto, disponibilizamos para o usuário as quantidades 1, 10, 50 e 100 números de blocos.

4.3. O número de passos (M)

O objetivo aqui é ver o comportamento do cálculo da energia para diferentes números de amostras M, mantendo fixos os parâmetros alfa, delta e o número de blocos. Na Fig. 4 apresentamos a energia estimada em função do número de amostras. Verificamos que quanto maior o número de amostras mais "estável" é o valor da energia já que a variância diminui com o aumento de amostras proporcionalmente a .

5. Resultados

A energia variacional do problema é obtida minimizando o valor da energia com respeito ao parâmetro a da Eq. (9). A maneira mais simples de realizar esta tarefa é através da comparação de resultados de simulações independentes para um intervalo de valores do parâmetro variacional a.

Adotando os valores mais adequados dos parâmetros de cálculo D e N que acabamos de discutir, realizamos simulações onde consideramos M igual a 5000 configurações. Os resultados que obtivemos estão ilustrados na Fig. 5. A partir desta figura podemos fazer diversas observações. À medida que nos aproximamos de a = 0.50 os erros diminuem. Em particular, para a = 0.50 a variância do resultado é nula. Para este valor do parâmetro a a função da Eq. (9) é realmente a própria solução para o estado fundamental do oscilador harmônico. A energia de todos os pontos amostrados coincide com o valor desta grandeza para este estado. Assim vemos de fato que o método aqui discutido não introduz aproximações além daqueles do método variacional, assim, os resultados são apenas afetados por incertezas estatísticas.

O exame dos resultados apresentados na Fig. 5 indica que o valor de a que minimiza a energia variacional deve estar no intervalo 0.47 < a < 0.53. Havendo necessidade de uma maior precisão para o valor do parâmetro a e sua correspondente energia, é suficiente prolongar as simulações apenas para a nesse intervalo.

O observador atento poderia ainda perguntar o que estaria ocorrendo com a energia para outros valores de a como, por exemplo, para a = 0.57. Se examinarmos seu erro e o que ocorre em sua vizinhança poderíamos supor que realmente houve uma flutuação estatística e que este valor não corresponde a um mínimo de energia. É importante, entretanto, observar que não há como garantir esta situação a priori, sendo, portanto, em última análise, necessário também prolongar a simulação para este valor de a. Esta é em verdade uma limitação do método. Não é possível determinar a priori quantas amostras são necessárias para atingirmos uma dada precisão. É importante também notar que o método de Monte Carlo, através de procedimentos mais elaborados como amostragem de importância ou re-pesagem das amostras possui recursos que podem minorar ou até mesmo eliminar estas dificuldades [5].

Em resumo, com os resultados aqui apresentados podemos dizer que o mínimo da energia variacional para o estado fundamental do problema é aproximadamente 0.5 em unidades tais que w = 1. Este valor corresponde a um mínimo para a = 0.5, centro simétrico do intervalo dos valores mais prováveis de a. A unidade de a é comprimento a menos dois. Convidamos vivamente o leitor a utilizar o programa que realiza estas simulações e que está disponível na rede para verificar e refinar os resultados aqui apresentados. Os mesmos estão em excelente acordo com aqueles esperados, Eq. (11) e discussão seguinte.

6. Conclusão

Na pesquisa dos sistemas de muitos corpos quânticos onde métodos variacionais são extensamente empregados, as funções tentativas utilizadas quase sempre levam a integrais sem soluções analíticas e, portanto, métodos numéricos devem ser empregados nos cálculos das quantidades de interesse. É neste contexto que o método variacional com Monte Carlo, aqui ilustrado, torna-se bastante importante.

Neste artigo aplicamos o método variacional com Monte Carlo ao oscilador harmônico quântico, um sistema físico bastante simples que dispensa a aplicação de métodos numéricos, mas que permite a comparação entre os nossos resultados e aqueles exatos. Além disso, uma aplicação computacional gráfica de fácil interatividade que realiza as simulações aqui descritas está disponível na rede de computadores e pode ser utilizada para que uma maior familiaridade com o método possa ser adquirida.

Recebido em 16/11/2005; Aceito em 23/1/2006

• [1] C. Cohen-Tannoudji, B. Dui and F. Laloe, Quantum Mechanics (Wiley-Interscience Publication, New York, 1977).
• [2] N. Metropolis A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys. 21, 1087 (1953).
• [3] H.M. Deitel e P.J. Deitel, Java, Como Programar(Bookman, Porto Alegre, 2003), traduçăo de Carlo Arthur Lang Lisboa, 4Ş ed.
• [4] R. Canabrava e S.A. Vitiello, Rev. Bras. Ens. Fis. 25, 35 (2003).
• [5] M.H. Kalos and P.A. Whitlock, Monte Carlo Methods: Basics (Wiley-Interscience Publication, New York, 1986), v. 1.

Datas de Publicação

• Publicação nesta coleção
  19 Jul 2006
• Data do Fascículo
  2006

Histórico

• Aceito
  23 Jan 2006
• Recebido
  16 Nov 2005