Computational Physics - Department of Physics

6.4 Linear Systems 177

whereLTis the upper matrix, implying that

LTi j=Lji.

The algorithm for the Cholesky decomposition is a special case of the general LU-decomposition
algorithm. The algorithm of this decomposition is as follows

• Calculate the diagonal elementLiiby setting up a loop fori= 0 toi=n− 1 (C++ indexing
  of matrices and vectors)

Lii=

(

Aii−

i− 1

∑

k= 0

L^2 ik

) 1 / 2

• within the loop overi, introduce a new loop which goes fromj=i+ 1 ton− 1 and calculate

Lji=

1

Lii

(

Ai j−

i− 1

∑

k= 0

Likljk

)

For the Cholesky algorithm we have always thatLii> 0 and the problem with exceedingly
large matrix elements does not appear and hence there is no need for pivoting.
To decide whether a matrix is positive definite or not needs some careful analysis. To find
criteria for positive definiteness, one needs two statements from matrix theory, see Golub
and Van Loan [28] for examples. First, the leading principalsubmatrices of a positive definite
matrix are positive definite and non-singular and secondly amatrix is positive definite if and
only if it has anLDLTfactorization with positive diagonal elements only in the diagonal matrix
D. A positive definite matrix has to be symmetric and have only positive eigenvalues.
The easiest way therefore to test whether a matrix is positive definite or not is to solve the
eigenvalue problemAx=λxand check that all eigenvalues are positive.

6.4.3 Solution of Linear Systems of Equations

With the LU decomposition it is rather simple to solve a system of linear equations

a 11 x 1 +a 12 x 2 +a 13 x 3 +a 14 x 4 =w 1

a 21 x 1 +a 22 x 2 +a 23 x 3 +a 24 x 4 =w 2

a 31 x 1 +a 32 x 2 +a 33 x 3 +a 34 x 4 =w 3

a 41 x 1 +a 42 x 2 +a 43 x 3 +a 44 x 4 =w 4.

This can be written in matrix form as
Ax=w.

whereAandware known and we have to solve forx. Using the LU dcomposition we write

Ax≡LUx=w. (6.21)

This equation can be calculated in two steps

Ly=w; Ux=y. (6.22)

To show that this is correct we use to the LU decomposition to rewrite our system of linear
equations as
LUx=w,