then, by definition,
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
If we substitute this into the difference equation, we get
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
Now we take the absolute value of this equation and obtain the stability requirement
we get two cases
which is always true. So this case does not pose any limitations.
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:
- Introducing the Heat Equation
- Discretization of the Heat Equation
- Stability Analysis and the Courant-Friedrichs-Lewy Criterion
- Programming the Heat Equation in C using Visual Studio
- MPICH2 and Visual Studio
- Parallel Computing in a Nutshell
- Measuring Parallel Performance
- Essential Ingredients of Message Passing
- MPI Warmup
- Parallelizing the Heat Equation