Équation de Schrodinger avec conditions aux limites périodiques

J'ai quelques questions concernant les points suivants:

J'essaie de résoudre l'équation de Schrodinger en 1D en utilisant la discrétisation de manivelle nicolson suivie de l'inversion de la matrice tridiagonale résultante. Mon problème a maintenant évolué vers un problème avec des conditions aux limites périodiques et j'ai donc modifié mon code pour utiliser l'algorithme Sherman Morrison.

Supposons que vmon RHS soit à chaque pas de temps lorsque je souhaite inverser la matrice tridiagonale. La taille de vest le nombre de points de grille que j'ai sur l'espace. Quand je me mets v[0]et v[-1]en termes les uns des autres comme cela est requis dans ma situation périodique, mon équation explose. Je ne peux pas dire pourquoi cela se produit. J'utilise python2.7 et le module de résolution intégré de scipy pour résoudre l'équation.

Cela m'amène à ma deuxième question: j'ai utilisé python parce que c'est le langage que je connais le mieux, mais je le trouve plutôt lent (même avec les optimisations offertes par numpy et scipy). J'ai essayé d'utiliser C ++ car je le connais assez bien. J'ai pensé utiliser le GSL qui serait optimisé pour BLAS, mais je n'ai trouvé aucune documentation pour créer des vecteurs complexes ou résoudre la matrice tridiagonale avec des vecteurs de valeur complexes.

Je voudrais des objets dans mon programme car je pense que ce serait le moyen le plus simple pour moi de généraliser plus tard pour inclure le couplage entre les fonctions d'onde, donc je m'en tiens à un langage orienté objet.

Je pouvais essayer d'écrire le solveur de matrice tridiagonale à la main, mais j'ai rencontré des problèmes quand je l'ai fait en python. Alors que j'évoluais sur de grandes périodes avec des pas de temps de plus en plus fins, l'erreur s'est accumulée et m'a donné un non-sens. Gardant cela à l'esprit, j'ai décidé d'utiliser les méthodes intégrées.

Tout conseil est fort apprécié.

EDIT: voici l'extrait de code correspondant. La notation est empruntée à la page de Wikipedia sur l'équation de la matrice tridiagonale (TDM). v est le RHS de l'algorithme crank nicolson à chaque pas de temps. Les vecteurs a, b et c sont les diagonales du TDM. L'algorithme corrigé pour le cas périodique provient du CFD Wiki . J'ai fait un petit changement de nom. Ce qu'ils ont appelé u, v J'ai appelé U, V (en majuscule). J'ai appelé q le complément, y la solution temporaire et la solution réelle self.currentState. L'affectation de v [0] et v [-1] est la cause du problème ici et a donc été commentée. Vous pouvez ignorer les facteurs de gamma. Ce sont des facteurs non linéaires utilisés pour modéliser les condensats 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

Deuxième question

Le code qui appelle Scipy / Numpy n'est généralement rapide que s'il peut être vectorisé; vous ne devriez rien avoir de "lent" à l'intérieur d'une boucle python. Même alors, il est à peu près inévitable qu'il soit au moins un peu plus lent que quelque chose utilisant une bibliothèque similaire dans un langage compilé.

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] + ...

C'est ce que je veux dire par "lent dans une boucle python". Python forest trop lent pour la plupart des applications numériques, et Scipy / Numpy n'affecte pas cela du tout. Si vous allez utiliser python, cette boucle interne doit être exprimée comme une ou deux fonctions Numpy / Scipy, que ces bibliothèques peuvent fournir ou non. S'ils ne fournissent pas quelque chose qui vous permet d'itérer sur des tableaux comme celui-ci et d'accéder à des éléments adjacents, python n'est pas le bon outil pour ce que vous voulez faire.

En outre, vous faites une inversion plutôt qu'une résolution matricielle. Une inversion matricielle, suivie d'une multiplication matrice-vecteur, est beaucoup plus lente qu'une résolution matricielle-vectorielle. C'est presque certainement la chose qui ralentit votre code plus que toute autre chose.

Si vous voulez utiliser C / C ++, GSL manque en quelque sorte d'algèbre linéaire complexe. Je recommanderais d'utiliser BLAS ou LAPACK directement, ou d'utiliser une bibliothèque comme PETSc ou Trilinos. Si vous avez installé MKL, vous pouvez également l'utiliser. Vous pouvez également consulter Fortran 2008, qui est orienté objet.

Vos matrices sont clairsemées, vous voudrez donc vous assurer d'utiliser des bibliothèques clairsemées.

Je dirais également que ce que vous faites ici semble suffisamment bas pour que l'orientation objet ne soit probablement pas votre principale préoccupation. Un tableau Fortran 90+ est probablement assez bon pour répondre à vos besoins, et les compilateurs F90 peuvent paralléliser automatiquement certaines boucles.

En outre, vous voudrez peut-être consulter Octave ou Matlab, qui ont la sparse()fonction. S'ils sont utilisés correctement, ils devraient pouvoir fonctionner assez rapidement.