Wavefunction for a 1D quantum harmonic oscillator

Robert Johansson ([email protected])

In [1]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
In [3]:
from scipy import *
from scipy import optimize
In [4]:
from wavefunction import *
from wavefunction.wavefunction1d import *

Problem parameters

In [6]:
h = 1
h_ = h/(2*pi)
e = 1.602e-19
cf = 1          # if cf = h_, use units where h_ = 1

mm = 1          # oscillator mass
omega = 2 * pi  # oscillator frequency in GHz
x0 = 0          # shift in oscillator potiential minimum

args = {'w': omega, 'm': mm, 'x0': x0}

k = -h_ ** 2 / (2 * mm)

x_min = -pi
x_max =  pi
N = 750 

Potential

In [7]:
def U_ho(x, args):
    """
    Harmonic oscillator potential
    """
    
    omega = args['w']
    m     = args['m']
    x0    = args['x0']
    
    u = 1/2.0 * m * (omega ** 2) * ((x - x0) ** 2)

    return u

Find the minimum point for the potential

In [9]:
x = linspace(x_min, x_max, N)
In [16]:
U = U_ho(x, args);

x_opt_min = optimize.fmin(U_ho, [0.0], (args,))
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 3
         Function evaluations: 6
In [18]:
print("Found potential minima at = %f" % x_opt_min)
Found potential minima at = 0.000000

Plot the potential

In [20]:
fig, ax = subplots()

ax.plot(x, U)
ax.plot(x_opt_min, U_ho(x_opt_min, args), 'o')

ax.set_ylim(-10, 80)
ax.set_xlabel(r'$x$', fontsize=18)
ax.set_ylabel(r'$U(x)$', fontsize=18);

Eigenstates and eigenvalues

In [23]:
u = assemble_u_potential(N, U_ho, x, args)

K = assemble_K(N, k, x_min, x_max)

V = assemble_V(N, u, x_min, x_max)

H = K + V

evals, evecs = solve_eigenproblem(H)
In [27]:
NN = 10

fig, axes = subplots(NN, 1, figsize=(10,NN*1), sharex=True, sharey=True)

for n in range(NN):
    Y = evecs[NN-n-1]
    axes[n].plot(x, Y.real)

axes[n].set_xlim(-1, 1)
axes[n].set_xlabel(r'$x$', fontsize=18)
#axes[n].set_ylabel(r'$\Psi_n(x)$', fontsize=18);
fig.tight_layout();
In [31]:
fig, ax = subplots(figsize=(12,8))

ax.plot(x, U, 'k')
for n in range(10):
    Y = evals[n] + evecs[n]

    mask = where(Y > U)    
    ax.plot(x[mask], evals[n].real * ones(shape(x))[mask], 'k--')

    mask = where(Y > U-2.0)
    ax.plot(x[mask], Y[mask].real)
    
ax.set_xlim(-1, 1)
ax.set_ylim(0, 10)
ax.set_xlabel(r'$x$', fontsize=18)
ax.set_ylabel(r'$U(x)$', fontsize=18);

Identify bound states

In [43]:
boundidx = arange(10)
nbound = len(boundidx)
            
print "Found bound states: ", boundidx
Found bound states:  [0 1 2 3 4 5 6 7 8 9]

Evaluate transition matrix elements

In [44]:
inner_prod = zeros((len(boundidx), len(boundidx))).astype(np.complex)
expect_pos = zeros((len(boundidx), len(boundidx))).astype(np.complex)
expect_kin = zeros((len(boundidx), len(boundidx))).astype(np.complex)

for i, l in enumerate(boundidx):
    for j, k in enumerate(boundidx):
    
        psi_l = wavefunction_normalize(x, evecs[l])
        psi_k = wavefunction_normalize(x, evecs[k])
    
        inner_prod[i,j] = inner_product(x, psi_l, psi_k)                 # <psi_l|psi_k>
        expect_pos[i,j] = inner_product(x, psi_l, x * psi_k)             # <psi_l|x|psi_k>
        expect_kin[i,j] = inner_product(x, psi_l, derivative(x, psi_k))  # <psi_l|p|psi_k>
In [54]:
print "eigenenergies = " 

evals[0:nbound].real
eigenenergies = 
Out[54]:
array([ 0.50058072,  1.50156845,  2.50220863,  3.50250109,  4.50244565,
        5.50204213,  6.50129033,  7.50019009,  8.49874122,  9.49694352])
In [55]:
print("Harmonicity:")

((evals[boundidx] - evals[boundidx[0]]) / (evals[boundidx[1]]-evals[boundidx[0]])).real
Harmonicity:
Out[55]:
array([ 0.        ,  1.        ,  1.9996528 ,  2.99895823,  3.99791609,
        4.99652621,  5.99478841,  6.99270251,  7.99026831,  8.98748564])
In [56]:
print "orthogonality:\n"
print_matrix(inner_prod)
orthogonality:

 1.000    0.000   -0.000   -0.000   -0.000    0.000    0.000    0.000   -0.000   -0.000  
 0.000    1.000    0.000    0.000    0.000    0.000    0.000    0.000   -0.000   -0.000  
-0.000    0.000    1.000    0.000   -0.000    0.000    0.000   -0.000    0.000   -0.000  
-0.000    0.000    0.000    1.000    0.000   -0.000    0.000    0.000   -0.000    0.000  
-0.000    0.000   -0.000    0.000    1.000    0.000   -0.000   -0.000    0.000   -0.000  
 0.000    0.000    0.000   -0.000    0.000    1.000    0.000   -0.000   -0.000    0.000  
 0.000    0.000    0.000    0.000   -0.000    0.000    1.000   -0.000   -0.000   -0.000  
 0.000    0.000   -0.000    0.000   -0.000   -0.000   -0.000    1.000    0.000   -0.000  
-0.000   -0.000    0.000   -0.000    0.000   -0.000   -0.000    0.000    1.000    0.000  
-0.000   -0.000   -0.000    0.000   -0.000    0.000   -0.000   -0.000    0.000    1.000  
In [57]:
print "expect_pos:\n"
print_matrix(expect_pos)
expect_pos:

-0.000   -0.113   -0.000   -0.000   -0.000   -0.000    0.000   -0.000   -0.000   -0.000  
-0.113   -0.000   -0.159   -0.000   -0.000    0.000   -0.000   -0.000   -0.000    0.000  
-0.000   -0.159    0.000   -0.195   -0.000   -0.000   -0.000   -0.000    0.000   -0.000  
-0.000   -0.000   -0.195   -0.000   -0.225    0.000   -0.000    0.000   -0.000   -0.000  
-0.000   -0.000   -0.000   -0.225   -0.000   -0.252   -0.000   -0.000   -0.000   -0.000  
-0.000    0.000   -0.000    0.000   -0.252    0.000   -0.276   -0.000   -0.000    0.000  
 0.000   -0.000   -0.000   -0.000   -0.000   -0.276    0.000   -0.298    0.000   -0.000  
-0.000   -0.000   -0.000    0.000   -0.000   -0.000   -0.298   -0.000   -0.318   -0.000  
-0.000   -0.000    0.000   -0.000   -0.000   -0.000    0.000   -0.318    0.000   -0.337  
-0.000    0.000   -0.000   -0.000   -0.000    0.000   -0.000   -0.000   -0.337   -0.000  
In [58]:
print "expect_kin:\n"
print_matrix(expect_kin)
expect_kin:

 0.000   -4.438   -0.000   -0.003   -0.000   -0.000    0.000   -0.000   -0.000   -0.000  
 4.438   -0.000   -6.272   -0.000   -0.006    0.000   -0.000   -0.000   -0.000    0.000  
 0.000    6.272    0.000   -7.678    0.000   -0.009   -0.000   -0.000    0.000   -0.000  
 0.003    0.000    7.678   -0.000   -8.861   -0.000   -0.013    0.000   -0.000   -0.000  
 0.000    0.006   -0.000    8.861    0.000   -9.902   -0.000   -0.017   -0.000   -0.000  
 0.000   -0.000    0.009    0.000    9.902   -0.000   -10.842   -0.000   -0.021    0.000  
-0.000    0.000    0.000    0.013    0.000    10.842    0.000   -11.704    0.000   -0.026  
 0.000    0.000    0.000   -0.000    0.017    0.000    11.704   -0.000   -12.506   -0.000  
 0.000    0.000   -0.000    0.000    0.000    0.021   -0.000    12.506    0.000   -13.257  
 0.000   -0.000    0.000    0.000    0.000   -0.000    0.026    0.000    13.257    0.000