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
fonte
Respostas:
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.
É 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.fonte