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