Homework #2

Due Friday, July 18, 2003

Walking The Straight Line

Answer all questions. Submit any pertinent hand calculations and/or Matlab code in a neat manner in order to help us give you partial credit. Readings: Van Loan 5.2.5, 6.

1) Consider the following matrix:

A = ( \begin{matrix}1 & 3 & 2 \\ 2 & 15 & 5 \\ −2 & 39 & 3 \end{matrix} ) .

(a) Calculate the LU decomposition of A without any pivoting. That is, calculate

M_{2}M_{1}A = U, L = M_{1}^{−1}M_{2}^{−1},

where M_{1} and M_{2} are elementary elimination matrices as given in section 6.3.2 of the text, and L and U are upper and lower triangular matrices, respectively. Show how you arrive at your M_{1}, M_{2}, U and L. (b) Calculate the LU decomposition with partial pivoting. That is, calculate

M_{2}P_{2}M_{1}P_{1}A = U

or, equivalently,

P A = LU.

Write out M_{1}, M_{2}, P_{1},P_{2}, P, L and U explicitly. Input A into Matlab and ask it to compute the LU decomposition using the built-in Matlab function lu: [L, U] = lu(A);. Compare your hand-caculations of L and U to these, if the results are different explain the differences. Now compute L and U again, this time using the following syntax: [L, U, P] = lu(A);. Does this result agree with your hand-calculations? Explain any differences. (c) The theory says that if you consider Ax = b with the equations reordered as given by P Ax = P b, where P = P_{2}P_{1} and P_{1} and P_{2}are from part (b), one can perform Guassian elimination without pivoting and still have all the elimination matrices with entries of magnitude less than or equal to 1. Show that this is true in this example.

2) Assume that A and C are given n − by − n matrices and that A is nonsingular. Assume that g and h are given n − by − 1 vectors. Consider the system

\begin{aligned}A^{T} y + Cz & = & g \\ Az & = & h \end{aligned} .

Show how to solve the system for z and y using two forward substitutions, two back substitutions and one LU decomposition. You may assume that pivoting is not needed anywhere.

3) Consider the n − by − n Hilbert matrix H with entries H_{ij} = 1/(i + j − 1), i, j = 1, 2, …, n. Note that Matlab has a function that constructs such matrices. Just type H = hilb(n). Suppose that we wish to solve Hx = b, where b = e_{1}, the first unit vector in R^{n}. Write a Matlab program that performs the following study:

First note that H is symmetric positive definite, which means that a Cholesky factorization exists. (a) Compute the Cholesky factorization H = U^{T} U, where U is upper triangular. In Matlab, just type U = chol(H). (b) Solve for x by forward and back substitution. Let’s call this solution x̃. Do this for n = 2, 3, …, 14. (c) For the purpose of camparison, calculate the true solution x = H^{−1}b. Note that Matlab calculates exactly the inverse of H, H^{−1}. To get it, type invhilb(n). (d) Plot the ∞-norm of the residual, ∥b − Ax̃∥_{∞}, versus matrix size, n, for n = 2, 3, …, 14. (e) Plot the relative error, \begin{matrix}∥x−x̃∥ & _{∞} \\ ∥x∥ & _{∞} \end{matrix}, versus matrix size, n, for n = 2, 3, …, 14. How large is n when the relative error is bigger than 1? (f) Plot the condition number of H versus matrix size, n, for n = 2, 3, …, 14. To get the condition number of a matrix A in Matlab type cond(A). Comment on the trend showed in this plot. (g) Note how big the condition number of H is when the relative error is bigger than 1.

Revision: Matlab may have trouble computing the Cholesky factorization of our Hilbert matrices when n gets bigger than a certain value and this is machine dependent; my machine went up to n = 14, while at least one student in our class could go up to only n = 12. To circumvent this problem, go up to the n when Matlab fails for you. In part (e), instead of making comparisons of the relative error against 1 compare against the largest observable relative error and denote this number by RE. Now in part (g) note how big the condition number of H is when the relative error is bigger than RE.