Question

In this exercise you will work with LU factorization of an matrix A. Theory: Any matrix A can be reduced to an echelon form by using only row replacement and row interchanging operations. Row interchanging is almost always necessary for a computer realization because it reduces the round off errors in calculations - this strategy in computer calculation is called partial pivoting, which refers to selecting for a pivot the largest by absolute value entry in a column. The MATLAB command [L, U] = lu(A) returns a permuted lower-triangular matrix L and an upper-triangular matrix U, such that, . The matrix U is an echelon form of a square matrix A or resembles an echelon form for a general matrix A; the rows of matrix L can be rearranged (permuted) to a lower-triangular matrix with 1’s on the main diagonal. To learn more about the LU factorization, please read section 2.5 of the textbook and the corresponding section in the “Study Guide” which you will find in MyMathLab under the “Tools for Success” link. Note: The MATLAB algorithm for LU factorization is somewhat different from the one presented in the textbook, Section 2.5. The computer algorithm is better adjusted to the computational purposes. Practical Applications of the LU Factorization for Square Matrices: Part I: When A is invertible, MATLAB uses LU factorization of , to find the inverse matrix by, first, inverting L and U and, then, computing . **Create a function in MATLAB that will begin with the commands: function [] = eluinv(A) [~,n]=size(A); [L,U] = lu(A) (Output matrices L and U) **Verify that A is equal to L*U. If it is the case, display the message: disp('Yes, I have got LU factorization') Note: You will need to use the function closetozeroroundoff for this check and for all other checks that follow. **Verify that U is an echelon form of A. If it is the case, output a message: disp('U is an echelon form of A') If it is not the case, the output message should be something like: disp('Something is wrong') **Next, check whether A is invertible or not. I suggest using the rank command. If A is not invertible, output: sprintf('A is not invertible') invA=[]; and terminate the program. **If A is invertible, calculate the inverses of the matrices L and U, invL and invU respectively, by applying the row-reduction algorithm to [L eye(n)] and [U eye(n)]. Due to the structure of the matrices L and U, these computations will not take too many arithmetic operations. **Calculate the matrix invA using the matrices invL and invU (see the theory above) and output the matrix invA. **Run the MATLAB command inv to compute the matrix P=inv(A), output P, and compare the matrices invA and P=inv(A). (Of course, you should use the function closetozeroroundoff of their difference.) If invA and P=inv(A) match, display the message: disp('Yes, LU factorization works for calculating the inverses') Otherwise, the message should be: disp('LU factorization does not work for me?!') **Type the function eluinv in your diary file. **Run the function eluinv(A) on the following matrices: (a) A=[1 1 4;0 -4 0;-5 -1 -8] (b) A=magic(4) (c) A=[2 1 -3 1;0 5 -3 5;-4 3 3 3;-2 5 1 3] (d) A=magic(3) Part II LU factorization can be successfully used in solving simultaneously a finite sequence of matrix equations with the same matrix A. This algorithm significantly reduces the number of arithmetic operations, or “flops” (floating point operations), provided L and U are calculated, in comparison with the other methods. Suppose that we have p matrix equations: Let . Solving simultaneously the equations above, we are looking for the solution in the form of a matrix (1) Assume that we have LU factorization of an invertible matrix A: . The logical chain below demonstrates the procedure of using LU factorization to find the solution X: (2) **Create a function in MATLAB: function [X1,X2] = msystem(A,B) [~,n]=size(A); [~,p]=size(B); [L,U] = lu(A); The inputs will be an invertible matrix A and matrix B. Line 4 of the code above will produce the matrices of LU factorization of A (do not output them here). First, we will calculate the solution according to the algorithm in (2). The output will be the matrix X2. Then, we will compute X1 according to the algorithm in (1). Calculating X2: **Find the solution (which is a matrix!) of the system using the row-reduction algorithm, which means that you will employ the command “rref” in your code. **Then, find the solution X2 of the equation (Y here is the output from the previous step) also using the row-reduction algorithm. The output matrix X2 must have p columns. (Remember, we are solving simultaneously p matrix equations!) Calculating X1: **Find the solution X1=inv(A)*B by using the MATLAB function inv. **Next, compare X1 and X2 by calculating closetozeroroundoff(X1-X2). If the last is a zero matrix, output the message: disp('The solutions are the same') else disp('There is a problem with my code!') **Type the function msystem in your diary file. **Run the function [X1,X2] = msystem(A,B) on the matrices: (a) A=[1 1 4;0 -4 0;-5 -1 -8], B=randi(10,3,4) (b) A=magic(3), B=[magic(3),eye(3)] (c) A=magic(5), B=randi(10,5,3)

Answer 1

SOLUTION

1.>Keep the following .m files in the same folder

2.>To run the commands save the .m files , and in command line type for example: - quer(ones(5,3))

quer.m

function [ ] = quer( A )

format, format compact
[m,n] = size(A);

[Q,R]=qr(A);
if (closetozeroroundoff( A - Q*R == 0))
disp(' the product of Q and R forms a QR decomposition of A');
else
disp(' No, it cannot be true ');
end

s=0;
t=0;
if( inv(Q) == transpose(Q))
s=1;
disp('Q is a unitary matrix');
end

if (istriu(R))
t=1;
disp('R is an upper-triangular matrix');
end

if( s==1 && t==1)
disp('Q*R forms an orthogonal-triangular decomposition of A');
end

if (m == n)
Anew=A;
count1=1;
P=closetozeroroundoff(Anew-triu(A));
while(P~=zeros(m,n))
Anew=R*Q;
[Q,R]=qr(Anew);
P=closetozeroroundoff(Anew-triu(A));
count1=count1+1;
end
disp(Anew); %matrix B
disp(count1); %no of iterations

for i=1:m
disp(Anew(i,i)); %diagonal matrix
end

e = eig(A);
disp(e) %eigen values

end

closetozeroroundoff.m

function [ B ] = closetozeroroundoff( A )
%UNTITLED3 Summary of this function goes here
% Detailed explanation goes here
[m,n] = size(A);
for i=1:m
for j=1:n
if abs(A(i,j))<10^(-7)
A(i,j)=0;
end
end
end
B=A;

end