CU Mechanical EngineeringNumerical Methods Matlab tutorials 

Writing Matlab Functions: GaussJordan elimination
In this example, we will write a function that solves a system of linear equations Ax = b using GaussJordan elimination with partial pivoting.
To begin, all functions within matlab must begin by declaring that mfile is a function. This is accomplished by having the first line of the mfile 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 GaussJordan 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:i1],[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 GaussJordan 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 mfile 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 (rightclick 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:i1],[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
Click here to go back to the tutorials home page.