1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
|
using Plots
using LinearAlgebra
m = 1
hbar = 1
xmin = -6.5
xmax = 6.5
N = 2000
x = collect(range(xmin, stop = xmax, length = N + 1)) # one extra, want 1000 spaces
dx = x[2] - x[1]
println("dx = ", dx)
p = 0.0
V0 = 2 * p^2 / (2.0 * m)
sig = 0.25
x0 = 2.0 # a bit away to not interfere with the walls
# setup the potential barrier
# V = 0 .* x
# # println("V = ", length(V))
# for i in 1:length(V)
# V[i] = abs(x[i])^n
# end
# make the initial wavefunction
psi0 = exp.(-((x[2:end-1] .+ x0) .^ 2) / (2 * sig^2)) .* exp.(1.0im * p * x[2:end-1])
C = sum(abs2.(psi0) .* dx)
psi0 = psi0 / sqrt(C)
# plot(x[2:end-1], abs.(psi0), xlabel = "x", ylabel = "Re(psi)", title = "Initial wavefunction", legend = :none)
# # make the Hamiltonian
# A = (hbar^2 / (2 * m * dx^2))
# # main diagonal
# d = A * ones(N - 1) + V[2:end-1] # add potential
# # upper and lower diagonals
# dl = (A / -2.0) * ones(N - 2)
# du = (A / -2.0) * ones(N - 2)
# H = Tridiagonal(dl, d, du)
# # diagonalize matrix
# E, psi = eigen(H)
# psi1 = psi[:, 1]
# print the energy levels
# println("Energy levels: ", E[1:10])
function get_energy_levels_for_n(n, pot_plot)
V = 0 .* x
for i in 1:length(V)
V[i] = abs(x[i])^n
end
# plot the potential
println("Plotting potential for n = ", n)
plot!(pot_plot, x, V, xlabel = "x", ylabel = "V(x)", title = "Potential V(x) = |x|^n", label = "n = $n")
println("Calculating energy levels for n = ", n)
# make the Hamiltonian
A = (hbar^2 / (2 * m * dx^2))
# main diagonal
d = A * ones(N - 1) + V[2:end-1] # add potential
# upper and lower diagonals
dl = (A / -2.0) * ones(N - 2)
du = (A / -2.0) * ones(N - 2)
H = Tridiagonal(dl, d, du)
println("Diagonalizing matrix for n = ", n)
# diagonalize matrix
E, psi = eigen(H)
println("Energy levels for n = ", n, ": ", E[1:10])
return E
end
po = plot!()
energy_levels = []
n_max = 19
for n in 1:n_max
E = get_energy_levels_for_n(n, po)
push!(energy_levels, E)
end
savefig(po, "hw8/10-3-potentials.png")
# plot the energy levels
pl = plot(title = "Energy levels for n = 1 to $n_max", xlabel = "n", ylabel = "Energy")
for i in 1:5
plot!(pl, 1:n_max, [E[i] for E in energy_levels], label = "Energy Level $i", marker = :circle)
end
savefig(pl, "hw8/10-3-energy-levels.png")
|
