9
$\begingroup$

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
$\endgroup$
  • 3
    $\begingroup$ It sounds (at first glance) like a bug in your periodic boundary conditions. Care to post a code snippet? $\endgroup$ – David Ketcheson Jan 19 '12 at 20:45
  • 2
    $\begingroup$ Welcome to Stack Exchange! In the future, if you have several questions, you may want to ask them separately. $\endgroup$ – Dan Jan 19 '12 at 22:49
  • $\begingroup$ Also: What exactly do you mean "set v[0] and v[-1] in terms of each other"? Are you setting the vector elements equal to each other after the solve, or are you using an off-tridiagonal element to couple them? $\endgroup$ – Dan Jan 19 '12 at 22:51
  • $\begingroup$ I have added my code above. If anything is unclear, please let me know. I shall remember to post separate questions next time. $\endgroup$ – WiFO215 Jan 20 '12 at 5:26
  • $\begingroup$ Thanks! It's a bit hard to read your code due to the formatting (very long lines). Also, commenting out the very part that you want people to pay attention to is confusing. Cod you write down the equations you are solving (with MathJax) using the same notation as your code? $\endgroup$ – David Ketcheson Jan 20 '12 at 6:18
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I shall certainly look into Fortran 2008. I have already got the 'almost tridiagonal' matrix. I mentioned above that I was using the Sherman Morrison algorithm. $\endgroup$ – WiFO215 Jan 20 '12 at 5:14
  • $\begingroup$ UPDATE: I've been trying to read up on ScaLAPACK because it looks very interesting. It allows for one to invert the matrices using a buzz word I've been hearing a lot "in parallel". All I know is that it uses all my processors and it therefore goes faster, but beyond that, I don't understand what it's about. Coming from a physics background, the only exposure to computing I have is with 101 courses in Python and C. Learning how to use this is going to take time. The documentation itself doesn't lend itself to clean reading. $\endgroup$ – WiFO215 Jan 20 '12 at 17:15
  • $\begingroup$ UPDATE 2: Man! This ScaLAPACK thing looks really complicated. I don't understand head or tail of what is on the website. I'm simply swimming in all the information. $\endgroup$ – WiFO215 Jan 20 '12 at 18:02
  • $\begingroup$ UPDATE 3: Alright, I've gone through the other recommendations PETSc and Trilinos. My final call is that I don't think I will be using these now as they look very complicated. That does not mean I will not read them. I will start reading them now, but by the time I comprehend and am capable of implementing them, months would have passed. I shall open a separate thread for my questions on the above as I am having difficulty with it.But that is for later. Now, I am back to tackling only question 1. $\endgroup$ – WiFO215 Jan 20 '12 at 18:47
  • $\begingroup$ I've updated my answer. $\endgroup$ – Dan Jan 20 '12 at 21:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.