diff options
Diffstat (limited to 'hw8/10-14.jl')
| -rw-r--r-- | hw8/10-14.jl | 68 |
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") |