In-Class Demonstration: LU-Decomposition

By Michael K. May, edited by Russell Blyth

restart: with(LinearAlgebra): with(plots): with(plottools):

We start with a random augmented matrix. We then repeat the matrix so we can keep the computations stable.

A := RandomMatrix(3,4,generator=rand(-5..5));

A := <<4 | -1 | -4 | -1>, <-1 | -4 | -2 | 3>, <4 | 4 | -3 | -3>>;

We augment the matrix by the identity matrix to record row reductions, then do normal step by step Gaussian elimination.

A1 := <A|IdentityMatrix(3)>;

A2 := RowOperation(A1,[2,1],1/4);

A3 := RowOperation(A2,[3,1],-1);

A4 := RowOperation(A3,[3,2],20/17);

We keep track of the coefficients used, namely 1/4, -1, 20/17.

We next take the appropriate columns of the coefficient matrix of A to find U and B .

L is the inverse of B.

U1 := DeleteColumn(A4,5..7); U := DeleteColumn(A4,4..7); B := DeleteColumn(A4,1..4);

L := B^(-1);

Note the the entries of L are simply the negatives of the coefficients used. Note also that L times U1 gives us the original matrix A.

L.U1; A;

Now we look at solving the system of equations represented by A.

C := Column(A,4);

S1 := LinearSolve(A);

We should get to the same solution if we solve the pair of triangular systems LY = C followed by UX = Y.

S2 := LinearSolve(L,C);

S3 := LinearSolve(U,S2);