diff options
author | sotech117 <michael_foiani@brown.edu> | 2024-04-27 04:25:23 -0400 |
---|---|---|
committer | sotech117 <michael_foiani@brown.edu> | 2024-04-27 04:25:23 -0400 |
commit | e650ed1e1e908e51c78c1b047bec0da7c4fea366 (patch) | |
tree | 1fe238de7ca199b7fdee9bc29395080b3c4790e7 /hw8/BoundStates.jl | |
parent | 02756d17bca6f2b3bafa3f7b9fb6e5af438e94a0 (diff) |
testing
Diffstat (limited to 'hw8/BoundStates.jl')
-rw-r--r-- | hw8/BoundStates.jl | 73 |
1 files changed, 73 insertions, 0 deletions
diff --git a/hw8/BoundStates.jl b/hw8/BoundStates.jl new file mode 100644 index 0000000..22e608c --- /dev/null +++ b/hw8/BoundStates.jl @@ -0,0 +1,73 @@ +#!/usr/bin/env julia + +"""Find eigenstates and eigenenergies of 1-D quantum problems""" + +using LinearAlgebra +using Plots + +N = 1000 # number of lattice points +L = 20.0 # x runs from -L/2 to L/2 +dx = L / N + +D = zeros(N, N) # discrete laplacian operator +V = zeros(N, N) # potential + +for i in 1:N + D[i, i] = -2.0 +end + +for i in 1:N-1 + D[i, i+1] = 1.0 + D[i+1, i] = 1.0 +end + +#println("\nLattice Laplacian operator") +#println(D) + +function potential(x) + """ The potential energy""" + #return 0.0 # particle in a box + #return 0.5 * x^2 # SHO with the spring constant k = 1 + #return -6.0 * x^2 + 8.0 * x^6 # potential with zero ground state energy + #return 0.1 * x^4 - 2.0 * x^2 + 0.0 * x # double-well potential + return 8 * x^6 - 8 * x^4 - 4 * x^2 + 1 # another double-well potential +end + +for i in 1:N + x = (i + 0.5) * dx - 0.5 * L # coordinates of lattice points + V[i, i] = potential(x) +end + +H = -0.5 * dx^(-2.0) * D + V # Hamiltonian. Here m = hbar = 1 + +#println("\nMatrix elements of Hamiltonian = ") +#println(H) + +e, v = eigen(H) # diagonalize Hamiltonian + +println("\nGround state energy = ", e[1]) +println("\n1st excited state energy = ", e[2]) +println("\n2nd excited state energy = ", e[3]) +println("\n3rd excited state energy = ", e[4]) +println("\n4th excited state energy = ", e[5]) + + +gs(x) = exp(-0.5 * x^2) # Gaussian that is exact ground state of SHO + + +plot(potential) + + +plot(v[:, 1]) +#plot(v[:,2]) +#plot(gs) + +#= +eList = zeros(0) +for i in 1:20 + push!(eList, e[i]) +end + + +bar(eList, orientation = :horizontal) +=# |