LU Decomposition

Matrices are used widely in Numerical Computing, many alogirithms inolve decomposing a matrix, A into an equivalent pair L and U, where L and U are lower and upper triangular matrices respectively and .

Consider the following system of equations,

here

so

and

Both L and U can be stored in the same matrix; since the diagonal of L is always unity it does not have to be store, the combined LU matrix can be stored as

Now, given that

then we solve for where

(recall that the diagonal of L will be all 1s) and then solve

The advantage of this method is that the triangular system of equations is very easy to slove. can be sloved by forward substitution and by backward substitution. is calculated by forward substitution as follows :

Once is known can easily be obtained by backwards substitution:

The LU decomposition algorithm is as follows:

For each row k from 1 to n-1 in order,

• column normlisation -- The elements in the colum below the diagonal element A(k,k) are divided by the diagonal element.
• submatrix modification -- The value of each element A(i,j) in the submatrix of A below and to the right of the diagonal element A(k,k), ie, A(k+1:n,k+1:n) is modified by subtracting A(i,k)*A(k,j)

The algorithm will not work if the A matrix is singular, that is, if any of the rows of A are all zero.

This algorithm omits the issue of pivoting so can be unstable - your program will work for the 3 by 3 matrix givenm here.

The an HPF program with subroutines to perform LU decomposition and foward and backwards substitution.

• Solution