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

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

codecogseqn6

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

codecogseqn12

which can be simplified by using the notation

CodeCogsEqn(16).gif

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

CodeCogsEqn(17).gif

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:

CodeCogsEqn(18).gif

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

CodeCogsEqn(26).png

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

CodeCogsEqn(19).gif

(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:

penetration_cpp

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

CodeCogsEqn(20).gif

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:

plot_34.jpgplot_142plot_289

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

CodeCogsEqn(21).gif

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

Save

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

CodeCogsEqn(5).gif ,

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

CodeCogsEqn(6).gif .

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

codecogseqn7 ,

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

CodeCogsEqn(8).gif .

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:

CodeCogsEqn(11).gif

codecogseqn10
With these substitutions and some rearrangement, the PDE becomes

CodeCogsEqn(12).gif

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

CodeCogsEqn(13).gif

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:

CodeCogsEqn(14).gif

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.

diffusion1

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:

unstable_explicit.jpg

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.

Save