# Example of Gaussian Elimination

 Consider the following system of 4 linear equations in in the 4 unknows x,y,z and w, ```2x - y + 3z + 4w = 9 x - 2z + 7w = 11 3x - 3y + z + 5w = 8 2x + y + 4z + 4w = 10 ``` The augmented matrix for this system is given by the matrix A, where,

> A := matrix(4,5,[2,-1,3,4,9,
> 1,0,-2,7,11,
> 3,-3,1,5,8,
> 2,1,4,4,10]);

```                             [2    -1     3    4     9]
[                        ]
[1     0    -2    7    11]
A := [                        ]
[3    -3     1    5     8]
[                        ]
[2     1     4    4    10]```
 To save typing let us alias the elementary row operations to P,Q and R,

`                                  I, P, Q, R`
> A1 := P(A,1,2);
```                             [1     0    -2    7    11]
[                        ]
[2    -1     3    4     9]
A1 := [                        ]
[3    -3     1    5     8]
[                        ]
[2     1     4    4    10]```
> A2 := R(A1,1,2,-2);
```                            [1     0    -2      7     11]
[                           ]
[0    -1     7    -10    -13]
A2 := [                           ]
[3    -3     1      5      8]
[                           ]
[2     1     4      4     10]```
> A3 := R(A2,1,3,-3);
```                            [1     0    -2      7     11]
[                           ]
[0    -1     7    -10    -13]
A3 := [                           ]
[0    -3     7    -16    -25]
[                           ]
[2     1     4      4     10]```
> A4 := R(A3,1,4,-2);
```                            [1     0    -2      7     11]
[                           ]
[0    -1     7    -10    -13]
A4 := [                           ]
[0    -3     7    -16    -25]
[                           ]
[0     1     8    -10    -12]```
> A5 := Q(A4,2,-1);
```                            [1     0    -2      7     11]
[                           ]
[0     1    -7     10     13]
A5 := [                           ]
[0    -3     7    -16    -25]
[                           ]
[0     1     8    -10    -12]```
> A6 := R(A5,2,3,3);
```                            [1    0     -2      7     11]
[                           ]
[0    1     -7     10     13]
A6 := [                           ]
[0    0    -14     14     14]
[                           ]
[0    1      8    -10    -12]```
> A7 := R(A6,2,4,-1);
```                            [1    0     -2      7     11]
[                           ]
[0    1     -7     10     13]
A7 := [                           ]
[0    0    -14     14     14]
[                           ]
[0    0     15    -20    -25]```
> A8 := Q(A7,3,-1/14);
```                            [1    0    -2      7     11]
[                          ]
[0    1    -7     10     13]
A8 := [                          ]
[0    0     1     -1     -1]
[                          ]
[0    0    15    -20    -25]```
> A9 := R(A8,3,4,-15);
```                             [1    0    -2     7     11]
[                         ]
[0    1    -7    10     13]
A9 := [                         ]
[0    0     1    -1     -1]
[                         ]
[0    0     0    -5    -10]```
> A10 := Q(A9,4,-1/5);
```                              [1    0    -2     7    11]
[                        ]
[0    1    -7    10    13]
A10 := [                        ]
[0    0     1    -1    -1]
[                        ]
[0    0     0     1     2]```
> SOL := backsub(A10);
`                             SOL := [-1, 0, 1, 2]`
 The previous elementary row operations can also be encoded with matrices. These matrices are the corresponding operation applied to the four by four identity matrix Id.

> Id := matrix(4,4,[1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1]);

```                                 [1    0    0    0]
[                ]
[0    1    0    0]
Id := [                ]
[0    0    1    0]
[                ]
[0    0    0    1]```
> E1 := P(Id,1,2);
```                                 [0    1    0    0]
[                ]
[1    0    0    0]
E1 := [                ]
[0    0    1    0]
[                ]
[0    0    0    1]```
> E2 := R(Id,1,2,-2);
```                                 [ 1    0    0    0]
[                 ]
[-2    1    0    0]
E2 := [                 ]
[ 0    0    1    0]
[                 ]
[ 0    0    0    1]```
> E3 := R(Id,1,3,-3): E4 := R(Id,1,4,-2): E5:=Q(Id,2,-1):E6:=R(Id,2,3,3):
> E7:=R(Id,2,4,-1):E8:=Q(Id,3,-1/14):E9:=R(Id,3,4,-15):E10:=Q(Id,4,-1/5):
 Let us now continue performing gaussian eliminations from A10 until we get the complete reduced row echelon form for the matrix.

> A11 := R(A10,3,2,7);

```                              [1    0    -2     7    11]
[                        ]
[0    1     0     3     6]
A11 := [                        ]
[0    0     1    -1    -1]
[                        ]
[0    0     0     1     2]```
> A12 := R(A11,3,1,2);
```                              [1    0    0     5     9]
[                       ]
[0    1    0     3     6]
A12 := [                       ]
[0    0    1    -1    -1]
[                       ]
[0    0    0     1     2]```
> A13 := R(A12,4,1,-5);
```                              [1    0    0     0    -1]
[                       ]
[0    1    0     3     6]
A13 := [                       ]
[0    0    1    -1    -1]
[                       ]
[0    0    0     1     2]```
> A14 := R(A13,4,2,-3);
```                              [1    0    0     0    -1]
[                       ]
[0    1    0     0     0]
A14 := [                       ]
[0    0    1    -1    -1]
[                       ]
[0    0    0     1     2]```
> A15 := R(A14,4,3,1);
```                               [1    0    0    0    -1]
[                      ]
[0    1    0    0     0]
A15 := [                      ]
[0    0    1    0     1]
[                      ]
[0    0    0    1     2]```
 The last column now contains the solution to the original system of equations. The above elementary operations are:

> E11:=R(Id,3,2,7):E12:=R(Id,3,1,2):E13:=R(Id,4,1,-5):E14:=R(Id,4,2,-3):
> E15 := R(Id,4,3,1):
 We now compute the product of all the matrices that represent the elementary row operations that we used to transform the augmented matrix A for the system into A15, the rref.

> Minv := E15 &* E14 &* E13 &* E12 &* E11 &* E10 &* E9 &*
> E8 &* E7 &* E6 &* E5 &* E4 &* E3 &* E2 &* E1:
> M := submatrix(A,1..4,1..4);

```                                [2    -1     3    4]
[                  ]
[1     0    -2    7]
M := [                  ]
[3    -3     1    5]
[                  ]
[2     1     4    4]```
 Let us check that Minv is in fact the inverse of M where A = [M,b],

> evalm(M &* Minv);

```                              [1    0    0    0]
[                ]
[0    1    0    0]
[                ]
[0    0    1    0]
[                ]
[0    0    0    1]```
> evalm(Minv);
```                         [-25     -3      13         ]
[---     --      --      1  ]
[14      14      14         ]
[                           ]
[-29                        ]
[---    1/35    1/7     3/5 ]
[35                         ]
[                           ]
[23      -2                 ]
[--      --     -2/7    -1/5]
[35      35                 ]
[                           ]
[31      11      -3         ]
[--      --      --     -1/5]
[70      70      14         ]```
 The inverse of M can also be obtained simply by the command

> inverse(M);

```                         [-25     -3      13         ]
[---     --      --      1  ]
[14      14      14         ]
[                           ]
[-29                        ]
[---    1/35    1/7     3/5 ]
[35                         ]
[                           ]
[23      -2                 ]
[--      --     -2/7    -1/5]
[35      35                 ]
[                           ]
[31      11      -3         ]
[--      --      --     -1/5]
[70      70      14         ]```
 The previous method however can, in principle, be done by hand and it is used to show important theorems about the existance of inverses and the nature of the solutions of a system of linear equations.

Link to the commands in this file
Carlos Rodriguez <carlos@math.albany.edu>