Solution of 1D Schrodinger using discretized hamiltonian

[beady]

Hi all,
I'm trying to solve the 1D Schrödinger equation for an arbritary potential, to calculate Franck Condon factors for absorption and emmision spectra. I can do this using iterative techniques (e.g. the Numerov method), but I can't seem to get it to work by discretrizing the hamiltonian, and then finding the eigen vectors/values. Here's what I'm doing, hopefully someone can help.

initial hamiltonian : $$H=\frac{1}{2 I}\frac{d^2}{dx^2} + W(x)$$

Split the space up into n points.

use the standard finite difference relations for a second derivative leads to:
$$\frac{1}{2I}\frac{v_{i-1}-2v_{i}+v_{i+1}}{\delta^2}+W_i v_i = E v_i$$

and the relevant forward and backward differences for the endpoints.

I then set up a matrix for the above linear equations, and then just solve the system using either Matlab or Mathematica eigenvalue and eigenvector routines.

The problem is that the resulting eigenvalues ( energies ) are completely wrong, and the wavefunctions ( individual eigenvectors ) are weird. It seems that the eigenvectors oscillate around around y=0. It's difficult to explain, but if the desired wavefunction is phi, then at even grid points (n=2,4,6...) , I get the correct phi(x) value, but for odd grid points (n=1,3,5,7...), I get -phi(x). This is easily corrected, but why is it happening?

Is it possible to calculate the wavefunction using this method? I may have made a coding error, but I'm interested to know if the method above should work, and if not, why not.

Thanks in advance

 

[inrainbows]

Hi,

did you ever have any luck solving this? I have the same problem using python (using the numpy.linalg package to perform the diagonalization). What's more it seems to persist even for a 5 point stencil for the second derivative.

 

[newbee]

What is the form of the potential for which you are getting the bad results.

 

[inrainbows]

Hi, actually I am just using a square well potential at the moment (for testing) - ie. V(x)=0. In my case, I suspect it is something to do with using the python linalg.eigh implementation, as it seems to work in Fortran with a call to the lapack function dstevx (which diagonalizes a tridiagonal matrix).

 

[newbee]

Ah. Well we use NumPy and SciPy around here so if there is a problem with it's linear algebra package we would certainly want to know about it. Please keep me informed. Thanks.

 

[f95toli]

Solving the 1D SE by calculatint eigenvectors/eigenvalues is fairly easy. The only "tricky" part is setting up the matrix (and relazing that it is is tri-diagonal).
However, there is one thing you need to watch out for: it is a good idea to put "walls" around the region of interest (i.e. set V(x)=a very large number) or otherwise the wavefunctions might "leak out" in non-physical way with some solvers (in principle you should also be able to simply use a very "deep" well but that is not very efficient).
You obviously have to make sure that you put the "walls" at an energy much higher than the eigenenergies you are interested in, or otherwise they will affect the result.