From bc515d3acdd94847b6e7aa6135bc234b46161db6 Mon Sep 17 00:00:00 2001 From: sotech117 Date: Tue, 7 May 2024 07:00:43 -0400 Subject: add hw9 and hw8 --- hw8/10-14-0.0001.png | Bin 0 -> 17051 bytes hw8/10-14-0.001.png | Bin 0 -> 60938 bytes hw8/10-14-0.01.png | Bin 0 -> 19005 bytes hw8/10-14-0.1.png | Bin 0 -> 17313 bytes hw8/10-14-E.png | Bin 0 -> 14507 bytes hw8/10-14-psi.mp4 | Bin 0 -> 3519376 bytes hw8/10-14-psi.png | Bin 0 -> 61965 bytes hw8/10-14.jl | 68 +++++++++++++++++++++++++------- hw8/10-17-1.png | Bin 0 -> 10080 bytes hw8/10-17-2.05.png | Bin 0 -> 24566 bytes hw8/10-17-2.png | Bin 0 -> 30436 bytes hw8/10-17-3.05.png | Bin 0 -> 28064 bytes hw8/10-17-3.png | Bin 0 -> 29569 bytes hw8/10-17-width-.05.png | Bin 0 -> 46547 bytes hw8/10-17-width-.1.png | Bin 0 -> 28921 bytes hw8/10-17-width-.2.png | Bin 0 -> 30538 bytes hw8/10-17.jl | 64 +++++++++++++++++++++--------- hw8/10-3-long.png | Bin 0 -> 27534 bytes hw8/10-3.jl | 39 +++++++++++++++--- hw8/10-3.png | Bin 0 -> 60650 bytes hw8/TimeDependentSchrodingerEquation.jl | 38 +++++++++--------- 21 files changed, 152 insertions(+), 57 deletions(-) create mode 100644 hw8/10-14-0.0001.png create mode 100644 hw8/10-14-0.001.png create mode 100644 hw8/10-14-0.01.png create mode 100644 hw8/10-14-0.1.png create mode 100644 hw8/10-14-E.png create mode 100644 hw8/10-14-psi.mp4 create mode 100644 hw8/10-14-psi.png create mode 100644 hw8/10-17-1.png create mode 100644 hw8/10-17-2.05.png create mode 100644 hw8/10-17-2.png create mode 100644 hw8/10-17-3.05.png create mode 100644 hw8/10-17-3.png create mode 100644 hw8/10-17-width-.05.png create mode 100644 hw8/10-17-width-.1.png create mode 100644 hw8/10-17-width-.2.png create mode 100644 hw8/10-3-long.png create mode 100644 hw8/10-3.png (limited to 'hw8') diff --git a/hw8/10-14-0.0001.png b/hw8/10-14-0.0001.png new file mode 100644 index 0000000..0dc656a Binary files /dev/null and b/hw8/10-14-0.0001.png differ diff --git a/hw8/10-14-0.001.png b/hw8/10-14-0.001.png new file mode 100644 index 0000000..fcf13d4 Binary files /dev/null and b/hw8/10-14-0.001.png differ diff --git a/hw8/10-14-0.01.png b/hw8/10-14-0.01.png new file mode 100644 index 0000000..8c6d987 Binary files /dev/null and b/hw8/10-14-0.01.png differ diff --git a/hw8/10-14-0.1.png b/hw8/10-14-0.1.png new file mode 100644 index 0000000..74f5045 Binary files /dev/null and b/hw8/10-14-0.1.png differ diff --git a/hw8/10-14-E.png b/hw8/10-14-E.png new file mode 100644 index 0000000..d13b721 Binary files /dev/null and b/hw8/10-14-E.png differ diff --git a/hw8/10-14-psi.mp4 b/hw8/10-14-psi.mp4 new file mode 100644 index 0000000..37c0d10 Binary files /dev/null and b/hw8/10-14-psi.mp4 differ diff --git a/hw8/10-14-psi.png b/hw8/10-14-psi.png new file mode 100644 index 0000000..7911c89 Binary files /dev/null and b/hw8/10-14-psi.png differ 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 Binary files /dev/null and b/hw8/10-17-1.png differ diff --git a/hw8/10-17-2.05.png b/hw8/10-17-2.05.png new file mode 100644 index 0000000..885a984 Binary files /dev/null and b/hw8/10-17-2.05.png differ diff --git a/hw8/10-17-2.png b/hw8/10-17-2.png new file mode 100644 index 0000000..7453918 Binary files /dev/null and b/hw8/10-17-2.png differ diff --git a/hw8/10-17-3.05.png b/hw8/10-17-3.05.png new file mode 100644 index 0000000..d18dc6d Binary files /dev/null and b/hw8/10-17-3.05.png differ diff --git a/hw8/10-17-3.png b/hw8/10-17-3.png new file mode 100644 index 0000000..e140e1c Binary files /dev/null and b/hw8/10-17-3.png differ diff --git a/hw8/10-17-width-.05.png b/hw8/10-17-width-.05.png new file mode 100644 index 0000000..61829cc Binary files /dev/null and b/hw8/10-17-width-.05.png differ diff --git a/hw8/10-17-width-.1.png b/hw8/10-17-width-.1.png new file mode 100644 index 0000000..b2e77cb Binary files /dev/null and b/hw8/10-17-width-.1.png differ diff --git a/hw8/10-17-width-.2.png b/hw8/10-17-width-.2.png new file mode 100644 index 0000000..ec49ca7 Binary files /dev/null and b/hw8/10-17-width-.2.png differ 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 Binary files /dev/null and b/hw8/10-3-long.png differ 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 Binary files /dev/null and b/hw8/10-3.png differ 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) -- cgit v1.2.3-70-g09d2