1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
|
#!/usr/bin/env julia
using Plots
using LinearAlgebra
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
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)
end
function normalization(psi, dx) # normalization of wavefunction
n2 = dot(psi, psi) * dx
return sqrt(n2)
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
dx = L / N
x = range(0.0, L, N) |> collect # lattice along x-axis
#println(x)
# initial wavefunction has position (x0), width (Delta), and momentum (k)
psi0 = initialWavefunction(x, 10.0, 1.0, 1000.0)
# normalize wavefunction
n = normalization(psi0, dx)
psi0 = psi0 / n
println("norm = ", normalization(psi0, dx))
# plot initial wavefunction and probability density
plot(prob(psi0))
# integrate forward in time
tf = 1.0
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)
psiMid = psi + 0.5 * dt * derivative(psi, dx)
psi = psi + dt * derivative(psiMid, dx)
end
return psi
end
# 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_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 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")
|
