Question

Question1 (50 pts): LU Factorization Code for Square Matrices without Row Exchange

I want you to write an LU decomposition program in Matlab for square matrices (n×n) where row exchange is not necessary (that is no pivot is 0). Here are some hints and requirements for your matlab code.

1. You should write comments for every procedure. Make sure that your code is well indexed (see below). Otherwise it will be hard for me to follow and bad programming practice for you.

2. This should be a Matlab Function called slu (don’t name it lu because matlab already has a program called lu):

3. Your program should be a function in Matlab. So your filename should be slu.m (.m is standard for matlab scripts or functions). In order to be a function you should start with

   function [L,U]=slu=slu(A)

   where A is the input square matrix. Always write what a program does using matlab

comment (%) so that it

4. First step is to determine the size of the matrix. Use size command to obtain the size of a matrix.

   [n,m]=size(A)

5. Then if n is not equal to m end the program. If n=m continue

6. So your code should look like this

   function [L,U]=slu=slu(A)
   %LU factorization of square matrices %with no row exchange
   [n,n]=size(A) %determine the size of A Tol=1.e-6 %tolerance level for zero

   for k=1:n %loop over pivots of A if(A(k,k) < tol) %check for pivot

   disp(‘cannot proceed without row exchange’) end %cannot proceed without row exchange L(k,k)=1;

for i=k+1:n %loop over each line to determine lik determine lik

for j=k:n start from row k+1 and col. k remove (pivot row×Lik) from each row of A

end

for j=k:n
write the pivot line to Ubecause pivot is fixed

end end

due 23/10/2019

Question2(50 pts): Solving matrix equation using LU

Next step is to solve the equation using substitution. So you will need to do the standard procedure for LU. First step is to solve Ly=b. And next is to solve Ux=y. For Ly=b you need to go from top to bottom, but for Ux=y you need to go from bottom xn to x1. This function will use slu function so you need to be on the same directory. Name your file “luslv.m”.

Function x=luslv(A,b)
%Solve Ax=b using L&U from slu(A) %No row exchanges

[L,U]=slu(A) %first decompose A to LU

First do a for loop from k=1:n and another for loop for j=1:k-1 and sum the contributions to b(k) from yk-1, yk-2...y1.

Then y(k)=b(k)-s %forward elim. To solve Ly=b

Then solve Ux=y. You need to start from n and go to 1. In matlab you can do it

for k=n:-1:1 %go backwards from n to 1
Sum up contributions from xk+1, xk+2...xn.

Remove the sum from y(k) and divide by the pivot

Answer 1

%%Matlab code for LU factorization and solving
clear all
close all
%Coefficient matrix A
A=[9 -2 6 1;-10.5 -10 100 -4;-9 -2 3 3;-3 2 250 -2];
b=[18;-2;2.5;2];

%Displaying all matrix
fprintf('The A matrix is \n')
disp(vpa(A,2))

fprintf('The b matrix is \n')
disp(vpa(b,2))

%Calling function lup
[L,U,P]=slu(A);
%loop fo finding all x
X=bsub(U,fsub(L,P*b(:,1)));

%Displaying solution matrix
fprintf('The solution matrix X is \n')
disp(vpa(X,12))

%error in solution
X2=A\b;

fprintf('Error in solution of x is %e \n',norm(X-X2))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%Matlab function for LU factorization using partial pivoting
function [L, U, P] = slu(A)
%A is the coefficient matrix
%L and U are lu factorized value of A
%P is parmutation matrix
%below is the algorithm for Gaussian elimination
    [m, n] = size(A); L=eye(n); P=eye(n); U=A;
    %initialization of L U and P
    %Loop for LU matrix formation for A matrix
    for k=1:m-1
        pivot=max(abs(U(k:m,k)));
        for j=k:m
            if(abs(U(j,k))==pivot)
                ind=j;
                break;
            end
        end
        %updated L U and P values
        U([k,ind],k:m)=U([ind,k],k:m);
        L([k,ind],1:k-1)=L([ind,k],1:k-1);
        P([k,ind],:)=P([ind,k],:);
        %Loop for L and U matrix formation
        for j=k+1:m
            L(j,k)=U(j,k)/U(k,k);
            U(j,k:m)=U(j,k:m)-L(j,k)*U(k,k:m);
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%Matlab function for backward subtitution
function x=bsub(U,y)
        %Code for Backward subtitution
        x=zeros(length(y),1);
        for i=length(y):-1:1
        s=0;
        for j=length(y):-1:1
            if i~=j
                s=s+U(i,j)*x(j);
            end
          
        end
        x(i)=(y(i)-s)/U(i,i);
        end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%Matlab function for forward subtitution
function y=fsub(L,b)
        %Code for forward subtitution
        y=zeros(length(b),1);
        for i=1:length(b)
        s=0;
        for j=1:length(b)
            if i~=j
                s=s+L(i,j)*y(j);
            end  
        end
        y(i)=(b(i)-s)/L(i,i);
        end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%