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

(Eq. 1)

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

(Eq. 2)

If we substitute this into the difference equation, we get

(Eq. 3)

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.

(Eq. 4)

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

(Eq. 5)

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

(Eq. 6)

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

(Eq. 7)

so that the general form of the round off error is

(Eq. 8)

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

(Eq. 9)

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

(Eq. 10)

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

(Eq. 11)

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

(Eq. 12)

to obtain

(Eq. 13)

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

(Eq. 14)

we get two cases

Case I:

(Eq. 15)

or

(Eq. 16)

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

Case II:

(Eq. 17)

or

(Eq. 18)

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

(Eq. 19)

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

(Eq. 20)

or

(Eq. 21)

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

(Eq. 22)

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.

REFERENCES

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. http://pleasemakeanote.blogspot.com/2008/06/parallel-computing-with-mpi-part-iii.html