The LU DecompositionWorksheet by Russell Blyth and Mike May, S.J.restart: with(LinearAlgebra): with(plots): with(plottools):

This worksheet investigates the LU decomposition of a matrix.The basic objectives are:1) Compute the LU decomposition, if possible.2) See how elementary matrices are involved in the process.3) Understand how the LU decomposition is useful in solving systems of equations.Note: If you need to re-execute any command in this worksheet you should re-execute the entire worksheet from the first command line above, since variable names are carried through the worksheet.

First define a random 3 x 4 matrix.A := RandomMatrix(3,4,generator=rand(-5..5));Let's keep track of the elementary matrices corresponding to the row operations used, as well as the cumulative product of the elementary matrices, which we can do by starting with the matrix A augmented by the appropriate identity matrix:Iden:=IdentityMatrix(3); AAug := <A | Iden>;Now we apply elementary row operations to the augmented matrix to get a row echelon form. We intend to use only elementary row operations of the form that add multiples of rows to lower rows (in order to get an LU-decomposition). We exhibit the elementary matrices at each step.(Note that the matrix A is not truly "random" - every time the worksheet is started, it produces the same matrix A. Thus we know in advance which row operations are required in this worksheet. This is handy for teaching purposes. If you want a truly "random" matrix, you need to precede the randmatrix command with a "randomize()", which generates a seed for the random functions in Maple, based on the computer clock.)type(AAug,Matrix);L1 := -AAug[2,1]/AAug[1,1]; AAugR := RowOperation(AAug,[2,1],L1); E1 := RowOperation(Iden,[2,1],L1);L2 := -AAugR[3,1]/AAugR[1,1]; AAugR := RowOperation(AAugR,[3,1],L2); E2 := RowOperation(Iden,[3,1],L2);L3 := -AAugR[3,2]/AAugR[2,2]; AAugR := RowOperation(AAugR,[3,2],L3); E3 := RowOperation(Iden,[3,2],L3);AAugR is in row echelon form. Note we don't bother to force leading ones. The matrix U can now be extracted:U := DeleteColumn(AAugR,5..7);The matrix B referred to in class is the product of the three elementary matrices used. This matrix B can also be extracted from AAugR:E3.E2.E1; B := DeleteColumn(AAugR,1..4);The matrix L we seek is just the inverse of this B:L := MatrixInverse(B);Now note the remarkable relationship between the entries of L, the elementary matrices, and the multipliers used in the row reductions! The matrix L can in fact be written down without actually computing the inverse of B.The off diagonal entries of L are the negatives of the multipliers used in the row operations.The off diagonal entries of L are also the off diagonal entries of the inverses of the elementary matrices corresponding to the row operations.[-L1, -L2, -L3]; [E1^(-1),E2^(-1),E3^(-1)];Check that A = LU(L . U); (A);

1) Repeat the above for a random 4 x 4 matrix AA (be careful not to reuse any of the variable names used previously, since some get reused below!)

Now let's see how to use the LU decomposition to solve a system AX = Y.We showed above that this is equivalent to a system L(UX) = Y. We can solve this system in two steps using only substitution on a triangular system in each step. Let Z=UX.Generate a RHS vector Y.Y := RandomMatrix(3,1,generator=rand(-5..5));Solve LZ = Y for Z; note this amounts to a simple forward-substitution (no further row reductions). We show the augmented matrix and have Maple compute.<L | Y>; Z:=LinearSolve(L,Y);Then solve UX = Z for X; note this amounts to a simple back-substitution (no further row reductions). We show the augmented matrix and again have Maple compute the solution for us.<U|Z>; X := LinearSolve(U,Z);(Maple does not always select the parameter the way we do.)Compare the result of solving the original system.LinearSolve(A,Y);

2) Generate a random 4 x 1 vector YY, and then repeat the steps above to solve (AA)(X) = YY.3) What advantage would the LU decomposition have for computer computation of solutions to linear systems?4) Use Maple's help command to find out about the LUdecomp command in the linalg package, Use the command to check your LU decomposition of AA in Exercise 1.