## Numerical solution of PDE:s, Part 6: Adiabatic approximation for quantum dynamics

Having solved the time-dependent Schrödinger equation both in real and imaginary time, we can move forward to investigate systems where the potential energy function $V$ has an explicit time dependence in it:

In this kind of systems, the expectation value of the Hamiltonian operator doesn’t have to stay constant.

Time-dependent perturbation theory is one method for finding approximate solutions for this kind of problems, but here I will handle a simpler example, which is called adiabatic approximation.

Suppose that the potential energy function $V(x,t)$ is known. Now, let’s say that we also know the solutions of the time-independent Schrödinger equation

for any value of $t$. I denote the solutions as $\psi_n (x;t)$, where it is understood that $x$ is a variable and $t$ is a parameter. Now, if the function $V(x,t)$ changes very slowly as a function of time, i.e. its partial derivative with respect to $t$ is small at all points of the domain, we can use the adiabatic approximation, which says that if the initial state $\Psi (x,0)$ is the ground state for the potential $V(x,0)$, then the state at time $t$ is the ground state for the potential $V(x,t)$.

So, we can change a ground state of one potential $V_1 (x)$ into the ground state of another potential $V_2 (x)&bg=ffffff&fg=000000$ by making a continuous change from $V_1 (x)$ to $V_2 (x)$ slowly enough.

Let’s test this by chooosing a function $V$ as

i.e. a Hookean potential that moves to the positive x-direction with constant speed. If we set the wavefunction at $t=0$ to

which is the ground state corresponding to $V(x,0)$, the time depelopment of the wavepacket should be like

which means that is moves with the same constant speed as the bottom of the potential $V$. This can be calculated with the R-Code below:

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

lx <- 6.0 #length of the computational domain
lt <- 15.0 #length of the simulation time interval
nx <- 100 #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 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(-(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 #off-diagonal elements
B[j,k] = kappa1
}
}
}

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

for(j in c(1:nx)) {
V[j] = as.complex(2*(j*dx-3-k*dt*0.05)*(j*dx-3-k*dt*0.05)) #time dependent potential
}

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

for(l in c(1:nx)) {
for(m in c(1:nx)) {
A[l,m]=0
B[l,m]=0
if(l==m) {
A[l,m] = 1 + 2*kappa1 + kappa2[m]
B[l,m] = 1 - 2*kappa1 - kappa2[m]
}
if((l==m+1) || (l==m-1)) {
A[l,m] = -kappa1
B[l,m] = kappa1
}
}
}
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) 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 in the following sequences of images you see that the approximation is quite good for $v = 0.05$

As an animation, this process looks like shown below:

By doing the same calculation again, but this time with $v = 3$, the image sequence looks like this:

where it is obvious that the approximation doesn’t work anymore.

## Numerical solution of PDE:s, Part 5: Schrödinger equation in imaginary time

In the last post, I mentioned that the solution of the time dependent 1D Schrödinger equation

can be written by expanding the initial state $\Psi (x,0)$ in the basis of the solutions of time-independent Schrödinger equation

and multiplying each term in the expansion by a complex-valued time dependent phase factor:

Now, assuming that the energies $E_n$ are all positive and are ordered so that $E_0$ is the smallest of them, we get an interesting result by replacing the time variable $t$ with an imaginary time variable $s = it$. The function $\Psi (x,s)$ is then

which is a sum of exponentially decaying terms, of which the first one, $c_0 \exp(-E_0 s)\psi_0 (x)$ decays slowest. So, in the limit of large $s$, the wave function is

i.e. after normalizing it, it’s approximately the same as the ground state $\psi_0$. Also, if $s$ is a large number, we have

or

So, here we have a way to use the TDSE to numerically calculate the ground state wavefunction and corresponding energy eigenvalue for any potential energy function $V(x)$. This is very useful, as the ground state can’t usually be solved analytically, except for some very simple potential energy functions such as the harmonic oscillator potential.

To test this imaginary time method, I will approximately calculate the ground state of an anharmonic oscillator, described by a Schrödinger equation

as an initial state, I will choose a Gaussian function

and the computational domain is defined by 0 < x < 6, 0 < t < 3, $\Delta x = 0.05$, $\Delta t = 0.01$.

An R-Code that performs this calculation is 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 <- 120 #number of discrete lattice points
nt <- 300 #number of timesteps
dx <- lx/nx #length of one discrete lattice cell
dt <- (-1i)*lt/nt #length of imaginary timestep

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

for(j in c(1:nx)) {
V[j] = as.complex(0.25*(j*dx-3)*(j*dx-3)+0.15*(j*dx-3)*(j*dx-3)*(j*dx-3)*(j*dx-3)) #anharmonic 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))) #Gasussian 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]
}

nrm = 0
for (l in c(1:nx)) nrm <- nrm + dx*abs(psi[l])*abs(psi[l])

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

And a set of plots of $|\Psi (x,t)|^2$ for several values of $t$ look like this:

Here the shape of the anharmonic potential has been plotted to the same images.

The problem with this method for finding a ground state is that if the system has more degrees of freedom than a single coordinate x, the matrices in the linear systems needed in the time stepping quickly become very large and the calculation becomes prohibitively slow. To make this worse, the Crank-Nicolson method of time stepping can’t be parallelized for multiple processors. However, there is another way to compute the evolution of a wave function in imaginary time, which is called Diffusion Monte Carlo, and that is easily parallelizable. DMC is one practical way for calculating ground state wave functions for multi-particle systems such as a helium or a lithium atom.

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

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:

## About the NASA’s TRAPPIST-1 findings

Yesterday, I was excitedly following the NASA press conference about the exoplanets found orbiting the star TRAPPIST-1 located 39 light years from earth. When I heard the day before that they’re going to announce something, the first idea that came to my mind was that they must have observed a combination of organic molecules and oxygen in some distant planet’s atmosphere, which would be a chemical nonequilibrium state that is difficult to explain without some kind of biology being present. Yes, of course you can’t have an atmosphere where there’s something like a mixture of methane and oxygen that is within flammability limits, but even very small concentrations of hydrocarbons and other organic compounds are apparently possible to notice from the IR spectrum of a distant object.

Actually, they didn’t have any spectral data yet, which could be used to make conclusions about the chemical composition of the atmospheres (that will be obtained later), but they had information about the masses, orbit radii and cross-sectional areas of the planets (these can be deduced from the oscillation of the star’s apparent position, see link “reduced mass“, and from the dimming of the light of the star when a planet temporarily moves past it), and the planet densities were too large for them to be gas-planets. Also, the temperatures of three of them seemed to be in the range where water can exist as a liquid at reasonable pressures.

Actually, a planet doesn’t necessarily need to even be in the so called habitable zone to possibly have liquid water – even on Mercury, the hottest planet on our Solar system, there are places that are always in shadow and can actually get quite cold because the rate of heat conduction through large masses of rock is slow compared to the rate of radiative heat loss to space. Therefore liquid water could exist in such places if contained inside some kind of “pocket” that would also protect it from low pressures. Also, the icy surface of Jupiter’s moon Europa can hide liquid water beneath it, because the tidal forces from Jupiter cause viscous (frictional) heating inside the moon, despite it being way too far from the Sun to be inside the definition of “habitable zone“.

So, more interesting stuff is probably coming, as the astronomers obtain infrared spectra from that star system.

Also, see the related discussion in Physics Forums here.

The featured image is from Wikimedia Commons.

## Numerical solution of PDE:s, Part 3: 2D diffusion problem

In the earlier posts related to PDE numerical integration, I showed how to discretize 1-dimensional diffusion or heat conduction equations either by explicit or implicit methods. The 1d model can work well in some situations where the symmetry of a physical system makes the concentration or temperature field practically depends on only one cartesian coordinate and is independent of the position along two other orthogonal coordinate axes.

The 2-dimensional version of the diffusion/heat equation is

assuming that the diffusion is isotropic, i.e. its rate does not depend on direction.

In a discretized description, the function C(x,y,t) would be replaced by a three-index object

assuming that one of the corners of the domain is at the origin. Practically the same can also be done with two indices, as in

where $N_x$ is the number of discrete points in x-direction.

Using the three-index version of the 2D diffusion equation and doing an explicit discretization, we get this kind of a discrete difference equation:

which can be simplified a bit if we have $\Delta x = \Delta y$ .

Below, I have written a sample R-Code program that calculates the evolution of a Gaussian concentration distribution by diffusion, in a domain where the x-interval is 0 < x < 6 and y-interval is 0 < y < 6. The boundary condition is that nothing diffuses through the boundaries of the domain, i.e. the two-dimensional integral of C(x,y,t) over the area $[0,6]\times [0,6]$ does not depend on t.

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

lx <- 6.0 #length of the computational domain in x-direction
ly <- 6.0 #length of the computational domain in y-direction
lt <- 6 #length of the simulation time interval
nx <- 30 #number of discrete lattice points in x-direction
ny <- 30 #number of discrete lattice points in y-direction
nt <- 600 #number of timesteps
dx <- lx/nx #length of one discrete lattice cell in x-direction
dy <- ly/ny #length of one discrete lattice cell in y-direction
dt <- lt/nt #length of timestep

D <- 1.0 #diffusion constant (assumed isotropic)

Conc2d = matrix(nrow=ny,ncol=nx)
DConc2d <- matrix(nrow=ny, ncol=nx) #a vector for the changes in concentration during a timestep
xaxis <- c(0:(nx-1))*dx #the x values corresponding to the discrete lattice points
yaxis <- c(0:(ny-1))*dy #the y values corresponding to the discrete lattice points

kappax <- D*dt/(dx*dx) #a parameter needed in the discretization
kappay <- D*dt/(dy*dy) #a parameter needed in the discretization

for (i in c(1:ny)) {
for (j in c(1:nx)) {
Conc2d[i,j] = exp(-(i*dy-3)*(i*dy-3)-(j*dx-3)*(j*dx-3)) #2D Gaussian initial concentration distribution
DConc2d[i,j] <- 0 #all initial values in DConc vector zeroed
}
}

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

for(k in c(1:nx))
{
Conc2d[1,k] <- Conc2d[k,2] #fluxes through the boundaries of the domain are forced to stay zero
Conc2d[nx,k] <- Conc2d[nx-1,k]
}

for(k in c(1:ny)) {
Conc2d[k,1] <- Conc2d[k,2]
Conc2d[k,nx] <- Conc2d[k,nx-1]
}
for (k in c(2:(ny-1))) {
for (l in c(2:(nx-1))) {
DConc2d[k,l] <- kappax*(Conc2d[k,l-1]-2*Conc2d[k,l]+Conc2d[k,l+1]) + kappay*(Conc2d[k-1,l]-2*Conc2d[k,l]+Conc2d[k+1,l]) #time stepping
}
}

for (k in c(2:(ny-1))) {
for (l in c(2:(nx-1))) {
Conc2d[k,l] <- Conc2d[k,l]+DConc2d[k,l] #add the changes to the vector Conc
}
}
k <- 0
l <- 0

if(j %% 3 == 1) { #make plots of C(x,y) on every third timestep
jpeg(file = paste("plot_",j,".jpg",sep=""))
persp(yaxis,xaxis,Conc2d,zlim=c(0,1))
title(paste("C(x,y) at t =",j*dt))
dev.off()
}
}

To plot the distribution C(x,y) with a color map instead of a 3D surface, you can change the line

persp(yaxis,xaxis,Conc2d,zlim=c(0,1))

to

image(yaxis,xaxis,Conc2d,zlim=c(0,0.3))

The two kinds of graphs are shown below for three different values of t.

Figure 1. Surface plots of the time development of a Gaussian mass or temperature distribution spreading by diffusion.

Figure 2. Time development of a Gaussian mass or temperature distribution spreading by diffusion, plotted with a red-orange-yellow color map.

To obtain an implicit differencing scheme, we need to write the discretized equation as a matrix-vector equation where $C_{i;j;k}$ is obtained from $C_{i;j;k+1}$ by backward time stepping. In the square matrix, there is one row (column) for every lattice point of the discrete coordinate system. Therefore, if there are N points in x-direction and N points in y-direction, then the matrix is an $N^2 \times N^2$ – matrix. From this it’s quite obvious that the computation time increases very quickly when the spatial resolution is increased.

Initially, one may think that the equation corresponding to diffusion in a really coarse $3 \times 3$ grid is the following one:

Where I’ve used the denotation $k_x = \frac{D\Delta t}{(\Delta x)^2}$ and $k_y = \frac{D\Delta t}{(\Delta x)^2}$ The problem with this is that there are unnecessary terms $-k_x$ in here, creating a cyclic boundary condition (mass that is diffusing through the boundary described by line x = L reappears from the boundary on the other side, x = 0. The correct algorithm for assigning the nonzero elements $A_{ij}$ of the matrix A is

1. $A_{ij} = 1 + 2k_x + 2k_y$ , when $i = j$
2. $A_{ij} = -k_x$ , when $j=i-1$ AND $i\neq 1$ (modulo $N_x$)
3. $A_{ij} = -k_x$ , when $j=i+1$ AND $i\neq 0$ (modulo $N_x$)
4. $A_{ij} = -k_y$ , when $j=i+N_x$ OR $j=i-N_x$

when you want to get a boundary condition which ensures that anything that diffuses through the boundaries is lost forever. Then the matrix-vector equation in the case of $3 \times 3$ lattice is

Unlike the linear system in the 1D diffusion time stepping, this is not a tridiagonal problem and consequently is slower to solve. An R-Code that produces a series of images of the diffusion process for a Gaussian concentration distribution in a $15 \times 15$ discrete lattice is given below.

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

lx <- 6.0 #length of the computational domain in x-direction
ly <- 6.0 #length of the computational domain in y-direction
lt <- 6 #length of the simulation time interval
nx <- 15 #number of discrete lattice points in x-direction
ny <- 15 #number of discrete lattice points in y-direction
nt <- 180 #number of timesteps
dx <- lx/nx #length of one discrete lattice cell in x-direction
dy <- ly/ny #length of one discrete lattice cell in y-direction
dt <- lt/nt #length of timestep

D <- 1.0 #diffusion constant

C = c(1:(nx*ny))
Cu = c(1:(nx*ny))
Conc2d = matrix(nrow=ny,ncol=nx)
xaxis <- c(0:(nx-1))*dx #the x values corresponding to the discrete lattice points
yaxis <- c(0:(ny-1))*dy #the y values corresponding to the discrete lattice points

kappax <- D*dt/(dx*dx) #a parameter needed in the discretization
kappay <- D*dt/(dy*dy) #a parameter needed in the discretization

A = matrix(nrow=(nx*ny),ncol=(nx*ny))

for(i in c(1:(nx*ny))) {
for(j in c(1:(nx*ny))) {
A[i,j] <- 0
if(i==j) A[i,j] <- 1+2*kappax+2*kappay
if(j==i+1 && (i%%nx != 0)) A[i,j] <- -kappax
if(j==i-1 && (i%%nx != 1)) A[i,j] <- -kappax
if(j==i+nx) A[i,j] <- -kappay
if(j==i-nx) A[i,j] <- -kappay
}
}

for (i in c(1:ny)) {
for (j in c(1:nx)) {
Conc2d[i,j] <- exp(-2*(i*dy-3)*(i*dy-3)-2*(j*dx-3)*(j*dx-3)) #2D Gaussian initial concentration distribution
}
}

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

for(k in c(1:ny)) {
for(l in c(1:nx)) {
C[(k-1)*nx+l] <- Conc2d[k,l]
}
}

Cu <- solve(A,C)

for(k in c(1:ny)) {
for(l in c(1:nx)) {
Conc2d[k,l] <- Cu[(k-1)*nx+l]
}
}

k <- 0
l <- 0

if(j %% 3 == 1) { #make plots of C(x,y) on every third timestep
jpeg(file = paste("plot_",j,".jpg",sep=""))
persp(y=xaxis,x=yaxis,z=Conc2d,zlim=c(0,1))
title(paste("C(x) at t =",j*dt))
dev.off()
}
}

Note that the function C(x,y) approaches a constant zero function as t increases, try modifying the code yourself to make a boundary condition that doesn’t let anything diffuse through the boundaries!

## Why does quantum mechanics need an infinite dimensional vector space?

When undergraduate physics students are trying to learn quantum mechanics, one of the most difficult obstacles is that abstract infinite dimensional vector spaces are needed, which makes it impossible for a person to mentally “see” the geometry in their minds in the way a three-dimensional space can be imagined. The fact that the coordinates of the points in that vector space are complex numbers certainly doesn’t make it any easier. When I was studying that myself in the early undergraduate studies, it took me over a year to get a hold of the concept that any self-adjoint operator acting on vectors in the quantum state space has a set of eigenvectors that form a basis of that space.

Sometimes, however, a quantum mechanics problem can be handled by approximating the state space with a finite-dimensional vector space. In fact, several physics Nobel prizes have been awarded for the study of two-state quantum systems. One example of a situation, where a quantum system is effectively a two-state system, is something where the ground state (n=0) and the first excited state (n=1) have energies that are close to each other, and the energy gap between the first and second (n = 2) excited states is unusually wide. Then, if the system is at low temperature, the Boltzmann factors of the lowest-lying two states,

and

are much larger numbers that the Boltzmann factors of any other states, and it is unlikely that a measurement of total energy will give a result other than one of the two lowest energy eigenvalues.

How, then, does one know that a quantum state space actually has to be infinite-dimensional? I once saw, in a mathematical physics textbook, a practice problem where this had to be proved, and I was proud for almost immediately finding the correct answer.

A cornerstone of nonrelativistic quantum mechanics is the position-momentum commutation relation

Now, if the quantum state space were finite-dimensional with dimension $N$, then these operators $\hat{x}$ and $\hat{p}$ would be $N\times N$ square matrices, and the constant $i\hbar$ on the RHS would be a diagonal matrix where all the diagonal elements have value $i\hbar$

Actually, this equation can’t hold in the finite-dimensional case. To see this, consider the calculation rules for taking the trace of a matrix sum or a matrix product:

Now, if you take the trace of both sides of the commutation relation equation, and use these rules, you will get the result

Which clearly isn’t possible. So, we have done a false assumption somewhere in the derivation of the contradictory equation above, and the false assumption must be the finite-dimensionality of the space where $x$ and $p$ are defined. An infinite dimensional operator doesn’t always have a trace at all, and this is the case with the operator $xp-px$.

An interested reader may try to form finite-dimensional matrix representations for the raising and lowering operators of a quantum harmonic oscillator, use those those to form matrices for operators $\hat{x}$ and $\hat{p}$, and see what happens when you calculate their commutator.

## A continuity curiosity

It’s sunday, and I don’t feel like writing lots of source code and posting it here today, so I’ll make a shorter post this time…

In physics and engineering, functions describing physical quantities are often classified as piecewise continuous, everywhere continuous, differentiable, continuously differentiable or complex analytical, which is a list of requirements of increasing strictness. In classical mechanics, most functions are continuous and “well-behaved”, excluding some fractal phase-space curves in chaotic systems.

In pure mathematics, there often happens some kind of play with badly discontinuous or in some other way “detached from reality” functions. Let me show one example of this kind of a function.

The continuity of a real function $f:\mathbb{R}\rightarrow\mathbb{R}$ at point $x_0$ is defined by

or in other words

which practically means that if we set any error bars $f(x_o)-\epsilon \leq f(x) \leq f(x_0)+\epsilon$ for $f(x)$, no matter how strict, there are always some error bars $x_0 - \delta \leq x \leq x_0 + \delta$ for $x$ that will keep the function value between those error bars.

So, here’s an example curious function that is special in the sense that it is continuous at only one point:

which is a real function that has value x when x is a rational number and has value -x when x is an irrational number (the notation $\mathbb{Q}^C$ means the complement of $\mathbb{Q}$).

Left as an exercise for the reader, is to explain how this function is continuous only at the point x=0.

Sometimes later, I’ll write about fractal curves, which are continuous everywhere but are nowhere differentiable and have no finite curve length.

## Numerical solution of PDE:s, Part 2: Implicit method

In the previous blog post, I showed how to solve the diffusion equation

using the explicit method, where the equation is converted to a discrete one

which can be simplified by using the notation

The problem with this approach is that if we want to have a high resolution, i.e. $\Delta x$ is very small, the timestep $\Delta t$ has to be made even much smaller to keep the procedure numerically stable. The stability of this kind of calculations can be investigated with Von Neumann stability analysis, and there you will find out that the instability acts by amplifying the short-wavelength Fourier components of the numerical discretization error.

The specific stability condition for the explicit solution of diffusion equation is

In the better method to solve the diffusion equation, the implicit method, we will not solve the numbers $f_{i;j+1}$ in a straightforward way from the numbers $f_{i;j}$. Instead, we use backward time stepping to write the $f_{i;j}$ as a function of the numbers $f_{i-1;j+1}$, $f_{i;j+1}$ and $f_{i+1;j+1}$, as in the equation below:

which represents a linear system of equations. More specifically, it is a tridiagonal system, and in matrix-vector form it reads

for the case of a relatively small mesh $n_x = 7$. So, now we have to solve this tridiagonal system on each timestep, and this take more computation time than the explicit method but has the advantage of making the calculation stable even when $\Delta t$ is not necessarily made much smaller than $\Delta x$.

A code for solving this kind of a linear system can be found from the book “Numerical Recipes for C” or from the corresponding book written for FORTRAN.

Now, let’s use the implicit method to solve a diffusion problem where the x-domain is

$x \in [0,6]$,

and the step sizes are

$\Delta x = 0.01$ and $\Delta t = 0.05$ .

The initial concentration profile is chosen to be

(don’t be confused by the fact that the function is now called C instead of f) and we use a boundary condition that forces the value of $C(x,t)$ at the left endpoint to be $C(0,t)=1$ and that at the right endpoint to be $C(6,t) = 0$. This means that the left boundary is an infinite source of concentration (or heat in the case of conduction problem) and the right boundary is an infinite “sink”. With physical intuition, it is easy to deduce that this kind of system evolves toward a steady state where the value of $C(x)$ decreases linearly from 1 to 0 on the interval $[0,6]$

A C++ code for doing this calculation is shown below.

// This program calculates the time development of a diffusion or heat conduction
// system with implicit finite differencing. The system described is the development of a concentration or temperature field
// between boundary points where it is constrained to stay at different constant values.
// Teemu Isojärvi, Feb 2017

#include
#include

using namespace std;

#define LX 6. // Length of spatial domain
#define NX 600 // Number of lattice points
#define LT 3. // Length of time interval
#define NT 60 // Number of timesteps

#define D 1. // Diffusion coefficient

int main(void)
{

double dx = (double)LX/(double)NX; // Lattice spacing
double dt = (double)LT/(double)NT; // Length of timestep

double c[NX]; // Concentration values at lattice points

double x; // Auxiliary position variable

double kappa = D*dt/(dx*dx); // Auxiliary variable for representing the linear system

double g[NX];
double b;
double q;

double u[NX]; // Vector for storing the solution on each timestep

for(int m = 0; m<NX; m++)
{
x = (double)m*dx;
c[m]=exp(-3*x*x); // Gaussian initial concentration
}

for(int n = 0; n<NT; n++)
{
c[0]=1;
c[NX-1]=0;

u[0]=(c[0]+kappa)/(2*kappa+1);
q = 2*kappa + 1;

for(int m = 1; m < NX; m++) { // First loop for solving the tridiagonal system g[m]=-kappa/q; q = (2*kappa + 1) + kappa*g[m]; u[m]=(c[m]+kappa*u[m-1])/q; } for(int m=(NX-2); m>=0; m--) u[m] -= g[m+1]*u[m+1]; // Second loop

for(int m=0; m<NX; m++) c[m] = u[m]; // Updating the concentration or temperature field

}

for(int m = 0; m<NX; m++)
{
x = (double)m*dx;
cout << x << " " << c[m] << "\n"; // Output with the results at time t = LT
}

return 0;
}

Running this program three times, with domain lenght parameters $L=0.5$, $L=1.0$ and $L=3.0$ and keeping the time step constant, we get data points that can be plotted in the same coordinate system with a graphing program, like below:

Figure 1. Time evolution of a concentration field C(x,t) in a system where the concentration is forced to stay at constant values at the endpoints of the domain.

The simulation seems to proceed as expected, approaching a linearly decreasing function $C(x)$

An equivalent code for FORTRAN is in the next box:

! Calculates the time development of a concentration distribution C(x,t) with implicit
! finite-differencing of a diffusion equation. The boundary condition is that the left boundary is an infinite source
! of solute/heat and the value of C(x) at x=0 stays at constant value 1. The value of C at the right boundary stays zero.
! Therefore, the function C(x) evolves towards a linearly decreasing function.
! Teemu Isojärvi, Feb 2017

PROGRAM MAIN

real :: DX, DT, LX, LT, D, KAPPA,B,Q ! Real variables for discretization and the diffusion constant
integer :: NT,NX ! Number of time and position steps

REAL :: C(600) ! Concentration field as an array of data points, dimension same as value of NX
REAL :: G(600)
REAL :: U(600)

REAL :: X ! Auxiliary variable

INTEGER :: M ! Integer looping variables
INTEGER :: N

LX = 6. ! Length of x-interval
LT = 3.0 ! Length of t-interval
NX = 600 ! Number of lattice points
NT = 300 ! Number of time steps
DX = LX/NX ! Distance between neighboring lattice points
DT = LT/NT ! Length of time step
D = 1. ! Diffusion/heat conduction coefficient
KAPPA = D*DT/(DX*DX)

do M = 1, NX ! Initial values of concentration
X = M*DX
C(M) = EXP(-3*X*X) ! Gaussian initial concentration distribution
end do

do N = 1, NT ! Time stepping loop

C(1)=1
C(NX)=0

U(1)=(C(1) + 2 * KAPPA)/(2 * KAPPA + 1)
Q = 2 * KAPPA + 1

do M = 2, NX ! First loop for solving the tridiagonal system
G(M) = -KAPPA/Q
Q = (2 * KAPPA + 1) + KAPPA * G(M)
U(M) = (C(M) + KAPPA * U(M-1))/Q
end do

do M = (NX-1), 1
U(M) = U(M) - G(M+1) * U(M+1) ! Second loop
end do

do M = 1, NX-1
C(M) = U(M) ! Updating the concentration or temperature field
end do

end do

do M = 1, NX
X = M*DX
print *,X,C(M) ! Print the x and concentration values at data points
end do

END

To test the implicit method with R-Code, let’s solve a problem where the length of the x-domain is 20, the time interval has length 3 and the initial distribution is

to ensure that we can set the boundary conditions $C(0,t)=C(L,t)=0$ (the Gaussian distribution doesn’t have time to spread all the way to the boundaries in a time interval of 3 units). The code for this calculation is shown below:

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

lx <- 20.0 #length of the computational domain
lt <- 3. #length of the simulation time interval
nx <- 4000 #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

D <- 1.0 #diffusion constant

Conc <- c(1:nx) #define a vector to put the x- and t-dependent concentration values in

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

kappa <- D*dt/(dx*dx) #a parameter needed in the discretization

offdiagonal <- rep(-kappa, times = nx-1)
ondiagonal <- rep(1+2*kappa, times = nx)

for (i in c(1:nx)) {
Conc[i] <- exp(-2*(i*dx-10)*(i*dx-10)) #a Gaussian initial concentration field

xaxis[i] <- i*dx #calculate the x coordinates of lattice points
}

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

sol <- Solve.tridiag(offdiagonal, ondiagonal, offdiagonal, Conc)

for (k in c(1:nx)) {
Conc[k] <- sol[k]
}

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

An animation of the results can be viewed in this link. If you test the code yourself, remember to install the limSolve package first, by writing the command

install.packages(“limSolve”)

in the R console. If you don’t want to load the video, here are the plots for 3 different values of t:

When solving PDE:s other than the ordinary diffusion equation, the implicit method is often even more necessary that it was in the examples here. For example, when solving the time-development of a quantum wave packet from the time-dependent Schroedinger equation, the solution function doesn’t stay normalized if one tries to do a simple explicit time stepping. The solution of Schroedinger equations with implicit method will be shown in the next post. The TDSE is a complex diffusion equation of the form

where the equation has been simplified by setting the Planck constant to value $2\pi$ and the particle mass to $m=1$.

The 1D diffusion and Schroedinger equations are simple in the sense that the linear system to be solved on each timestep is tridiagonal, with makes the solution considerably faster than is the case with a general linear system. Systems that are not tridiagonal will appear when one wants to solve equations with more than one space coordinate (i.e. x and y instead of just x), or when the equation contains higher than 2nd order derivatives with respect to the position coordinate (this can be demonstrated with the thin-film equation, which is fourth-order with respect to x).

## Numerical solution of PDE:s, Part 1: Explicit method

In this post and a few ones following it, I will show how to write C++, Fortran or R-Code programs that solve partial differential equations like diffusion equation or Schroedinger equation numerically, either by explicit or implicit method. An explicit method is a very simple direct calculation, while the implicit (Crank-Nicolson) method involves solving a linear system of equations on every discrete timestep.

A simple example of a PDE is the Poisson equation

,

which describes things like the electric potential field produced by a charge distribution which is the function $\rho$ multiplied by a constant, and can be solved if the function $\rho$ and the boundary values of the function $f$ on some closed surface are known.

In the case of these posts, the equations that I will handle are time-evolution equations, where we are solving a function $f(x,t)$, where the variables are 1D position and time, and the boundary conditions are the function $f(x,0)$ (function $f$ at an initial moment $t=0$) and some conditions for $f(0,t)$ and $f(L,t)$, where $x=0$ and $x=L$ are the left and right boundaries of a computational domain.

The equation that describes both diffusion (spreading of a dissolved substance in solution, or spreading of an unevenly distributed component in gas phase) and heat conduction in a 1-dimensional domain is of the form

.

Here $f(x,t)$ is a function that gives solute concentration or temperature at point $x$ at time $t$, and $D$ is a constant that tells how quickly the diffusion or conduction takes place in medium.

To solve this kind of an equation numerically, we have to choose some domains in time and space

,

and define the object $f_{i;j}$, where the $i$ and $j$ are discrete indices and there is an approximate correspondence

.

Here $\Delta x$ is the spatial discretization step and $\Delta t$ is the timestep. The initial condition means that we know the numbers $f_{i;0}$ for all values of $i$. This is equivalent to knowing the function $f(x,0)$.

The most simple way of converting the PDE into a difference equation for the object $f_{i;j}$ is to use straightforward finite differencing for the derivatives in $x$ and $t$:

With these substitutions and some rearrangement, the PDE becomes

So, now we know what kind of numbers to add to the values $f_{i;j}$ for $i = 0,1,2,\dots ,n_x - 1$ at timestep $j$ to get the numbers $f_{i;j+1}$ for $i = 0,1,2,\dots ,n_x - 1$. Note that I now call the total number of discrete x-axis points $n_x$.

Next, as an example I will solve this discretized equation with the domains and step sizes

$L = 6$
$n_x = 30$
$t_{end} = 6$
$n_t = 600$
$D = 1$

and an initial condition

i.e. a Gaussian temperature/concentration profile. I will also set a boundary condition for the position derivative of $f(x,t)$ at the domain boundaries:

which basically means that the endpoints are impenetrable walls that don’t let anything to diffuse/conduct through them. In other words, the system can be described as “closed” or “adiabatic” depending on whether it’s a diffusion or conduction system.

An example C++ code to do this numerical calculation and to print the values $f_{i;j}$ at the final time $j = n_t$ is shown below:

// This program calculates the time development of a diffusion or heat conduction
// system with explicit finite differencing.
// Teemu Isojärvi, Feb 2017

#include <math.h>;
#include <iostream>;

using namespace std;

#define LX 6. // Length of spatial domain
#define NX 30 // Number of lattice points
#define LT 6 // Length of time interval
#define NT 600 // Number of timesteps

#define D 1. // Diffusion coefficient

int main(void)
{

double dx = (double)LX/(double)NX; // Lattice spacing
double dt = (double)LT/(double)NT; // Length of timestep

double c[NX]; // Concentration values at lattice points
double dc[NX]; // Change of concentration during timestep

double x; // Auxiliary position variable

for(int m = 0; m<NX; m++)
{
x = ((double)m+0.5)*dx;
c[m]=exp(-2*(x-3)*(x-3)); // Gaussian initial concentration
}

for(int n = 0; n<NT; n++)
{
c[0]=c[1]; // Derivative zeroed at left boundary point
c[NX-1] = c[NX-2]; // Derivative zeroed at right boundary point

for(int m = 0; m<NX; m++)
{
dc[m] = (D*dt/(dx*dx))*(c[m+1]-2*c[m]+c[m-1]); // Main finite differencing
}

for(int m = 1; m<NX-1; m++) c[m] += dc[m]; // Updating the concentration
}

for(int m = 0; m<NX; m++)
{
x = (double)m*dx;
cout << x << " " << c[m] << "\n"; // Output with the results at time t = LT
}

return 0;
}

The code can easily be edited to use different initial or boundary conditions, or different discretization step sizes or domain intervals. Doing the calculation three times with time interval lenghths L=0.5, L=1, and L=6, and plotting the output data in a single graph with Grace or similar free graphing program, the result is Fig 1.

Figure 1. Diffusive time evolution of a Gaussian temperature/concentration distribution

To make this kind of graphs, you need to run the program from command line and direct the standard output into a file – in Linux this is done with a command like

./diffusion > out.txt

and in Windows command line it’s done as

diffusion.exe > out.txt

The same code is also easy to write with FORTRAN. This code is shown below.

! Calculates the time development of an initially uneven concentration distribution C(x,t) with explicit
! finite-differencing of a diffusion equation. The boundary condition is that no solute diffuses through the boundaries
! Teemu Isojärvi, Feb 2017

PROGRAM MAIN

real :: DX, DT, LX, LT, D ! Real variables for discretization and the diffusion constant
integer :: NT,NX ! Number of time and position steps

REAL :: C(30) ! Concentration field as an array of data points, dimension same as value of NX
REAL :: DC(30) ! Change of the above quantity during a time step

REAL :: X ! Auxiliary variable

INTEGER :: M ! Integer looping variables
INTEGER :: N

LX = 6. ! Length of x-interval
LT = 6. ! Length of t-interval
NX = 30 ! Number of lattice points
NT = 600 ! Number of time steps
DX = LX/NX ! Distance between neighboring lattice points
DT = LT/NT ! Length of time step
D = 1 ! Diffusion/heat conduction coefficient

do M = 1, NX ! Initial values of concentration
X = M*DX
C(M) = EXP(-2*(X-3)*(X-3)) ! Gaussian initial concentration distribution
end do

do N = 1, NT ! Time stepping loop

C(1)=C(2)
C(NX)=C(NX-1)

do M = 2, NX-1
DC(M) = (D*DT/(DX*DX))*(C(M+1)-2*C(M)+C(M-1)) ! Main time stepping calculation
end do

do M = 1, NX
C(M) = C(M)+DC(M) ! Update the data points
end do

end do

do M = 1, NX
X = M*DX
print *,X,C(M) ! Print the x and concentration values at data points
end do

END

Finally, I will show a way to write, with R-Code programming language, a program that will do the same calculation and automatically produce output files with plots of the function $f(x,t)$ on every one timestep out of 3:

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

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

D <- 1.0 #diffusion constant

Conc <- c(1:nx) #define a vector to put the x- and t-dependent concentration values in
DConc <- c(1:nx) #a vector for the changes in concentration during a timestep
xaxis <- c(1:nx) #the x values corresponding to the discrete lattice points

kappa <- D*dt/(dx*dx) #a parameter needed in the discretization

for (i in c(1:nx)) {
Conc[i] = exp(-2*(i*dx-3)*(i*dx-3)) #Gaussian initial concentration distribution
DConc[i] <- 0 #all initial values in DConc vector zeroed
xaxis[i] <- i*dx #calculate the x coordinates of lattice points
}

for (j in c(1:nt)) { #main time stepping loop
Conc[1] <- Conc[2] #derivatives of C(x) at both ends of domain forced to stay zero
Conc[nx] <- Conc[nx-1]

for (k in c(2:(nx-1))) {
DConc[k] <- kappa*(Conc[k-1]-2*Conc[k]+Conc[k+1]) #calculate the changes in concentration
}

for (l in c(1:nx)) {
Conc[l] <- Conc[l]+DConc[l] #add the changes to the vector Conc
}
k <- 0
l <- 0
if(j %% 3 == 0) { #make plots of C(x) on every third timestep
jpeg(file = paste("plot_",j,".jpg",sep=""))
plot(xaxis,Conc,xlab="position (x)", ylab="concentration",ylim=c(0:1))
title(paste("C(x) at t =",j*dt))
lines(xaxis,Conc)
dev.off()
}


So, the output image files are named as plot_1.jpg, plot_2.jpg, and so on. Now we can use the free ImageJ program (or something else with the same functionality) to import this sequence of images and to save it as an AVI video file. The video can be viewed here.

As an another example, I will give an R-Code that calculates the diffusion or conduction through a layer of intermediate material, given that the function $f(x,t)$ is set to have a constant nonzero value at $x=0$ and the boundary condition that its derivative with respect to x is zero at the domain endpoints:

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

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

D <- 1.0 #diffusion constant

Conc <- c(1:nx) #define a vector to put the x- and t-dependent concentration values in
DConc <- c(1:nx) #a vector for the changes in concentration during a timestep
xaxis <- c(1:nx) #the x values corresponding to the discrete lattice points

kappa <- D*dt/(dx*dx) #a parameter needed in the discretization

for (i in c(1:nx)) {
Conc[i] <- 0.5*(1+tanh(5-5*i*dx)) #a hyperbolic tangent initial concentration field
DConc[i] <- 0 #all elements zeroed initially
xaxis[i] <- i*dx #calculate the x coordinates of lattice points
}

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

Conc[1] <- 1 #force the C to stay at value 1 at left boundary
Conc[2] <- 1 #force the derivative of C to zero at left boundary
Conc[nx] <- Conc[nx-1] #force the derivative of C to zero at right boundary

for (k in c(2:(nx-1))) {
DConc[k] <- kappa*(Conc[k-1]-2*Conc[k]+Conc[k+1]) #calculate the changes in concentration
}

for (l in c(2:(nx-1))) {
Conc[l] <- Conc[l]+DConc[l] #add the changes to the vector Conc
}

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

The animation made of this simulation can be viewed here

So, here were the first examples of PDE solving with discretization. The problem with the explicit finite differencing used here, is that if one decides to use a very small spatial step size $\Delta x$, the timestep $\Delta t$ may have to be made prohibitively small to keep the simulation numerically stable. Suppose we do the simulation with parameters $L = 6$, $n_x = 300$, $t_{end} = 1$, $n_t = 100$. Now the end result plotted with Grace looks like this:

Figure 2. Unstable behavior of explicit method discretization of diffusion equation.

So obviously the numerical errors done in the discretization accumulate exponentially and “blow up” when trying to use too short a $\Delta x$ compared to $\Delta t$. In later blog posts I will show how to use these programming languages to write an implicit PDE solver that corrects this problem for the diffusion/conduction equation, as well as makes it possible to solve Schroedinger equation (a complex-valued version of diffusion equation, which can have wave-like solutions with reflection, interference, etc.). Also, I will extend the method to problems with more than one spatial dimension, and show how to solve a nonlinear PDE called thin-film equation.

## Linear algebra for balancing chemical reaction equations

The balancing of chemical reaction equations is familiar for many, either from high school chemistry or introductory university chemistry. A simple example of such a problem is the reaction equation for the combustion of methane in oxygen:

$aCH_4 + bO_2 \rightarrow cCO_2 + dH_2 O$,

for which we are supposed to find coefficients a,b,c and d that make the reaction possible in terms of conservation of the number of atoms of any element that is present. In this case, it is rather easy to see that the number of $CO_2$ molecules must be equal to the number of $CH_4$ molecules, and that the number of water molecules must be twice that. Therefore the possible equations are

$CH_4 + 2O_2 \rightarrow CO_2 + 2H_2 O$ , or

$2CH_4 + 4O_2 \rightarrow 2CO_2 + 4H_2 O$, or

$3CH_4 + 6O_2 \rightarrow 3CO_2 + 6H_2 O$ ,

or any other “multiple” of the first equation (obtained by multiplying all the coefficients in the equation with the same number).

Usually this kind of balancing is done mainly by intuition and trying to guess the right coefficients, but this can be difficult if the reactants and product contain many different elements, each of which are present in more than one kind of molecule. One example of a more complicated (but still quite simple) reaction is the reaction of chlorine dioxide ($ClO_2$) with alkaline water to form chlorate ($ClO_3^-$), chlorite ($ClO_2^-$) and water:

$aClO_2 + bOH^- \rightarrow cClO_2^- + dClO_3^- + eH_2 O$ .

Now, an easily applicable way to solve a set of acceptable coefficients a,b,c,d and e is to explicitly write the balancing condition equations for each element that appears in the equation. These balancing conditions are:

$a = c + d$                                 (balancing of chlorine)
$2a + b = 2c + 3d + e$             (balancing of oxygen)
$b = 2e$                                   (balancing of hydrogen)
$b = c + d$                                (balancing of electric charge)

This can be solved with many kinds of mathematical software, either free (there’s lots of free math programs especially for Linux) or commercial (like Wolfram Mathematica).

But there’s one another thing. Because the number of unknowns in the system of equations above is larger than the number of equations by one, we need to decide the value of one of the coefficients $a, \dots , b$ arbitrarily. So let’s say that $a = 1$. Now the set of equations is:

$a = 1$                                   (arbitrary choice of one coeffiecient)
$a = c + d$                            (balancing of chlorine)
$2a + b = 2c + 3d + e$          (balancing of oxygen)
$b = 2e$                              (balancing of hydrogen)
$b = c + d$                           (balancing of electric charge)

Using the column vector

the system of equations can be written in matrix form

This can be solved in Mathematica by writing the command

Solve[{a == 1, a == c + d, 2a + b == 2*c + 3*d + e, b == 2e, b == c + d}]

and the result is

$a = 2$,

$b = 2$,

$c = 1$,

$d = 1$,

$e = 1$,
which leads to the reaction equation

$2ClO_2 + 2OH^- \rightarrow ClO_2^- + ClO_3^- + H_2 O$

which doesn’t even need to be multiplied with any constant to make all the coefficients integers (usually that is necessary).
IMAGE SOURCE: The feature image of an Erlenmeyer flask on this blog post is a public domain image from Wikipedia, address https://en.wikipedia.org/wiki/Erlenmeyer_flask#/media/File:Duran_erlenmeyer_flask_narrow_neck_50ml.jpg .