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

Tags:
1. Dec 20, 2017

### Muizz

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:

with $t = \frac{\hbar^2}{2ma^2}$.

Now, my specific problem is this:

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. Dec 20, 2017

### Ray Vickson

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.

3. Dec 20, 2017

### 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.

4. Dec 20, 2017

### Staff: Mentor

It's the action of numpy.linalg.eigh.

5. Dec 20, 2017

### Muizz

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