Schrodinger equation with periodic boundary conditions

Question

I have a couple of questions regarding the following:

I am trying to solve the Schrodinger equation in 1D using the crank nicolson discretization followed by inverting the resulting tridiagonal matrix. My problem has now evolved into a problem with periodic boundary conditions and so I have modified my code to use the Sherman Morrison algorithm.

Suppose v is my RHS in at each time step when I wish to invert the tridiagonal matrix. The size of v is the number of grid points I have over space. When I set v[0] and v[-1] in terms of each other as is required in my periodic situation, my equation blows up. I cannot tell why this is happening. I am using python2.7 and scipy's inbuilt solve_banded to solve the equation.

This leads me to my second question: I used python because it is the language I know best, but I find it rather slow (even with the optimizations offered by numpy and scipy). I have tried using C++ as I am reasonably familiar with it. I thought I'd use the GSL which would be BLAS optimized, but found no documentation to create complex vectors or solve the tridiagonal matrix with such complex valued vectors.

I would like objects in my program as I feel it would be the easiest way for me to generalize later to include coupling between wavefunctions thus I am sticking to an object oriented language.

I could try writing the tridiagonal matrix solver by hand, but I ran into problems when I did so in python. As I evolved over large times with finer and finer time steps, the error accumulated and gave me nonsense. Keeping this in mind, I decided to use the in-built methods.

Any advice is much appreciated.

EDIT: Here is the relevant code snippet. The notation is borrowed from Wikipedia's page on the tridiagonal matrix (TDM) equation. v is the RHS of the crank nicolson algorithm at each time step. The vectors a, b and c are the diagonals of the TDM. The corrected algorithm for the periodic case is from the CFD Wiki. I have done a little renaming. What they have called u, v I have called U, V (capitalized). I have called q the complement, y the temporary solution and the actual solution self.currentState. The assignment of v[0] and v[-1] is what is causing the problem here and thus has been commented out. You may ignore the factors of gamma. They are non-linear factors used to model Bose Einstein Condensates.

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

Answer 1

Second question

Code that calls Scipy/Numpy is usually only fast if it can be vectorized; you shouldn't have anything "slow" inside of a python loop. Even then, it's pretty much unavoidable that it will be at least a little slower than something using a similar library in a compiled language.

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

This is what I mean by "slow in a python loop". Python's for is unacceptably slow for most numerical applications, and Scipy/Numpy do not affect this at all. If you're going to use python, this inner loop should be expressed as one or two Numpy/Scipy functions, which those libraries may or may not provide. If they don't provide something that allows you to iterate over arrays like this and access adjacent elements, python is the wrong tool for what you want to do.

Also, you're doing an inversion rather than a matrix-vector solve. A matrix inversion, followed by a matrix-vector multiply, is much slower than a matrix-vector solve. This is almost certainly the thing slowing down your code more than anything else.

If you want to use C/C++, GSL is kind of lacking when it comes to complex linear algebra. I'd recommend either using BLAS or LAPACK directly, or using a library like PETSc or Trilinos. If you have MKL installed, you can use that too. You might also want to check out Fortran 2008, which is object-oriented.

Your matrices are sparse, so you'll want to make sure you use sparse libraries.

I would also say that what you are doing here seems low-level enough that object-orientation should probably not be your primary concern. A Fortran 90+ array is probably a pretty good match to what you need, and F90 compilers can auto-parallelize some loops.

Also, you may want to check out Octave or Matlab, which have the sparse() function. If used properly, these should be able to run pretty quickly.