Another service from Omega

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,

> alias(P=swaprow,Q=mulrow,R=addrow);

                                  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>
Last modified: Wed May 27 12:03:46 EDT 1998