Gauss Elimination

Gauss Elimination

Matrix Problem Statement

As per hypothesis, LU Factorization should be more efficient than Gauss Elimination if you are solving many linear systems of equations ([A][X] = [B]), where the coefficient matrix ([A]) remains the same but the constants matrix ([B]) changes.

This project is meant to test if this hypothesis is true and how size effects efficiency.

Solution 

GaussLUTime.m 

function [GaussTime,LUTime] = GaussLUTime(A)

if size(A,1) ~= size(A,2)

disp(‘Input should be a square matrix.’)

return

end

B_Matrix = -10 + (10+10)*rand(size(A,1),200);

%%Using Gaussian elimination

A_gauss = A;

tic;

n = size(A_gauss,1);

for index=1:200

B = B_Matrix(:,index);

%create upper triangular matrix

s=0;

for j=1:n-1

ifA_gauss(j,j)==0

k=j;

for k=k+1:n

ifA_gauss(k,j)==0

continue

end

break

end

B=A_gauss(j,:); C=B(j);

A_gauss(j,:)=A_gauss(k,:); B(j)=B(k);

A_gauss(k,:)=B; B(k)=C;

end

for i=1+s:n-1

L=A_gauss(i+1,j)/A_gauss(j,j);

A_gauss(i+1,:)=A_gauss(i+1,:)-L*A_gauss(j,:);

B(i+1)=B(i+1)-L*B(j);

end

s=s+1;

end

%—————————————————————–

%Solution of equations

x(n)=B(n)/A_gauss(n,n);

for i=n-1:-1:1

sum=0;

for j=i+1:n

sum=sum+A_gauss(i,j)*x(j);

end

x(i)=(1/A_gauss(i,i))*(B(i)-sum);

end

end

GaussTime = toc;

%%Using LU decomposition

A_LU = A;

tic;

for index=1:200

B = B_Matrix(:,index);

x = getXbyLU(A_LU,B);

end

LUTime = toc;

end 

getXbyLU.m 

function x=getXbyLU(A,B)

[L,U,P] = LUdecomposition(A);

% Matrix Forward Substitution

B=P*B;

m=length(B);

x=zeros(m,1);

x(1)=B(1)/L(1,1);

for j=2:m

x(j)=(B(j) -sum(L(j,1:j-1)’.*x(1:j-1)))/L(j,j);

end

%Matrix Backward Substitution

B=x;

m=length(B);

x=zeros(m,1);

x(m)=B(m)/U(m,m);

for j=(m-1):-1:1

x(j)=(B(j) -sum(U(j,j+1:end)’.*x(j+1:end)))/U(j,j);

end

end 

LUdecomposition.m 

function [L,U,P] = LUdecomposition(A)   %  LU factorization

% [L,U,P] = Lu(A) returns unit lower triangular matrix L, upper

% triangular matrix U, and permutation matrix P so that P*A = L*U.

m=length(A);

U=A;

L=zeros(m,m);

PV=(0:m-1)’;

for j=1:m

% Pivot Voting (Max value in this column first)

[~,ind]=max(abs(U(j:m,j)));

ind=ind+(j-1);

t=PV(j); PV(j)=PV(ind); PV(ind)=t;

t=L(j,1:j-1); L(j,1:j-1)=L(ind,1:j-1); L(ind,1:j-1)=t;

t=U(j,j:end); U(j,j:end)=U(ind,j:end); U(ind,j:end)=t;

% LU

L(j,j)=1;

for i=(1+j):size(U,1)

c= U(i,j)/U(j,j);

U(i,j:m)=U(i,j:m)-U(j,j:m)*c;

L(i,j)=c;

end

end

P=zeros(m,m);

P(PV(:)*m+(1:m)’)=1;

end