xxxxxxxxxx
 
html"<style>
main {
    max-width:75%;
    padding: unset;
    margin-right: unset !important;
    margin-bottom: 20px;
    align-self: center !important;
}
pluto-helpbox {
    visibility: hidden;
}
</style>"

5.1 ms

xxxxxxxxxx
 
begin
    using PlutoUI
    using Plots
    using SparseArrays
    using LinearAlgebra
    using Statistics
end

13.5 s

Simulating the 1-D TDSE

xxxxxxxxxx
 
md"""
# Simulating the 1-D TDSE
"""

8.8 μs

$i ℏ \frac{\partial ψ (x, t)}{\partial t} = (- \frac{ℏ^{2}}{2 m} \frac{\partial^{2}}{\partial x^{2}} + V (x)) ψ (x, t)$

xxxxxxxxxx
 
md"""
$$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)$$ 
​
"""

6.8 μs

Crank-Nicholson scheme

xxxxxxxxxx
 
md"### Crank-Nicholson scheme"

8.3 μs

We can write a forward and backward euler expansion for the time derivative (setting $ℏ = 1$ ):

$ψ (x, t + Δ t) = ψ (x, t) - i \cdot \hat{H} ψ (x, t) \cdot Δ t$

$ψ (x, t + Δ t) = ψ (x, t) - i \cdot \hat{H} ψ (x, t + Δ t) \cdot Δ t$

Taking the mean value of these two expressions, we obtain the crank nicholson scheme:

$(1 + \frac{1}{2} i \hat{H} Δ t) ψ (x, t + Δ t) = (1 - \frac{1}{2} i \hat{H} Δ t) ψ (x, t)$

This gives us a unitary evolution operator, which will preserve the normalization. The second derivative in $\hat{H}$ can be approximated using finite differences:

$\frac{\partial^{2} ψ_{j}^{n}}{\partial x^{2}} = \frac{ψ_{j + 1}^{n} - 2 ψ_{j}^{n} + ψ_{j - 1}^{n}}{(Δ x)^{2}}$

This can then be written in matrix form like so;

$M_{1} ψ^{n + 1} = M_{2} ψ^{n} ⟹ ψ^{n + 1} = M_{1}^{- 1} \cdot M_{2} ψ^{n}$

xxxxxxxxxx
 
md"""
We can write a forward and backward euler expansion for the time derivative (setting $$\hbar = 1$$):
​
$$\psi(x, t + \Delta t) = \psi(x, t) - i \cdot \hat{H}\psi(x,t) \cdot \Delta t$$
​
$$\psi(x, t + \Delta t) = \psi(x, t) - i \cdot \hat{H}\psi(x,t + \Delta t) \cdot \Delta t$$
​
Taking the mean value of these two expressions, we obtain the crank nicholson scheme:
​
$$\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)$$
​
This gives us a unitary evolution operator, which will preserve the normalization. The second derivative in $$\hat{H}$$ can be approximated using finite differences:
​
$$\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}$$
​
This can then be written in matrix form like so; 
​
$$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$$
"""

14.4 μs

$M_{1} = (\begin{matrix} ξ & - α & 0 & . . . & . . . & . . . & 0 \\ - α & ξ & - α & 0 & . & . & . \\ 0 & - α & ξ & - α & 0 & . & . \\ ⋮ & ⋱ & ⋱ & ⋱ & ⋱ & ⋱ & ⋮ \\ . & . & 0 & - α & ξ & - α & 0 \\ . & . & . & 0 & - α & ξ & - α \\ 0 & . . . & . . . & . . . & 0 & - α & ξ \end{matrix}) M_{2} = (\begin{matrix} γ & α & 0 & . . . & . . . & . . . & 0 \\ α & γ & α & 0 & . & . & . \\ 0 & α & γ & α & 0 & . & . \\ ⋮ & ⋱ & ⋱ & ⋱ & ⋱ & ⋱ & ⋮ \\ . & . & 0 & α & γ & α & 0 \\ . & . & . & 0 & α & γ & α \\ 0 & . . . & . . . & . . . & 0 & α & γ \end{matrix})$

xxxxxxxxxx
 
md"""
$$M_1 = \begin{pmatrix}
\xi & -\alpha & 0 & ... & ... & ... & 0\\
-\alpha & \xi & -\alpha & 0 & . &. & .\\
0 & -\alpha & \xi & -\alpha & 0 & . & .\\
\vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots\\
. & . & 0 & -\alpha & \xi & -\alpha & 0\\
. & . & . & 0 & -\alpha & \xi & -\alpha\\
0 & ... & ... & ...& 0 & -\alpha & \xi &
\end{pmatrix}
​
\hspace{2cm}
​
M_2 = \begin{pmatrix}
\gamma & \alpha & 0 & ... & ... & ... & 0\\
\alpha & \gamma & \alpha & 0 & . &. & .\\
0 & \alpha & \gamma & \alpha & 0 & . & .\\
\vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots\\
. & . & 0 & \alpha & \gamma & \alpha & 0\\
. & . & . & 0 & \alpha & \gamma & \alpha\\
0 & ... & ... & ...& 0 & \alpha & \gamma &
\end{pmatrix}$$
"""

7.2 μs

$α = \frac{i Δ t}{4 m (Δ x)^{2}} ξ = 1 + \frac{i Δ t}{2} \cdot (\frac{1}{m (Δ x)^{2}} + V (x)) γ = 1 - \frac{i Δ t}{2} \cdot (\frac{1}{m (Δ x)^{2}} + V (x))$

xxxxxxxxxx
 
md"""
$$\boxed{\alpha = \frac{i\Delta t}{4m(\Delta x)^2}}
\hspace{2cm}
\boxed{\xi = 1 + \frac{i\Delta t}{2} \cdot \left ( \frac{1}{m(\Delta x)^2} + V(x)\right)}
\hspace{2cm}
\boxed{\gamma = 1 - \frac{i\Delta t}{2} \cdot \left ( \frac{1}{m(\Delta x)^2} + V(x)\right)}$$
​
"""

8.4 μs

Boundary Condition

By virtue of this specific construction of the evolution matrices, the system has fixed boundary conditions ( $ψ (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.

xxxxxxxxxx
 
md""" #### Boundary Condition
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.
"""

29.9 μs

Initial Condition

xxxxxxxxxx
 
md"#### Initial Condition"

12.6 μs

initCond (generic function with 1 method)

xxxxxxxxxx
 
function initCond(xgrid, sType, μ, σ, p)
    
    if(sType == 1) #moving gaussian
       return @. Complex(1/(σ*sqrt(2*π))^0.5 * exp(p*im*xgrid - ((xgrid - μ)^2/(4*σ^2))))
    end
end

55.6 μs

Crank-Nicholson evolution

xxxxxxxxxx
 
md"#### Crank-Nicholson evolution"

8.4 μs

schrodingerEquation (generic function with 1 method)

xxxxxxxxxx
 
function schrodingerEquation(dx::Float64, xi::Float64, xf::Float64, m::Int64, dt::Float64, mass::Float64, wavepacket::Tuple{Float64, Float64, Float64}, potential)
    
    xgrid = collect(xi:dx:xf) #discretized spatial grid
    n = length(xgrid)
    
    V = potential.(xgrid)
    
    #CONSTRUCTING THE EVOLUTION MATRIX
    diag = ones(Complex, n)
    α = im*dt/(4*mass*dx^2) * diag[2:end]
    ξ = 1 .+ (im * dt/2) * (1/(mass*dx^2) .+ V)
    γ = 1 .- (im * dt/2) * (1/(mass*dx^2) .+ V)
    
    M₁ = LinearAlgebra.SymTridiagonal(ξ,-α)
    M₂ = LinearAlgebra.SymTridiagonal(γ, α)
    
    #INITIAL CONDITIONS
    ψ = zeros(Complex, (n, m))
    ψ[:, 1] = initCond(xgrid, 1, wavepacket...)
    
    #TIME EVOLUTION LOOP
    for j in 1:(m-1)
        ψ[:, j+1] = M₁\(M₂*ψ[:, j])
    end
    
    return (xgrid, V, ψ)
end

516 μs

Potential function

xxxxxxxxxx
 
md"#### Potential function"

6.4 μs

V (generic function with 1 method)

xxxxxxxxxx
 
function V(x)
    return (x > -5.0 && x < 5.0) ? 2.5 : 0.0
    #return 0
    #return 0.05*x^2
end

45.0 μs

Plotting

xxxxxxxxxx
 
md"#### Plotting"

6.6 μs

plotProbabilityDensity (generic function with 1 method)

xxxxxxxxxx
 
function plotProbabilityDensity(xgrid, V, ψ, m)
    #GIF CREATION
    anim = @gif for j = 1:10:m
        plot(xgrid, abs2.(ψ[:, j]), ylim = (-0.05, 0.5), label = "")
        plot!(xgrid, V, color = :black, label = "")
    end
end

168 μs

plotRealImg (generic function with 1 method)

xxxxxxxxxx
 
function plotRealImg(xgrid, V, ψ, m)
    anim = @gif for j = 1:10:m
        plot(xgrid, V, color = :black, label = "", ylim = (-0.5, 0.5))
        plot!(xgrid, abs2.(ψ[:, j]), color = :black, label = "|ψ|^2")
        plot!(xgrid, real.(ψ[:, j]), color = :red, label = "Re(ψ)")
        plot!(xgrid, imag.(ψ[:, j]), color = :blue, label = "Im(ψ)")
    end
end

80.8 μs

Simulation results

xxxxxxxxxx
 
md"## Simulation results"

5.9 μs

nSteps

xxxxxxxxxx
 
nSteps = 2200

176 ns

xxxxxxxxxx
 
data = schrodingerEquation(0.01, -40.0, 40.0, nSteps, 0.025, 5.0, (-25.0, 2.0, 6.0), V);

20.2 s