Table of Contents
Done
ToBeDone
Start with \[ E\left[ \psi \right] = \frac{\int dX\psi^{*}(X)H\psi(X)}{\int dX\psi^{*}(X)\psi(X)} = \frac{ \left\langle \psi \right| H \left| \psi \right\rangle}{\left\langle \psi | \psi \right\rangle}\]
Through the variational method, we can obtain approximately the correct energy of a Hamiltonian using any basis wave fucntions that satisfy the boundary conditions.
\[ \mathbf{H} = \left\langle \psi_m \right| H \left| \psi_n \right\rangle \\ \mathbf{S} = \left\langle \psi_m \mid \psi_n \right\rangle \\ \mathbf{C} = \ \]
The generalized eigenvalue eqaution \[ \beta_k \mathbf{A}\mathbf{x}_K = \alpha_k \mathbf{B}\mathbf{x}_K \\ k = 0,1,2,...\]
\(\mathbf{HC} = E\mathbf{SC}\) where \(\mathbf{C}\) is the matrix of basis vectors
eye(N)
eig( M )
inv( M )
M[1:,i]
quadgk( f, a, b)
#this is just here to make sure that I can plot every time I restart the ipython notebook
using PyPlot
whos()
#globals used in the Hydrogen 4Basis
H4B_ALPHA = [13.00773, 1.962079, 0.444529, 0.1219492]
pie3_2 = sqrt( pi^3 )
#compose two functions
function myCompose(f,g)
#given two functions of 1 variable each
return (x,y)->(f(x)*g(y))
end
function compose(f,g)
return (x)->(f(x)*g(x))
end
# generates the nmth element of the S matrix from CP3.2
function S_inf(n,m)
if( (n+m)%2 == 0 )
return ( (2.0/(n+m+5)) - (4.0/(n+m+3)) + (2.0/(n+m+1)))
else
return 0
end
end
# generates the nmth element of the H matrix from CP3.2
function H_inf(n,m)
if((n+m)%2 == 0)
return ( -8.0 * ( (1-m-n-(2m*n)) / ((m+n+3)*(m+n+1)*(m+n-1)) ) )
else
return 0
end
end
function printMatElem( f, n, m)
for i in [1:n]
for j in [1:m]
temp = f(i,j)
print("$(temp), ")
end
println()
end
end
function fillMat(n,m,f)
mat = Array(Float64,n,m)
for i in [1:n]
for j in [1:m]
mat[i,j] = f((i-1),(j-1))
end
end
return mat
end
function infWellEig(n)
return eig( inv(fillMat(n,n,S_inf))*fillMat(n,n,H_inf) )
end
function infWellStatesandEig(n)
mat =
(vals, vecs) = eig(mat)
return (mat, vals, vecs)
end
# used to clobber terms in a matrix with NaN or less than the tolerance tol
function squishMatrix(mat; tol = 1e-12)
sz = size(mat)
for i = [1:sz[1]]
for j = [1:sz[2]]
if( (mat[i,j] < tol) || (mat[i,j] = NaN) )
mat[i,j] = 0
end
end
end
return matm2_dif = squishMat(m2)
end
function infWellEnergies(n)
(energies,vecs) = infWellEig(n)
return sort(energies)
end
function infWellBasisFuncs(x,n)
return (x^(n+2))-(x^n)
end
#draws the wave functions
function infWellWavenVec(N; steps = 1000, x_start = -1.0, x_stop = 1.0, eig_vecs=nothing, n = 0)
#has the user passed any eigenvectors?
if (eig_vecs == nothing)
ev = eye(N)
else
ev = eig_vecs
end
x = linspace(x_start,x_stop,steps)
wave = 0*Array(Float64,steps)
for i = [1:steps]
for j = [1:N]
bfj_x = infWellBasisFuncs(x[i], (j-1))
if (bfj_x != NaN) #if the basis function j at x is NaN, skip it or it will create garbage
wave[i] += ev[j, (n+1)] * bfj_x #j-1 because the eigenfunctions starts an N=0,1,2,..
end
end
end
return (x,wave)
end
# infWellPlotn
function infWellPlotn(N; n = 0, steps = 1000, x_start = -1.0, x_stop = 1.0, eig_vecs=nothing)
(x,wave) = infWellWavenVec(N,n=n,steps=steps, x_start=x_start, x_stop=x_stop, eig_vecs=eig_vecs)
plot(x,wave)
end
# quantumWaveKVec
# description:
# generates a vector of the probability density based on the projection into an eigen base
# N - the number of basis functions, or the order of the approx Matrix
# domain
# basisFunction is a function of two values, position and wave number
# projectionEigenVector
function quantumWaveKVec(N, domain, basisFunctions, projectionEigenVector)
sz = length(domain)
wave = Array(Float64,sz)
bfk_x = Array(Float64,N)
for i = [1:sz]
for k = [1:N]
bfk_tmp = basisFunctions(domain[i], (k-1))
if (bfk_tmp != NaN)
wave[i] += ev[j, (n+1)] * bfj_x #j-1 because the eigenfunctions starts an N=0,1,2,..
end
end
end
return wave
end
# print the energies and Plot the prob densities for an inifite square well with N basis functions
function infSqrWellApproxPlotv1(N; steps = 1000, width = 2.0, center = 0.0, showWhich=-1)
#step on generate the eigen vectors and values
(vals, vecs) = infWellEig(N)
#find the order of sorting
sortOrder = sortperm(vals)
# create the domain and range for the wave functions
domain = linspace( ((-width/2.0)-center), ((width/2.0)+center), steps)
projectionEnergy = 0*Array(Float64,N)
projectionRange = 0*Array(Float64,steps,N) # an Nx1000 matrix for simpelist use
for i in [1:N]
projectionEnergy[i] = vals[sortOrder[i]]
projectionRange[1:,i] = quantumWaveKVec(N,domain, infWellBasisFuncs, vecs[1:,sortOrder[i]])
end
if (showWhich == -2)
for i in [1:N]
plot(domain,projectionRange[1:,i])
end
elseif (showWhich == -1) #show the first five (default)
for i in [1:5]
plot(domain,projectionRange[1:,i], label="Wave $(i-1) with Energy $(round(projectionEnergy[i],4))")
end
else
plot(domain,projectionRange[1:,showWhich])
end
if ( sign(projectionRange[ int(steps/2), 1]) == -1)
legend(loc="lower right")
else
legend(loc="upper right")
end
end
# infWellBasisPlot plots the first 5 orders of the basis wave functions
function infWellBasisPlot(;steps = 1000)
domain = linspace(-1.0,1.0,steps)
pRange = 0*Array(Float64,steps,5)
for i = [1:steps]
for j = [1:5]
pRange[i,j] = infWellBasisFuncs( domain[i], (j-1) )
end
end
for j = [1:5]
plot(domain, pRange[1:,j], label="Eqn order $(j)")
end
legend()
end
# calc the prob densities for an inifite square well with N basis functions
function infSqrWellApprox(N; steps = 1000, width = 2.0, center = 0.0)
#step on generate the eigen vectors and values
(vals, vecs) = infWellEig(N)
#find the order of sorting
sortOrder = sortperm(vals)
# create the domain and range for the wave functions
domain = linspace( ((-width/2.0)-center), ((width/2.0)+center), steps)
projectionEnergy = 0*Array(Float64,N)
projectionRange = 0*Array(Float64,steps,N) # an Nx1000 matrix for simpelist use
for i in [1:N]
projectionEnergy[i] = vals[sortOrder[i]]
projectionRange[1:,i] = quantumWaveKVec(N,domain, infWellBasisFuncs, vecs[1:,sortOrder[i]])
end
return (domain, projectionRange, projectionEnergy)
end
# plot the infSqrWellWaveApprox
function infSqrWellApproxPlot(N)
(d, pR, pE) = infSqrWellApprox(N)
for i in [1:5]
plot(domain,projectionRange[1:,i], label="Wave $(i-1) with Energy $(round(projectionEnergy[i],4))")
end
mid = int(size(pR)[1] /2.0)
if ( sign(pR[ int(steps/2), 1]) == -1)
legend(loc="lower right")
else
legend(loc="upper right")
end
end
function hydrogen4Basis(r,a)
return exp( -1 * H4B_ALPHA[a] * r^2)
end
# the matrix element generation functions have their indicies increased by one to work with fillMat funciton
function hydrogen4BOverlap(p,q)
return sqrt( (pi / ( H4B_ALPHA[p+1] + H4B_ALPHA[q+1] ))^3 )
end
function hydrogen4BKinetic(p,q)
return (3.0 * (H4B_ALPHA[q+1]*H4B_ALPHA[p+1]*pie3_2) / sqrt( (H4B_ALPHA[p+1]+H4B_ALPHA[q+1])^5 ))
end
function hydrogen4BCoulomb(p,q)
return (-2pi / (H4B_ALPHA[p+1]+H4B_ALPHA[q+1]))
end
function hydrogen4BHamiltonian(p,q)
return hydrogen4BKinetic(p,q) + hydrogen4BCoulomb(p,q)
end
#return the hamilton and the overlap matricies
function hydrogen_4BH(N)
return fillMat(N,N,hydrogen4BHamiltonian)
end
function hydrogen_4BS(N)
return fillMat(N,N,hydrogen4BOverlap)
end
function hydrogen_4BEigen(;N=4)
return eig( inv(hydrogen_4BS(N))*hydrogen_4BH(N) )
end
ffo = myCompose(x->x, x->x)
println(ffo(2,2))
println(ffo(2,4))
println(quadgk(sin,0,pi))
println(quadgk(compose(sin,cos),0,pi))
# test that the H_inf returns values for H matrix
printMatElem( H_inf, 5, 5)
# test that the S_inf returns values for H matrix
printMatElem( S_inf, 5, 5)
# test fillMat for the S matrix for the infitite sqaure well
S_temp = fillMat(5,5,S_inf)
# test fillMat for the H matrix for the infitite sqaure well
H_temp = fillMat(5,5,H_inf)
# do the eigen values match those in the book in table 3.1
eig(inv(S_temp)*H_temp )
# YES THEY DO! but they are out of order
(m1,v1,e1) = infWellStatesandEig(5)
infWellEnergies(8)
infWellEnergies(12)[1:5]
infWellEnergies(16)[1:5]
(p,q) = infWellEig(8)
# test ploting and
x_tmp = linspace(-1.0,1.0,100)
w_tmp = 0*Array(Float64,100)
ev_tmp = eye(5)
for i in [1:100]
for j in [1:5]
bf = infWellBasisFuncs(x_tmp[i],(j-1))
if( bf != NaN )
w_tmp[i] += (ev_tmp[j,1] * bf)
end
end
end
plot(x_tmp,w_tmp)
infWellWavenVec(5, n=0)
infWellPlotn(5)
infSqrWellApproxPlotv1(5)
infSqrWellApproxPlotv1(8)
infSqrWellApproxPlotv1(12)
infSqrWellApproxPlotv1(16)
infSqrWellApproxPlotv1(20)
20 is the last value before the the values start becomming imaginary.
infWellBasisPlot()
hfoo = hydrogen_4BEigen()
dot(hfoo[2][1:,1],hfoo[2][1:,2])
heig = hfoo[2]
hS = hydrogen_4BS(4)
hEV = eigvecs(hS)
hs = squishMat(conj(transpose(hEV))*hS*hEV)
hV = hEV* hs^(-1/2)
hId = squishMat( transpose(hV)*hS*hV )
hH = hydrogen_4BH(4)
# hEO_1 is the value of Energies and Waves ("orbitals") of hydrogen atom using the method in sec 3.3
hEO_1 = eig( transpose(hV)*hH*hV )
# this is my method for finding the Enegies and Orbitals by finding the eigens of the inverse of S multiplied by H
hEO_2 = eig( inv(hydrogen_4BS(4)) * hydrogen_4BH(4) )
#now find the difference between the two methods
hDiff = hEO_2[2]-hEO_1[2]
cos(-pi/4)^2