我有几个关于以下问题的问题：

我正在尝试使用曲柄尼科尔森离散化求解一维薛定谔方程，然后反转得到的三对角矩阵。我的问题现在已经演变成周期性边界条件的问题，因此我修改了我的代码以使用 Sherman Morrison 算法。

假设v当我希望反转三对角矩阵时，我的 RHS 在每个时间步。的大小v是我在空间上的网格点数。当我按照周期性情况的要求设置v[0]和v[-1]相互关联时，我的等式就会爆炸。我不知道为什么会这样。我正在使用 python2.7 和 scipy 的内置 solve_banded 来求解方程。

这引出了我的第二个问题：我使用 python 是因为它是我最了解的语言，但我发现它相当慢（即使使用 numpy 和 scipy 提供的优化）。我已经尝试过使用 C++，因为我相当熟悉它。我以为我会使用经过 BLAS 优化的 GSL，但没有找到文档来创建复数向量或用这种复数向量求解三对角矩阵。

我想要我的程序中的对象，因为我觉得这将是我以后概括包括波函数之间的耦合的最简单方法，因此我坚持使用面向对象的语言。

我可以尝试手动编写三对角矩阵求解器，但是当我在 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