Schrödinger-Gleichung mit periodischen Randbedingungen

9

Ich habe einige Fragen zu folgenden Themen:

Ich versuche, die Schrödinger-Gleichung in 1D unter Verwendung der Kurbel-Nicolson-Diskretisierung zu lösen, gefolgt von der Invertierung der resultierenden tridiagonalen Matrix. Mein Problem hat sich nun zu einem Problem mit periodischen Randbedingungen entwickelt, und deshalb habe ich meinen Code geändert, um den Sherman Morrison-Algorithmus zu verwenden.

Angenommen, vmeine RHS befindet sich zu jedem Zeitpunkt, zu dem ich die tridiagonale Matrix invertieren möchte. Die Größe von vist die Anzahl der Gitterpunkte, die ich über dem Raum habe. Wenn ich setze v[0]und v[-1]in Bezug aufeinander, wie es in meiner periodischen Situation erforderlich ist, explodiert meine Gleichung. Ich kann nicht sagen, warum das passiert. Ich benutze python2.7 und scipy's eingebautes Solve_Banded, um die Gleichung zu lösen.

Dies führt mich zu meiner zweiten Frage: Ich habe Python verwendet, weil es die Sprache ist, die ich am besten kenne, aber ich finde es ziemlich langsam (selbst mit den Optimierungen, die von numpy und scipy angeboten werden). Ich habe versucht, C ++ zu verwenden, da ich damit einigermaßen vertraut bin. Ich dachte, ich würde die GSL verwenden, die BLAS-optimiert wäre, fand aber keine Dokumentation, um komplexe Vektoren zu erstellen oder die tridiagonale Matrix mit solch komplexen Vektoren zu lösen.

Ich möchte Objekte in meinem Programm haben, da ich der Meinung bin, dass dies für mich der einfachste Weg ist, später zu verallgemeinern, um die Kopplung zwischen Wellenfunktionen einzuschließen, sodass ich mich an eine objektorientierte Sprache halte.

Ich könnte versuchen, den tridiagonalen Matrixlöser von Hand zu schreiben, aber ich hatte Probleme, als ich dies in Python tat. Als ich mich über große Zeiten mit immer feineren Zeitschritten weiterentwickelte, häufte sich der Fehler und gab mir Unsinn. Vor diesem Hintergrund habe ich mich für die integrierten Methoden entschieden.

Jeder Rat wird sehr geschätzt.

EDIT: Hier ist das relevante Code-Snippet. Die Notation stammt von der Wikipedia-Seite zur TDM-Gleichung (Tridiagonal Matrix). v ist die RHS des Kurbel-Nicolson-Algorithmus zu jedem Zeitschritt. Die Vektoren a, b und c sind die Diagonalen des TDM. Der korrigierte Algorithmus für den periodischen Fall stammt aus dem CFD-Wiki . Ich habe ein wenig umbenannt. Was sie u, v genannt haben Ich habe U, V (groß geschrieben) genannt. Ich habe q das Komplement, y die temporäre Lösung und die eigentliche Lösung self.currentState genannt. Die Zuordnung von v [0] und v [-1] verursacht hier das Problem und wurde daher auskommentiert. Sie können die Faktoren von Gamma ignorieren. Sie sind nichtlineare Faktoren, die zur Modellierung von Bose-Einstein-Kondensaten verwendet werden.

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
quelle
3
Es klingt (auf den ersten Blick) wie ein Fehler in Ihren periodischen Randbedingungen. Möchtest du ein Code-Snippet posten?
David Ketcheson
2
Willkommen bei Stack Exchange! Wenn Sie in Zukunft mehrere Fragen haben, möchten Sie diese möglicherweise separat stellen.
Dan
Außerdem: Was genau meinst du mit "setze v [0] und v [-1] in Bezug aufeinander"? Setzen Sie die Vektorelemente nach dem Lösen gleich oder verwenden Sie ein nicht tridiagonales Element, um sie zu koppeln?
Dan
Ich habe meinen Code oben hinzugefügt. Wenn etwas unklar ist, lassen Sie es mich bitte wissen. Ich werde daran denken, das nächste Mal separate Fragen zu stellen.
WiFO215
Vielen Dank! Aufgrund der Formatierung (sehr lange Zeilen) ist es etwas schwierig, Ihren Code zu lesen. Es ist auch verwirrend, genau den Teil zu kommentieren, auf den die Leute achten sollen. Codieren Sie die Gleichungen, die Sie lösen (mit MathJax), mit derselben Notation wie Ihr Code?
David Ketcheson

Antworten:

2

Zweite Frage

Code, der Scipy / Numpy aufruft, ist normalerweise nur dann schnell, wenn er vektorisiert werden kann. Sie sollten nichts "Langsames" in einer Python-Schleife haben. Selbst dann ist es ziemlich unvermeidlich, dass es zumindest etwas langsamer ist als etwas, das eine ähnliche Bibliothek in einer kompilierten Sprache verwendet.

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

Das meine ich mit "langsam in einer Python-Schleife". Pythons forist für die meisten numerischen Anwendungen inakzeptabel langsam, und Scipy / Numpy haben keinen Einfluss darauf. Wenn Sie Python verwenden, sollte diese innere Schleife als eine oder zwei Numpy / Scipy-Funktionen ausgedrückt werden, die diese Bibliotheken möglicherweise bereitstellen oder nicht. Wenn sie nichts bieten, mit dem Sie über Arrays wie dieses iterieren und auf benachbarte Elemente zugreifen können, ist Python das falsche Werkzeug für das, was Sie tun möchten.

Außerdem führen Sie eher eine Inversion als eine Matrixvektorlösung durch. Eine Matrixinversion, gefolgt von einer Matrixvektormultiplikation, ist viel langsamer als eine Matrixvektorlösung. Dies ist mit ziemlicher Sicherheit die Sache, die Ihren Code mehr als alles andere verlangsamt.

Wenn Sie C / C ++ verwenden möchten, fehlt GSL, wenn es um komplexe lineare Algebra geht. Ich würde empfehlen, entweder BLAS oder LAPACK direkt zu verwenden oder eine Bibliothek wie PETSc oder Trilinos zu verwenden. Wenn Sie MKL installiert haben, können Sie dies auch verwenden. Vielleicht möchten Sie auch Fortran 2008 ausprobieren, das objektorientiert ist.

Ihre Matrizen sind spärlich, daher sollten Sie sicherstellen, dass Sie spärliche Bibliotheken verwenden.

Ich würde auch sagen, dass das, was Sie hier tun, so niedrig erscheint, dass die Objektorientierung wahrscheinlich nicht Ihr Hauptanliegen sein sollte. Ein Fortran 90+ -Array passt wahrscheinlich ziemlich gut zu Ihren Anforderungen, und F90-Compiler können einige Loops automatisch parallelisieren.

Vielleicht möchten Sie auch Octave oder Matlab ausprobieren, die diese sparse()Funktion haben. Bei richtiger Verwendung sollten diese ziemlich schnell ausgeführt werden können.

Dan
quelle
Ich werde sicherlich auf Fortran 2008 eingehen. Ich habe bereits die "fast tridiagonale" Matrix. Ich habe oben erwähnt, dass ich den Sherman Morrison-Algorithmus verwendet habe.
WiFO215
UPDATE: Ich habe versucht, ScaLAPACK nachzulesen, weil es sehr interessant aussieht. Es ermöglicht einem, die Matrizen mit einem Modewort zu invertieren, das ich viel "parallel" gehört habe. Ich weiß nur, dass es alle meine Prozessoren verwendet und daher schneller geht, aber darüber hinaus verstehe ich nicht, worum es geht. Ich komme aus der Physik und bin mit 101 Kursen in Python und C nur mit Computern vertraut. Das Erlernen der Verwendung wird einige Zeit in Anspruch nehmen. Die Dokumentation selbst eignet sich nicht zum sauberen Lesen.
WiFO215
UPDATE 2: Mann! Diese ScaLAPACK-Sache sieht wirklich kompliziert aus. Ich verstehe weder Kopf noch Schwanz von dem, was auf der Website steht. Ich schwimme einfach in allen Informationen.
WiFO215
UPDATE 3: Okay, ich habe die anderen Empfehlungen PETSc und Trilinos durchgesehen. Mein letzter Anruf ist, dass ich nicht glaube, dass ich diese jetzt verwenden werde, da sie sehr kompliziert aussehen. Das heißt nicht, dass ich sie nicht lesen werde. Ich werde jetzt anfangen, sie zu lesen, aber bis ich sie verstehe und in der Lage bin, sie umzusetzen, wären Monate vergangen. Ich werde einen separaten Thread für meine Fragen zu den oben genannten Themen eröffnen, da ich Schwierigkeiten damit habe. Aber das ist für später. Jetzt bin ich zurück, um nur Frage 1
anzugehen.
Ich habe meine Antwort aktualisiert.
Dan