Tengo un par de preguntas sobre lo siguiente:
Estoy tratando de resolver la ecuación de Schrodinger en 1D usando la discretización de manivela nicolson seguida de invertir la matriz tridiagonal resultante. Mi problema ahora se ha convertido en un problema con condiciones de límite periódicas, por lo que he modificado mi código para usar el algoritmo Sherman Morrison.
Supongamos vque mi RHS está en cada paso de tiempo cuando deseo invertir la matriz tridiagonal. El tamaño de ves el número de puntos de cuadrícula que tengo sobre el espacio. Cuando establezco v[0]y v[-1]en términos mutuos como se requiere en mi situación periódica, mi ecuación explota. No puedo decir por qué sucede esto. Estoy usando python2.7 y scipy incorporado solve_banded para resolver la ecuación.
Esto me lleva a mi segunda pregunta: utilicé Python porque es el lenguaje que mejor conozco, pero lo encuentro bastante lento (incluso con las optimizaciones ofrecidas por numpy y scipy). He intentado usar C ++ ya que estoy bastante familiarizado con él. Pensé que usaría el GSL que estaría optimizado para BLAS, pero no encontré documentación para crear vectores complejos o resolver la matriz tridiagonal con vectores de valores tan complejos.
Me gustaría tener objetos en mi programa, ya que creo que sería la forma más fácil de generalizar más tarde para incluir el acoplamiento entre las funciones de onda, por lo que me apego a un lenguaje orientado a objetos.
Podría intentar escribir el solucionador de matriz tridiagonal a mano, pero me encontré con problemas cuando lo hice en Python. A medida que evolucioné durante grandes tiempos con pasos de tiempo cada vez más finos, el error se acumuló y me dio tonterías. Teniendo esto en cuenta, decidí usar los métodos incorporados.
Cualquier consejo es muy apreciado.
EDITAR: Aquí está el fragmento de código relevante. La notación está tomada de la página de Wikipedia en la ecuación de la matriz tridiagonal (TDM). v es el RHS del algoritmo de manivela nicolson en cada paso de tiempo. Los vectores a, byc son las diagonales del TDM. El algoritmo corregido para el caso periódico es de CFD Wiki . He cambiado un poco el nombre. Lo que han llamado u, v He llamado U, V (en mayúscula). He llamado q el complemento, y la solución temporal y la solución real self.currentState. La asignación de v [0] y v [-1] es lo que está causando el problema aquí y, por lo tanto, se ha comentado. Puede ignorar los factores de gamma. Son factores no lineales utilizados para modelar 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
fuente
Respuestas:
Segunda pregunta
El código que llama a Scipy / Numpy generalmente solo es rápido si se puede vectorizar; No deberías tener nada "lento" dentro de un bucle de Python. Incluso entonces, es casi inevitable que sea al menos un poco más lento que algo que use una biblioteca similar en un lenguaje compilado.
Esto es lo que quiero decir con "lento en un bucle de Python". Python
fores inaceptablemente lento para la mayoría de las aplicaciones numéricas, y Scipy / Numpy no afectan esto en absoluto. Si va a usar Python, este ciclo interno debe expresarse como una o dos funciones Numpy / Scipy, que esas bibliotecas pueden proporcionar o no. Si no proporcionan algo que le permita iterar sobre matrices como esta y acceder a elementos adyacentes, python es la herramienta incorrecta para lo que desea hacer.Además, está haciendo una inversión en lugar de una resolución de matriz-vector. Una inversión de matriz, seguida de una multiplicación de matriz-vector, es mucho más lenta que una resolución de matriz-vector. Esto es casi seguro que ralentiza su código más que cualquier otra cosa.
Si desea usar C / C ++, GSL es un poco escaso cuando se trata de álgebra lineal compleja. Recomiendo usar BLAS o LAPACK directamente, o usar una biblioteca como PETSc o Trilinos. Si tiene instalado MKL, también puede usarlo. También es posible que desee consultar Fortran 2008, que está orientado a objetos.
Sus matrices son escasas, por lo que querrá asegurarse de utilizar bibliotecas dispersas.
También diría que lo que está haciendo aquí parece lo suficientemente bajo como para que la orientación a objetos probablemente no sea su principal preocupación. Una matriz Fortran 90+ es probablemente una muy buena combinación con lo que necesita, y los compiladores F90 pueden auto paralelizar algunos bucles.
Además, es posible que desee consultar Octave o Matlab, que tienen la
sparse()función. Si se usa correctamente, estos deberían poder ejecutarse bastante rápido.fuente