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