Equação de Schrodinger com condições de contorno periódicas

9

Eu tenho algumas perguntas sobre o seguinte:

Estou tentando resolver a equação de Schrodinger em 1D usando a discretização de crank nicolson seguida pela inversão da matriz tridiagonal resultante. Meu problema agora evoluiu para um problema com condições de contorno periódicas e, portanto, modifiquei meu código para usar o algoritmo Sherman Morrison.

Suponha que vmeu RHS esteja em cada etapa do tempo em que desejo inverter a matriz tridiagonal. O tamanho de vé o número de pontos de grade que tenho no espaço. Quando eu defino v[0]e v[-1]em termos um do outro, conforme é necessário em minha situação periódica, minha equação explode. Não sei dizer por que isso está acontecendo. Eu estou usando python2.7 e resolve_banded do scipy para resolver a equação.

Isso me leva à minha segunda pergunta: usei python porque é a linguagem que melhor conheço, mas acho bastante lenta (mesmo com as otimizações oferecidas por numpy e scipy). Eu tentei usar C ++ porque estou razoavelmente familiarizado com ele. Eu pensei em usar o GSL que seria otimizado para o BLAS, mas não encontrei documentação para criar vetores complexos ou resolver a matriz tridiagonal com esses vetores valiosos complexos.

Gostaria de objetos no meu programa, pois acho que seria a maneira mais fácil de generalizar mais tarde para incluir o acoplamento entre as funções de onda, por isso estou aderindo a uma linguagem orientada a objetos.

Eu poderia tentar escrever o solucionador de matriz tridiagonal manualmente, mas tive problemas quando o fiz em python. À medida que evoluía ao longo do tempo, com etapas cada vez mais finas, o erro se acumulava e me dava bobagem. Tendo isso em mente, decidi usar os métodos internos.

Qualquer conselho é muito apreciado.

EDIT: Aqui está o trecho de código relevante. A notação é emprestada da página da Wikipedia na equação da matriz tridiagonal (TDM). v é o RHS do algoritmo crank nicolson a cada etapa do tempo. Os vetores a, bec são as diagonais do TDM. O algoritmo corrigido para o caso periódico é do CFD Wiki . Eu fiz um pouco de renomeação. O que eles chamaram de u, v Eu chamei de U, V (maiúsculo). Chamei q de complemento, y a solução temporária e a solução real self.currentState. A atribuição de v [0] e v [-1] é o que está causando o problema aqui e, portanto, foi comentada. Você pode ignorar os fatores de gama. São fatores não lineares usados ​​para modelar os condensados ​​de Bose Einstein.

for T in np.arange(self.timeArraySize):
        for i in np.arange(0,self.spaceArraySize-1):
            v[i] = Y*self.currentState[i+1] + (1-2*Y)*self.currentState[i] + Y*self.currentState[i-1] - 1j*0.5*self.timeStep*potential[i]*self.currentState[i] - self.gamma*1j*0.5*self.timeStep*(abs(self.currentState[i])**2)*self.currentState[i]
            b[i] = 1+2*Y + 1j*0.5*self.timeStep*potential[i] + self.gamma*self.timeStep*1j*0.5*(abs(self.currentState[i])**2)

        #v[0] = Y*self.currentState[1] + (1-2*Y)*self.currentState[0] + Y*self.currentState[-1] - 1j*0.5*self.timeStep*potential[0]*self.currentState[0]# - self.gamma*1j*0.5*self.timeStep*(abs(self.currentState[0])**2)*self.currentState[0]
        #v[-1] = Y*self.currentState[0] + (1-2*Y)*self.currentState[-1] + Y*self.currentState[-2] - 1j*0.5*self.timeStep*potential[-1]*self.currentState[-1]# - self.gamma*1j*0.5*self.timeStep*(abs(self.currentState[-1])**2)*self.currentState[-1]
        b[0] = 1+2*Y + 1j*0.5*self.timeStep*potential[0] + self.gamma*self.timeStep*1j*0.5*(abs(self.currentState[0])**2)
        b[-1] = 1+2*Y + 1j*0.5*self.timeStep*potential[-1] + self.gamma*self.timeStep*1j*0.5*(abs(self.currentState[-1])**2)

        diagCorrection[0], diagCorrection[-1] = - b[0], - c[-1]*a[0]/b[0]

        tridiag = np.matrix([
            c,
            b - diagCorrection,
            a,
        ])

        temp = solve_banded((1,1), tridiag, v)

        U = np.zeros(self.spaceArraySize, dtype=np.complex64)
        U[0], U[-1] = -b[0], c[-1]

        V = np.zeros(self.spaceArraySize, dtype=np.complex64)
        V[0], V[-1] = 1, -a[0]/b[0]

        complement = solve_banded((1,1), tridiag, U)

        num = np.dot(V, temp)
        den = 1 + np.dot(V, complement)

        self.currentState = temp  - (num/den)*complement
WiFO215
fonte
3
Parece (à primeira vista) como um bug em suas condições de contorno periódicas. Gostaria de publicar um trecho de código?
David Ketcheson
2
Bem-vindo ao Stack Exchange! No futuro, se você tiver várias perguntas, faça-as separadamente.
Dan
Além disso: o que exatamente você quer dizer "defina v [0] e v [-1] em termos um do outro"? Você está definindo os elementos vetoriais iguais entre si após a resolução ou está usando um elemento fora do tridiagonal para acoplá-los?
Dan
Eu adicionei meu código acima. Se alguma coisa não estiver esclarecida, comunique-me por favor. Lembrarei de postar perguntas separadas da próxima vez.
WiFO215 20/01/12
Obrigado! É um pouco difícil de ler seu código devido à formatação (linhas muito longas). Além disso, comentar a parte em que você quer que as pessoas prestem atenção é confuso. Cod você escreve as equações que está resolvendo (com MathJax) usando a mesma notação que o seu código?
David Ketcheson

Respostas:

2

Segunda questão

O código que chama Scipy / Numpy geralmente é rápido apenas se puder ser vetorizado; você não deve ter nada "lento" dentro de um loop python. Mesmo assim, é praticamente inevitável que seja pelo menos um pouco mais lento do que algo usando uma biblioteca semelhante em uma linguagem compilada.

for i in np.arange(0,self.spaceArraySize-1):
            v[i] = Y*self.currentState[i+1] + (1-2*Y)*self.currentState[i]   ...
            b[i] = 1+2*Y + 1j*0.5*self.timeStep*potential[i] + ...

É isso que quero dizer com "lento em um loop python". O Python foré inaceitavelmente lento para a maioria dos aplicativos numéricos, e o Scipy / Numpy não afeta isso. Se você for usar python, esse loop interno deve ser expresso como uma ou duas funções Numpy / Scipy, que essas bibliotecas podem ou não fornecer. Se eles não fornecerem algo que permita iterar sobre matrizes como esta e acessar elementos adjacentes, python é a ferramenta errada para o que você deseja fazer.

Além disso, você está fazendo uma inversão em vez de uma solução de vetor de matriz. Uma inversão de matriz, seguida por uma multiplicação de vetor de matriz, é muito mais lenta que uma solução de vetor de matriz. Isso é quase certamente a coisa que torna seu código mais lento do que qualquer outra coisa.

Se você deseja usar o C / C ++, o GSL está faltando no que diz respeito à álgebra linear complexa. Eu recomendaria usar BLAS ou LAPACK diretamente, ou usar uma biblioteca como PETSc ou Trilinos. Se você possui o MKL instalado, também pode usá-lo. Você também pode querer conferir o Fortran 2008, que é orientado a objetos.

Suas matrizes são esparsas, portanto, certifique-se de usar bibliotecas esparsas.

Eu diria também que o que você está fazendo aqui parece baixo o suficiente para que a orientação a objetos provavelmente não seja sua principal preocupação. Uma matriz Fortran 90+ é provavelmente uma boa combinação para o que você precisa, e os compiladores F90 podem paralelizar automaticamente alguns loops.

Além disso, você pode querer conferir o Octave ou o Matlab, que têm a sparse()função. Se usado corretamente, eles devem poder ser executados rapidamente.

Dan
fonte
Certamente examinarei o Fortran 2008. Já tenho a matriz 'quase tridiagonal'. Mencionei acima que estava usando o algoritmo Sherman Morrison.
WiFO215 20/01/12
ATUALIZAÇÃO: Eu tenho tentado ler sobre o ScaLAPACK porque parece muito interessante. Ele permite inverter as matrizes usando uma palavra que estou ouvindo muito "em paralelo". Tudo o que sei é que ele usa todos os meus processadores e, portanto, vai mais rápido, mas, além disso, não entendo do que se trata. Vindo de uma formação em física, a única exposição à computação que tenho é com 101 cursos em Python e C. Aprender a usar isso levará tempo. A documentação em si não se presta a uma leitura limpa.
WiFO215 20/01/12
ATUALIZAÇÃO 2: Cara! Essa coisa do ScaLAPACK parece realmente complicada. Eu não entendo a cabeça ou o rabo do que está no site. Estou simplesmente nadando em todas as informações.
WiFO215
ATUALIZAÇÃO 3: Tudo bem, eu já passei pelas outras recomendações PETSc e Trilinos. Minha ligação final é que acho que não vou usá-las agora, pois elas parecem muito complicadas. Isso não significa que não vou lê-los. Começarei a lê-los agora, mas quando compreender e sou capaz de implementá-los, os meses já teriam passado. Vou abrir um tópico separado para minhas perguntas acima, pois estou tendo dificuldades com isso. Mas isso é para mais tarde. Agora, estou de volta para combater única questão 1.
WiFO215
Eu atualizei minha resposta.
Dan