General Notation

For large systems of equations, an iterative solution scheme for the unknown vector can always be written in the form

where is the iteration matrix, is a constant vector and p is an iteration counter. Convergence of this scheme is guaranteed if the largest eigenvalue of the iteration matrix is less that unity, where

= spectral radius =

Therefore, if < 1 the iterative scheme will converge. If << 1, the iterative scheme converges very rapidly. If but less than unity, the scheme will be slowly converging. The iteration algorithm will diverge if the spectral radius is greater than unity.

Convergence is tested during the iterative process by computing the largest relative change from one iteration to the next, and comparing the absolute value of this result with some desired tolerance. If the maximum relative change is less than the desired accuracy, then the process is terminated. If this condition is not satisfied, then another iteration is performed.

The Gauss Seidel Method

Let’s take the original system of equations given by and convert it into the classical Gauss Seidel iterative scheme. To do this, let’s break the original matrix into three specific components, or

where the three matrices on the right hand side, in sequence, are strictly lower triangular, diagonal, and strictly upper triangular matrices.

Now, substituting this into the original balance expression gives

or

If we premultiply by and notice that the solution vector appears on both sides of the equation, we can write the equation in an iterative form as

Clearly this is in the standard form for iterative solutions as defined in eqn. (6.9), where the iteration matrix is given by

and the constant vector is written as

This form of the iteration strategy is useful for the study of the convergence properties of model problems. It is, however, not particularly useful as a program algorithm for code implementation.

For actual implementation on the computer, one writes these equations differently, never having to formally take the inverse as indicated above. In practice, eqn. (6.12) is written in iterative form as

or

This specific form is somewhat odd at first glance, since appears on both sides of the equation. This is justified because of the special form of the strictly lower triangular matrix, . This can be seen more clearly if the matrix equations are written using discrete notation.

In discrete form eqn. (6.16) can be expanded as

where the diagonal elements of are simply and the limits associated with the summations account for the special structure of the matrices.

As a specific example, consider the 3x3 system shown below:

The iterative equations for this case can be written individually as

Thus, if the equations are taken in sequence, all the terms on the right hand side are known, thereby allowing computation of a new estimate for the full solution vector in terms of the very latest information available for the calculation. This is the basic idea behind the Gauss Seidel method.

The Successive Over-Relaxation (SOR) Method

To improve the rate of convergence, one might consider using a weighted average of the results of the two most recent estimates to obtain the next best guess of the solution. If the solution is converging, this might help extrapolate to the real solution more quickly. This idea is the basis of the SOR method.

In particular, let be some weight factor with a value between 0 and 2. Now, let’s compute the next value of to use in the Gauss Seidel method as a linear combination of the current value, , and the previous solution, , as follows:

Note that if is unity, we simply get the standard Gauss Seidel method (or whatever base iterative scheme is in use). When is greater that unity, the system is said to be over-relaxed, indicating that the latest value, , is being weighted more heavily (weight for is negative). If, however, is less than one, the system is under-relaxed, this time indicating that the previous solution, , is more heavily weighted (positive weight values).

The idea, of course, is to choose the relaxation parameter to improve convergence (reduce the spectral radius). This is most often done in a trial and error fashion for certain classes of problems (experience helps here). Some more advanced codes do try to estimate this quantity as part of the iterative calculation, although this is not particularly easy.

A variety of schemes for improving convergence have been developed over the years, with many taking advantage of the particular structure of the algebraic equations or some characteristic of the physical system under study. The specifics of these algorithms are beyond the scope of this course, but the interested student is encouraged to see a good text on advanced numerical methods for further details as desired.