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
=