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

Homework Help: Numerically find the energy of the delta-well's bound state

  1. Dec 20, 2017 #1
    I'm working on an assignment where I'm required to numerically find the energy of a delta-potential's bound state. To do this, we've converted the time-independent schrödinger equation to an eigenvalue problem with E the eigen value, ψ the eigen vector and H a matrix as follows:
    391fcef12d.png
    with ##t = \frac{\hbar^2}{2ma^2}##.

    Now, my specific problem is this:
    3e839211e6.png

    I've already done all the theoretical calculations, and am just left with figuring out the code. I've tried a bunch of things, but I haven't been able to get the bound state to appear. My code is as follows:


    Code (Python):
    import numpy as np
    import matplotlib
    matplotlib.use("TkAgg")
    import matplotlib.pyplot as plt

    #################################################
    # Parameter input
    #################################################

    hbar = 1.05e-34
    h    = 6.626e-34
    m    = 1
    N    = 400
    a    = 1
    L    = N * a
    t    = hbar**2 / ( 2 * m * a**2)

    #################################################
    # Parameter sanity check
    #################################################
    assert N%2 == 0, 'N must be an even number'
    assert m    > 0, 'Mass must be positive'
    assert a    > 0, 'Distance between points must be positive'
    assert t   >= 0, 'T must be non-negative. Make sure hbar is positive.'

    #################################################
    # Base code
    #################################################
    for i in range(0,100,10):
        # Building the hamiltonian operator matrix and calculating its eigen
        # values (eval) and eigen vectors (evec)
        n  = int(N/2)
        a1 = np.full((N-1,), -t)
        a2 = np.full((N,), 2*t)
        Hamiltonian = np.diag(a1, -1) +  np.diag(a2) + np.diag(a1, 1)
        Hamiltonian[n,n] = Hamiltonian[n,n]- 1*10**i  # Adds potential in middle.
        eval, evec = np.linalg.eigh(Hamiltonian) # use .eig() if Hamiltonian is not symmetric

        plt.plot(evec)
        plt.show()
    I'm quite desperate, hence the post here. Thanks in advance for all your help!

    For those amongst you who'd prefer a nicely formatted image of the assignment: https://puu.sh/yKp4D/b668471918.png
     
    Last edited: Dec 20, 2017
  2. jcsd
  3. Dec 20, 2017 #2

    Ray Vickson

    User Avatar
    Science Advisor
    Homework Helper

    I do not see any place where you are computing an eigenvalue or an eigenvector, unless is is in a procedure that is called somewhere.

    Anyway, the type of approach you attempt is quite possibly bound to fail, because of numerical scaling problems and severe roundoff error problems. Your matrix H has a possibly small entry ##t##. It is much better to set ##H = t A## where ##A## is the matrix you would get if you put ##t = 1##. Now, if ##\lambda## is an eigenvalue of ##A## and ##v## is an eigenvector, then ##t \lambda## is an eigenvalue of ##H = tA## and ##v## is an eigenvector (of both ##A## and ##H##). All the numerical work of finding eigenvalues/eigenvectors of ##A## is divorced from the numerical value of ##t##. However, you can just multiply your answer by ##t## as the very last step of you solution procedure, after all the hard work has been done.
     
  4. Dec 20, 2017 #3

    DrClaude

    User Avatar

    Staff: Mentor

    According to your calculations, at what is the energy of the bound state?

    By the way, it looks strange that are are using SI units for ħ but take the mass and a to be 1. I would use units such that ħ = 1.
     
  5. Dec 20, 2017 #4

    DrClaude

    User Avatar

    Staff: Mentor

    It's the action of numpy.linalg.eigh.
     
  6. Dec 20, 2017 #5
    Right. That makes a lot of sense (and explains why they're telling my to vary ##V_0/t##.. Changing this little thing gets a much nicer curve, I'll still need to verify if it's the correct curve, but I call this progress :) Thanks!
    Edit: Turns out I'm ~2 orders of magnitude off now.. Progress indeed.

    That's my doing. I figured I knew ##\hbar## but not the rest so I only kept it in SI units.

    According to my calculations, the bound state occurs at ##E = -\frac{m}{2\hbar^2}##. Note how I leave out a factor (as I do in my code) because I don't know how to compensate for ##\lambda\delta(x)##.
     
    Last edited: Dec 20, 2017
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted