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 multiplied by a constant, and can be solved if the function and the boundary values of the function 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 , where the variables are 1D position and time, and the boundary conditions are the function (function at an initial moment ) and some conditions for and , where and 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 is a function that gives solute concentration or temperature at point at time , and 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 , where the and are discrete indices and there is an approximate correspondence

.

Here is the spatial discretization step and is the timestep. The initial condition means that we know the numbers for all values of . This is equivalent to knowing the function .

The most simple way of converting the PDE into a difference equation for the object is to use straightforward finite differencing for the derivatives in and :

With these substitutions and some rearrangement, the PDE becomes

So, now we know what kind of numbers to add to the values for at timestep to get the numbers for . Note that I now call the total number of discrete x-axis points .

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

and an initial condition

i.e. a Gaussian temperature/concentration profile. I will also set a boundary condition for the position derivative of 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 at the final time is shown below:

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.

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 on every one timestep out of 3:

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 is set to have a constant nonzero value at and the boundary condition that its derivative with respect to x is zero at the domain endpoints:

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 , the timestep may have to be made prohibitively small to keep the simulation numerically stable. Suppose we do the simulation with parameters , , , . Now the end result plotted with Grace looks like this:

So obviously the numerical errors done in the discretization accumulate exponentially and “blow up” when trying to use too short a compared to . 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.