add hw9 and hw8

1 files changed, 54 insertions, 14 deletions

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