Wednesday, June 18, 2008

Parallel Computing with MPI - Part III: Stability Analysis of the Discretized Heat Equation

In today's part of this series, we'll go over a basic stability analysis for the discretized heat equation that was derived in the previous post.

Here, I will follow John Anderson's neat exposition of the stability analysis for difference equations. We start by defining a couple of terms first. Let,

D = exact solution of the difference equation (infinite accuracy computer)
N = numerical solution of the difference equation (finite accuracy computer)
ε = Round off error

then, by definition,

ε = N - D

which is the error between the exact solution and the finite accuracy numerical solution of the difference equation. Our difference equation can therefore be represented by N since we are most likely to solve it on a finite accuracy computer, right? To make headway, we replace T in our difference equation by N. at the outset, we get

but, the numerical solution can be expressed in terms of the round off error and the exact difference equation solution as

If we substitute this into the difference equation, we get
however, the exact solution D to the difference equation already satisfies the discretized equation. Then, D simply cancels out from the above equation, and therefore, the round off error also satisfies the difference equation, i.e.

Since we are solving using a real computer, then round off errors are inevitable. Therefore, we want to study how these errors evolve from one timestep to another. If ε gets smaller or at worst stays the same, then the solution is stable. Otherwise, it gets larger, then the solution is unstable. Mathematically, for a numerical solution to be stable, this can be expressed as

Now, we wish to seek a general form of the round off error. Thankfully, once upon a time, a mathematician called Joseph Fourier showed that any function can be expressed as an infinite series of sines and cosines, or more conveniently, exponentials. Therefore, it is reasonable to assume the following general form for the round off error

where i is the imaginary unit, kn is a wave number, and An(t) reflects the time dependence of the round off error. One may also assume a general form for An(t). Since the error will either grow or decay with time, then one can set

so that the general form of the round off error is

Now, we can substitute this form into the discretized equation accordingly. For example,

Upon substitution and simplification, i.e. dividing by ε(x,y,t), we get

nice! Notice that the infinite sum is no longer needed. Also, the exponential terms inside the parentheses can be collapsed into cosines

This form, albeit quite simple, does not help in the subsequent steps. We need to invoke one extra trigonometric identity

to obtain

Now we take the absolute value of this equation and obtain the stability requirement

we get two cases

Case I:


which is always true. So this case does not pose any limitations.

Case II:


but, since the square of the sine is always less than or equal to 1, then

going back to the original inequality, for stability, we must satisfy


This is known as the Courant-Friedrichs-Lewy or the CFL condition. Although not in the context in which it was originally derived, a CFL condition may be obtained for almost any numerical analysis of differential equation.

In general, if the grid spacing was not uniform, we would have

Note that in Horak and Gruber's paper (Parallel Numerical Solution of 2-D Heat Equation), there is a typo in the stability condition equation on page 49 (first equation on that page). The "2" in the denominator should actually be a "4" as was shown in this post. The CFL condition that the authors refer to actually corresponds to the one dimensional heat equation.

In my next post, I'll post a sequential program for solving the heat equation in C using visual studio.


1. ANDERSON, J. D. 1995 Computational fluid dynamics : the basics with applications. McGraw-Hill.

2. TANNEHILL, J. C., ANDERSON, D. A. & PLETCHER, R. H. 1997 Computational fluid mechanics and heat transfer. 2nd edn. Taylor & Francis.

Table of Contents:

  1. Introducing the Heat Equation
  2. Discretization of the Heat Equation
  3. Stability Analysis and the Courant-Friedrichs-Lewy Criterion
  4. Programming the Heat Equation in C using Visual Studio
  5. MPICH2 and Visual Studio
  6. Parallel Computing in a Nutshell
  7. Measuring Parallel Performance
  8. Essential Ingredients of Message Passing
  9. MPI Warmup
  10. Parallelizing the Heat Equation

Cite as:
Saad, T. "Parallel Computing with MPI - Part III: Stability Analysis of the Discretized Heat Equation". Weblog entry from Please Make A Note.


  1. i thought you might want to know you have a couple of typos
    keep up the good work

  2. you followed John Anderson's neat exposition of the stability analysis for difference equations.
    Could you provide some references


  3. I thought the stability condition for a one dimensional transient heat conduction problem using the explicit method was r <= 0.5, not r <= 0.25. Is 0.25 correct?