6.4 Linear Systems 175
ui j=ai j−i− 1
∑
k= 1likuk j.- Then calculate the diagonal elementuj j
uj j=aj j−j− 1
∑
k= 1ljkuk j. (6.19)- Finally, calculate the elementsli j,i>j
li j=1
uj j(
ai j−i− 1
∑
k= 1likuk j)
, (6.20)
The algorithm is known as Doolittle’s algorithm since the diagonal matrix elements ofLare 1.
For the case where the diagonal elements ofUare 1 , we have what is called Crout’s algorithm.
For the case whereU=LTso thatuii=liifor 1 ≤i≤nwe can use what is called the Cholesky
factorization algorithm. In this case the matrixAhas to fulfill several features; namely, it
should be real, symmetric and positive definite. A matrix is positive definite if the quadratic
formxTAx> 0. Establishing this feature is not easy since it implies the use of an arbitrary
vectorx 6 = 0. If the matrix is positive definite and symmetric, its eigenvalues are always real
and positive. We discuss the Cholesky factorization below.
A crucial point in the LU decomposition is obviously the casewhereuj jis close to or equals
zero, a case which can lead to serious problems. Consider thefollowing simple 2 × 2 example
taken from Ref. [30]
A=(
0 1
1 1
)
.
The algorithm discussed above fails immediately, the first step simple states thatu 11 = 0. We
could change slightly the above matrix by replacing 0 with 10 −^20 resulting in
A=(
10 −^201
1 1
)
,
yielding
u 11 = 10 −^20
l 21 = 1020
andu 12 = 1 and
u 22 =a 11 −l 21 = 1 − 1020 ,
we obtain
L=(
1 0
10201
)
,
and
U=(
10 −^201
0 1 − 1020
)
,
With the change from 0 to a small number like 10 −^20 we see that the LU decomposition is now
stable, but it is not backward stable. What do we mean by that?First we note that the matrix
Uhas an elementu 22 = 1 − 1020. Numerically, since we do have a limited precision, which for
double precision is approximatelyεM∼ 10 −^16 it means that this number is approximated in
the machine asu 22 ∼− 1020 resulting in a machine representation of the matrix as