Équation de Schrodinger avec conditions aux limites périodiques

9

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
WiFO215
la source
3
Cela ressemble (à première vue) à un bug dans vos conditions aux limites périodiques. Vous souhaitez publier un extrait de code?
David Ketcheson
2
Bienvenue sur Stack Exchange! À l'avenir, si vous avez plusieurs questions, vous voudrez peut-être les poser séparément.
Dan
Aussi: Que voulez-vous dire exactement "définir v [0] et v [-1] en termes l'un par rapport à l'autre"? Définissez-vous les éléments vectoriels égaux les uns aux autres après la résolution, ou utilisez-vous un élément hors tridiagonal pour les coupler?
Dan
J'ai ajouté mon code ci-dessus. Si quelque chose n'est pas clair, faites-le moi savoir. Je me souviendrai de poster des questions distinctes la prochaine fois.
WiFO215
Merci! La lecture de votre code est un peu difficile à cause du formatage (très longues lignes). De plus, commenter la partie à laquelle vous voulez que les gens prêtent attention est déroutant. Cod vous écrivez les équations que vous résolvez (avec MathJax) en utilisant la même notation que votre code?
David Ketcheson

Réponses:

2

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.

Dan
la source
Je vais certainement regarder Fortran 2008. J'ai déjà la matrice «presque tridiagonale». J'ai mentionné ci-dessus que j'utilisais l'algorithme Sherman Morrison.
WiFO215
MISE À JOUR: J'ai essayé de lire sur ScaLAPACK parce que ça a l'air très intéressant. Il permet d'inverser les matrices en utilisant un mot à la mode que j'entends beaucoup "en parallèle". Tout ce que je sais, c'est qu'il utilise tous mes processeurs et qu'il va donc plus vite, mais au-delà, je ne comprends pas de quoi il s'agit. Issu d'un milieu physique, la seule exposition à l'informatique que j'ai est avec 101 cours en Python et C. Apprendre à utiliser cela va prendre du temps. La documentation elle-même ne se prête pas à une lecture propre.
WiFO215
MISE À JOUR 2: Mec! Cette chose ScaLAPACK semble vraiment compliquée. Je ne comprends pas ce qui se trouve sur le site Web. Je nage simplement dans toutes les informations.
WiFO215
MISE À JOUR 3: D'accord, j'ai parcouru les autres recommandations PETSc et Trilinos. Mon dernier appel est que je ne pense pas que je vais les utiliser maintenant car ils ont l'air très compliqués. Cela ne veut pas dire que je ne les lirai pas. Je vais commencer à les lire maintenant, mais au moment où je les comprendrai et serai capable de les mettre en œuvre, des mois se seraient écoulés. Je vais ouvrir un fil séparé pour mes questions sur ce qui précède car j'ai du mal à le faire, mais c'est pour plus tard. Maintenant, je suis de retour à aborder uniquement la question 1.
WiFO215
J'ai mis à jour ma réponse.
Dan