## 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!

## 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.