Gaussian Elimination Method With Partial Pivoting in Matlab

1. Using the Gaussian Elimination Method in Matlab
2. Using Pivoting as an Assisting Function in Matlab

The article will help the reader understand how to use Gaussian Elimination Method in Matlab.

Using the Gaussian Elimination Method in Matlab

The Gauss method is a classical method for solving linear algebraic equations (SLA) systems.

This is a method of sequential exclusion of variables; when using elementary transformations, the system of equations is reduced to an equivalent system of triangular form.

All the system variables are sequentially located, starting from the last (by number).

The Gauss-Jordan algorithm divides the first row by coefficient A[1,1] at the first step.

Then the algorithm adds the first row to the rest of the rows with such coefficients that their coefficients in the first column turn to zeros - for this, when adding the first row to the i-th, it is necessary to multiply it by -A[i,1].

For each operation with matrix A (division by a number and addition to one row by another), the corresponding operations are performed with vector b; in a sense, it behaves as if it were the m+1st column of matrix A.

As a result, at the end of the first step, the first column of matrix A will become single (i.e., it will contain one in the first row and zeros in the rest).

Similarly, the second step of the algorithm is performed; only now the second column and the second row are considered. First, the second row is divided by A[2,2], then subtracted from all other rows with coefficients such as resetting the second column of the matrix A.

And so on until we process all rows or columns of matrix A. If n= m, then it is obvious from the construction of the algorithm that the matrix A will turn out to be singular, which is what we needed.

Using Pivoting as an Assisting Function in Matlab

Of course, the scheme described above is incomplete. It works only if at each i-th step the element A[i, i] is different from zero - otherwise, we will not be able to zero the remaining coefficients in the current column by adding the i-th row to them.

There is a process of selecting a reference element (in one word, “pivoting”) to make the algorithm work in such cases. It consists of the fact that the matrix’s rows and columns are rearranged so that the desired element A[i, i] has a non-zero number.

Note that the rearrangement of rows is much easier to implement on a computer than the rearrangement of columns: after all, when swapping some two columns, you need to remember that these two variables swapped places, so that when restoring the answer, correctly restore which answer belongs to which variable.

When rearranging the rows, no such additional actions need to be performed.

In other words, before performing the i-th phase of the Gauss-Jordan algorithm with the partial pivoting heuristic, it is necessary to find the maximum modulo in the i-th column among the elements with indices from i to n, and exchange the row with this element with the i-th row.

function x = pivot(~, ~)
% Program to Solve Ax = b using Gaussian Elimination Method with partial-
% pivoting (Row-exchange)
%==========================================================================

% OUTPUT:
%==========================================================================

% Solution vector, x

%==========================================================================
A=[1 2 3; 4 5 6];
b=[1 -2 5];
[m,n] = size(A); % checking the size of matrix


%==========================================================================
% Initialization
x = zeros(m,1);
l = zeros(m,m-1);

%==========================================================================
% Main Program
%==========================================================================

% Reducing Matrix A to upper triangular form

for k = 1:m-1
    % =========Performing Partial-pivoting=================================
        for p = k+1:m
            if (abs(A(k,k)) < abs(A(p,k)))
                A([k p],:) = A([p k],:);
                  b([k p]) = b([p k]);
            end
        end
    % =====================================================================
% =========Setting Matrix A to triangular form=================================
    for i = k+1:m
        l(i,k) = A(i,k)/A(k,k);
        for j = k+1:n
            A(i,j) = A(i,j)-l(i,k)*A(k,j);
        end
        b(i) = b(i)-l(i,k)*b(k);
    end
end

for k = 1:m-1
    for i = k+1:m
        A(i,k) = 0;
    end
end

%==========================================================================
% Backward substitution of Matrix A
%==========================================================================
x(m) = b(m)/A(m,m);

for i = m-1:-1:1
    sum = 0;
    for j = i+1:m
        sum = sum + A(i,j)*x(j);
    end
    x(i) = (b(i)- sum)/A(i,i);
end

%==========================================================================
end

Output:

ans =

    -3
     2