diff options
Diffstat (limited to 'hw8/10-3.jl')
-rw-r--r-- | hw8/10-3.jl | 92 |
1 files changed, 92 insertions, 0 deletions
diff --git a/hw8/10-3.jl b/hw8/10-3.jl new file mode 100644 index 0000000..8616a69 --- /dev/null +++ b/hw8/10-3.jl @@ -0,0 +1,92 @@ +#!/usr/bin/env julia + +"""Find eigenstates and eigenenergies of 1-D quantum problems""" + +using LinearAlgebra +using Plots + +#println("\nLattice Laplacian operator") +#println(D) + +function potential(x, n) + """ 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 + return abs(x)^n +end + +#println("\nMatrix elements of Hamiltonian = ") +#println(H) + + +function map_n_to_energies(n) + 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 + + + for i in 1:N + x = (i + 0.5) * dx - 0.5 * L # coordinates of lattice points + V[i, i] = potential(x, n) + end + + H = -0.5 * dx^(-2.0) * D + V # Hamiltonian. Here m = hbar = 1 + + + e, v = eigen(H) # diagonalize Hamiltonian + + println("\n n = ", n) + println("Ground state energy = ", e[1]) + println("1st excited state energy = ", e[2]) + println("2nd excited state energy = ", e[3]) + println("3rd excited state energy = ", e[4]) + println("4th excited state energy = ", e[5], "\n") + + return e +end + +n_max = 18 +n_to_e = [map_n_to_energies(n) for n in 1:n_max] + +# plot e[0] for all N +eList = zeros(0) +for i in 1:n_max + push!(eList, n_to_e[i][1]) +end + +plot(eList) + +# 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) +=# |