add hw9 and hw8

21 files changed, 152 insertions, 57 deletions

diff --git a/hw8/10-14-0.0001.png b/hw8/10-14-0.0001.png
new file mode 100644
index 0000000..0dc656a
--- /dev/null
+++ b/hw8/10-14-0.0001.png
diff --git a/hw8/10-14-0.001.png b/hw8/10-14-0.001.png
new file mode 100644
index 0000000..fcf13d4
--- /dev/null
+++ b/hw8/10-14-0.001.png
diff --git a/hw8/10-14-0.01.png b/hw8/10-14-0.01.png
new file mode 100644
index 0000000..8c6d987
--- /dev/null
+++ b/hw8/10-14-0.01.png
diff --git a/hw8/10-14-0.1.png b/hw8/10-14-0.1.png
new file mode 100644
index 0000000..74f5045
--- /dev/null
+++ b/hw8/10-14-0.1.png
diff --git a/hw8/10-14-E.png b/hw8/10-14-E.png
new file mode 100644
index 0000000..d13b721
--- /dev/null
+++ b/hw8/10-14-E.png
diff --git a/hw8/10-14-psi.mp4 b/hw8/10-14-psi.mp4
new file mode 100644
index 0000000..37c0d10
--- /dev/null
+++ b/hw8/10-14-psi.mp4
diff --git a/hw8/10-14-psi.png b/hw8/10-14-psi.png
new file mode 100644
index 0000000..7911c89
--- /dev/null
+++ b/hw8/10-14-psi.png
diff --git a/hw8/10-14.jl b/hw8/10-14.jl
index 27d1b5c..8a6307a 100644
--- a/hw8/10-14.jl
+++ b/hw8/10-14.jl
@@ -52,9 +52,12 @@ plot(prob(psi0))
 
 # integrate forward in time
 tf = 1.0
-dt = 0.1
+dt = 0.0007
 tspan = (0.0, tf)
 
+println("dx = ", dx)
+println("dt = ", dt)
+
 function timeEvolve(psi0, tf, dt) # second order Runge-Kutta algorithm
 	psi = psi0
 	for t in range(0, stop = tf, step = dt)
@@ -64,19 +67,56 @@ function timeEvolve(psi0, tf, dt) # second order Runge-Kutta algorithm
 	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
+# get energy of psi
+function energy(psi, dx)
+	e_complex = dot(psi, H(psi, dx))
+	return real(e_complex)
+end
 
 # compare initial and final wavefunction probabilities
-psi = timeEvolve(psi0, tf, dt)
-# psi = sol[:, end]
-times = sol.t
-
+psi_1 = timeEvolve(psi0, tf, dt)
+psi_2 = timeEvolve(psi0, tf, 0.00088)
+psi_3 = timeEvolve(psi0, tf, 0.00089)
+psi_4 = timeEvolve(psi0, tf, 0.0009)
 # 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")
-
-
+println("norm euler dt=.0007 ->", normalization(psi_1, dx))
+println("norm euler dt=.00088 ->", normalization(psi_2, dx))
+println("norm euler dt=.00089 -> ", normalization(psi_3, dx))
+println("norm euler dt=.0009 -> ", normalization(psi_4, dx))
+
+
+tendency(psi, dx, t) = derivative(psi, dx) # use ODE solver
+problem = ODEProblem(tendency, psi0, tspan, dx) # specify ODE
+sol = solve(problem, Tsit5(), reltol = 1e-13, abstol = 1e-13) # solve using Tsit5 algorithm to specified accuracy
+psi_5 = sol[:, end]
+times_2 = sol.t
+println("norm ode :", normalization(psi_5, dx))
+
+# cluclate the energy at each timestep
+# energies_over_time = [energy(sol[:, i], dx) for i in 1:length(times_2)]
+# plot(times_2, energies_over_time, title = "Energy of Wave Packet over Time", xlabel = "t", ylabel = "E", lw = 1.5, ylim = (85, 90), legend = false)
+# savefig("hw8/10-14-E.png")
+
+# make an animation of the wave packet
+# anim = @animate for i in 1:length(times_2)
+# 	plot(x, prob(sol[:, i]), title = "Propagation of Wave Packet", xlabel = "x", ylabel = "|ψ(x)|²", lw = 1.5, ylim = (0, 1), label = "t = $(round(times_2[i], digits = 3))", legend = :topright)
+# end
+# mp4(anim, "hw8/10-14-psi.mp4", fps = 75)
+
+#  plot([prob(psi0), prob(psi_1), prob(psi_2)], label = ["Initial" "Final (Euler dt = $dt) @ t=$tf" "Final (ODE abserr = 10^(-12)) @ t=$tf"], lw = 1.5, title = "Propagation of Wave Packet", xlabel = "x", ylabel = "|ψ(x)|²")
+# plot(
+# 	[prob(psi0), prob(psi_1), prob(psi_2), prob(psi_3), prob(psi_4)],
+# 	label = ["Initial" "Final (Euler dt = $dt) @ t=$tf" "Final (Euler dt = 0.00088) @ t=$tf" "Final (Euler dt = 0.00089) @ t=$tf" "Final (ODE abserr = 10^(-12)) @ t=$tf"],
+# 	lw = 1.5,
+# 	title = "Errors on Propagation of Wave Packet k_0 = 1000.0",
+# 	xlabel = "x",
+# )
+plot(
+	[prob(psi0), prob(psi_1), prob(psi_2), prob(psi_3), prob(psi_4), prob(psi_5)],
+	label = ["Initial" "Final (Euler dt = $dt) @ t=$tf" "Final (Euler dt = 0.00088) @ t=$tf" "Final (Euler dt = 0.00089) @ t=$tf" "Final (Euler dt = 0.0009) @ t=$tf" "Final (ODE abserr = 10^(-12)) @ t=$tf"],
+	lw = 1.5,
+	title = "Propagation of Wave Packet (k_0 = 1000.0)",
+	xlabel = "x",
+	ylabel = "|ψ(x)|²",
+)
+savefig("hw8/10-14-psi.png")
diff --git a/hw8/10-17-1.png b/hw8/10-17-1.png
new file mode 100644
index 0000000..5646660
--- /dev/null
+++ b/hw8/10-17-1.png
diff --git a/hw8/10-17-2.05.png b/hw8/10-17-2.05.png
new file mode 100644
index 0000000..885a984
--- /dev/null
+++ b/hw8/10-17-2.05.png
diff --git a/hw8/10-17-2.png b/hw8/10-17-2.png
new file mode 100644
index 0000000..7453918
--- /dev/null
+++ b/hw8/10-17-2.png
diff --git a/hw8/10-17-3.05.png b/hw8/10-17-3.05.png
new file mode 100644
index 0000000..d18dc6d
--- /dev/null
+++ b/hw8/10-17-3.05.png
diff --git a/hw8/10-17-3.png b/hw8/10-17-3.png
new file mode 100644
index 0000000..e140e1c
--- /dev/null
+++ b/hw8/10-17-3.png
diff --git a/hw8/10-17-width-.05.png b/hw8/10-17-width-.05.png
new file mode 100644
index 0000000..61829cc
--- /dev/null
+++ b/hw8/10-17-width-.05.png
diff --git a/hw8/10-17-width-.1.png b/hw8/10-17-width-.1.png
new file mode 100644
index 0000000..b2e77cb
--- /dev/null
+++ b/hw8/10-17-width-.1.png
diff --git a/hw8/10-17-width-.2.png b/hw8/10-17-width-.2.png
new file mode 100644
index 0000000..ec49ca7
--- /dev/null
+++ b/hw8/10-17-width-.2.png
diff --git a/hw8/10-17.jl b/hw8/10-17.jl
index 6dfff2b..d1d18b7 100644
--- a/hw8/10-17.jl
+++ b/hw8/10-17.jl
@@ -4,20 +4,29 @@ using Plots
 using LinearAlgebra
 using DifferentialEquations
 
+k_0 = 700.0
 # useful functions:
 
-function H(psi, dx) # action of Hamiltonian on wavefunction
+function H(psi, dx, k = 0, width = 0.1) # 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)
 
+	# add the potential function V_0 = 2 * k_0^2 for specific region
+	for i in 1:length(psi)
+		x = i * dx
+		if x >= 30.0 && x <= 30.0 + width
+			Hpsi[i] = 2 * k^2 * psi[i]
+		end
+	end
+
 	# 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)
+derivative(psi, dx) = -1.0im * H(psi, dx, k_0)
 
 function initialWavefunction(x::Vector{Float64}, x0 = 10.0, Delta = 1.0, k = 4.0)
 	Delta2 = Delta^2
@@ -32,15 +41,16 @@ 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
+N = 5000  # number of lattice points
+L = 40.0 # x runs from 0 to L
 dx = L / N
 
 x = range(0.0, L, N) |> collect # lattice along x-axis
-#println(x)
+println("dx = ", dx)
+
 
 # initial wavefunction has position (x0), width (Delta), and momentum (k)
-psi0 = initialWavefunction(x, 10.0, 0.05, 700.0)
+psi0 = initialWavefunction(x, 10.0, 0.05, k_0)
 
 # normalize wavefunction
 n = normalization(psi0, dx)
@@ -51,31 +61,49 @@ println("norm = ", normalization(psi0, dx))
 plot(prob(psi0))
 
 # integrate forward in time
-tf = 1.0
+tf = 3.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
+# 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)
+
+# 		# println("t = ", t, " norm = ", real(normalization(psi, dx)))
+
+# 		if real(normalization(psi, dx)) < 0.999 || real(normalization(psi, dx)) > 1.010
+# 			println("dt = ", dt, " t = ", t, " norm = ", normalization(psi, dx))
+# 			println("Normalization deviated too far from 1.0")
+# 			break
+# 		end
+# 	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
+sol = solve(problem, Tsit5(), reltol = 1e-14, abstol = 1e-14) # solve using Tsit5 algorithm to specified accuracy
 
 # compare initial and final wavefunction probabilities
-#psi = timeEvolve(psi0, tf, dt)
+# 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)])
+plot(prob(psi0))
+savefig("10-17-1.png")
+# plot a vertical line where the barrier is
+# barrier_on_lattice = 6.5 / dx
+barrier_on_lattice = 30.0 / dx
+plot(prob(psi), label = "final (t=$tf)", xlabel = "x", ylabel = "|ψ|^2", title = "Initial and final probability densities k_0=$k_0", lw = 1.5)
+plot!([barrier_on_lattice, barrier_on_lattice], [0.0, 0.3], color = :red, label = "barrier (width .2)", lw = 1.5, linestyle = :dash)
+savefig("10-17-width-.1.png")
+# plot([prob(psi0), prob(psi)], label = ["initial" "final (t=$tf)"], xlabel = "x", ylabel = "probability density", title = "Initial and final probability densities k_0=$k_0", lw = 1.5)
+# plot!([barrier_on_lattice, barrier_on_lattice], [0.0, 2.0], color = :red, label = "barrier (width .2)", lw = 1.5, linestyle = :dash)
+# savefig("10-17-width-.2.png")
 
 
diff --git a/hw8/10-3-long.png b/hw8/10-3-long.png
new file mode 100644
index 0000000..8d782a4
--- /dev/null
+++ b/hw8/10-3-long.png
diff --git a/hw8/10-3.jl b/hw8/10-3.jl
index 8616a69..2f4a690 100644
--- a/hw8/10-3.jl
+++ b/hw8/10-3.jl
@@ -60,16 +60,43 @@ function map_n_to_energies(n)
 	return e
 end
 
-n_max = 18
-n_to_e = [map_n_to_energies(n) for n in 1:n_max]
+n_s = collect(0:1:18)
+n_to_e = [map_n_to_energies(n) for n in n_s]
 
 # plot e[0] for all N
-eList = zeros(0)
-for i in 1:n_max
-	push!(eList, n_to_e[i][1])
+ground_state = []
+excited_1 = []
+excited_2 = []
+excited_3 = []
+for i in 1:length(n_to_e)
+	push!(ground_state, n_to_e[i][1])
+	push!(excited_1, n_to_e[i][2])
+	push!(excited_2, n_to_e[i][3])
+	push!(excited_3, n_to_e[i][4])
 end
 
-plot(eList)
+plot(ground_state, label = "groud state energy for n", xlabel = "n (level)", ylabel = "energy", title = "excited energy levels for V(n) = abs(x)^n", marker = :circle)
+plot!(excited_1, label = "1st excited state", marker = :circle)
+plot!(excited_2, label = "2nd excited state", marker = :circle)
+plot!(excited_3, label = "3rd excited state", marker = :circle)
+
+# plot the energies for an inifinite square well as a horizontial line
+# function excited_state_to_energy_inf_square_well(n)
+# 	return n^2 * pi^2 / 2
+# end
+
+# ground_state_inf_square_well = [excited_state_to_energy_inf_square_well(1) for i in 1:length(n_s)]
+# excited_1_inf_square_well = [excited_state_to_energy_inf_square_well(2) for i in 1:length(n_s)]
+# excited_2_inf_square_well = [excited_state_to_energy_inf_square_well(3) for i in 1:length(n_s)]
+# excited_3_inf_square_well = [excited_state_to_energy_inf_square_well(4) for i in 1:length(n_s)]
+
+# plot!(ground_state_inf_square_well, label = "ground state energy for infinite square well")
+# plot!(excited_1_inf_square_well, label = "1st excited state for infinite square well")
+# plot!(excited_2_inf_square_well, label = "2nd excited state for infinite square well")
+# plot!(excited_3_inf_square_well, label = "3rd excited state for infinite square well")
+
+
+savefig("hw8/10-3.png")
 
 # gs(x) = exp(-0.5 * x^2) # Gaussian that is exact ground state of SHO
 
diff --git a/hw8/10-3.png b/hw8/10-3.png
new file mode 100644
index 0000000..23c15a9
--- /dev/null
+++ b/hw8/10-3.png
diff --git a/hw8/TimeDependentSchrodingerEquation.jl b/hw8/TimeDependentSchrodingerEquation.jl
index dacbced..1d4f315 100644
--- a/hw8/TimeDependentSchrodingerEquation.jl
+++ b/hw8/TimeDependentSchrodingerEquation.jl
@@ -7,26 +7,26 @@ 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
+	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)
+	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)
+	n2 = dot(psi, psi) * dx
+	return sqrt(n2)
 end
 
 prob(psi) = real(psi .* conj(psi))
@@ -56,17 +56,17 @@ 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
+	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
+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)