Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Solution of 1D Schrodinger using discretized hamiltonian

  1. Sep 11, 2006 #1
    Hi all,
    I'm trying to solve the 1D Schrodinger 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 : [tex]H=\frac{1}{2 I}\frac{d^2}{dx^2} + W(x) [/tex]

    Split the space up into n points.

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

    and the relevent 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
    Last edited: Sep 11, 2006
  2. jcsd
  3. Oct 29, 2008 #2

    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.
  4. Oct 29, 2008 #3
    What is the form of the potential for which you are getting the bad results.
  5. Oct 29, 2008 #4
    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).
  6. Oct 29, 2008 #5
    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.
  7. Oct 29, 2008 #6


    User Avatar
    Science Advisor
    Gold Member

    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.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook