Download Numerical Computing Problem Set 2 Solutions: Forward Substitution, LU Factorization - Prof and more Assignments Mathematics in PDF only on Docsity!
D. Schwendeman
Numerical Computing Problem Set 2
Due: Thursday, 9/22/
- Text exercise 2.4, page 97.
- Text exercise 2.7, page 97.
- Consider the system Ax = b, where A is an n Ć n lower triangular matrix and b is an n-vector. (a) Describe an algorithm to solve this system for x using forward substitution. What is the operation count of your algorithm (i.e., the total number of floating point additions, multiplications, and divisions required)? How does this count behave as n becomes large? (b) Write a Matlab program that computes the solution x for a given vector b and lower tri- angular matrix A. Run your program for a 5Ć5 system generated using A=tril(rand(5,5)), x=rand(5,1) and b=A*x. Compare the solution computed using your program with that generated using Matlabās magic backslash.
- Let
A =
ļ£® ļ£Æļ£°^52 14
3 2 7
ļ£¹ ļ£ŗļ£»
(a) Find elimination matrices M 1 and M 2 such that M 2 M 1 A = U , where U is an upper triangular matrix. (b) Find the matrix L in the LU-factorization of A.
- Write a Matlab program that computes L and U in the LU factorization of an n Ć n matrix A (without pivoting). Your code should check for zero pivots and give an error message if one is encountered. Run your code on the matrix A=rand(7,7), a 7 Ć 7 matrix with random elements. Check the results of your code by comparing the product LU with the original matrix A. Run the built-in program [l,u]=lu(A). Typically, Matlabās results and yours (even if correct) will not agree (even if exact arithmetic is used). Give an explanation why you think this might be so.
Problem 1: Exercise 2.4, Page 97
Part (a) A matrix is singular if and only if its determinant is zero (there are, of course, other equivalent statements; see page 51 of the text for a discussion). To verify that det(A) = 0, we use the following Matlab code, which returns a result of 0.
1 det([1,1,0;1,2,1;1,3,2])
Part (b) We notice immediately that the vector b is just twice the second column of A; hence the vector x = [0, 2 , 0]T^ is a solution to the linear system Ax = b. Since A is singular, any vector of the form x + Ī³z is also a solution to the linear system Ax = b, where A(z) = 0 and Ī³ ā R (see page 51 for a discussion). Hence there are infinitely many solutions.
Problem 2: Exercise 2.7, Page 97
Part (a) Using the familiar form of a 2 Ć 2 determinant, we have det(A) = 1 ā (1 + Īµ)(1 ā Īµ) = 1 ā (1 ā Īµ^2 ) = Īµ^2.
Part (b) If we compute the determinant using the formula det(A) = Īµ^2 , then we compute this number as zero if it is less than the underflow value: |Īµ^2 | <UFL, or |Īµ| < āUFL. If we tell a program like Matlab to compute det(A) for us, it will need to compute the values 1 Ā± Īµ; in this case, if |Īµ| < Īµmach, then f l(1 Ā± Īµ) = 1, so det(A) will be computed as 0 (assuming the computer follows the normal order of operations).
Part (c)
A = LU =
1 ā Īµ 1
ļ£Æļ£° 1 1 +^ Īµ 0 Īµ^2
Part (d) Again, there are two options: we can either choose to use the simplification det(U ) = Īµ^2 , or we can allow a program like Matlab to compute the value from the matrix U. The problematic values of Īµ for these cases are the same as in part (b).
1 A = tril(rand(5,5)); b = rand(5,1); abs(fwdSubst(A,b) - A \ b) If we get small values (I happened to get a vector of 0ās when I tested it), weāve probably written our code correctly.
Problem 4
For this problem, use as a reference sections 2.4.3 and 2.4.4 from the text.
Part (a)
M 1 =
, M 2 =
To obtain these matrices, we look at the top of page 67; think of M 1 as acting on the vector a that is the first column of the matrix A. Similarly, think of M 2 as acting on the second column of A. In general, we obtain Mk by thinking of how it acts on the kth column of A; we use the process described in the book to obtain the kth column of Mk in this way, and we fill in the rest of Mk as though it were the appropriate identity matrix.
Part (b)
L =
, U =
If M 2 M 1 A = U = M 2 M 1 (LU ), then we must have M 2 M 1 L = I, as explained in the book. So we obtain L from the equation L = M 1 ā 1 M 2 ā 1. (You can verify that LU = A.)
Problem 5
Algorithm 2.3 from the textbook computes the āin-placeā LU factorization of A (provided you follow the textās comment above the algorithm), and is implemented below.
1 function [L,U] = lufact(A)
2 n = size(A,1); 3 for k=1:n- 4 if A(k,k)== 5 disp('pivot is zero'); 6 return 7 else 8 A(k+1:n,k)=A(k+1:n,k)/A(k,k); 9 A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); 10 end 11 end 12 L = tril(A,-1)+eye(n); 13 U = triu(A); To test this code, enter the following into Matlabās command line.
1 A = rand(7,7); [L1,U1] = lu(A); [L2,U2] = lufact(A); 2 norm(L1 - L2) 3 norm(U1 - U2)
My resultant norms were on the order of 10. There isnāt an exact match because Matlabās built-in code is highly optimized, and includes pivoting (note that the action of pivoting is different from the elements of the matrix identified as being the pivots; see page 70 of the textbook for an explanation of pivoting).
