10 files changed, 492 insertions, 0 deletions
diff --git a/hw8/10-14-old.png b/hw8/10-14-old.png
new file mode 100644
index 0000000..c1d5618
--- /dev/null
+++ b/hw8/10-14-old.png
diff --git a/hw8/10-14.jl b/hw8/10-14.jl
new file mode 100644
index 0000000..27d1b5c
--- /dev/null
+++ b/hw8/10-14.jl
@@ -0,0 +1,82 @@
+#!/usr/bin/env julia
+
+using Plots
+using LinearAlgebra
+using DifferentialEquations
+
+# useful functions:
+
+function H(psi, dx) # action of Hamiltonian on wavefunction
+	Hpsi = zeros(ComplexF64, size(psi))
+	# -(1/2) * laplacian(psi) (m = hbar = 1)
+	Hpsi[2:end-1] = 0.5 * (2.0 * psi[2:end-1] - psi[3:end] - psi[1:end-2]) / (dx * dx)
+
+	# periodic boundary conditions
+	#Hpsi[1] = 0.5 * (2.0*psi[1] - psi[2] - psi[end])/(dx*dx)
+	#Hpsi[end] = 0.5 * (2.0*psi[end] - psi[end-1] - psi[1])/(dx*dx)
+	return Hpsi
+end
+
+derivative(psi, dx) = -1.0im * H(psi, dx)
+
+function initialWavefunction(x::Vector{Float64}, x0 = 10.0, Delta = 1.0, k = 4.0)
+	Delta2 = Delta^2
+	return exp.(-(x .- x0) .^ 2 / Delta2) .* exp.(1.0im * k * x)
+end
+
+function normalization(psi, dx) # normalization of wavefunction
+	n2 = dot(psi, psi) * dx
+	return sqrt(n2)
+end
+
+prob(psi) = real(psi .* conj(psi))
+
+# The actual simulation
+N = 400  # number of lattice points
+L = 20.0 # x runs from 0 to L
+dx = L / N
+
+x = range(0.0, L, N) |> collect # lattice along x-axis
+#println(x)
+
+# initial wavefunction has position (x0), width (Delta), and momentum (k)
+psi0 = initialWavefunction(x, 10.0, 1.0, 1000.0)
+
+# normalize wavefunction
+n = normalization(psi0, dx)
+psi0 = psi0 / n
+println("norm = ", normalization(psi0, dx))
+
+# plot initial wavefunction and probability density
+plot(prob(psi0))
+
+# integrate forward in time
+tf = 1.0
+dt = 0.1
+tspan = (0.0, tf)
+
+function timeEvolve(psi0, tf, dt) # second order Runge-Kutta algorithm
+	psi = psi0
+	for t in range(0, stop = tf, step = dt)
+		psiMid = psi + 0.5 * dt * derivative(psi, dx)
+		psi = psi + dt * derivative(psiMid, dx)
+	end
+	return psi
+end
+
+# tendency(psi, dx, t) = derivative(psi, dx) # use ODE solver
+# problem = ODEProblem(tendency, psi0, tspan, dx) # specify ODE
+# sol = solve(problem, Tsit5(), reltol = 1e-12, abstol = 1e-12) # solve using Tsit5 algorithm to specified accuracy
+
+# compare initial and final wavefunction probabilities
+psi = timeEvolve(psi0, tf, dt)
+# psi = sol[:, end]
+times = sol.t
+
+# check that normalization hasn't deviated too far from 1.0
+println("norm = ", normalization(psi, dx))
+
+plot([prob(psi0), prob(psi)])
+savefig("10-14.png")
+
+
diff --git a/hw8/10-14.png b/hw8/10-14.png
new file mode 100644
index 0000000..3e18e5b
--- /dev/null
+++ b/hw8/10-14.png
diff --git a/hw8/10-14t.png b/hw8/10-14t.png
new file mode 100644
index 0000000..3e18e5b
--- /dev/null
+++ b/hw8/10-14t.png
diff --git a/hw8/10-17.jl b/hw8/10-17.jl
new file mode 100644
index 0000000..6dfff2b
--- /dev/null
+++ b/hw8/10-17.jl
@@ -0,0 +1,81 @@
+#!/usr/bin/env julia
+
+using Plots
+using LinearAlgebra
+using DifferentialEquations
+
+# useful functions:
+
+function H(psi, dx) # action of Hamiltonian on wavefunction
+	Hpsi = zeros(ComplexF64, size(psi))
+	# -(1/2) * laplacian(psi) (m = hbar = 1)
+	Hpsi[2:end-1] = 0.5 * (2.0 * psi[2:end-1] - psi[3:end] - psi[1:end-2]) / (dx * dx)
+
+	# periodic boundary conditions
+	#Hpsi[1] = 0.5 * (2.0*psi[1] - psi[2] - psi[end])/(dx*dx)
+	#Hpsi[end] = 0.5 * (2.0*psi[end] - psi[end-1] - psi[1])/(dx*dx)
+	return Hpsi
+end
+
+derivative(psi, dx) = -1.0im * H(psi, dx)
+
+function initialWavefunction(x::Vector{Float64}, x0 = 10.0, Delta = 1.0, k = 4.0)
+	Delta2 = Delta^2
+	return exp.(-(x .- x0) .^ 2 / Delta2) .* exp.(1.0im * k * x)
+end
+
+function normalization(psi, dx) # normalization of wavefunction
+	n2 = dot(psi, psi) * dx
+	return sqrt(n2)
+end
+
+prob(psi) = real(psi .* conj(psi))
+
+# The actual simulation
+N = 400  # number of lattice points
+L = 20.0 # x runs from 0 to L
+dx = L / N
+
+x = range(0.0, L, N) |> collect # lattice along x-axis
+#println(x)
+
+# initial wavefunction has position (x0), width (Delta), and momentum (k)
+psi0 = initialWavefunction(x, 10.0, 0.05, 700.0)
+
+# normalize wavefunction
+n = normalization(psi0, dx)
+psi0 = psi0 / n
+println("norm = ", normalization(psi0, dx))
+
+# plot initial wavefunction and probability density
+plot(prob(psi0))
+
+# integrate forward in time
+tf = 1.0
+dt = 5e-7
+tspan = (0.0, tf)
+
+function timeEvolve(psi0, tf, dt) # second order Runge-Kutta algorithm
+	psi = psi0
+	for t in range(0, stop = tf, step = dt)
+		psiMid = psi + 0.5 * dt * derivative(psi, dx)
+		psi = psi + dt * derivative(psiMid, dx)
+	end
+	return psi
+end
+
+tendency(psi, dx, t) = derivative(psi, dx) # use ODE solver
+problem = ODEProblem(tendency, psi0, tspan, dx) # specify ODE
+sol = solve(problem, Tsit5(), reltol = 1e-12, abstol = 1e-12) # solve using Tsit5 algorithm to specified accuracy
+
+# compare initial and final wavefunction probabilities
+#psi = timeEvolve(psi0, tf, dt)
+psi = sol[:, end]
+times = sol.t
+
+# check that normalization hasn't deviated too far from 1.0
+println("norm = ", normalization(psi, dx))
+
+plot([prob(psi0), prob(psi)])
+
+
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)
+=#
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)
+=#
diff --git a/hw8/RadialBoundStates.jl b/hw8/RadialBoundStates.jl
new file mode 100644
index 0000000..ed7bcc9
--- /dev/null
+++ b/hw8/RadialBoundStates.jl
@@ -0,0 +1,59 @@
+#!/usr/bin/env julia
+
+"""Find eigenstates and eigenenergies of central potential problems"""
+
+using LinearAlgebra
+using Plots
+
+N = 5000  # number of lattice points
+L = 20.0 # r runs from 0 to L
+dr = L / N
+
+D = zeros(N, N) # discrete radial 2nd derivative 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(r, ℓ = 0)
+	""" The potential energy"""
+	#return 0.5 * ell * (ℓ+1.0) * pow(r, -2.0) # V=0: Free particle in spherical coordinates	
+
+	return -1.0 / r + 0.5 * ℓ * (ℓ + 1.0) * r^(-2.0) # Hydrogen atom
+
+	#return -r^(-1.1) + 0.5 * ℓ * (ℓ+1.0) * r^(-2.0) # modified Coulomb potential
+end
+
+for i in 1:N
+	r = (i + 0.5) * dr # radial coordinates of lattice points
+	V[i, i] = potential(r, 0)
+end
+
+H = -0.5 * dr^(-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])
+
+
+plot(potential)
+
+
+plot(v[:, 1])
+#plot(v[:,2])
diff --git a/hw8/TimeDependentSchrodingerEquation.jl b/hw8/TimeDependentSchrodingerEquation.jl
new file mode 100644
index 0000000..dacbced
--- /dev/null
+++ b/hw8/TimeDependentSchrodingerEquation.jl
@@ -0,0 +1,81 @@
+#!/usr/bin/env julia
+
+using Plots
+using LinearAlgebra
+using DifferentialEquations
+
+# useful functions:
+
+function H(psi, dx) # action of Hamiltonian on wavefunction
+    Hpsi = zeros(ComplexF64, size(psi))
+    # -(1/2) * laplacian(psi) (m = hbar = 1)
+    Hpsi[2:end-1] = 0.5 * (2.0*psi[2:end-1] - psi[3:end] - psi[1:end-2])/(dx*dx)
+
+    # periodic boundary conditions
+    #Hpsi[1] = 0.5 * (2.0*psi[1] - psi[2] - psi[end])/(dx*dx)
+    #Hpsi[end] = 0.5 * (2.0*psi[end] - psi[end-1] - psi[1])/(dx*dx)
+    return Hpsi
+end
+
+derivative(psi, dx) = -1.0im * H(psi, dx)
+
+function initialWavefunction(x::Vector{Float64}, x0 = 10.0, Delta = 1.0, k = 4.0)
+    Delta2 = Delta^2
+    return exp.(- (x .- x0).^2 / Delta2 ) .* exp.(1.0im * k * x)
+end
+
+function normalization(psi, dx) # normalization of wavefunction
+    n2 = dot(psi, psi) * dx
+    return sqrt(n2)
+end
+
+prob(psi) = real(psi .* conj(psi))
+
+# The actual simulation
+N = 400  # number of lattice points
+L = 20.0 # x runs from 0 to L
+dx = L / N
+
+x = range(0.0, L, N) |> collect # lattice along x-axis
+#println(x)
+
+# initial wavefunction has position (x0), width (Delta), and momentum (k)
+psi0 = initialWavefunction(x, 10.0, 1.0, 0.0)
+
+# normalize wavefunction
+n = normalization(psi0, dx)
+psi0 = psi0 / n
+println("norm = ", normalization(psi0, dx))
+
+# plot initial wavefunction and probability density
+plot(prob(psi0))
+
+# integrate forward in time
+tf = 1.0
+dt = 0.0001
+tspan = (0.0, tf)
+
+function timeEvolve(psi0, tf, dt) # second order Runge-Kutta algorithm
+    psi = psi0
+    for t in range(0, stop=tf, step=dt)
+        psiMid = psi + 0.5 * dt * derivative(psi, dx)
+        psi = psi + dt * derivative(psiMid, dx)
+    end
+    return psi
+end
+
+tendency(psi, dx, t) = derivative(psi, dx) # use ODE solver
+problem = ODEProblem(tendency, psi0, tspan, dx) # specify ODE
+sol = solve(problem, Tsit5(), reltol=1e-12, abstol=1e-12) # solve using Tsit5 algorithm to specified accuracy
+
+# compare initial and final wavefunction probabilities
+#psi = timeEvolve(psi0, tf, dt)
+psi = sol[:, end]
+times = sol.t
+
+# check that normalization hasn't deviated too far from 1.0
+println("norm = ", normalization(psi, dx))
+
+plot([prob(psi0), prob(psi)])
+
+
diff --git a/hw8/tightbinding.jl b/hw8/tightbinding.jl
new file mode 100644
index 0000000..0d1cc8c
--- /dev/null
+++ b/hw8/tightbinding.jl
@@ -0,0 +1,24 @@
+using LinearAlgebra
+
+# Define the tight binding Hamiltonian
+function tight_binding_hamiltonian(n_sites::Int64, t::Float64)
+    H = zeros(n_sites, n_sites)
+    for i in 1:(n_sites-1)
+        H[i, i+1] = t
+        H[i+1, i] = t
+    end
+    H
+end
+
+# Find the eigenenergies of the Hamiltonian
+function eigenenergies(H::Array{Float64, 2})
+    eigvals(H)
+end
+
+n_sites = 10 # number of sites in the chain
+t = 1.0 # hopping parameter
+H = tight_binding_hamiltonian(n_sites, t)
+energies = eigenenergies(H)
+println(energies)
+
+