Gauss Elimination

Here are the M-files to implement naive Gauss Elimination and Gauss Elimination with partial pivoting.

Naive Gauss Elimination M-file

function x = GaussNaive(A,b)

% GaussNaive: naive Gauss elimination

%   x = GaussNaive(A,b): Gauss elimination without pivoting.

% input:

%   A = coefficient matrix

%   b = right hand side vector

% output:

%   x = solution vector

[m,n] = size(A);

% m = no of rows, n = no of columns, size of matrix A is m x n

if m~=n, error('Matrix A must be square'); end

% if m not equal to n, error message is displayed as 'Matrix A must be square'

nb = n+1;

% no of columns for augmented matrix = no of columns of the matrix A + 1 column for the vector B

Aug = [A b];

% Coefficient matrix A and right-hand-side vector b are combined in an augmented matrix

% forward elimination

for k = 1:n-1

    % k represents column 1 to (n-1)th column

  for i = k+1:n

      % i represents row 2 to nth row

    factor = Aug(i,k)/Aug(k,k);

    % to compute the factor in order to obtain the upper triangular matrix

    Aug(i,k:nb) = Aug(i,k:nb)-factor*Aug(k,k:nb);

    % to compute the new equations until transform the matrix into upper triangular form

    % the new ith row kth column elements is equal to the old ith row kth column elements 

    % minus the multiplication of factor and kth row kth column elements 

    % since using augmented matrix, need to use nb (n+1) because of additional 1 column vector b

  end

end

% back substitution

x = zeros(n,1);

% solution vector x is a matrix of zero coefficients of n x 1 dimension

% n =  no of equations

x(n) = Aug(n,nb)/Aug(n,n);

% to calculate the last unknown x in nth equation

% by dividing the constant b in last equation with the coefficient of x(n)

% formula is same as Eq (9.12) in textbook

for i = n-1:-1:1

    % i represents no of equations (n) minus 1 with increment of -1 until 1

  x(i) = (Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);

  % to calculate the remaining unknown x's in (n-1)th equation and so on until the first equation

  % formula is same as Eq (9.13) in textbook

end

Gauss Elimination with Partial Pivoting M-file

function x = GaussPivot(A,b)

% GaussPivot: Gauss elimination pivoting

%   x = GaussPivot(A,b): Gauss elimination with pivoting.

% input:

%   A = coefficient matrix

%   b = right hand side vector

% output:

%   x = solution vector

[m,n]=size(A);

% m = no of rows, n = no of columns,size of matrix A is m x n 

if m~=n, error('Matrix A must be square'); end

% if m not equal to n, error message will be displayed as 'Matrix A must be square'

nb=n+1;

% no of columns of augmented matrix = no of columns of matrix A + 1 column of vector b 

Aug=[A b];

% forward elimination

for k = 1:n-1

    % k represents column 1 to (n-1)th column

  % partial pivoting

  [big,i]=max(abs(Aug(k:n,k)));

  % to determine the coefficient with the largest absolute value in the column below the pivot element

  % big = largest element

  % i = index corresponding to the largest element (no of row)

  ipr=i+k-1;

  % new row (with largest element below the pivot element)= old row (with pivot element)+ k - 1

  if ipr~=k

      % if new row doesn't equal to k (row with largest element isn't row 1)

    Aug([k,ipr],:)=Aug([ipr,k],:);

    % to switch the row with largest pivot element to the first row

  end

  for i = k+1:n

      % i represents row 2 to nth row

    factor=Aug(i,k)/Aug(k,k);

    % to compute the factor in order to obtain the upper triangular matrix

    Aug(i,k:nb)=Aug(i,k:nb)-factor*Aug(k,k:nb);

    % to compute the new equations until transform the matrix into upper triangular form

    % the new ith row kth column elements is equal to the old ith row kth column elements 

    % minus the multiplication of factor and kth row kth column elements 

    % since using augmented matrix, need to use nb (n+1) because of additional 1 column vector b

  end

end

% back substitution

x=zeros(n,1);

% solution vector x is a matrix of zero coefficients of n x 1 dimension

% n =  no of equations

x(n)=Aug(n,nb)/Aug(n,n);

% to calculate the last unknown x in nth equation

% by dividing the constant b in last equation with the coefficient of x(n)

% formula is same as Eq (9.12) in textbook

for i = n-1:-1:1

    % i represents no of equations (n) minus 1 with increment of -1 until 1

  x(i)=(Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);

  % to calculate the remaining unknown x's in (n-1)th equation and so on until the first equation

  % formula is same as Eq (9.13) in textbook

end