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