# Numerical solution of PDE:s, Part 4: Schrödinger equation

Posted by

In the earlier posts, I showed how to numerically solve a 1D or 2D diffusion or heat conduction problem using either explicit or implicit finite differencing. In the 1D example, the relevant equation for diffusion was

and an important property of the solution was the conservation of mass,

i.e. the integral of the concentration field over whole domain stays constant.

Next, I will show how to integrate the 1D time-dependent Schrödinger equation, which in a nondimensional form where we set $\hbar = 1$ and $m = 1$ reads:

Here $i$ is the imaginary unit and $V(x)$ is the potential energy as a function of $x$. The solutions of this equation must obey a conservation law similar to the mass conservation in the diffusion equation, the conservation of norm:

where the quantity $|\Psi (x,t)|$ is the modulus of the complex-valued function $\Psi (x,t)$ . This property of the solution is also called unitarity of the time evolution.

Apart from the TDSE, another way to represent the time development of this system is to find the normalized solutions $\psi_0 (x)$, $\psi_1 (x)$, $\psi_2 (x) \dots$ of the time-independent Schrödinger equation

and write the initial state $\Psi (x,0)$ as a linear combination of those basis functions:

This is possible because the solutions of the time-independent equation form a basis for the set of acceptable wave functions $\psi (x)$. Then, every term in that eigenfunction expansion is multiplied by a time dependent phase factor $\exp(-iE_n t)$:

The numbers $E_n$ are the eigenvalues corresponding to the solutions $\psi_n (x)$ and the function $\psi_0 (x)$ is called the ground state corresponding to potential $V(x)$, while the functions $\psi_1 (x)$ is the first excited state and so on.

The Schrödinger equation can’t be discretized by using either the explicit or implicit method that we used when solving the diffusion equation. The method is either numerically unstable or doesn’t conserve the normalization of the wave function (or both) if you try to do that. The correct way to discretize the Schrödinger equation is to replace the wave function with a discrete equivalent

and the potential energy function $V(x)$ with $V_{i;j}$ (or $V_i$ in the case of time-independent potential), and write an equation that basically tells that propagating the state $\Psi_{i;j}$ forward by half a time step gives the same result as propagating the state $\Psi_{i;j+1}$ backwards by half a time step:

Here we have

and

This kind of discretization is called the Crank-Nicolson method. As boundary conditions, we usually set that at the boundaries of the computational domain the wavefunction stays at value zero: $\Psi (0,t) = \Psi (L,t) = 0$ for any value of $t$. In the diffusion problem, this kind of a BC corresponded to infinite sinks at the boundaries, that annihilated anything that diffused through them. In the Schrödinger equation problem, which is a complex diffusion equation, the equivalent condition makes the boundaries impenetrable walls that deflect elastically anything that collides with them.

An R-Code that calculates the time evolution of a Gaussian initial wavefunction

in an area of zero potential:

for a domain 0 < x < 6, a lattice spacing $\Delta x = 0.05$, time interval 0 < t < 2 and time step $\Delta t = 0.01$, is given below:

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

lx <- 6.0 #length of the computational domain
lt <- 2.0 #length of the simulation time interval
nx <- 120 #number of discrete lattice points
nt <- 200 #number of timesteps
dx <- lx/nx #length of one discrete lattice cell
dt <- lt/nt #length of timestep

V = c(1:nx) #potential energies at discrete points

for(j in c(1:nx)) {
V[j] = 0 #zero potential
}

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

for(j in c(1:nx)) {
kappa2[j] <- as.complex(kappa1*2*dx*dx*V[j])
}

psi = as.complex(c(1:nx)) #array for the wave function values

for(j in c(1:nx)) {
psi[j] = as.complex(exp(-2*(j*dx-3)*(j*dx-3))) #Gaussian initial wavefunction
}

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

A = matrix(nrow=nx,ncol=nx) #matrix for forward time evolution
B = matrix(nrow=nx,ncol=nx) #matrix for backward time evolution

for(j in c(1:nx)) {
for(k in c(1:nx)) {
A[j,k]=0
B[j,k]=0
if(j==k) {
A[j,k] = 1 + 2*kappa1 + kappa2[j]
B[j,k] = 1 - 2*kappa1 - kappa2[j]
}
if((j==k+1) || (j==k-1)) {
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)) {
psi[l] <- sol[l]
}

if(k %% 3 == 1) { #make plots of |psi(x)|^2 on every third timestep
jpeg(file = paste("plot_",k,".jpg",sep=""))
plot(xaxis,abs(psi)^2,xlab="position (x)", ylab="Abs(Psi)^2",ylim=c(0,2))
title(paste("|psi(x,t)|^2 at t =",k*dt))
lines(xaxis,abs(psi)^2)
dev.off()
}
}

The output files are plots of the absolute squares of the wavefunction, and a few of them are shown below.

In the next simulation, I set the domain and discrete step sizes the same as above, but the initial state is:

Which is a Gaussian wave packet that has a nonzero momentum in the positive x-direction. This is done by changing the line

for(j in c(1:nx)) {
psi[j] = as.complex(exp(-2*(j*dx-3)*(j*dx-3))) #Gaussian initial wavefunction
}+(1i)*j*dx

into

for(j in c(1:nx)) {
psi[j] = as.complex(exp(-2*(j*dx-3)*(j*dx-3)+(1i)*j*dx)) #Gaussian initial wavefunction
}

The plots of $|\Psi (x,t)|^2$ for several values of $t$ are shown below

and there you can see how the wave packet collides with the right boundary of the domain and bounces back.

In the last simulation, I will set the potential function to be

which is a harmonic oscillator potential, and with the nondimensional mass $m =1$ and Planck constant $\hbar = 1$ the ground state $\psi _0 (x)$ of this system is

If I’d set the initial state to be $\Psi (x,0) = \psi_0 (x)$, or any other solution of the time-independent SE, the modulus of the wavefunction would not change at all. To get something interesting to happen, I instead set an initial state that is a displaced version of the ground state:

The solution can be obtained with the code shown below:

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

lx <- 6.0 #length of the computational domain
lt <- 3.0 #length of the simulation time interval
nx <- 360 #number of discrete lattice points
nt <- 300 #number of timesteps
dx <- lx/nx #length of one discrete lattice cell
dt <- lt/nt #length of timestep

V = c(1:nx) #potential energies at discrete points

for(j in c(1:nx)) {
V[j] = as.complex(2*(j*dx-3)*(j*dx-3)) #Harmonic oscillator potential with k=4
}

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

for(j in c(1:nx)) {
kappa2[j] <- as.complex(kappa1*2*dx*dx*V[j])
}

psi = as.complex(c(1:nx)) #array for the wave function values

for(j in c(1:nx)) {
psi[j] = as.complex(exp(-(j*dx-2)*(j*dx-2))) #Gaussian initial wavefunction, displaced from equilibrium
}

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

A = matrix(nrow=nx,ncol=nx) #matrix for forward time evolution
B = matrix(nrow=nx,ncol=nx) #matrix for backward time evolution

for(j in c(1:nx)) {
for(k in c(1:nx)) {
A[j,k]=0
B[j,k]=0
if(j==k) {
A[j,k] = 1 + 2*kappa1 + kappa2[j]
B[j,k] = 1 - 2*kappa1 - kappa2[j]
}
if((j==k+1) || (j==k-1)) {
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)) {
psi[l] <- sol[l]
}

if(k %% 3 == 1) { #make plots of Abs(psi(x))^2 on every third timestep
jpeg(file = paste("plot_",k,".jpg",sep=""))
plot(xaxis,abs(psi)^2, xlab="position (x)", ylab="Abs(Psi)^2",ylim=c(0,2))
title(paste("|psi(x,t)|^2 at t =",k*dt))
lines(xaxis,abs(psi)^2)
lines(xaxis,V)
dev.off()
}
}

and the solution at different values of $t$ look like this (images and video):

Here the shape of the Hookean potential energy is plotted in the same images. So, here you see how the center of the Gaussian wavefunction oscillates around the point $x = 3$, just like a classical mechanical harmonic oscillator does when set free from a position that is displaced from equilibrium.

By changing the code that produces the output images, we can also get a sequence of plots of the imaginary part of the wavefunction:

if(k %% 3 == 1) { #make plots of Im(psi(x)) on every third timestep
jpeg(file = paste("plot_",k,".jpg",sep=""))
plot(xaxis,Im(psi), xlab="position (x)", ylab="Im(Psi)",ylim=c(-1.5,1.5))
title(paste("Im(psi(x,t)) at t =",k*dt))
lines(xaxis,Im(psi))
lines(xaxis,V)
dev.off()
}

and the resulting plots look like this: