Wavefunction of the 1D Schrodinger Equation using Julia
One of the cornerstones of quantum mechanics is that of wave-particle duality.
Down at the microscope level, things behave extremely counterintuitively. And things quite literally, appear and disappear at the blink of an eye.
And matter, behaves more like waves, fluctuating back and forth.
The governing equations of these waves behave according to the Schrodinger equations, discovered by the Austrian Physicist, Erwin Schrodinger.
For this example, we shall be dealing with one of the simplest cases. But even then, it requires non-trivial mathematical treatment.
A particle confined to a 1 dimensional box.
It may exist anywhere within the box, but not outside.
A tiny little pocket universe, if you will.
It’s behaviour is governed by a 2nd order differential equation as shown below.
Where h_bar is the plank constant, and m is the mass of the particle.
The general solution
2nd order differential equations admit general solutions of the form:
Utilising Euler’s formula for complex numbers
We can rewrite the solution into the following form:
Boundary conditions
Since the particle is confined to exist within the box, it is mandated to be unable to exist at the boundaries of the box.
At zero and at L
Hence,
for x = 0, C sin(k(0)) + D cos(k(0)) = 0, > D = 0
for x = L, C sin(kL) + D cos(kL) = 0
C sin(kL) + 0 = 0Csin(kL) = 0
kL = nπ , where n = 1, 2, 3,......=> k = nπ/L
C is the normalisation constant,
√(2/L) 
The wave-function
Hence we arrive the wave-function given below.
It governs the probability of locating this particle within the well.
Julia code
I actually chanced upon the Python code snippets, provided in the above link.
And decided to re-write the code in Julia.
As no special library is required, we just need the plotting library Plots.
using Plots
plotly()
Plots.PlotlyBackend()#Declare plotly backend
#At the time of publication, default "GR" backend is glitchy
Create the wavefunction and feed it the relevant inputs
#Wave function
function Ψ(x,n,L)
√(2.0/L)*sin(Float16(n*π*x)/L)
end#n = energy level (n = 1,2,3, .... ∞ )
#L = box width in angstromL = 20
n = 3
x = 0:0.1:L
y = Ψ.(x,n,L)
Creating the plots to visualise the positional probability of the particle.
The wave function provides the probability amplitude of locating the particle.
We need to square it, in order to obtain the probability distribution.
p1=plot(x,y, w=3, label=”probability amplitude: Ψ(x,$n,$L)”,
xlabel=”x”, ylabel=”Ψ”, ylim=(-1,1))p2=plot(x,y.^2, w=3, label=”probability distribution: Ψ²(x,$n,$L)”, xlabel=”x”, ylabel=”Ψ²”, line=:red)plot(p1, p2, layout=(1,2), legend=:bottomright, legendfontsize=7)
And as seen, from the red graph below, the are 3 places inside the potential well that the particle is the most likely to be located.
And finally, as a sanity check to ensure that our results are sound, let’s integrate the expression provided by psi_square to ensure that we are correct.
We shall integrate under expression, to see if we get the value of 1.
As the particle is confined to exist within the box, the likelihood of it being within the box, has to be 100%.
And we use the following Gauss Konrod integration algorithm
using QuadGK
g(x) = (√(2.0/L)*sin(Float64(n*π*x)/L))^2integral, err = quadgk(x -> g(x), 0, 20, rtol=1e-5)# rtol here refers to the error tolerance.
# Often in numerical methods, we cannot obtain an exact solution, and we terminate the solver once it get close enough. .
.
.< (1.0, 0.0) >
And as expected, we obtain the value of one.
Since the expression is simple enough to solve, the relative error is zero.
The entirety of the code is consolidated below for the reader’s convenience.
