Geoscience Reference
In-Depth Information
Figure 4.29 Sketch of LU decomposition of matrix M .
In order to make up the difference between A and M , one choice is to let the matrix
C contain just the two non-zero diagonals of M that correspond to the zero diago-
nals of A . This is the standard incomplete LU decomposition method. However, this
method converges slowly. A better choice, which was proposed by Stone (1968), is
to allow C to have non-zero elements on the diagonals corresponding to all seven
non-zero diagonals of LU . C must contain the 'two extra' diagonals of M , and the ele-
ments on the remaining diagonals of C are chosen to ensure C
0. The resulting
equations relating L and U with A are (see Ferziger and Peric, 1995)
L W , l =
A W , l /(
1
+ α
s U N , l m j )
L S , l
=
A S , l
/(
1
+ α
s U E , l 1
)
L P , l =
A P , l + α
(
L W , l U N , l m j +
L S , l U E , l 1 )
L W , l U E , l m j
L S , l U N , l 1
s
U N , l
= (
A N , l
α
s L W , l U N , l m j )/
L P , l
U E , l = (
A E , l α
s L S , l U E , l 1 )/
L P , l
(4.227)
where L W , l , L S , l , and L P , l are the coefficients of matrix L ; U N , l and U E , l are the
coefficients of matrix U (the coefficients in the main diagonal are set to be 1 to get the
unique solution in LU decomposition);
α s is a coefficient less than 1; l is the index of
points in the matrix A , related to i and j by l
=
j
+ (
i
1
)
m j ; and m j is the number of
points on the j line.
The set of equations (4.227) can be solved in a sequential order beginning at the
southeast corner of the grid. For points next to boundaries, any matrix element that
carries the index of a boundary point is set to be zero.
After L and U are determined, the approximate formulation for Eq. (4.225) is
obtained as
LU
=
R
(4.228)
By defining U
=
, Eq. (4.228) can be split as
L
=
R
(4.229)
U
=
Search WWH ::




Custom Search