add hw9 and hw8

1 files changed, 46 insertions, 18 deletions

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")