Lecture 1

The central concern of numerical modeling in fluid-related disciplines is this: "How best can one approximate space-time continuum by discrete numbers?" This is not solely a question of numerical resolution: even if you have a supercomputer that affords world's finest resolution, you can end up with a totally useless result if you discretize the governing equation in a wrong way.

Consider a simple linear damping problem:

for which we know the exact solution: . Let us approximate the left-hand side of the above PDE by finite difference, representing time t with a series of discrete time steps , where is integer and positive constant. There are a number of ways to do this: it may not appear obvious at first what difference the following three schemes bring about:




It turns out that the above three approximations yield very distinct solution behaviors. The first or the forward scheme performs variably depending on the magnitude of (Figure 1). As long as is small, the approximate solution stays close to the exact solution, but beyond unity the solution develops spurious oscillation, and for > 2 the solution diverges with time (numerically unstable). Certainly the usefulness of this scheme hinges on the smallness of .

The second or the leapfrog scheme gives rise to two branches of solutions, viz. physical mode and computational mode. The physical mode gives a fairly accurate approximation to the exact solution provided is small (Figure 2). However, the computational mode is unstable (Figure 3; note different vertical scale). Ironically, the growth rate is greatest when the time resolution is strongest (i.e., small ). The excitation of the unstable computational mode may be delayed by a judicious choice of initial condition (specification of q at t = 0 and t = ), but in practice, the roundoff error of the computer seeds the computational mode and thus it is very difficult to suppress it completely. Because of this, as far as this problem is concerned, the leapfrog is least desirable of the three schemes.

The third or the backward scheme is devoid of numerical instability regardless of the magnitude of (absolutely stable). It is also as accurate as the physical mode of the leapfrog scheme (Figure 4). For this particular problem, this scheme appears to be the winner.


How do we measure accuracy and stability of a scheme? Typically, accuracy is measured by the size of the truncation error, that is, the difference between the left-hand side of the above equations and true time derivative. A finite difference scheme is said to be first-order accurate when the leading order of truncation error is linear with respect to the inclement of differencing (), and second-order accurate when it is quadratic. In the above examples, (1) and (3) are first-order accurate whereas (2) is second-order accurate:



Stability concerns whether the solution is bounded for a duration of time. As such it is a global measure of the solution. Nonetheless, for linear problems with constant coefficients that permit a modal solution of the form , it suffices to examine the amplification factor

to determine the numerical stability: a scheme is numerically stable if

This is the Von Neumann stability analysis. (Actually, when the modulus of the amplification factor is equal to unity in the leapfrog scheme, there can be a non-exponential--algebraic--growth. However, this is a subtle point and has little practical significance.)

It is important to note that accuracy and stability of a particular scheme is problem-specific. For example, consider an oscillatory equation

We can apply similar finite difference approximations to the above:

The Von Neumann stability analysis shows that the amplification factor for this problem is a complex number. That means that the finite difference approximation affects both amplitude and phase of the oscillatory solution. The effects of the three finite difference approximations on the damped and oscillatory solutions are summarized in the following table.

To HomeTo Lecture 2