EXERCISE 4 Enter the matrix A and the vector b in MATLAB 6 A

EXERCISE 4 Enter the matrix A and the vector b in MATLAB 6 A=1-4-3-7-8 3 The exact solution to the system Axb is the vector x(-1,8,1, -2)7. (a) Enter [L,U,P] -lu(A) to find the LU decomposition of the matrix PA and verify that PA LU (b) Use the LU decomposition you found in part (a) to solve the system Ax -b. Call the computed solution x lu. (c) Enter the vector x and compare your solution x lu from part (b) with the exact solution x by computing norm(x lu - x) (the norm function gives the magnitude of the vector, that is, for a vector aa a2,...,an)7 the norm of a is detined as: norm (a) = vas + + . . . + an ) EXERCISE 5 In this question we will compare the speed of two methods for solving the equation Ax -b when A is an invertible square matrix. We will use the MATLAB tic and toc command to measure the computation times. Enter: A. rand(500) ; ones (500,1); b#A*x; x - Important: Be sure to use semicolon after each command so that matrices and vectors are not dis- played. Do not print or include these large matrices and vectors in your lab write-up. (a) Solve Ax-b using the reduced row echelon form and store the solution in x rref tic; R-rref([A, b)); x-rref -R(: ,end); toc (b) Solve Ax - b using the LU decomposition as you did in EXERCISE 4(a)(b), and calculate the elapsed time using the tic toc function. Store the solution in x_lu. (Make sure you use semicolon; the only output should be the elapsed ti. Which method is faster? (c) Compare the solutions from parts (a) and (b) with the exact solution x by computing norm(x rref x) and norm(x.lu - x) . How accurate are the solutions from parts (a) and (b)? 2011 Stefania Tracogna. SoMSS. ASU

Solution

function [L, U, P] = lu_decomposition_pivot(A) n = size(A,1); Ak = A; L = zeros(n); U = zeros(n); P = eye(n); for k = 1:n-1 for i = k+1:n [~,r] = max(abs(Ak(:,k))); Ak([k r],:) = Ak([r k],:); P([k r],:) = P([r k],:); L(i,k) = Ak(i,k) / Ak(k,k); for j = k+1:n U(k,j-1) = Ak(k,j-1); Ak(i,j) = Ak(i,j) - L(i,k)*Ak(k,j); end end end L(1:n+1:end) = 1; U(:,end) = Ak(:,end); return