Numerical solution of PDE:s, Part 7: 2D Schrödinger equation

Haven’t been posting for a while, but here’s something new… Earlier I showed how to solve the 1D Schrödinger equation numerically in different situations. Now I’m going to show how to calculate the evolution of a 2D wavepacket in a potential energy field that has been constructed to mimic the classical “two-slit experiment” which shows how the mechanics of low-mass particles like electrons can exhibit interference similar to the mechanics of classical waves (sound, light, water surface, and so on).

A 2D Schrödinger equation for a single particle in a time-independent background potential V(x,y) is


Where the particle mass has been set to 1 and the Planck’s constant to 2\pi.

To solve this numerically, we need the Crank-Nicolson method, as was the case when solving the 1D problem. More specifically, the linear system to be solved is




where the wavefunction now has two position indices and one time index, and the potential energy has only two position indices.

To form a model of the two-slit experiment, we choose a domain 0 < x < 6; 0 < y < 6 and make a potential energy function defined by

IF (x < 2.2 OR x > 3.8 OR (x > 2.7 AND x < 3.3)) THEN IF (3.7 < y < 4) THEN V(x,y) = 30

IF (x < 0.5 OR x > 5.5 OR y < 0.5 OR y > 5.5) THEN V(x,y) = 30

Otherwise V(x,y) = 0.

which corresponds to having hard walls surrounding the domain and a barrier with two holes around the line y = 3.85

For an initial condition, we choose a Gaussian wavepacket that has a nonzero expectation value of the momentum in y-direction:


An R-Code that solves this problem for a time interval 0 < t < 1 is

library(graphics) #load the graphics library needed for plotting

lx <- 6.0 #length of the computational domain (both x and y)
lt <- 1.0 #length of the simulation time interval
nx <- 47 #number of discrete lattice points in x and y directions
nt <- 60 #number of timesteps
dx <- lx/nx #length of one discrete lattice cell, same in x and y directions
dt <- lt/nt #length of timestep

V = matrix(nrow=nx,ncol=nx) #potential energies at discrete xy points

for(j in c(1:nx)) { #construct the potential field with a double-slit barrier
for(k in c(1:nx)) {
V[j,k] = 0+0i
if((j*dx<2.2)||(j*dx>3.8)||((j*dx>2.7) && (j*dx<3.3))) {
if((k*dx>3.7) && (k*dx<4.0)) {
V[j,k] = 30+0i #No significant density is going to go through these barriers
if((j*dx>5.5) || (j*dx<0.5) || (k*dx>5.5) || (k*dx<0.5)) {
V[j,k] = 30+0i

kappa1 = (1i)*dt/(2*dx*dx) #an element needed for the matrices
kappa2 <- c(1:(nx*nx)) #another element

for(j in c(1:nx)) {
for(k in c(1:nx)) {
kappa2[(j-1)*nx+k] <- kappa1*2*dx*dx*V[j,k]

psi = c(1:(nx*nx)) #array for the wave function values

for(j in c(1:nx)) {
for(k in c(1:nx)) {
psi[(j-1)*nx+k] = exp(-2*(k*dx-1.7)*(k*dx-1.7)-2*(j*dx-3)*(j*dx-3)+(2i)*k*dx) #Gaussian initial wavefunction
if((j*dx > 5.5)||(j*dx < 0.5)||(k*dx > 5.5)||(k*dx < 0.5)) {
psi[(j-1)*nx+k] = as.complex(0)

xaxis <- c(1:nx)*dx #the x values corresponding to the discrete lattice points
yaxis <- c(1:nx)*dx #the y values

A = matrix(nrow=nx*nx,ncol=nx*nx) #matrix for forward time evolution
B = matrix(nrow=nx*nx,ncol=nx*nx) #matrix for backward time evolution
P = matrix(nrow=nx,ncol=nx) #matrix for the solution after time stepping
IP = matrix(nrow=4*nx,ncol=4*nx) #matrix for the higher resolution image obtained by interpolation

for(j in c(1:(nx*nx))) { #Set the values for time evolution matrix elements
for(k in c(1:(nx*nx))) {
if(j==k) {
A[j,k] = 1 + 4*kappa1 + kappa2[j]
B[j,k] = 1 – 4*kappa1 – kappa2[j]
if((k==j+1) || (k==j-1)) {
A[j,k] = -kappa1
B[j,k] = kappa1
if((k==j+nx)||(k==j-nx)) {
A[j,k] = -kappa1
B[j,k] = kappa1

for(k in c(1:nt)) { #main time stepping loop

sol <- solve(A,B%*%psi) #solve the system of equations

for(l in c(1:(nx*nx))) {
psi[l] <- sol[l]

for(l in c(1:nx)) {
for(m in c(1:nx)) {
P[l,m] = abs(psi[(l-1)*nx + m])*abs(psi[(l-1)*nx + m]) #square of the absolute value of wave function
if(abs(V[l,m]) > 5) P[l,m] = 2
for(l in c(1:(nx-1))) {
for(m in c(1:(nx-1))) { #make a bitmap with 4 times more pixels, using linear interpolation
IP[4*l-3,4*m-3] = P[l,m]
IP[4*l-2,4*m-3] = P[l,m]+0.25*(P[l+1,m]-P[l,m])
IP[4*l-1,4*m-3] = P[l,m]+0.5*(P[l+1,m]-P[l,m])
IP[4*l,4*m-3] = P[l,m]+0.75*(P[l+1,m]-P[l,m])

for(l in c(1:(4*nx))) {
for(m in c(1:(nx-1))) {
IP[l,4*m-2] = IP[l,4*m-3]+0.25*(IP[l,4*m+1]-IP[l,4*m-3])
IP[l,4*m-1] = IP[l,4*m-3]+0.5*(IP[l,4*m+1]-IP[l,4*m-3])
IP[l,4*m] = IP[l,4*m-3]+0.75*(IP[l,4*m+1]-IP[l,4*m-3])

jpeg(file = paste(“plot_abs_”,k,”.jpg”,sep=””)) #save the image
image(IP, zlim = c(0,0.15))


The code produces a sequence of image files, where the probability density is plotted with colors, as an output. Some representative images from this sequence (converted to grayscale) is shown below:

A video of the time evolution is shown below:

The treshold for maximum white color has been chosen to be quite low, to make the small amount of probability density that crosses the barrier visible.

The discrete grid of points has been made quite coarse here to keep the computation time reasonable, and the resolution has been increased artificially by using linear interpolation between the discrete points.

So, now we’ve seen how to solve the motion of 2D wavepackets moving around obstacles. In the next numerical methods post, I’ll go through the numerical solution of a nonlinear PDE.