ECE 5340/6340 Lecture 5:            Gaussian Elimination

 

 

METHODS OF SOLVING MATRIX EQUATIONS:        

a)   Direct

·       Gaussian Elimination  <<< We will study

·       LU Decomposition

b)  Iterative

·       SOR:  Successive Over-Relaxation <<

·       Conjugate Gradient Method

 

GAUSSIAN ELIMINATION

 

Example:  Solve

 

                   3x1    + 2x2          + 4x3          = 19

                   2x1    + 6x2          + 5x3          = 29

                   x1      + x2            + x3            = 6

 

Exact solution: x1 = 1,  x2 =2,  x3 = 3

 

1) Write as an augmented matrix:   A|b

 

3 2 4 | 19

2 6 5 | 29

1 1 1 | 6

 

 

 

2)  Eliminate all but x1 from first equation:

 

  (3    2       4     |           19  )* -2/3

+(2    6       5     |           29  )

   0    14/3 7/3  | 49/3  << New second row

 

  (3    2       4     |           19  )* -1/3

+(1    1       1     |           6   )

                      0    1/3    -1/3 |           -1/3  << New third row

 

                   3       2       4     |           19     << New matrix

0       14/3 7/3  | 49/3

0       1/3    -1/3 |           -1/3

 

3) Eliminate x2 from all but second row:

    (0   14/3 7/3  | 49/3) * (-14/3)/(1/3)

+ (0   1/3    -1/3 |           -1/3)

    0    0     -1/2  |         -3/2   << New third row

 

3       2       4     |           19     << New matrix

0       14/3 7/3  | 49/3

0        0     -1/2  |         -3/2

 

4) Solve using Back Substitution starting at bottom:

 

-1/2 x3 = -3/2  ŕ x3 = 3

 

14/3 x2 + 7/3 x3 = 49/3 ŕ x2 = 2

 

3 x1 + 2 x2 + 4 x3 = 19 ŕ x1 = 1

 

 

 

 

5)  Check your solution in original matrix:

 

| 3 2 4 |  |1| = |19|

| 2 6 5 |  |2|    |29|

| 1 1 1 |  |6|    | 6 |

 

WHAT CAN GO WRONG (Besides programming error)

·       Pivoting

·       Scaling

·       Ill-conditioning

 

NUMERICAL ERROR and Floating point arithmetic:

          Reading (Handout)

 

Computer represents real numbers by truncated real numbers.  This causes numerical errors, which are generally small.  IF, however, you are working in a range where these errors are significant, then they can cause problems.

 

Example:  Add 1/3 + 1/3 + 1/3 using 4 bit floating point arithmetic.  (Computers are generally 8- or 16-bit in single precision).  Exact answer:  1

 

          0.3333  = a

          0.3333  = b

      + 0.3333 = c

         0.9999   

Now, suppose you had used this sum as a variable in your program… IF ( (a+b+c) = 1) then_do_something, your program will malfunction because of the floating point arithmetic.

0.9999 does NOT equal integer 1.  It is ZERO!

 

Computers TRUNCATE.  They do not round. 

 

Many numerical errors are CUMULATIVE (they get worse as the problem is larger, or as an iterative solution is allowed to run longer.)

 

PIVOTING:

 

         

The pivot element should be the largest MAGNITUDE element in the row.

 

Example:  Solve (upside down step 2 above)

 

0       1/3    -1/3 |           -1/3

0       14/3 7/3  | 49/3

3       2       4     |           19 

 

TO reduce row 3, you would multiply row 1 by 3/0, which blows up.   Pivoting AT EACH STEP will prevent this.  If element is not zero, but is small, numerical error increases.

 

PARTIAL PIVOTING algorithm:

Interchange rows until largest magnitude element is the pivot element.

         

                   3       2       4     |           19      <<< This is OK for step 1.

0       1/3    -1/3 |           -1/3

0       14/3 7/3  | 49/3

 

Now pivot for step 2:

 

3       2       4     |           19     

0       14/3 7/3  | 49/3

0       1/3    -1/3 |           -1/3

 

SCALING: 

 

Recall that a matrix represents a set of vectors.  If one of the vectors is significantly longer than the others, this does not affect the solution IN THEORY, but it causes numerical errors which affect the solution in practice.

 

Example:  Solve this with 4-bit floating point arithmeetic:

 

          x1 +  1012 x2 = 1 + 1012

          x1 - x2 =  0

 

Answer: x1 = x2 =  1

Reduce:

          x1 +  1012 x2 = 1 + 1012 = 1000000000001

         

            x1 +  1012 x2 = 1012  << (4 decimal places accurate)

          -(x1 -           x2 =  0)

                    1012 x2 = 1012  << (4 decimal places accurate)

 

                             ŕ x2 = 1

          Substitute into second equation: x1 = 1

          Substitute into first equation (just as legit): x1 = 0

 

          This solution is flakey!

 

How do you tell when a solution is scaled poorly, and what do you do about it?

 

          The b vector should all be similar order of magnitude.  One easy way of getting this is to divide each equation through by its bi value.  Then the b vector = 1.

 

         

 

          PIVOTING ALGORITHM

 

At each step of gaussian elimination, put the largest

element in the column on the diagonal.

 

 

DO  L = 1,M-1                                                            ! which step of the elimination you are on

 

c --- Find pivot element and location --

            pivot = 0                                              ! initialize pivot element

            ipivot = 0                                              ! initialize pivot row location

            DO I = L , M                                       ! find pivot element

                        IF ( | a(I,L) | > pivot) THEN      

                                    pivot = a(I,L)                ! store new pivot element

                                    ipivot = I                      ! store location of new pivot element

                        ENDIF

            ENDDO

 

c --- Exchange Lth row and pivot row to put pivot element on top ---

            DO J = L , N                                                   ! for each non-zero element in the row

                        holder = a(L,J)                         ! hold the value currently in the top row

                        a(L,J) = a(ipivot,J)                    ! move the element in ipivot row to top row

                        a(ipivot,J) = holder                   ! put top row element into ipivot row

            ENDDO

 

ENDDO

 

 

SCALING ALGORITHM (One of many methods)

 

Before your start elimination, find the magnitude of each vector (row in the array), and make it a unit vector.  This makes all the vectors the same size.

 

DO I = 1 , M                                                   ! for Each row

 

C – Find length of each vector (matrix row)

DO J = 1 , N

            vector_length = vector_length +  a(I ,J ) 2

ENDDO

vector_length = sqrt(vector_length)

 

--Scale each vector  to a unit length (=1.0) –

DO J = 1 , N

            a(I , J ) = a(I , J) / vector_length

ENDDO

 

ENDDO