Równanie Schrodingera z okresowymi warunkami brzegowymi

Mam kilka pytań dotyczących następujących kwestii:

Próbuję rozwiązać równanie Schrodingera w 1D, używając dyskretyzacji korbowego nicolsonu, a następnie odwracając otrzymaną macierz tridiagonalną. Mój problem przekształcił się teraz w problem z okresowymi warunkami brzegowymi, dlatego zmodyfikowałem swój kod, aby używać algorytmu Shermana Morrisona.

Załóżmy, że vjest mój RHS na każdym etapie, kiedy chcę odwrócić matrycę tridiagonal. Rozmiar vto liczba punktów siatki, które mam nad przestrzenią. Kiedy ustawiam v[0]i v[-1]pod względem siebie nawzajem, jak jest to wymagane w mojej okresowej sytuacji, moje równanie wybuchnie. Nie umiem powiedzieć, dlaczego tak się dzieje. Korzystam z Python2.7 i wbudowanego rozwiązania Scipy Soln_banded, aby rozwiązać równanie.

To prowadzi mnie do drugiego pytania: użyłem Pythona, ponieważ jest to język, który znam najlepiej, ale uważam, że jest on raczej powolny (nawet w przypadku optymalizacji oferowanych przez Numpy i Scipy). Próbowałem używać C ++, ponieważ znam go dość dobrze. Pomyślałem, że użyję GSL, który byłby zoptymalizowany za pomocą BLAS, ale nie znalazłem dokumentacji do tworzenia złożonych wektorów lub rozwiązywania macierzy tridiagonalnej za pomocą wektorów o złożonych wartościach.

Chciałbym mieć obiekty w moim programie, ponieważ uważam, że byłby to dla mnie najłatwiejszy sposób na uogólnienie w celu włączenia sprzężenia między funkcjami falowymi, dlatego trzymam się języka zorientowanego obiektowo.

Mógłbym spróbować napisać trójwymiarowy solver macierzy ręcznie, ale napotkałem problemy, kiedy to zrobiłem w Pythonie. Gdy ewoluowałem przez długi czas z coraz to krótszymi krokami czasowymi, błąd narastał i dał mi bzdury. Mając to na uwadze, postanowiłem skorzystać z wbudowanych metod.

Wszelkie porady są mile widziane.

EDYCJA: Oto odpowiedni fragment kodu. Notacja została zapożyczona ze strony Wikipedii na równaniu macierzy tridiagonal (TDM). v jest RHS algorytmu korbowego nicolsona na każdym etapie. Wektory a, b i c są przekątnymi TDM. Poprawiony algorytm dla przypadku okresowego pochodzi z Wiki CFD . Zmieniłem trochę nazwę. To, co nazwali u, v Nazwałem U, V (wielką literą). Nazwałem q uzupełnienie, y rozwiązanie tymczasowe i rzeczywiste rozwiązanie self.currentState. Przypisanie v [0] i v [-1] jest przyczyną problemu i dlatego zostało skomentowane. Możesz zignorować czynniki gamma. Są to czynniki nieliniowe stosowane do modelowania kondensatów Bosego Einsteina.

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

Drugie Pytanie

Kod wywołujący Scipy / Numpy jest zwykle szybki tylko wtedy, gdy można go wektoryzować; nie powinieneś mieć nic „powolnego” w pętli pytona. Nawet wtedy jest prawie nieuniknione, że będzie co najmniej trochę wolniejsze niż coś używającego podobnej biblioteki w skompilowanym języku.

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

To mam na myśli przez „powolne w pętli pytona”. Python forjest niedopuszczalnie wolny dla większości aplikacji numerycznych, a Scipy / Numpy w ogóle na to nie wpływają. Jeśli zamierzasz używać Pythona, ta wewnętrzna pętla powinna być wyrażona jako jedna lub dwie funkcje Numpy / Scipy, które te biblioteki mogą, ale nie muszą, zapewniać. Jeśli nie zapewniają czegoś, co pozwala na iterację po takich tablicach i dostęp do sąsiednich elementów, python jest niewłaściwym narzędziem do tego, co chcesz zrobić.

Robisz też inwersję zamiast rozwiązania wektor macierzowy. Inwersja macierzy, po której następuje mnożenie macierz-wektor, jest znacznie wolniejsza niż rozwiązanie macierz-wektor. To prawie na pewno spowalnia Twój kod bardziej niż cokolwiek innego.

Jeśli chcesz używać C / C ++, GSL jest trochę brakuje, jeśli chodzi o złożoną algebrę liniową. Polecam albo użyć bezpośrednio BLAS lub LAPACK, albo biblioteki takiej jak PETSc lub Trilinos. Jeśli masz zainstalowany MKL, możesz go również użyć. Możesz także sprawdzić Fortran 2008, który jest zorientowany obiektowo.

Twoje macierze są rzadkie, więc upewnij się, że używasz rzadkich bibliotek.

Powiedziałbym również, że to, co tu robisz, wydaje się na tyle niskie, że orientacja obiektowa prawdopodobnie nie powinna być twoim głównym zmartwieniem. Macierz Fortran 90+ jest prawdopodobnie całkiem dobrym dopasowaniem do tego, czego potrzebujesz, a kompilatory F90 mogą automatycznie zrównoleglać niektóre pętle.

Możesz także wypróbować Octave lub Matlab, które mają tę sparse()funkcję. Przy prawidłowym użyciu powinny one działać dość szybko.