112
5. Solution of Linear Equation Systems
both 2D and 3D versions) and the 3D CGSTAB solver are available via Internet; see appendix for details. A biconjugate gradient solver with incomplete
Cholesky preconditioning for nine-point schemes in 2D is also provided.
5.3.8 Multigrid Methods
The final method for solving linear systems to be discussed here is the multigrid method. The basis for the multigrid concept is an observation about iterative methods. Their rate of convergence depends on the eigenvalues of the
iteration matrix associated with the method. In particular, the eigenvalue(s)
with largest magnitude (the spectral radius of the matrix) determines how
rapidly the solution is reached; see Sect. 5.3.2. The eigenvector(s) associated
with this eigenvalue(s) determines the spatial distribution of the iteration
error and varies considerably from method to method. Let us briefly review
the behavior of these entities for some of the methods presented above. The
properties are given for Laplace equation; most of them generalize to other
elliptic partial differential equations.
For Laplace equation, the two largest eigenvalues of the Jacobi method
are real and of opposite sign. One eigenvector represents a smooth function of
the spatial coordinates, the other, a rapidly oscillating function. The iteration
error for the Jacobi method is thus a mixture of very smooth and very rough
components; this makes acceleration difficult. On the other hand, the GaussSeidel method has a single real positive largest eigenvalue with an eigenvector
that makes the iteration error a smooth function of the spatial coordinates.
The largest eigenvalues of the SOR method with optimum over-relaxation
factor lie on a circle in the complex plane and there are a number of them;
consequently, the error behaves in a very complicated manner. In ADI, the
nature of the error depends on the parameter but tends to be rather complicated. Finally, SIP has relatively smooth iteration errors.
Some of these methods produce errors that are smooth functions of the
spatial coordinates. Let us consider one of these methods. The iteration error
8' and residual pn after the nth iteration are related by Eq. (5.15). In the
Gauss-Seidel and SIP methods, after a few iterations, the rapidly varying
components of the error have been removed and the error becomes a smooth
function of the spatial coordinates. If the error is smooth, the update (an
approximation to the iteration error) can be computed on a coarser grid.
On a grid twice as coarse as the original one in two dimensions, iterations
cost 114 as much; in three dimensions, the cost is 118 the fine grid cost.
Furthermore, iterative methods converge much faster on coarser grids. GaussSeidel converges four times as fast on a grid twice as coarse; for SIP the ratio
is less favorable but still substantial.
This suggests that much of the work can be done on a coarser grid. To
do this, we need t o define: the relationship between the two grids, the finite
difference operator on the coarse grid, a method of smoothing (restricting)
the residual from the fine grid to the coarse one and a method of interpolating
5. Solution of Linear Equation Systems
both 2D and 3D versions) and the 3D CGSTAB solver are available via Internet; see appendix for details. A biconjugate gradient solver with incomplete
Cholesky preconditioning for nine-point schemes in 2D is also provided.
5.3.8 Multigrid Methods
The final method for solving linear systems to be discussed here is the multigrid method. The basis for the multigrid concept is an observation about iterative methods. Their rate of convergence depends on the eigenvalues of the
iteration matrix associated with the method. In particular, the eigenvalue(s)
with largest magnitude (the spectral radius of the matrix) determines how
rapidly the solution is reached; see Sect. 5.3.2. The eigenvector(s) associated
with this eigenvalue(s) determines the spatial distribution of the iteration
error and varies considerably from method to method. Let us briefly review
the behavior of these entities for some of the methods presented above. The
properties are given for Laplace equation; most of them generalize to other
elliptic partial differential equations.
For Laplace equation, the two largest eigenvalues of the Jacobi method
are real and of opposite sign. One eigenvector represents a smooth function of
the spatial coordinates, the other, a rapidly oscillating function. The iteration
error for the Jacobi method is thus a mixture of very smooth and very rough
components; this makes acceleration difficult. On the other hand, the GaussSeidel method has a single real positive largest eigenvalue with an eigenvector
that makes the iteration error a smooth function of the spatial coordinates.
The largest eigenvalues of the SOR method with optimum over-relaxation
factor lie on a circle in the complex plane and there are a number of them;
consequently, the error behaves in a very complicated manner. In ADI, the
nature of the error depends on the parameter but tends to be rather complicated. Finally, SIP has relatively smooth iteration errors.
Some of these methods produce errors that are smooth functions of the
spatial coordinates. Let us consider one of these methods. The iteration error
8' and residual pn after the nth iteration are related by Eq. (5.15). In the
Gauss-Seidel and SIP methods, after a few iterations, the rapidly varying
components of the error have been removed and the error becomes a smooth
function of the spatial coordinates. If the error is smooth, the update (an
approximation to the iteration error) can be computed on a coarser grid.
On a grid twice as coarse as the original one in two dimensions, iterations
cost 114 as much; in three dimensions, the cost is 118 the fine grid cost.
Furthermore, iterative methods converge much faster on coarser grids. GaussSeidel converges four times as fast on a grid twice as coarse; for SIP the ratio
is less favorable but still substantial.
This suggests that much of the work can be done on a coarser grid. To
do this, we need t o define: the relationship between the two grids, the finite
difference operator on the coarse grid, a method of smoothing (restricting)
the residual from the fine grid to the coarse one and a method of interpolating