The Multigrid Workbench: Linear Iterations

Iterative methods for a linear system of equations

	Au = f,		(1)

like the Gauß-Seidel method, over-relaxation (SOR), or the Gauß-Jacobi method can all be written in the general form
 

	Buk+1 = f - ( A - B ) uk		(2)

or
 

	uk+1 = uk + B-1 ( f - Auk )		(3)

where the upper index indicates the number of the iterate. B is a matrix characteristic for the method.

Clearly, the exact solution u^*=A^-1f is a fixed point of iteration (2). To check whether the iteration converges u^k=u^*+e^k and u^k+1=u^*+e^k+1 may be inserted into (3). This simplifies to
 

	ek+1 = ( I - B-1A ) ek.		(4)

Convergence is obtained, when the error norm is reduced. Thus
 

	||ek+1 || || I - B-1A || || ek ||.		(5)

assures convergence when || I - B^-1A || < 1. This statement is a precise way of saying that for (3) to converge fast, B must approximate A. Clearly, when B=A this norm vanishes and we trivially get a method which converges in one step. However, this iteration would require the inversion of B, or the solution of a system, where B is the system matrix. This would essentially mean the original system is solved directly, and nothing would be won.

The difficulty is to find the right compromise of having B easily invertible, and being a good approximation to A.

Relaxation methods

Relaxation methods use particularly simple choices of B. For the Gauß-Jacobi method B=diag(A), for the Gauß-Seidel method, B is taken as the lower triangle of A. These B are certainly easily invertible, however, they do not make very good approximations to A, at least when A is large. However, they provide good smoothers for the multigrid method.

Coarse grid correction

Here we now study an alternative method, the so-called coarse grid correction method.

Our idea is to construct B^-1 form a similar grid problem as A, however, a simpler one by taking fewer gridlines. Let us assume that A_C is constructed just like A, but from a grid with only (N-1)x(N-1) gridlines. The dimension of the matrix is therefore less than ¼ of the original dimension. It is therefore obviously easier to invert A_C than A. Unfortunately, however, equation (3) cannot be used any more, because now A_C^-1 does not have a dimension compatible with f, A, and u. To fix this, we introduce the interpolation (prolongation) operator P, and the restriction operator R. R and P are mappings between these two spaces:

	R: |RNxN -> |R(N-1)/2x(N-1)/2		(6)

	P: |R(N-1)/2x(N-1)/2 -> |RNxN		(7)

and define

	BC = P AC-1 R		(8)

Now we replace B^-1 in (3) by B_C. Note that since B_C has rank smaller than A, therefore there exists no formulation of the iteration equivalent to (2).

For the same reason, the product B_CA is singular, and

	|| I - BCA || 1.

Thus the coarse grid correction cannot be a convergent method. However, it is a good iteration with fast convergence for the low frequency solution modes.

Multigrid method

The multigrid method alternates between smoothers and coarse grid corrections. The combination of the two results in a new iterative algorithm, where the different convergence properties of the two basic iteration combine favorably. The two grid algorithm resulting from this idea is extended to a multigrid algorithm by extending this idea recursively to the solution of A^-1.

Other iterative methods

Besides the simple iterative schemes many more possibilities exist. We refer to the chapter on numerical linear algebra in the Computational Science Education Project.

Ulrich Ruede , Thu Feb 2 21:06:09 MEZ 1995
Updated by Craig C. Douglas