Gauss elimination to solve A*x = b linear system (MATLAB)

Question

I'm trying to make a code that solves A*x = b, linear systems.

I made the code below using the gauss elimination process, and it works everytime if A doesn't have any 0's in it. If A has zeros in it, then sometimes it works, sometimes it doesn't. Basically I'm trying an alternative to the "A\b" in MATLAB.

Is there a better/simpler way of doing this?

A = randn(5,5);
b = randn(5,1);

nn = size(A);
n = nn(1,1);
U = A;
u = b;

for c = 1:1:n
    k = U(:,c);
    for r = n:-1:c
        if k(r,1) == 0
            continue;
        else
            U(r,:) = U(r,:)/k(r,1);
            u(r,1) = u(r,1)/k(r,1);
        end
    end
    for r = n:-1:(c+1)
        if k(r,1) == 0
            continue;
        else
            U(r,:) = U(r,:) - U(r-1,:);
            u(r,1) = u(r,1) - u(r-1,1);
        end
    end
end
x = zeros(size(b));
for r = n:-1:1
    if r == n
        x(r,1) = u(r,1);
    else
        x(r,1) = u(r,1);
        x(r,1) = x(r,1) - U(r,r+1:n)*x(r+1:n,1);
    end
end
error = A*x - b;
for i = 1:1:n
    if abs(error(i)) > 0.001
        disp('ERROR!');
        break;
    else
        continue;
    end
end
disp('x:');
disp(x);

Working example with 0's:

A = [1, 3, 1, 3;
    3, 4, 4, 1;
    3, 0, 3, 9;
    0, 4, 0, 1];

b = [3;
    4;
    5;
    6];

Example that fails (A*x-b isn't [0])

A = [1, 3, 1, 3;
    3, 4, 4, 1;
    0, 0, 3, 9;
    0, 4, 0, 1];
b = [3;
    4;
    5;
    6];

Explanation of my algorithm: Lets say I have the following A matrix:

|4, 1, 9|
|3, 4, 5|
|1, 3, 5|

For the first column, I divide each line by the first number in the row, so every row starts with 1

|1, 1/4, 9/4|
|1, 4/3, 5/3|
|1,   3,   5|

Then I subtract the last row with the one above it, and then I'll do the same for the row above and so on.

|1,     1/4,     9/4|
|0, 4/3-1/4, 5/3-9/4|
|0,   3-4/3,   5-5/3|

|1,    0.25,   2.250|
|0,   1.083, -0.5833|
|0,   1.667,   3.333|

Then I repeat the same for the rest of the columns.

|1,    0.25,   2.250|
|0,       1, -0.5385|
|0,       1,   1.999|

|1,    0.25,    2.250|
|0,       1,  -0.5385|
|0,       0,  -8.7700|

|1,    0.25,    2.250|
|0,       1,  -0.5385|
|0,       0,        1|

The same operations I do in A I do in b so the system stays equivalent.

re-UPDATE:

I added this right after "for c = 1:1:n"

So before doing anything it sorts the rows of A (and b) in order to make the "c" column have decrescent entries (0's will be left on the bottom rows of A). Right now it seems to work for any invertible square matrix, although I'm not sure it will.

r = c;
a = r + 1;
while r <= n
    if r == n
        r = r + 1;
    elseif a <= n
        while a <= n
            if abs(U(r,c)) < abs(U(a,c))
                UU = U(r,:);
                U(r,:) = U(a,:);
                U(a,:) = UU;
                uu = u(r,1);
                u(r,1) = u(a,1);
                u(a,1) = uu;
            else
                a = a+1;
            end
        end
    else
        r = r+1;
        a = r+1;
    end
end

Answer 1

Try this:

Ab = [A,b] % Extended matrix of the system of equations
rref(Ab)   % Result of applying the Gauss-Jordan elimination to the extended matrix

See rref documentation for more details and examples.

Answer 2

Gaussian elimination with pivoting is as following.

function [L,U,P] = my_lu_piv(A)
n = size(A,1);
I = eye(n);
O = zeros(n);
L = I;
U = O;
P = I;
    function change_rows(k,p)
        x = P(k,:); P(k,:) = P(p,:); P(p,:) = x;
        x = A(k,:); A(k,:) = A(p,:); A(p,:) = x;
        x = v(k); v(k) = v(p); v(p) = x;
    end

    function change_L(k,p)
        x = L(k,1:k-1); L(k,1:k-1) = L(p,1:k-1);
        L(p,1:k-1) = x;
    end
for k = 1:n
    if k == 1, v(k:n) = A(k:n,k);
    else
        z = L(1:k-1,1:k -1)\ A(1:k-1,k);
        U(1:k-1,k) = z;
        v(k:n) = A(k:n,k)-L(k:n,1:k-1)*z;
    end
    if k<n
        x = v(k:n); p = (k-1)+find(abs(x) == max(abs(x))); % find index p
        change_rows(k,p);
        L(k+1:n,k) = v(k+1:n)/v(k);
        if k > 1, change_L(k,p); end
    end
    U(k,k) = v(k);
end
end

In order to solve the system..

% Ax = b (1) original system % LU = PA
(2) factorization of PA or A(p,:) into the product LU % PAx = Pb (3) multiply both sides of (1) by P % LUx = Pb
(4) substitute (2) into (3) % let y = Ux (5) define y as Ux % let c = Pb (6) define c as Pb % Ly = c
(7) subsitute (5) and (6) into (4) % U*x = y (8) a rewrite of (5)

To do this..

% [L U p] = lu (A) ; % factorize % y = L \ (P*b) ; % forward solve of (7), a lower triangular system % x = U \ y ; % backsolve of (8), an upper triangular system

Answer 3

Gaussian algorithm assumes that the matrix is converted to an upper triangular matrix. This does not happen in your example. The result of your algorithm is

A =

   1   3   1   3
   3   4   4   1
   0   0   3   9
   0   4   0   1

U =

   1.00000   3.00000   1.00000   3.00000
  -0.00000   1.00000  -0.20000   1.60000
   0.00000   0.00000   1.00000   3.00000
   0.00000   4.00000  -0.00000   1.00000

As you can see, it's not upper triangular. You are skipping rows, if the pivot element is zero. That does not work. To fix this you need to swap columns in the matrix and rows in the vector if the pivot element is zero. At the end you have to swap back rows in your result b resp. u.

Gaussian algorithm is:

1 Set n = 1
2 Take pivot element (n, n)
3 If (n, n) == 0, swap column n with column m, so that m > n and (n, m) != 0 (swap row m and n in vector b)
4 Divide n-th row by pivot element (divide n-th row in vector b)
5 For each m > n
6   If (m, n) != 0
7      Divide row m by m and subtract element-wise row n (same for vector b)
8 n = n + 1
9 If n <= number of rows, go to line 2

In terms of numerical stability it would be best to use the maximum of each row as pivot element. Also you can use the maximum of the matrix as pivot element by swapping columns and rows. But remember to swap in b and to swap back in your solution.