Уравнение Шредингера с периодическими граничными условиями

У меня есть пара вопросов относительно следующего:

Я пытаюсь решить уравнение Шредингера в 1D, используя дискретизацию кривошипа Николсона с последующим инвертированием получающейся трехдиагональной матрицы. Моя проблема теперь превратилась в проблему с периодическими граничными условиями, и поэтому я изменил свой код для использования алгоритма Шермана Моррисона.

Предположим v, моя RHS находится на каждом временном шаге, когда я хочу инвертировать трехдиагональную матрицу. Размер v- это количество точек сетки, которые у меня есть в пространстве. Когда я установил v[0]и v[-1]с точки зрения друг друга , как это требуется в моей периодической ситуации, мое уравнение взрывает. Я не могу сказать, почему это происходит. Я использую Python2.7 и встроенный в Scipy решениеущество для решения уравнения.

Это подводит меня ко второму вопросу: я использовал python, потому что это язык, который я знаю лучше всего, но я нахожу его довольно медленным (даже с оптимизацией, предлагаемой numpy и scipy). Я пытался использовать C ++, так как я достаточно знаком с ним. Я думал, что буду использовать GSL, который будет оптимизирован для BLAS, но не нашел документации для создания сложных векторов или решения трехдиагональной матрицы с такими комплексными векторами.

Я хотел бы, чтобы объекты в моей программе были, так как я чувствую, что для меня было бы проще всего обобщить позже, чтобы включить связь между волновыми функциями, поэтому я придерживаюсь объектно-ориентированного языка.

Я мог бы попытаться написать решатель трехдиагональной матрицы вручную, но у меня возникли проблемы, когда я делал это в Python. Когда я эволюционировал в течение длительного времени с более мелкими и более мелкими временными шагами, ошибка накапливалась и давала мне глупость. Помня об этом, я решил использовать встроенные методы.

Любой совет высоко ценится.

РЕДАКТИРОВАТЬ: Вот соответствующий фрагмент кода. Запись заимствована из страницы Википедии об уравнении трехдиагональной матрицы (TDM). v - RHS алгоритма Никольсона для каждого временного шага. Векторы a, b и c являются диагоналями TDM. Исправленный алгоритм для периодического случая взят из CFD Wiki . Я сделал небольшое переименование. То, что они назвали u, v Я назвал U, V (с большой буквы). Я назвал q дополнение, y временное решение и фактическое решение self.currentState. Назначение v [0] и v [-1] является причиной проблемы и поэтому было закомментировано. Вы можете игнорировать факторы гамма. Это нелинейные факторы, используемые для моделирования бозе-эйнштейновских конденсатов.

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

Второй вопрос

Код, который вызывает Scipy / Numpy, обычно быстр, только если его можно векторизовать; у вас не должно быть ничего "медленного" внутри цикла питона. Даже в этом случае неизбежно, что это будет, по крайней мере, немного медленнее, чем что-либо, использующее подобную библиотеку на скомпилированном языке.

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

Это то, что я подразумеваю под «медленно в цикле питона». Python forнедопустимо медленен для большинства числовых приложений, и Scipy / Numpy не влияют на это вообще. Если вы собираетесь использовать python, этот внутренний цикл должен быть выражен в виде одной или двух функций Numpy / Scipy, которые эти библиотеки могут предоставлять или не предоставлять. Если они не предоставляют что-то, что позволяет вам перебирать массивы, подобные этой, и получать доступ к смежным элементам, python - неправильный инструмент для того, что вы хотите сделать.

Кроме того, вы делаете инверсию, а не матрично-векторное решение. Инверсия матрицы, за которой следует умножение матрицы на вектор, происходит намного медленнее, чем матрица на вектор. Это почти наверняка замедляет ваш код больше, чем что-либо еще.

Если вы хотите использовать C / C ++, GSL не хватает, когда речь идет о сложной линейной алгебре. Я бы порекомендовал либо использовать BLAS или LAPACK напрямую, либо использовать такие библиотеки, как PETSc или Trilinos. Если у вас установлен MKL, вы тоже можете его использовать. Вы также можете попробовать Fortran 2008, который является объектно-ориентированным.

Ваши матрицы разрежены, поэтому вы должны убедиться, что используете разреженные библиотеки.

Я также сказал бы, что то, что вы делаете здесь, кажется достаточно низкоуровневым, поэтому объектная ориентация, вероятно, не должна быть вашей главной задачей. Массив Fortran 90+, вероятно, очень хорошо соответствует тому, что вам нужно, и компиляторы F90 могут автоматически распараллеливать некоторые циклы.

Кроме того, вы можете проверить Octave или Matlab, которые имеют sparse()функцию. При правильном использовании они должны работать довольно быстро.