Gaussian Elimination not working

Question

I tried to write a program that does LU split using Gaussian elimination. My program splits Matrix A into A=LR where L,R are triangular matrices. This works perfectly well. Let Ax=b be a system of equation. My input to the program is (A,b) and I want the operations that are done to get the upper triangular matrix R to be applied to b like we do it in school to solve systems using gaussian elimination. That part seems not to be working.

Can someone give me a hint why its not working ?

function [ x ] = Gauss( A,b )
n=length(b);
for j=1:n-1
    if A(j,j)==0
        break;
    end;

    for i=j+1:n
       A(i,j)=A(i,j)/A(j,j);
       b(j)=b(i)-A(i,j)*b(j)

       for k=j+1:n
          A(i,k)=A(i,k)-A(i,j)*A(j,k);
       end
    end
end
 b
end

Answer 1

You say that your code creates the upper triangular matrix of A correctly, but it doesn't. Let me show you an example.

Let A and b be

A =

 3     2     3     4
 4     3     2     1
 1     0     4     0
 0     5     0     3

b =

 2
 4
 6
 7

If we run your code as it is and look at A and b we get

A =

3.0000    2.0000    3.0000    4.0000
1.3333    0.3333   -2.0000   -4.3333
0.3333   -2.0000   -1.0000  -10.0000
     0   15.0000  -30.0000 -232.0000


b =

7.0000
-203.0000
187.0000
7.0000

Which is neither a triangular matrix, nor the b that we expected. But if we modify slightly your program to be:

function [ x ] = Gauss( A,b )
n=length(b);
for j=1:n-1
    if A(j,j)==0
        break;
    end;

    for i=j+1:n
       f=A(i,j)/A(j,j);  %%Save the proportion between the rows in a
                         %%different variable outside the matrix, or
                         %%you will loose the value that was originally there


       b(i)=b(i)-f*b(j);   %%The operation has to be done in the row you are currently working

       for k=1:n    %%You have to make the operation in the full row,
                    %%not only in the remaining columns, also you can 
                    %%make this without a for loop using `:`
                    %%indexing, but if you dont know about it, 
                    %%leave as it is, it works
          A(i,k)=A(i,k)-f*A(j,k);
       end
    end
end
A
b
end

You get this result

A =

3.0000    2.0000    3.0000    4.0000
     0    0.3333   -2.0000   -4.3333
     0         0   -1.0000  -10.0000
     0         0         0 -232.0000


b =

2.0000
1.3333
8.0000
227.0000

Which is a upper triangular matrix, and the b we want. Hopefully you can take it from here, just as a reference the next steps should look like this

A =

1.0000    0.6667    1.0000    1.3333
     0    1.0000   -6.0000  -13.0000
     0         0    1.0000   10.0000
     0         0         0    1.0000


b =

0.6667
4.0000
-8.0000
-0.9784

then

 A =

 1     0     0     0
 0     1     0     0
 0     0     1     0
 0     0     0     1


b =

-1.1379
1.9871
1.7845
-0.9784

In which A is already an identity matrix, which mean that b is already our answer, that we can corroborate doing

A\b

ans=

-1.1379
1.9871
1.7845
-0.9784