xxxxxxxxxx
12
1
html"<style>
2
main {
3
max-width:75%;
4
padding: unset;
5
margin-right: unset !important;
6
margin-bottom: 20px;
7
align-self: center !important;
8
}
9
pluto-helpbox {
10
visibility: hidden;
11
}
12
</style>"
xxxxxxxxxx
7
1
begin
2
using PlutoUI
3
using Plots
4
using SparseArrays
5
using LinearAlgebra
6
using Statistics
7
end
Simulating the 1-D TDSE
xxxxxxxxxx
3
1
md"""
2
# Simulating the 1-D TDSE
3
"""
xxxxxxxxxx
4
1
md"""
2
$$i\hbar\frac{\partial \psi(x, t)}{\partial t} = \left(-\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V(x) \right) \psi(x, t)$$
3
4
"""
Crank-Nicholson scheme
xxxxxxxxxx
1
1
md"### Crank-Nicholson scheme"
We can write a forward and backward euler expansion for the time derivative (setting
Taking the mean value of these two expressions, we obtain the crank nicholson scheme:
This gives us a unitary evolution operator, which will preserve the normalization. The second derivative in
This can then be written in matrix form like so;
xxxxxxxxxx
19
1
md"""
2
We can write a forward and backward euler expansion for the time derivative (setting $$\hbar = 1$$):
3
4
$$\psi(x, t + \Delta t) = \psi(x, t) - i \cdot \hat{H}\psi(x,t) \cdot \Delta t$$
5
6
$$\psi(x, t + \Delta t) = \psi(x, t) - i \cdot \hat{H}\psi(x,t + \Delta t) \cdot \Delta t$$
7
8
Taking the mean value of these two expressions, we obtain the crank nicholson scheme:
9
10
$$\left(1 + \frac{1}{2} i\hat{H}\Delta t\right)\psi(x, t + \Delta t) = \left(1 - \frac{1}{2} i\hat{H}\Delta t\right)\psi(x, t)$$
11
12
This gives us a unitary evolution operator, which will preserve the normalization. The second derivative in $$\hat{H}$$ can be approximated using finite differences:
13
14
$$\frac{\partial^2 \psi_j^n}{\partial x^2} = \frac{\psi_{j+1}^n - 2\psi_j^n + \psi_{j-1}^n}{(\Delta x)^2}$$
15
16
This can then be written in matrix form like so;
17
18
$$M_1 \psi^{n+1} = M_2 \psi^n \hspace{1cm} \implies \hspace{1cm} \psi^{n+1} = M_1 ^{-1} \cdot M_2 \psi^n$$
19
"""
xxxxxxxxxx
23
1
md"""
2
$$M_1 = \begin{pmatrix}
3
\xi & -\alpha & 0 & ... & ... & ... & 0\\
4
-\alpha & \xi & -\alpha & 0 & . &. & .\\
5
0 & -\alpha & \xi & -\alpha & 0 & . & .\\
6
\vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots\\
7
. & . & 0 & -\alpha & \xi & -\alpha & 0\\
8
. & . & . & 0 & -\alpha & \xi & -\alpha\\
9
0 & ... & ... & ...& 0 & -\alpha & \xi &
10
\end{pmatrix}
11
12
\hspace{2cm}
13
14
M_2 = \begin{pmatrix}
15
\gamma & \alpha & 0 & ... & ... & ... & 0\\
16
\alpha & \gamma & \alpha & 0 & . &. & .\\
17
0 & \alpha & \gamma & \alpha & 0 & . & .\\
18
\vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots\\
19
. & . & 0 & \alpha & \gamma & \alpha & 0\\
20
. & . & . & 0 & \alpha & \gamma & \alpha\\
21
0 & ... & ... & ...& 0 & \alpha & \gamma &
22
\end{pmatrix}$$
23
"""
xxxxxxxxxx
8
1
md"""
2
$$\boxed{\alpha = \frac{i\Delta t}{4m(\Delta x)^2}}
3
\hspace{2cm}
4
\boxed{\xi = 1 + \frac{i\Delta t}{2} \cdot \left ( \frac{1}{m(\Delta x)^2} + V(x)\right)}
5
\hspace{2cm}
6
\boxed{\gamma = 1 - \frac{i\Delta t}{2} \cdot \left ( \frac{1}{m(\Delta x)^2} + V(x)\right)}$$
7
8
"""
Boundary Condition
By virtue of this specific construction of the evolution matrices, the system has fixed boundary conditions (
xxxxxxxxxx
3
1
md""" #### Boundary Condition
2
By virtue of this specific construction of the evolution matrices, the system has fixed boundary conditions ($\psi(x) = 0$ at the boundaries). This effectively emulates an infinite wall on either sides, resulting in a particle-in-a-box system in the absence of any other potential inside.
3
"""
Initial Condition
xxxxxxxxxx
1
1
md"#### Initial Condition"
initCond (generic function with 1 method)
xxxxxxxxxx
6
1
function initCond(xgrid, sType, μ, σ, p)
2
3
if(sType == 1) #moving gaussian
4
return @. Complex(1/(σ*sqrt(2*π))^0.5 * exp(p*im*xgrid - ((xgrid - μ)^2/(4*σ^2))))
5
end
6
end
Crank-Nicholson evolution
xxxxxxxxxx
1
1
md"#### Crank-Nicholson evolution"
schrodingerEquation (generic function with 1 method)
xxxxxxxxxx
27
1
function schrodingerEquation(dx::Float64, xi::Float64, xf::Float64, m::Int64, dt::Float64, mass::Float64, wavepacket::Tuple{Float64, Float64, Float64}, potential)
2
3
xgrid = collect(xi:dx:xf) #discretized spatial grid
4
n = length(xgrid)
5
6
V = potential.(xgrid)
7
8
#CONSTRUCTING THE EVOLUTION MATRIX
9
diag = ones(Complex, n)
10
α = im*dt/(4*mass*dx^2) * diag[2:end]
11
ξ = 1 .+ (im * dt/2) * (1/(mass*dx^2) .+ V)
12
γ = 1 .- (im * dt/2) * (1/(mass*dx^2) .+ V)
13
14
M₁ = LinearAlgebra.SymTridiagonal(ξ,-α)
15
M₂ = LinearAlgebra.SymTridiagonal(γ, α)
16
17
#INITIAL CONDITIONS
18
ψ = zeros(Complex, (n, m))
19
ψ[:, 1] = initCond(xgrid, 1, wavepacket...)
20
21
#TIME EVOLUTION LOOP
22
for j in 1:(m-1)
23
ψ[:, j+1] = M₁\(M₂*ψ[:, j])
24
end
25
26
return (xgrid, V, ψ)
27
end
Potential function
xxxxxxxxxx
1
1
md"#### Potential function"
V (generic function with 1 method)
xxxxxxxxxx
5
1
function V(x)
2
return (x > -5.0 && x < 5.0) ? 2.5 : 0.0
3
#return 0
4
#return 0.05*x^2
5
end
Plotting
xxxxxxxxxx
1
1
md"#### Plotting"
plotProbabilityDensity (generic function with 1 method)
xxxxxxxxxx
7
1
function plotProbabilityDensity(xgrid, V, ψ, m)
2
#GIF CREATION
3
anim = for j = 1:10:m
4
plot(xgrid, abs2.(ψ[:, j]), ylim = (-0.05, 0.5), label = "")
5
plot!(xgrid, V, color = :black, label = "")
6
end
7
end
plotRealImg (generic function with 1 method)
xxxxxxxxxx
8
1
function plotRealImg(xgrid, V, ψ, m)
2
anim = for j = 1:10:m
3
plot(xgrid, V, color = :black, label = "", ylim = (-0.5, 0.5))
4
plot!(xgrid, abs2.(ψ[:, j]), color = :black, label = "|ψ|^2")
5
plot!(xgrid, real.(ψ[:, j]), color = :red, label = "Re(ψ)")
6
plot!(xgrid, imag.(ψ[:, j]), color = :blue, label = "Im(ψ)")
7
end
8
end
Simulation results
xxxxxxxxxx
1
1
md"## Simulation results"
2200
xxxxxxxxxx
1
1
nSteps = 2200
xxxxxxxxxx
1
1
data = schrodingerEquation(0.01, -40.0, 40.0, nSteps, 0.025, 5.0, (-25.0, 2.0, 6.0), V);