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

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 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 is obtained from 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 – 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 grid is the following one:

Where I’ve used the denotation and The problem with this is that there are unnecessary terms 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 of the matrix A is

1. , when

2. , when AND (modulo )

3. , when AND (modulo )

4. , when OR

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