104
5. Solution of Linear Equation Systems
The coefficients must be calculated in this order. For nodes next to boundaries, any matrix element that carries the index of a boundary node is understood to be zero. Thus, along the west boundary (2 = 2), elements with
index 1 - N j are zero; along the south boundary ( j = 2), elements with index
1 - 1 are zero; along the north boundary ( j = N j - I ) , elements with index
1 + 1 are zero; finally, along the east boundary (i = Ni - I ) , elements with
index 1 + N j are zero.
We now turn to solving the system of equations with the aid of this
approximate factorization. The equation relating the update to the residual
is (see Eq. (5.19)):
The equations are solved as in in generic LU decomposition. Multiplication
of the above equation by L-' leads to:
This equation is to be solved by marching in the order of increasing 1. When
the computation of R is complete, we need to solve Eq. (5.43):
61 = R~ - ~ ; @ + 1 - u,$@+Nj
(5.45)
in order of decreasing index 1.
In the SIP method, the elements of the matrices L and U need be calculated only once, prior to the first iteration. On subsequent iterations, we
need calculate only the residual, then R and finally 6, by solving the two
triangular systems.
Stone's method usually converges in a small number of iterations. The rate
of convergence can be improved by varying a from iteration to iteration (and
point to point). These methods converge in fewer iterations but they require
the factorization to be redone each time a: is changed. Since computing L and
U is as expensive as an iteration with a given decomposition, it is usually
more efficient overall to keep a: fixed.
Stone's method can be generalized to yield an efficient solver for the ninediagonal matrices that arise when compact difference approximations are
applied in two dimensions and for the seven-diagonal matrices that arise when
central differences are used in three dimensions. A 3D (7-point) vectorized
version is given by Leister and PeriC (1994); two 9-point versions for 2D
problems are described by Schneider and Zedan (1981) and PeriC (1987).
Computer codes for five-diagonal (2D) and seven-diagonal (3D) matrices are
5. Solution of Linear Equation Systems
The coefficients must be calculated in this order. For nodes next to boundaries, any matrix element that carries the index of a boundary node is understood to be zero. Thus, along the west boundary (2 = 2), elements with
index 1 - N j are zero; along the south boundary ( j = 2), elements with index
1 - 1 are zero; along the north boundary ( j = N j - I ) , elements with index
1 + 1 are zero; finally, along the east boundary (i = Ni - I ) , elements with
index 1 + N j are zero.
We now turn to solving the system of equations with the aid of this
approximate factorization. The equation relating the update to the residual
is (see Eq. (5.19)):
The equations are solved as in in generic LU decomposition. Multiplication
of the above equation by L-' leads to:
This equation is to be solved by marching in the order of increasing 1. When
the computation of R is complete, we need to solve Eq. (5.43):
61 = R~ - ~ ; @ + 1 - u,$@+Nj
(5.45)
in order of decreasing index 1.
In the SIP method, the elements of the matrices L and U need be calculated only once, prior to the first iteration. On subsequent iterations, we
need calculate only the residual, then R and finally 6, by solving the two
triangular systems.
Stone's method usually converges in a small number of iterations. The rate
of convergence can be improved by varying a from iteration to iteration (and
point to point). These methods converge in fewer iterations but they require
the factorization to be redone each time a: is changed. Since computing L and
U is as expensive as an iteration with a given decomposition, it is usually
more efficient overall to keep a: fixed.
Stone's method can be generalized to yield an efficient solver for the ninediagonal matrices that arise when compact difference approximations are
applied in two dimensions and for the seven-diagonal matrices that arise when
central differences are used in three dimensions. A 3D (7-point) vectorized
version is given by Leister and PeriC (1994); two 9-point versions for 2D
problems are described by Schneider and Zedan (1981) and PeriC (1987).
Computer codes for five-diagonal (2D) and seven-diagonal (3D) matrices are