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. is very small, the timestep 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 in a straightforward way from the numbers . Instead, we use backward time stepping to write the as a function of the numbers , and , 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 . 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 is not necessarily made much smaller than .

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

,

and the step sizes are

and .

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 at the left endpoint to be and that at the right endpoint to be . 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 decreases linearly from 1 to 0 on the interval

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 , and 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

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 (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 and the particle mass to .

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