Question

Make the program(such as matlab) for solving the linear system

Ax = b

using the Naive-Gaussian elimination & also the Gaussian elimination with partial pivoting.

Answer 1

Main program files :

#############################################################################################

function[x]=solve_gaussElimin(A,b)

n=length(b);% number of elements in b vector

% We start LU decomposition
for k = 1:n-1 % At every stage (column)
  
% Compute the corresponding multipliers of L
% and store it in A(To save the memory)
J=k+1:n;
A(J,k) = A(J,k) / A(k,k);
  
% update submatrix by outer product
A(J,J) = A(J,J) - A(J,k) * A(k,J);
end
% Here we get the factorisation PA=LU.
L=tril(A,-1)+eye(n);
U=triu(A);

y=forward_sub(L,b);

x=backward_sub(U,y);

################################################################################

function [x]=solve_Gauss_partial_pivot(A,b)

% A = [1 1 1 1 1;
% -2 -1 0 1 2;
% 4 1 0 1 4;
% -8 -1 0 1 8;
% 16 1 0 1 16]; % The A matrix
%   
% b = [2 0 8/3 0 32/5]'; %b vector

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

% Ax=b => PAx=Pb => LUx=Pb. This can be solved in 2 steps

% Step 1 : solve y for Ly=Pb
% Step 2 : solve x for Ux=y

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

% Step 1

% permute b according to P
b1=P*b;
y=forward_sub(L,b1);

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

% Step 2

x=backward_sub(U,y);

##########################################################################################

Necessary files :

##########################################################################################

function[x]=forward_sub(L,b)
% forward substitution
n=size(L,1);
x=zeros(n,1);

x(1) = b(1)/L(1,1);
for k = 2:n
x(k) = (b(k) - L(k,1:k-1) * x(1:k-1))/ L(k,k);
end

###########################################################################################

function[x]=backward_sub(U,b)
% forward substitution
n=size(U,1);

x=zeros(n,1); % Initiating
% backward substitution
x(n) = b(n) / U(n,n);
for k = n-1:-1:1
x(k) = ( b(k) - U(k,k+1:n)*x(k+1:n) ) / U(k,k);
end

############################################################################################

function [P,L,U] =Partial_pivoting(A)

n = size(A,1);

% Order in which the Identity matrix is at the begining.
p = 1:n;
P=eye(n);
% We start LU decomposition with partial pivoting
for k = 1:n-1 % At every stage (column)
  
% find the row index such that it has the absolute maximum in that column
[val,q] = max(abs(A(k:n,k)));
q = q + k-1;
  
% swap the rows "k" and "q" and do the same permutation in "p"
A([k,q],:)=A([q,k],:);
p([k,q])=p([q,k]);
  
% Compute the corresponding multipliers of L and store it in A(To save the memory)
J=k+1:n;
A(J,k) = A(J,k) / A(k,k);
  
% update submatrix by outer product
A(J,J) = A(J,J) - A(J,k) * A(k,J);
end
% Here we get the factorisation PA=LU.
L=tril(A,-1)+eye(n);
U=triu(A);
P((1:n),:)=P(p,:);
end