University of Colorado ::: Mechanical Engineering

CU Mechanical Engineering

Numerical Methods Matlab tutorials

Writing Matlab Functions: Gauss-Jordan elimination

In this example, we will write a function that solves a system of linear equations Ax = b using Gauss-Jordan elimination with partial pivoting.

To begin, all functions within matlab must begin by declaring that m-file is a function. This is accomplished by having the first line of the m-file having the word 'function' followed by the name of the variable that will be returned. This is consequently set equal to the name of the name of the file with the quantities that are passed to the function in parenthesis. For the purpose of this exercise, the name of the function will be 'gauss_jordan_elim' and the first line of the function is as follows:

function [x,err]=gauss_jordan_elim(A,b)

As with all functions in Matlab, the first line of code must begin by defining the function. This function solves a set of linear first order equations where the leading coefficients can be represented in a matrix 'A' such that the system becomes one of the form 'Ax=b'. This function thus requires two inputs: the coefficient matrix 'A' and the solution vector 'b'.

The first segment of code runs a series of tests to ensure that the inputs meet the required format.

[n,m]=size(A); % size of matrix A
err =0;
x=zeros(n,1);
if n ~= m
    disp('error: n~=m');
    err = 1;
end

The above section of code begins by using the 'size' function which can take any matrix (scalar, vector, etc.) and returns the dimensions. In this case, it returns the rows ('n') and columns ('m') of the input matrix 'A'. Next, we define a variable 'err' to track and monitors errors which might be encountered. The following line generates a column vector of zeros. The following conditional statement tests to see if the input matrix is square. If not, the error variable is assigned a value of '1'.

The next piece of code firsts tests to make sure the vector is of the right dimension using the not equals '~=' operator and takes the transpose if necessary to make the matrix dimensions agree.

if length(b) ~= n 
    disp('error: wrong size of b');
    err = 2;
else
    if size(b,2) ~= 1
        b=b';
    end
    if size(b,2) ~= 1
        disp('error: b is a matrix');
        err = 3;
    end
end

As you can see, in each of the above sections, tests are carried out to ensure the inputs fufill all the requirements of general matrix operations. If not, a specific error is code is specified and an error message is displayed on the screen.

If there are no erros, the code proceeds.

if err == 0
    Aa=[A,b];
    for i=1:n 
        [Aa(i:n,i:n+1),err]=gauss_pivot(Aa(i:n,i:n+1));
        if err == 0
            Aa(1:n,i:n+1)=gauss_jordan_step(Aa(1:n,i:n+1),i);
        end
    end
    x=Aa(:,n+1);
end
A=0;

Above, a new matrix 'Aa' is defined by concatenating the matrix 'A' and the column vector 'b'. The next step defines a 'for' loop to systematically perform the Gauss-Jordan elimination through the use of two different functions. This code uitilizes inline functions for this task.

It is necessary to define two parts of the above block of code. First is the gauss_pivot part:

function [A1,err]=gauss_pivot(A)

[n,m]=size(A); % size of matrix A
A1=A;
err = 0; % error flag
if A1(1,1) == 0
    check = logical(1); % logical(1) - TRUE
    i = 1;
    while check
        i = i + 1;
        if i > n  
            disp('error: matrix is singular');
            err = 1;
            check = logical(0);
        else
            if A(i,1) ~= 0 & check
                check = logical(0);
                b=A1(i,:);      % next three operations exchange rows 1 and i
                A1(i,:)=A1(1,:);
                A1(1,:)=b;
            end
        end
    end
end

The Gauss Pivot function begins by defining a conditional variable to monitor the progress inside the following 'while' loop which will execute as long as the conditional statement is true. This means that as long as 'check' is equal to one then the loop will continue to execute. The code continues to increment the index variable 'i' and perform the pivot. It should again be noted that each conditional or loop statement must be terminated by an 'end' statement.

The second is the gauss_jordan_step function:

function A1=gauss_jordan_step(A,i)

[n,m]=size(A); % size of matrix A
A1=A;
s=A1(i,1);
A1(i,:) = A(i,:)/s;
k=[[1:i-1],[i+1:n]];
for j=k
    s=A1(j,1);
    A1(j,:)=A1(j,:)-A1(i,:)*s;
end

This small function performs a series of simple matrix manipulations to complete the Gauss-Jordan elimination. It should be noted that the indexing variable 'i' must be passed in the function call so that the manipulations happen at the right place. This function copies the matrix 'A' into a new matrix 'A1' and then selects an element 's' from the matrix.

The m-file can now be saved as 'gauss_jordan_elim.m' The following are several sample outputs from the code illustrating the ways to call the function with different types of return arguments.

a =

     0     1    -1
     1     3     2
    -1     2     1
	
b =

     4
     4
     8
	 
>> gauss_jordan_elim(a,b)

ans =

    -3
     3
    -1
	
a =

     8     3     7     4
     6     3     3     7
     8     3     7     4
     7     5     6     4

b =

     9
     6
     5
     9
	 
>> gauss_jordan_elim(a,b)
error: matrix is singular

ans =

    0.5000
    0.5000
    0.5000
   -4.0000

The entire matlab script can be seen below or downloaded here (right-click to save):

function [x,err]=gauss_jordan_elim(A,b)

[n,m]=size(A); % size of matrix A
err =0;
x=zeros(n,1);
if n ~= m
    disp('error: n~=m');
    err = 1;
end
if length(b) ~= n 
    disp('error: wrong size of b');
    err = 2;
else
    if size(b,2) ~= 1
        b=b';
    end
    if size(b,2) ~= 1
        disp('error: b is a matrix');
        err = 3;
    end
end
if err == 0
    Aa=[A,b];
    for i=1:n 
        [Aa(i:n,i:n+1),err]=gauss_pivot(Aa(i:n,i:n+1));
        if err == 0
            Aa(1:n,i:n+1)=gauss_jordan_step(Aa(1:n,i:n+1),i);
        end
    end
    x=Aa(:,n+1);
end
A=0;

%----------------------------------------------------

function A1=gauss_jordan_step(A,i)

[n,m]=size(A); % size of matrix A
A1=A;
s=A1(i,1);
A1(i,:) = A(i,:)/s;
k=[[1:i-1],[i+1:n]];
for j=k
    s=A1(j,1);
    A1(j,:)=A1(j,:)-A1(i,:)*s;
end

%----------------------------------------------------------

function [A1,err]=gauss_pivot(A)

[n,m]=size(A); % size of matrix A
A1=A;
err = 0; % error flag
if A1(1,1) == 0
    check = logical(1); % logical(1) - TRUE
    i = 1;
    while check
        i = i + 1;
        if i > n  
            disp('error: matrix is singular');
            err = 1;
            check = logical(0);
        else
            if A(i,1) ~= 0 & check
                check = logical(0);
                b=A1(i,:);      % next three operations exchange rows 1 and i
                A1(i,:)=A1(1,:);
                A1(1,:)=b;
            end
        end
    end
end