1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Numerically find the energy of the delta-well's bound state
Loading...