In []:

In [28]:

HERMITE6 = [ (x)->(1), # Hermite_0 polynomial
             (x)->(2x), 
             (x)->(4x^2-2),
             (x)->(8x^3 - 12x),
             (x)->(16x^4 - 48x^2 + 12),
             (x)->(32x^5 - 160 - 120x) ]
SQRPI = 1.7724538509055159 # square root of pi

Out[28]:

1.7724538509055159

In [4]:

function getHermite(n)
    return HERMITE6[n+1] # because Jl start indicies at 1
end

Out[4]:

getHermite (generic function with 1 method)

In [5]:

function computeHermite(x,n)
    f = getHermite(n)
    return f(x)
end

Out[5]:

computeHermite (generic function with 1 method)

In [67]:

function qho_WaveFact(n)
    return (1 / sqrt( 2^n * factorial(n) ))
end

function qho_WaveNatX(x,n)
    norm = 0.10627834566699294
    WF = qho_WaveFact(n)
    H_n = computeHermite(x,n)
    
    return ( norm* WF * H_n * exp( -0.5 * x^2 ) )
end

Out[67]:

qho_WaveNatX (generic function with 1 method)

In [15]:

function gaussianPert(x;alpha=4)
    return ( (1/alpha)*exp(-0.5x^2) )
end

Out[15]:

gaussianPert (generic function with 1 method)

In [42]:

function qho_EnergyN(n)
    norm4 = 0.10627834566699294^4
    hBar = 1.05457e-34
    return ( (n + 0.5) * norm4 * hBar^2 * pi )
end

Out[42]:

qho_EnergyN (generic function with 1 method)

In [69]:

wsp = linspace(-10,10,1000)
psi1 = zeros(1000)
phi1 = zeros(1000)
psi2 = zeros(1000)
phi2 = zeros(1000)
gPert1 = zeros(1000)
for i in [1:1000]
    psi1[i] = qho_WaveNatX( wsp[i], 0 )
    phi1[i] = qho_WaveNatX( wsp[i], 0 )
    psi2[i] = qho_WaveNatX( wsp[i], 1 )
    phi2[i] = qho_WaveNatX( wsp[i], 1 )
    gPert1[i] = gaussianPert( wsp[i] )
end

println(dot( psi1, phi1 ))
println(dot( psi2, phi2 ))
println(dot( psi1, phi2 ))
println(dot( psi2, phi1 ))

1.0
0.9999999999999998
-2.035545151696134e-17
-2.035545151696134e-17

In [40]:

for i in [1:1000]
    phi1[i] = gPert1[i] * phi1[i]
end
dot(phi1, psi1)

Out[40]:

0.2041241452319315

In [37]:

1/sqrt(88.53406985273054)

Out[37]:

0.10627834566699294

In [43]:

qho_EnergyN(0)

Out[43]:

2.2286909229549866e-72

In [44]:

16 * 4

Out[44]:

In [45]:

64 - 3 - 3

Out[45]:

In [46]:

58 / 5

Out[46]:

11.6

In [47]:

58 / 7

Out[47]:

8.285714285714286

In [55]:

acos(2/3)*180/pi

Out[55]:

48.189685104221404

In [91]:

function fillWaveMat(n, N, dVec)
    m = zeros(n,N)
    for i = [1:N]
        for j = [1:n]
            m[j,i] = qho_WaveNatX( dVec[j], i-1 )
        end
    end
    
    return m
end

Out[91]:

fillWaveMat (generic function with 2 methods)

In [90]:

function fillPsiNCol(n, N; left=-4, right=4)
    return fillWaveMat(n, N, linspace(left,right,n) )
end

Out[90]:

fillPsiNCol (generic function with 2 methods)

In [123]:

function fillPsiNRow(n, N; left=-4, right=4)
    return transpose( fillPsiNCol(n, N) )
end

Out[123]:

fillPsiNRow (generic function with 2 methods)

In [92]:

p_mat = fillPsiNCol( 100, 4 )

Out[92]:

100x4 Array{Float64,2}:
 3.56524e-5   -0.000201681  0.000781512  -0.00238774
 4.90963e-5   -0.00027212   0.00103177   -0.0030795 
 6.71697e-5   -0.000364617  0.00135204   -0.00393963
 9.1298e-5    -0.000485159  0.00175846   -0.00499892
 0.000123286  -0.000641054  0.00226983   -0.00629077
 0.000165398  -0.000841124  0.00290769   -0.00785047
 0.00022045   -0.0010959    0.00369636   -0.00971416
 0.000291914  -0.0014178    0.00466278   -0.0119174 
 0.000384029  -0.0018213    0.00583624   -0.0144934 
 0.000501922  -0.00232306   0.00724784   -0.0174707 
 0.000651738  -0.00294198   0.00892972   -0.0208704 
 0.000840764  -0.00369917   0.010914     -0.0247036 
 0.00107755   -0.00461785   0.0132315    -0.0289674 
 ⋮                                                  
 0.000840764   0.00369917   0.010914      0.0247036 
 0.000651738   0.00294198   0.00892972    0.0208704 
 0.000501922   0.00232306   0.00724784    0.0174707 
 0.000384029   0.0018213    0.00583624    0.0144934 
 0.000291914   0.0014178    0.00466278    0.0119174 
 0.00022045    0.0010959    0.00369636    0.00971416
 0.000165398   0.000841124  0.00290769    0.00785047
 0.000123286   0.000641054  0.00226983    0.00629077
 9.1298e-5     0.000485159  0.00175846    0.00499892
 6.71697e-5    0.000364617  0.00135204    0.00393963
 4.90963e-5    0.00027212   0.00103177    0.0030795 
 3.56524e-5    0.000201681  0.000781512   0.00238774

In [97]:

p_tmat = fillPsiNRow( 100, 4)

Out[97]:

4x100 Array{Float64,2}:
  3.56524e-5    4.90963e-5   6.71697e-5   …  4.90963e-5  3.56524e-5 
 -0.000201681  -0.00027212  -0.000364617     0.00027212  0.000201681
  0.000781512   0.00103177   0.00135204      0.00103177  0.000781512
 -0.00238774   -0.0030795   -0.00393963      0.0030795   0.00238774

In [117]:

function qho_HKet( h, ket )
    for i = [1:size(ket)[1]]
        for j = [1:size(ket)[2]]
               ket[i,j] = ket[i,j] * h[i]
        end
    end 
    
    return ket
end

Out[117]:

qho_HKet (generic function with 1 method)

In [118]:

function qho_PertKet()
    ket = fillPsiNCol( 100, 4)
    h = linspace(-4,4,100)
    for i = [1:100]
        h[i] = gaussianPert(h[i])
    end
    return qho_HKet( h, ket )
end

Out[118]:

qho_PertKet (generic function with 2 methods)

In [119]:

qho_PertKet()

Out[119]:

100x4 Array{Float64,2}:
 2.99001e-9  -1.69141e-8  6.5542e-8   -2.00249e-7
 5.67013e-9  -3.14271e-8  1.1916e-7   -3.55651e-7
 1.06131e-8  -5.76109e-8  2.13628e-7  -6.22478e-7
 1.96073e-8  -1.04194e-7  3.77651e-7  -1.07358e-6
 3.57538e-8  -1.8591e-7   6.58268e-7  -1.82437e-6
 6.4351e-8   -3.27254e-7  1.13129e-6  -3.05436e-6
 1.14318e-7  -5.68297e-7  1.91681e-6  -5.03745e-6
 2.0045e-7   -9.73563e-7  3.20181e-6  -8.18338e-6
 3.46915e-7  -1.64528e-6  5.27221e-6  -1.30927e-5
 5.92608e-7  -2.74279e-6  8.55737e-6  -2.06273e-5
 9.99175e-7  -4.51033e-6  1.36901e-5  -3.19963e-5
 1.66281e-6  -7.316e-6    2.15851e-5  -4.88572e-5
 2.73132e-6  -1.17051e-5  3.35386e-5  -7.34251e-5
 ⋮                                               
 1.66281e-6   7.316e-6    2.15851e-5   4.88572e-5
 9.99175e-7   4.51033e-6  1.36901e-5   3.19963e-5
 5.92608e-7   2.74279e-6  8.55737e-6   2.06273e-5
 3.46915e-7   1.64528e-6  5.27221e-6   1.30927e-5
 2.0045e-7    9.73563e-7  3.20181e-6   8.18338e-6
 1.14318e-7   5.68297e-7  1.91681e-6   5.03745e-6
 6.4351e-8    3.27254e-7  1.13129e-6   3.05436e-6
 3.57538e-8   1.8591e-7   6.58268e-7   1.82437e-6
 1.96073e-8   1.04194e-7  3.77651e-7   1.07358e-6
 1.06131e-8   5.76109e-8  2.13628e-7   6.22478e-7
 5.67013e-9   3.14271e-8  1.1916e-7    3.55651e-7
 2.99001e-9   1.69141e-8  6.5542e-8    2.00249e-7

In [124]:

HT = fillPsiNRow(100,4)*qho_PertKet()

Out[124]:

4x4 Array{Float64,2}:
  0.0505713     1.09804e-18  -0.0119198    -5.25389e-20
  6.74619e-19   0.0337142    -9.90611e-20  -0.0137638  
 -0.0119198     2.2551e-19    0.0252856    -5.20546e-19
  2.32798e-19  -0.0137638    -9.27017e-20   0.0206031

In [127]:

hte = eigvecs(HT)

Out[127]:

4x4 Array{Float64,2}:
  0.92941       0.369048     -2.44289e-17  -9.83259e-17
  0.0           0.0           0.533851      0.845579   
 -0.369048      0.92941      -5.261e-18     5.804e-17  
  7.20767e-18  -1.81518e-17   0.845579     -0.533851

Warning: Possible conflict in library symbol dsyevr_
Warning: Possible conflict in library symbol dgeev_

In [139]:

qho_e = transpose(hte)*HT*hte

Out[139]:

4x4 Array{Float64,2}:
  0.0553044     6.93889e-18  -4.54135e-19  -5.68567e-18
  8.67362e-18   0.0205526    -5.65905e-19   1.56229e-18
 -5.09548e-19  -3.85088e-19   0.0119135     8.67362e-19
 -5.97443e-18   9.06585e-19   3.46945e-18   0.0424039

In [136]:

# used to clobber terms in a matrix with NaN
function squishMat(mat; tol = 1e-16)
    sz = size(mat)
    for i = [1:sz[1]]
        for j = [1:sz[2]]
            if ((mat[i,j] == NaN) || (mat[i,j] <= tol))
                mat[i,j] = 0
            end
        end
    end
    return mat
end

Out[136]:

squishMat (generic function with 1 method)

In [140]:

squishMat(qho_e)

Out[140]:

4x4 Array{Float64,2}:
 0.0553044  0.0        0.0        0.0      
 0.0        0.0205526  0.0        0.0      
 0.0        0.0        0.0119135  0.0      
 0.0        0.0        0.0        0.0424039

In [164]:

function diffMatN(N)
    m = zeros(N,N)
    m[1,1] = -1
    m[1,2] = 1
    
    for i = [2:(N-1)]
        m[i,(i-1)] = -1
        m[i,i] = 0
        m[i,(i+1)] = 1
    end
    
    m[N,(N-1)] = -1
    m[N,N] = 1

    return m
end

Out[164]:

diffMatN (generic function with 2 methods)

In [165]:

d1 = diffMatN(10)

Out[165]:

10x10 Array{Float64,2}:
 -1.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  0.0
 -1.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0  0.0
  0.0  -1.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0  -1.0   0.0   1.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0  -1.0   0.0   1.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0  -1.0   0.0   1.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0  -1.0   0.0   1.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0  -1.0   0.0   1.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0   0.0  -1.0   0.0  1.0
  0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -1.0  1.0

In [160]:

d1*d1*d1*d1

Out[160]:

10x10 Array{Float64,2}:
  0.0      2.0e-8  -2.0e-8  -1.0e-8  …   0.0      0.0      0.0      0.0   
 -2.0e-8   4.0e-8   1.0e-8  -4.0e-8      0.0      0.0      0.0      0.0   
 -2.0e-8  -1.0e-8   6.0e-8   0.0         1.0e-8   0.0      0.0      0.0   
  1.0e-8  -4.0e-8   0.0      6.0e-8      0.0      1.0e-8   0.0      0.0   
  1.0e-8   0.0     -4.0e-8   0.0        -4.0e-8   0.0      1.0e-8   0.0   
  0.0      1.0e-8   0.0     -4.0e-8  …   0.0     -4.0e-8   0.0      1.0e-8
  0.0      0.0      1.0e-8   0.0         6.0e-8   0.0     -4.0e-8   1.0e-8
  0.0      0.0      0.0      1.0e-8      0.0      6.0e-8  -1.0e-8  -2.0e-8
  0.0      0.0      0.0      0.0        -4.0e-8   1.0e-8   4.0e-8  -2.0e-8
  0.0      0.0      0.0      0.0        -1.0e-8  -2.0e-8   2.0e-8   0.0

In []: