Generally, let L_k be a lower triangular matrix which is ones on the diagonal, zeros everywhere else, but possibly nonzero in the k_th column.

Prove that the inverse of L_k is a matrix with identical structure except the entries below the diagonal have their signs switched. (Do this proof with matrix notation, keeping track of indices.)

Let L be a lower triangular matrix which is ones, zeros everywhere else, but the k+1st, k+2nd, ... nth columns can have non-zero entries.

Prove that L_k * L is a lower triangular matrix which is ones, zeros everywhere else, but the kth, k+1st, k+2nd, ... nth columns have non-zero entries which are exactly those non-zero entries directly inherited from L and L_k.

Write a *.m file which given a number n returns a vector "icount". The program will first find 100 random n x n matrices, 100 random vectors b, and will find the resulting 100 solutions. You then make icount(1) = flops. Reset the flop counter. Now find a single n x n matrix. Perform its LU decomposition. Find 100 random vectors b, and find the solution x by first solving L y = b using low_solve.m and then solving U x = y using up_solve.m. Then make icount(2) = flops.

I've provided the routine up_solve.m, you have to write the routine low_solve.m (Let me know if you fail in copying and cannibalizing the routine up_solve.m to do this.)

You now have a *.m file which given n will return icount. Using this routine, create a long vector as follows:

>> [icount] = expt(2);

>> y(1,1) = icount(1); y(2,1) = icount(2);

>> [icount] = expt(3);

>> y(1,2) = icount(1); y(2,2) = icount(2);

>> [icount] = expt(4);

>> y(1,3) = icount(1); y(2,3) = icount(2);

>> [icount] = expt(5);

>> y(1,4) = icount(1); y(2,4) = icount(2);

Repeat this until you've performed the routine for the 10 x 10 matrices. (Hint: you might want to do this with another *.m file, if you're quicker with an editor than with typing.)

Let x = 2:10. Plot y(1,:) against x. >> hold on. Plot y(2,:) against x using a different line-type.

Explain your results.

1) Verify the above statements. (This is just a formal manipulation.)

2) How do we know that the new matrices L'_1 and L'_2 are actually lower triangular? If you can't prove it, at least go through a 3 x 3 example by hand and verify that they really are.

3) Matlab will give us these matrices: [L,U,P] = LU(A) where PA = LU . To solve A x = b, this means it's good enough if we solve PA x = Pb = b'. So first multiply b by P to find a new right hand side. We are now on familiar territory. PA x = b' ==> LU x = b'. First solve L y = b', then solve U x = y.

4) Now, take the matrix A from exercise 6 on page 242. Take b = [0.3097; 0.2838; 0.4742; 0.4716]. First use my LU solver:

>> [L,U] = LU(A);

>> [y] = low_solve(L,b);

>> [x] = up_solve(U,y);

>> max(A*x-b)

Then use the matlab solver with pivoting:

>> [L,U,P] = lu(A);

>> bp = P*b;

>> [y] = low_solve(L,bp);

>> [x] = up_solve(U,y);

>> max(A*x-b)

What did you find? What did you learn? Could you have used the matlab version of the lu docomposition routine? (What happens if you'd typed "[L,U] = lu(A)"?)